Two-Fluid RANS-RSTM-PDF Model for Turbulent Particulate Flows

Abstract

A novel three-dimensional (3D) model based on Reynolds turbulence stress model (RSTM) closure of equations of carrier and particulate phases was elaborated for channel turbulent flows. The essence of the model is the direct calculation of normal and shear components of the Reynolds stresses for the particulate phase similar to the carrier fluid. The model is based on the Eulerian approach, which is applied for the 3D RANS modeling of the carrier flow and the particulate phase and the statistical probability dense function (PDF) approach focusing on the mathematical description of the second moments of the particulate phase.

1. Introduction

Turbulent channel particulate flows have numerous engineering applications ranging from pneumatic conveying systems to coal gasifiers, chemical reactor design and are one of the most thoroughly investigated subject in the area of the particulate flows. These flows are very complex and influenced by various physical phenomena, such as particle-turbulence and particle-particle interactions, deposition, gravitational and viscous drag forces, particle rotation and lift forces, turbulent dispersion, etc.

One of the most applicable approaches in computational fluid dynamics (CFD) is the Reynolds-averaged Navier-stokes (RANS) equation approach, which is found in many industrial implementations and is very likely to be claimed and applied in the foreseeable future.

The procedure of analysis of the predicted fluid parameters becomes much more complicated in case of the three-dimensional (3D) dense particulate flows, with additional inclusion of two or/and four coupling phenomena of inter-particle collisions. Within a frame of the RANS modeling, one of the challenging and advantaged theoretical approaches is the Reynolds turbulence stress model (RSTM), since k-ε closure model cannot describe properly the flow of complex geometry.

For the proper modeling of the particulate flows, which may include multifold processes, e.g., the particle-turbulence, or/and particle-particle, and particle-wall interactions and other relevant effects, one can apply the probability dense function (PDF) approach, which gives reasonable formalism for the closure of the governing mass and momentum equations of the particulate phase. Within the PDF approach, the closure of the governing equations of the particulate phase is based on a solution of the differential transport equations written for each particle velocity covariance, taking into account possible mechanisms of the particle-turbulence and particle-particle interactions. Such procedure is similar to the RSTM closure.

Currently, the probability dense function (PDF) approach is widely applied for the numerical modeling of the particulate flows. The PDF models, e.g., [1, 2], contain more complete differential transport equations, which are written for various velocity correlations and consider both the turbulence augmentation and attenuation due to the presence of particles.

As opposed to the round channel flows, the rectangular or/and square channel flows, even in case of unladen flows, are considerably anisotropic with respect to the components of the turbulence energy, that is more expressed near the channel walls and corners being notable for the secondary flows. In addition, the presence of particles enhances such anisotropy. Such particulate flows are studied by RSTM, which are based on the transport equations written for all components of the Reynolds stress tensor-associated with the particulate phase.

The RSTM approach allows to completely analyze the influence of particles on axial, transverse, and spanwise components of the turbulence kinetic energy, including also possible modifications of the cross-correlation velocity moments.

A number of studies based on the RSTM approach showed its good performance and capability for simulation of the complicated flows [3], as well for the turbulent subsonic [4] and supersonic flows [5] and viscoelastic flows [6].

Taulbee [7] used the RSTM approach to calculate the particle-laden shear flow by applying the direct numerical simulation (DNS). The flow Reynolds number was low (Re = 952). Therefore, the method by [7] cannot be applied to the real turbulent flows characterized by considerably higher Reynolds numbers, unlike the present numerical simulation which handles with the flow Reynolds numbers of 56,000 and 140,000.

The RSTM approach was also applied in the model [8]. The only difference was in the fact that the closure of equations of motion of the particulate phase was based on the Boussinesq hypothesis, where the turbulent viscosity of the particulate phase was introduced and calculated using the algebraic expressions obtained in the PDF approach [9]. However, the model [8] was inherently eclectic, since it applied the RSTM closure for a carrier flow and at the same time the Boussinesq hypothesis for the particulate phase.

Recently, Mukin and Zaichik [10] have proposed the nonlinear algebraic Reynolds stress model for the gas flow laden with small heavy particles based on the PDF approach. The original equations written for each component of Reynolds stress were reduced to their general form in terms of the turbulence energy and its dissipation rate with additional effect of the particulate phase. However, the model [10] does not enable a direct solution of the differential transport equations and it applies the k-ε solution.

The study of the particle dispersion that occurs in the velocity uniform shear turbulent flow assumes knowing the internal structure, general relationships, and methods of the flow generation.

Based on the analysis of results of numerous investigations of turbulent dispersion of finite inertia particles, it should be singled out three effects that are of high importance and should be considered: (i) the inertia effect implying that the dispersion of solid particles might exceed the dispersion of the fluid particles in the absence of a body force [11, 12]; (ii) the crossing trajectory effect [13] meaning that in the presence of a drift velocity a finite inertia particle will disperse less than a fluid particle; and (iii) the continuity effect [14] where the dispersion in the direction of the drift velocity exceeds the dispersion in other two directions.

In the mathematical description of the particle turbulent dispersion, there are a number of models and numerical simulations, which can be classified into the Lagrangian and Eulerian (two-fluid) approximations. They relate to the specific flow structures, such as confined flows in pipes and channels or free jets, wakes, and wall boundary flows.

Taylor [15] made the Lagrangian analysis of the dispersion of a particle in a stationary homogeneous turbulence, which showed that turbulent dispersion varies in time and derived the asymptotic expressions of dispersion for short and large times. Later on, the Lagrangian approach for the turbulent dispersion of particles was further developed in [16–21].

One of the most comprehensive numerical research of the particle dispersion in the uniform shear flow, based on the Lagrangian approach, was carried out in [22]. Ahmed and Elghobashi [22] applied DNS to the investigation of the particle inertia effect, effect of direction of gravitation as well as magnitude of a shear number. In particular, Ahmed and Elghobashi [22] have revealed that gravity reduces the particle turbulent dispersion and diffusivity in all directions due to the drift velocity effects.

On the contrary, there are very less studies on the particle turbulent dispersion, which are based on the Eulerian approach. The most significant are the numerical investigations [9, 23], where the statistical PDF models were applied to the particle behavior in the turbulent flows.

The 3D RSTM approach presented here applies the closure of equations of motion of the particulate phase, which is carried out similarly to the closure of the carrier flow, i.e., the equations are written for the normal and shear components of the Reynolds stress. The Reynolds stress equations are derived from the PDF model [9] and presented in a general case. The advantage of the given model is in use of the same closure for both the carrier flow and particulate phase, namely, the Reynolds differential equations.

The given 3D RSTM model has been applied for the turbulent dispersion of solid particles in a turbulent horizontal channel flow imposed to uniform shear.

The obtained numerical results have been verified and validated by comparison with the experimental data.

2. Numerical model

The present numerical model being proposed for the evaluation of the particle dispersion is the Eulerian approach, which applies both the 3D RANS modeling of the carrier flow and the particulate phase [24] and the statistical PDF approach focusing on the mathematical description of the second moments of the particulate phase [9]. Within the Eulerian approach, the particulate phase is considered as the diluted medium; therefore, the effect of the particle collision is negligible that means the application of the one-way coupling.

The numerical simulation considered the turbulent dispersion of solid particles in horizontal channel uniform shear turbulent flow for two different cases: i) shear of the mean flow velocity is along the direction of gravity (Figure 1a) and ii) shear of the mean flow velocity is directed normally to gravity (Figure 1b). Here u is the mean axial velocity of gas.

Figure 1.

Shear of the mean flow velocity in horizontal channel flow. (a) shear of the mean flow velocity is along the direction of gravity; (b) shear of the mean flow velocity is directed normally to gravity.

The particles were brought into the uniform shear gas flow, which has been preliminarily computed to obtain the velocity flow field.

The system of the momentum and closure equations of the gas phase are identical for the unladen flows, while the particle-laden flows are under the impact of the viscous drag force. The Cartesian coordinates are used here.

2.1. Governing equations for the particulate phase

The 3D governing equations for the particulate phase are written as follows:

The particle mass conservation equation:

$\frac{\partial \left(\alpha u_{\text{s}}\right)}{\partial x}+\frac{\partial \left(\alpha v_{\text{s}}\right)}{\partial y}+\frac{\partial \left(\alpha w_{\text{s}}\right)}{\partial z}=\frac{\partial }{\partial x}{D}_{\text{s}}\frac{\partial \alpha }{\partial x}+\frac{\partial }{\partial y}{D}_{\text{s}}\frac{\partial \alpha }{\partial y}+\frac{\partial }{\partial z}{D}_{\text{s}}\frac{\partial \alpha }{\partial z}$E1

x-component of the momentum equation:

$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}^2-{\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}\right]+\frac{\partial }{\partial y}\alpha \left[u_{\text{s}}{v}_{\text{s}}-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}\right]+\frac{\partial }{\partial z}\alpha \left[u_{\text{s}}{w}_{\text{s}}-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}\right]=\alpha {C^{\prime }}_{\text{D}}\frac{\left(u-u_{\text{s}}\right)}{{\tau }_{\text{p}}}\hfill \\ +\alpha \text{ sgn}\left(g_x\right)\left(1-\frac{\rho }{{\rho }_{\text{p}}}\right)\hfill \end{array}$E2

y-component of the momentum equation:

$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{v}_{\text{s}}-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}\right]+\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}^2-{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}\right]+\frac{\partial }{\partial z}\alpha \left[v_{\text{s}}{w}_{\text{s}}-{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}\right]=\alpha {C^{\prime }}_{\text{D}}\frac{\left(v-v_{\text{s}}\right)}{{\tau }_{\text{p}}}\hfill \\ +\alpha \mathrm{sgn}\left(g_y\right)\left(1-\frac{\rho }{{\rho }_{\text{p}}}\right)\hfill \end{array}$E3

z-component of the momentum equation:

$\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{w}_{\text{s}}-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}\right]+\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{w}_{\text{s}}-{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}\right]+\frac{\partial }{\partial z}\alpha \left[w_{\text{s}}^2-{\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}\right]=\alpha {C^{\prime }}_{\text{D}}\frac{\left(w-w_{\text{s}}\right)}{{\tau }_{\text{p}}}$E4

where u, v, and w are the axial, transverse, and spanwise time-averaged velocity components of gas, respectively; u_s, v_s and w_s are the axial, transverse, and spanwise time-averaged velocity components of particulate phase, respectively; ρ is the material density of gas; ρ_p is the material density of particles; α is the particle mass concentration; g_y is the y-component of gravity.

The relative friction coefficient $C_{\text{D}}^{\text{'}}$ is calculated according to [25].

The closure model for the transport equations of the particulate phase was applied to the PDF model [26], where $D_{\text{s}}$ is the coefficient of the turbulent diffusion of the particulate phase.

The equations for the second-order moments of the fluctuating velocity (turbulent stresses) of the particulate phase are written based on the PDF approach in [9]. These equations describe convective and diffusive transfer, generation of particle velocity fluctuations due to the velocity gradient, generation of fluctuations resulting from particle entrainment into the fluctuating motion of carrier gas flow, and dissipation of turbulent stresses of the particulate phase caused by interfacial forces:

Equation of the x-normal component of the Reynolds stress:

$\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}-{\tau }_{\text{p}}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}\right]+\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)}{3}\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}\right]$

$\begin{array}{l}+\frac{\partial }{\partial z}\alpha [w_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{w^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}]\\ =\frac{\partial }{\partial x}\left\{\alpha {\tau }_{\text{p}}\left[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}\right]\right\}\end{array}$

$+\frac{\partial }{\partial y}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}$

$\begin{array}{l}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}]\}\\ +\frac{\partial }{\partial z}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}\end{array}$

$+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}$

$+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}]\}$E5

Equation of the y-normal component of the Reynolds stress:

$\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}\right]+\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}-{\tau }_{\text{p}}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}\right]$

$\begin{array}{l}+\frac{\partial }{\partial z}\alpha [w_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{w^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}]\\ =\frac{\partial }{\partial x}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}\end{array}$

$+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}]\}$

$+\frac{\partial }{\partial y}\{\alpha {\tau }_{\text{p}}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}]\}$

$+\frac{\partial }{\partial z}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}$

$+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}]\}$E6

Equation of the z-normal component of the Reynolds stress:

$\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}\right]+\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}\right]$

$+\frac{\partial }{\partial z}\alpha [w_{\text{s}}{\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}-{\tau }_{\text{p}}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{w^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}]=\frac{\partial }{\partial x}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}$

$\begin{array}{l}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_s}}}{\partial y}\\ +2\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}]\}\end{array}$

$\begin{array}{l}+\frac{\partial }{\partial y}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\\ +2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}+2\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}\end{array}$

$+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}]\}$

$+\frac{\partial }{\partial z}\left\{\alpha {\tau }_{\text{p}}\left[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}\right]\right\}$E7

Equation of the xy shear stress component of the Reynolds stress:

$$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}\right]\\ +\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}\right]+\frac{\partial }{\partial z}\alpha [w_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}\end{array}$$E400

$$-\frac{{\tau }_{\text{p}}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)}{3}\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}]$$

$$=\frac{\partial }{\partial x}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}+\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}$$

$$\begin{array}{l}\times \left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{w}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}]\}\\ +\frac{\partial }{\partial y}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}+\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}\end{array}$$

$$\begin{array}{l}\times \left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}\\ +2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}]\}\end{array}$$

$$\begin{array}{l}+\frac{\partial }{\partial z}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}\\ +\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\end{array}$$

$\begin{array}{l}\times \frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}\\ +\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}]\}\end{array}$E8

Equation of the xz shear stress component of the Reynolds stress:

$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}\right]\\ +\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}\right]\end{array}$E100

$\begin{array}{l}+\frac{\partial }{\partial z}\alpha \left[w_{\text{s}}{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}\right]\\ =\frac{\partial }{\partial x}\{\frac{\alpha {\tau }_{\text{p}}}{3}[2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}+2\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}\end{array}$E101

$\begin{array}{l}\times \left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}\\ +\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}]\}\end{array}$E102

$\begin{array}{l}+\frac{\partial }{\partial y}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}\\ +\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\end{array}$E103

$\begin{array}{l}\times \frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}+\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial z}\\ +\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\end{array}$E104

$\begin{array}{l}\times \frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}]\}+\frac{\partial }{\partial z}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}\\ +\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}\end{array}$E105

$+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}]\}$E9

Equation of the yz shear stress component of the Reynolds stress:

$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{{\tau }_p}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}\right]\\ +\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}\right]\end{array}$

$\begin{array}{l}\frac{\partial }{\partial x}\alpha \left[u_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{{\tau }_p}{3}\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^l{\displaystyle \overline{{u^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}\right]\\ +\frac{\partial }{\partial y}\alpha \left[v_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{2{\tau }_{\text{p}}}{3}\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}\right]\end{array}$

$\begin{array}{l}+\frac{\partial }{\partial z}\alpha \left[w_{\text{s}}{\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}-\frac{2{\tau }_p}{3}\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}\right]\\ =\frac{\partial }{\partial x}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{{v^{\prime }}^2}}\right)\end{array}$

$\begin{array}{l}\times \frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}\\ +\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial y}+\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\end{array}$

$\begin{array}{l}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial z}+\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}\\ +\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^k{\displaystyle \overline{{w^{\prime }}^2}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}]\}\end{array}$

$\begin{array}{l}+\frac{\partial }{\partial z}\{\frac{\alpha {\tau }_{\text{p}}}{3}[\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}+g_{\text{u}}^l{\displaystyle \overline{u^{\prime }{v}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}\\ +\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}+\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}\end{array}$

$+2\left({\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^k{\displaystyle \overline{u^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}+2\left({\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}+g_{\text{u}}^n{\displaystyle \overline{v^{\prime }{w}^{\prime }}}\right)\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial y}]\}$E10

where $g_{\text{u}}^l,g_{\text{u}}^n,\text{and}{g}_{\text{u}}^k$ are the coefficients characterizing the entrainment of particles into the fluctuating motion of the flow [9] for x, y, and z directions, respectively.

2.2. Boundary conditions

The wall conditions are set for the gas at the side walls of the channel based on the control volume method by [27, 28] in a similar way as in the case of the coincidence of shear of the mean flow velocity and gravity.

The numerical simulation considers the turbulent dispersion of solid particles in horizontal channel uniform shear turbulent flow for two different cases: i) shear of the mean flow velocity is along the direction of gravity (Figure 1a) and ii) shear of the mean flow velocity is directed normally to gravity (Figure 1b). Therefore, two sets of the boundary conditions are used for the calculations.

The boundary conditions for the particulate phase are set at the flow axis as follows:

Case 1 for z = 0:

$\frac{\partial u_{\text{s}}}{\partial z}=v_{\text{s}}=w_{\text{s}}=\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial z}=\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial z}=\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial z}={\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}={\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}={\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}=\frac{\partial \alpha }{\partial z}=0$E11

Case 2 for y = 0:

$\frac{\partial u_{\text{s}}}{\partial y}=v_{\text{s}}=w_{\text{s}}=\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial y}=\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial y}=\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial y}={\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}={\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}={\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}=\frac{\partial \alpha }{\partial y}=0$E12

The boundary conditions for the particulate phase are set at the channel walls according to [9]:

Case 1 for y = 0.5h_y:

$\frac{\partial u_{\text{s}}}{\partial y}=2\left(\frac{1-\chi e_x}{1+\chi e_x}-\frac{1-\chi }{1+\chi }\right){\left(\frac{2}{\pi {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}\right)}^{0.5}\frac{u_{\text{s}}}{{\tau }_{\text{p}}}$E13

$v_{\text{s}}=w_{\text{s}}=\frac{\partial \alpha }{\partial y}=0$E14

$P_{\text{uus}}=-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}\frac{\partial u_{\text{s}}}{\partial y}=P_{\text{uus}}={\eta }_x{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}\frac{\partial u_{\text{s}}}{\partial y}$E15

and applying the expression ${\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}=-{\eta }_x{\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}$, where ${\eta }_x$ is the coefficient of friction between the particles and the wall,

$P_{\text{vvs}}=P_{\text{wws}}=P_{\text{uvs}}=P_{\text{uws}}=P_{\text{vws}}=0$

Case 2 for z = 0.5h_z:

$\frac{\partial u_{\text{s}}}{\partial y}=2\left(\frac{1-\chi e_x}{1+\chi e_x}-\frac{1-\chi }{1+\chi }\right){\left(\frac{2}{\pi {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}\right)}^{0.5}\frac{u_{\text{s}}}{{\tau }_{\text{p}}}$E16

$v_{\text{s}}=w_{\text{s}}=\frac{\partial \alpha }{\partial z}=0$E17

$P_{\text{uus}}=-{\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}\frac{\partial u_{\text{s}}}{\partial z}={\eta }_x{\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}\frac{\partial u_{\text{s}}}{\partial z}$E18

$P_{\text{vvs}}=P_{\text{wws}}=P_{\text{uvs}}=P_{\text{uws}}=P_{\text{vws}}=0$E19

e_x is the coefficient of restitution in the axial direction, which is modeled as:

$e_x=\{\begin{array}{l}1-{\xi }_x\text{, 0}\le {\xi }_x\le \frac{2}{7}\hfill \\ \frac{\text{5}}{\text{7}}\text{, }{\xi }_x>\frac{2}{7}\hfill \end{array}$E20

Here, the parameter ${\xi }_x={\eta }_x\left(1+e_{\text{r}}\right)tan{\theta }_x$, where $e_{\text{r}}$ is the coefficient of restitution of the particle velocity normal to the wall; ${\theta }_x={\tan }^{-1}\left(v_{\text{s}}/u_{\text{s}}\right)$ is the angle of attack between the trajectory of the particle and the wall; $\chi$ is the reflection coefficient, which is the probability of the particles recoiling off the boundaries and back to the flow. The coefficient of restitution reflects the loss of the particle momentum as the particle hits the walls. In the given model, $\chi$= 1/3, $e_{\text{r}}$= 1 and ${\eta }_x$= 0.39 [29].

The conditions for the transverse and spanwise components of the gas velocity are set at the channel walls in terms of impenetrability and no-slip.

The set of boundary conditions for gas and particulate phase at the exit of the channel is written, respectively, as follows:

$\frac{\partial u_{\text{s}}}{\partial x}=\frac{\partial v_{\text{s}}}{\partial x}=\frac{\partial w_{\text{s}}}{\partial x}=\frac{\partial \alpha }{\partial x}=\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}}{\partial x}=\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}^2}}}{\partial x}=\frac{\partial {\displaystyle \overline{{w^{\prime }}_{\text{s}}^2}}}{\partial x}=\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{v^{\prime }}_{\text{s}}}}}{\partial x}=\frac{\partial {\displaystyle \overline{{u^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}=\frac{\partial {\displaystyle \overline{{v^{\prime }}_{\text{s}}{w^{\prime }}_{\text{s}}}}}{\partial x}=0$E21

2.3. Computational method

The control volume method was applied to solve the 3D partial differential equations written for the unladen flow and the particulate phase (Eqs. (1)–(11)), taking into account the boundary conditions (Eqs. (12)–(21)). The governing equations were solved using the implicit lower and upper (ILU) matrix decomposition method with the flux-blending-differed correction and upwind-differencing schemes by [27]. This method is utilized for the calculations of the particulate turbulent flows in channels of the rectangular and square cross sections. The calculations were performed in the dimensional form for all the flow conditions. The number of the control volumes was 1120000.

3. Laboratory experiments

The obtained numerical results have been verified and validated in comparison with the data obtained by the experimental facility of Tallinn University of Technology.

The experimental method for the determination of the particle dispersion was based on recording the particle trajectories by means of a high-speed video camera on separate regions of a flow that locate at various distances from a point source of particles, and the subsequent processing of the frames [30].

The experimental setup for the investigations of particle dispersion (Figure 2) allowed to generate the shear flow similarly to [31] by means of flat plates installed with a varied pitch. The test section was 2 m long with 400 × 200 mm cross section.

Figure 2.

Experimental setup.

Two cases of spatial orientation of shear of the mean flow velocity were investigated. Figure 2 shows the top view of the setup for the case when shear is along the direction of gravity (Figure 1a). For investigations of the particle dispersion when shear is directed normally to gravity (Figure 1b), the setup was turned sideways as a whole at an angle of 90° around the axis of the flow.

The mean flow velocity was 5.1 m/s. Glass spherical particles (physical density of 2500 kg/m³) with an average diameter of 55 μm were used in the experiment runs. The root-mean-square deviation of the diameter of particles did not exceed 0.1. The particles were entered into the flow through the source point which was the L-shaped tubule of 200 μm inner diameter.

All measurements and data processing were carried out at the flow location x = 1212 mm.

The data processing technique [30] was applied to determine the particle spatial displacement along the y-axis, namely D_y, which characterizes quantitatively the particle turbulent dispersion. D_y is calculated as the axial displacement of the maximum value of distribution of the particle mass concentration determined at the location x = 1212 mm relative to the initial flow location that disposes near the exit of the source point.

4. Results and discussions

The numerical results presented below have been obtained at two locations of the flow: initial location signed “ini” and disposed at the exit of the particle source point and the location 2x/h_y = 12.63 from the exit of the particle point source. The turbulent dispersion of 55-μm glass spherical particles was examined. The flow mass loading was about 10⁻⁶ kg dust/kg air.

Figures 3–15 show the numerical data obtained by the presented model for two cases of spatial orientation of shear of the mean flow velocity: shear is along the direction of gravity (case 1), and shear is directed normally to gravity (case 2).

Figure 3.

Transverse distributions of axial velocities of gas and particles, case 1. Here and below u₀ is the mean flow velocity; u₀ = 5.1 m/s.

Figure 4.

Spanwise distributions of axial velocities of gas and particles, case 2.

Figure 5.

Transverse distributions of a turbulence kinetic energy; cases 1 and 2.

Figure 6.

Spanwise distributions of a turbulence kinetic energy, cases 1 and 2.

Figure 7.

Transverse distributions of xy shear stress component of the Reynolds stress of gas and particles, case 1. Here ${\displaystyle \overline{u^{\prime }{v}^{\prime }}}={\displaystyle \overline{u^{\prime }{v}^{\prime }}}/u_0^2$.

Figure 8.

Spanwise distributions of xz shear stress component of the Reynolds stress of gas and particles, case 2. Here ${\displaystyle \overline{u^{\prime }{w}^{\prime }}}={\displaystyle \overline{u^{\prime }{w}^{\prime }}}/u_0^2$.

Figure 9.

Transverse distributions of x-normal components of the Reynolds stress of gas and particles, case 1. Here ${\displaystyle \overline{{u^{\prime }}_{}^2}}={\displaystyle \overline{{u^{\prime }}_{}^2}}/u_0^2$ and ${\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}={\displaystyle \overline{{u^{\prime }}_{\text{s}}^2}}/u_0^2$.

Figure 10.

Transverse distributions of particles mass concentration, case 1.

Figure 11.

Spanwise distributions of axial velocities of gas and particles, case 1, location 2x/h_y = 12.63.

Figure 12.

Spanwise distribution of particles mass concentration, case 1.

Figure 13.

Transverse distribution of particles mass concentration, case 2.

Figure 14.

Spanwise distribution of particles mass concentration, case 2.

Figure 15.

Transverse distributions of axial velocities of gas and particles, cases 1 and 2, location 2x/h_y = 12.63.

Figure 3 shows the transverse distributions of axial velocities of gas and particles for case 1. It is evident that the linear profiles of the averaged axial velocity components of gas and particulate phase across the flow are almost preserved starting from the initial cross section till the pipe exit. Besides, they occupy almost the whole turbulent core of the flow with slight increase of the values in the turbulent core and decrease near the walls due to the effect of a viscous dissipation. The similar profiles are observed with respect of distribution of the same averaged axial velocity components for gas and particulate phase along the spanwise direction (Figure 4).

Since the axial velocity increases toward the bottom wall, the profiles of a turbulence kinetic energy have their higher values near the bottom wall area (Figure 5). However, along the spanwise direction, the profiles of the turbulence kinetic energy are symmetrical, since there is no change of the axial velocity along this direction (Figure 6).

The profiles of the Reynolds shear stresses of gas and particulate phase are shown in Figures 7 and 8. Here it is evident that there is some kind of plateau in the turbulent core. This confirms that we deal with the shear flow; hence, it must be the constant value of the Reynolds shear stresses observed for cases 1 and 2, i.e., for the xy-plane (case 1) and xz-plane (case 2) Reynolds shear stress components. Here, the linear distributions of the averaged axial velocity components across the flow take place along the spanwise direction.

Figure 9 show the transverse distributions of x-normal components of the Reynolds stress of gas and particulate phase obtained for case 1. It can be seen that unlike ${\displaystyle \overline{{u^{\prime }}_{}^2}}$, the maximum value of ${\displaystyle \overline{{u^{\prime }}_s^2}}$ distribution located near the channel top wall is larger than the one near the bottom wall. This is due to the effect of particle inertia and their crosswise motion that cause different axial particle accelerations near the top and bottom walls (Figure 3).

Figures 10–13 present the transverse and spanwise distributions of the particle mass concentration c/c₀ across the flow at the initial location and the location 2x/h_y = 12.63 for both the cases of spatial orientation of shear of the mean flow velocity. Here c₀ is the value of the particle mass concentration at the initial location at the flow axis. These distributions reflect the character of the particle turbulent dispersion that occurs in the given channel shear flow. It is obvious that a) due to gravity the particles go down, and thus the mass concentration profile shifts toward the bottom wall (case 1) and b) the profiles become wider relative to their initial distributions due to the particle turbulent dispersion (Figure 10).

Since in case 1 there is symmetrical distribution of parameters along the spanwise direction (Figures 6 and 11), the symmetrical distribution of the mass concentration along this direction (Figure 12) can be observed, both at the initial and exit cross sections.

A similar situation is observed for case 2, when the linear change of the axial velocity takes place along the spanwise direction. Here the particles go down due to gravity (see Figure 13), and simultaneously there is no shift of the distribution of the mass concentration along the spanwise direction (Figure 14).

Table 1 presents the values of the particle spatial displacement D_y obtained experimentally and numerically for two cases of spatial orientation of shear of the mean flow velocity. This displacement characterized quantitatively the particle turbulent dispersion. It is evident that the numerical values of displacement fit satisfactory with the experimental ones that validate the reliability of the presented model.

	Case 1		Case 2
	Experiment	Modeling	Experiment	Modeling
Dy, mm	46.0	43.7	48.2	50.0

Table 1.

Particle displacement.

Table 1 shows that the particle dispersion in case 1 is smaller than in case 2. This fact can be explained by the particle axial velocity taking place in case 2 is smaller than the one for case 1 in the same y location (Figures 10, 13, and 15).

5. Conclusions

The 3D Reynolds stress turbulence model (RSTM) based on the 3D RANS and statistical PDF approaches has been elaborated for the turbulent dispersion of solid particles in particulate horizontal channel shear flow domain.

The main distinctive feature of the given model is in use of the same closure for both the carrier flow and particulate phase, namely the Reynolds differential equation.

The presented model has several important advantages over the Lagrangian approach:

1. direct simulation of the particle concentration;

2. direct simulation of the particles influence on a carrier flow;

3. there is no basic limit for the parameters of a particulate flow, namely the flow Reynolds number and value of the particle concentration.

Based on the given model, two cases of spatial orientation of shear of the mean flow velocity have been examined. It has been obtained that the effect of orientation of shear appears through decrease of the particle dispersion in case of directional coincidence between shear and gravity as compared with the case of their mutual perpendicularity.

The validity of the elaborated model has been confirmed by experimental investigations of effect of shear of the mean flow velocity on the turbulent particle dispersion.

Acknowledgments

The authors are grateful for the technical support of Texas Advanced Computing Center (TACC) in Austin, USA. The authors are grateful for the fulfilled research to the Estonia-Norway project EMP230 support. This study is related to the activity of the European network action COST MP1106 “Smart and green interfaces—from single bubbles and drops to industrial, environmental, and biomedical applications.”