Numerical simulation of grid-generated turbulent particulate flow by three-dimensional Reynolds stress

A three-dimensional (3D) Reynolds stress turbulence model based on 3D Reynolds-averaged Navier–Stokes equations has been elaborated for grid-generated turbulence in particulate downward flow arranged in the channel domain of the square cross section. The model presented considers both the enhancement and attenuation of turbulence by means of the additional terms of the transport equations of the normal Reynolds stress components. It allows us to carry out calculations covering the long distance of the channel length without using algebraic assumptions for various components of the Reynolds stress. The results obtained show the effects of particles and mesh size of the turbulence generating grids on turbulence modification. In particular, the presence of solid particles at the initial period of turbulence decay results in the pronounced enhancement of turbulence that diminishes appreciably downwards in the area of typical channel turbulent flow. As the results show, the character of modification of all three normal components of the Reynolds stress taking place at the initial period of turbulence decay are uniform almost all over the channel cross sections. The increase in the grid mesh size slows down the rate of the turbulence enhancement which is caused by particles.


INTRODUCTION
Turbulent gas-solid particle flows in channels have numerous engineering applications ranging from pneumatic conveying systems to coal gasifiers, chemical reactor design, and are one of the most thoroughly investigated subjects in the area of 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, etc.
The mutual effect of particles and a flow turbulence has been the subject of numerous theoretical studies during several decades.These studies have reported about the influence of gas turbulence on particles (one-way coupling) and/or particles on the turbulence of a carrier gas flow (two-way coupling) in case of high flow mass loading (four-way coupling).The influence of particles on gas turbulence consists in turbulence attenuation or augmentation depending on the relation between the parameters of gas and particles.
There are different approaches and numerical models that describe the mutual effect of gas turbulence and particles.The k ε − models, earlier elaborated for turbulent particulate flows, e.g., Elghobashi and Abou-Arab (1983), Pourahmadi and Humphrey (1983), Rizk and Elghobashi (1989), Simonin (1990), Deutsch and Simonin (1991), considered turbulence attenuation only by the additional terms of the equations of the turbulence kinetic energy and its dissipation rate.The results obtained by these models were validated by the experimental data on turbulent particulate free-surface flows (Shraiber et al., 1990).
The models by Crowe and Gillandt (1998) and Crowe (2000) considered both the turbulence augmentation and attenuation in pipe particulate flows depending on the flow mass loading and the Stokes number.Later on these models have been expanded for free-surface flows.As opposed to the k ε − models, Crowe and Gillandt (1998) and Crowe (2000) considered both the turbulence augmentation caused by the velocity slip between gas and particles and the turbulence attenuation due to the change of the turbulence macroscale which occurred in the particulate flow as compared to the unladen flow.The given approach has been successfully tested for various pipe and channel particulate flows.
Currently, the probability dense function (PDF) approach is widely applied to the numerical modelling of particulate flows.The PDF models, for example, by Reeks (1991Reeks ( , 1992)), Zaichik and Vinberg (1991), Zaichik et al. (2007), Kartushinsky et al. (2009a), contain more complete differential transport equations, which are written for various velocity correlations and consider both the turbulence augmentation and attenuation due to particles.
As opposed to pipe flows, rectangular and square channel flows, even in case of unladen flows, are considerably anisotropic with respect to the components of the turbulence energy, which is vividly expressed near the channel walls and corners being notable as for the secondary flows.In addition, the presence of particles aggravates such anisotropy.Such flows are studied by the Reynolds stress turbulence models (RSTM), which are based on the transport equations for all components of the Reynolds stress tensor and the turbulence dissipation rate.The RSTM approach allows us to completely analyse the influence of particles on longitudinal, radial, and azimuthal components of the turbulence kinetic energy, including also possible modifications of the cross-correlation velocity moments.A few studies based on the RSTM approach have shown its good performance and capability for simulation of complicated flows (Gerolymos and Vallet, 2003), as well for turbulent particulate flows (Taulbee et al., 1999).
Recently Mukin and Zaichik (2012) have proposed a 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 the 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.Eventually, the model by Mukin and Zaichik (2012) operated with the k ε − solution and did not allow analysis of the effect of particles on each component of the Reynolds stress.
The presented 3D RSTM model has been applied to the downward grid-generated turbulent particulate flow in the frame of the channel domain.The influence of solid particles on turbulence modification was under consideration for two different values of the length scale of the initial turbulence, which were set by the mesh size of the turbulence generating grid.Particularly, the modification of three normal components of the Reynolds stress which occurred in the axial direction and in the flow cross section was analysed.The advantage of the given model consists in that it considers the particulate phase, following Crowe (2000), taking into account both the enhancement and attenuation of turbulence.
The model, discussed here, has been verified and validated by comparison of the numerical results with the experimental data by Hussainov et al. (2005) for the turbulent downward vertical channel flow gridgenerated by gas-solid particles.

Model description
The sketch of the computational flow domain is shown in Fig. 1.Here u and s u are the velocities of gas and particles, respectively.

Governing equations for the Reynolds stress turbulence model
The numerical simulation of the stationary incompressible 3D turbulent particulate flow in the square cross section channel was performed by the 3D RANS model with applying the 3D Reynolds stress turbulence model to the closure of the governing equations of gas.The particulate phase was modelled in the frame of the 3D Euler approach with the equations closed by the two-way coupling model by Crowe (2000) and the eddy-viscosity concept.
The particles were brought into the developed isotropic turbulent flow set-up in the channel domain which had been preliminarily computed to obtain the flow velocity 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.Therefore, here we present only the system of equations of the gas phase written for the case of the particle-laden flow in the Cartesian coordinates.
The 3D governing equations for the stationary gas phase of the laden flow are written together with the closure equations as follows: 1. continuity equation where , u , v and w are the axial, transverse, and spanwise time-averaged velocity components of the gas phase, respectively; 5. the transport equation of the x-normal component of the Reynolds stress 8. the transport equation of the xy shear stress component of the Reynolds stress ; 9. the transport equation of the xz shear stress component of the Reynolds stress ; 10. the transport equation of the yz shear stress component of the Reynolds stress 11. the transport equation of the dissipation rate of the turbulence kinetic energy The given system of transport equations Eqs (1)-( 11) is based on the model by Launder et al. (1975) with using the numerical constants taken from Pope (2008) 2)-( 7) pertain to the presence of particles in the flow and contain the particle mass concentration .
α The influence of particles on gas is considered by the aerodynamic drag force in the momentum equations (the last term of the right-hand sides of Eqs ( 2)-( 4), and by the turbulence generation and attenuation effects contained in the transport equations of components of the Reynolds stress (the penultimate and last terms of the right-hand sides of Eqs ( 5)-( 7), respectively).The given model employs the two-way coupling approach by Crowe (2000), where the turbulence generation terms are proportional to the squared slip velocity, and the turbulence attenuation terms are expressed via the hybrid length scale h L and the hybrid dissipation rate h ε of the particle-laden flow, where h L is calculated as the harmonic average of the integral length scale of the unladen flow 0 L and the interparticle distance .
u v w The production terms P are determined according to Pope (2008) as follows: The diffusive or second-order partial differentiation over Cartesian coordinates, i.e. the first three terms in Eqs ( 5)-( 11) are given, e.g., by (Pope, 2008).The anisotropy terms R of the normal and shear components of the Reynolds stress 2 , u′ 2 , v′ 2 , w′ , u v ′ ′ , u w ′ ′ v w ′ ′ are defined by various pressure-rateof-strain models of the isotropic turbulence written in terms of variation of constants R C and 2 C (Pope, 2008) as follows: The relative friction coefficient D C′ is expressed as The closure model for the transport equations of the particulate phase was applied to the PDF model by Zaichik and Alipchenkov (2005).There the turbulent kinetic energy of the dispersed phase, the coefficients of turbulent viscosity, and turbulent diffusion of the particulate phase are respectively determined as follows: 0 1 exp , where t ν is turbulent viscosity, 2 0 0 0.09 , is the particle response time with respect to correction of the particle motion to the non-Stokesian regime.

Boundary conditions for the Reynolds stress turbulence model
The flow discussed here is vertical, and symmetrical with respect to the vertical axis.Therefore, the symmetry conditions are set at the flow axis, and the wall conditions are set at the wall.
The axisymmetric conditions are written as follows: The wall conditions are written as follows: for 0.5 : where h is the channel width; * v is the friction velocity of gas; the empirical constant 0.41; ae = the wall coordinates y + and z + correspond to the transverse and spanwise directions, respectively: The friction velocity of gas * v can be determined according to Perić and Scheuerer (1989) as follows: where c µ is the numerical constant of the k ε − model, 0.09.c µ = For the normal and shear stresses and dissipation rate of the unladen flow calculated at the wall, the boundary conditions are set based on the "wall-function" according to Pope (2008) with the following relationships for the production and dissipation terms: for 0.5 : The boundary conditions for the particulate phase are set at the wall as follows: for 0.5 : for 0.5 : At the exit of the channel the following boundary conditions are set: Additionally, the initial boundary conditions are set for three specific cases: 1. the low level of the initial intensity of turbulence that usually occurs at the axis of the channel turbulent flow; 2. the high level of the initial turbulence generated by two different grids: a) small grid of the mesh size 4.8 mm; M = b) large grid with the mesh size of 10 mm.M =

Numerical method
The control volume method was applied to solve the 3D partial differential equations written for the unladen flow (Eqs (1)-( 11)) and the particulate phase (Eqs ( 26)-( 29)), respectively, with taking the boundary conditions (Eqs (30)-( 41)) into account.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 Perić and Scheuerer (1989).This method is utilized for the calculation of particulate turbulent flows in channels of 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 1 120 000.

RESULTS AND DISCUSSION
The presented RSTM model has been verified and validated by comparison of the numerical results with the experimental data obtained by Hussainov et al. (2005) for the grid-generated turbulent downward vertical channel flow of the 200-mm square cross section loaded with 700-µm glass beads of the physical density 3 2500 kg/m .p ρ = The mean flow velocity U was 9.5 m/s, the flow mass loading * m was 0.14 kg dust/kg air.The grids of the square mesh size 4.8 M = and 10 mm were used for generating the initial turbulence length scale of the flow.
The validity criterion is based on the satisfactory agreement of the axial turbulence decay curves occurring behind different grids in the unladen and particle-laden flows obtained by the given RSTM model and by the experiments of Hussainov et al. (2005).Figure 2 demonstrates such agreement for the grid 4.8 M = mm.Figure 3 shows the decay curves calculated by the presented RSTM model for the grids 4.8 M = and 10 mm.As follows from Figs 2 and 3, the pronounced turbulence enhancement by particles is observed for both grids.The character of the turbulence attenuation occurring along the flow axis agrees with the behaviour of the decay curves in the grid-generated turbulent flows described by Hinze (1975).
One can see that the turbulence enhancement occupies over 75% of the half-width of the channel, and it takes place at the initial period of the turbulence decay of the particle-laden flow as compared to the unladen flow.The distributions of ,    3).This means that the grid-generated turbulence of the particulate flow decays downstream, causing a decrease in the rate of turbulence enhancement due to the particles occurring beyond the initial period of turbulence decay.As a result, the turbulence is attenuated, which is expressed in terms of decrease in u ∆ towards the pipe wall (see Fig. 7).Such a tendency has been shown qualitatively by Kartushinsky et al. (2009b).
The defined increase in ,   The analysis of Fig. 11 shows that the increase in the grid mesh size results in a weaker contribution of particles to turbulence enhancement and dissipation of the kinetic energy taking place over the cross section for the initial period of turbulence decay.This can be explained by a higher rate of particle involvement into the turbulent motion due to the longer residence time that comes from the larger size of the eddies.

CONCLUSIONS
The 3D Reynolds stress turbulence model (RSTM) based on the 3D RANS numerical approach has been elaborated for the grid-generated turbulence in the particulate downward channel flow domain.
The influence of the particulate phase on turbulence modification, namely the modification of three normal components of the Reynolds stress that occurred in the axial direction and in the flow cross section, was examined for two different values of the length scale of the initial turbulence, which was set by the mesh size of the turbulence generating grid.
The obtained numerical results allow us to draw the following conclusions: 1.The presence of the considered solid particles at the initial period of turbulence decay results in a pronounced turbulence enhancement for both turbulence generating grids.2. The character of modification of the normal components of the Reynolds stress which has taken place at the initial period of turbulence decay is uniform almost all over the channel cross section.As the turbulence level diminishes downstream, the extent of this uniformity (plateau width) and the effect of particles on turbulence decrease.3. The rate of turbulence enhancement beyond the initial period of turbulence decay reduces.The character of turbulence modification changes from initially uniform distributions of the normal components of the Reynolds stress to those inherent in the turbulent channel flow.4. The increase in the grid mesh size slows down the rate of turbulence enhancement caused by particles.
The model presented considers both the enhancement and attenuation of turbulence by means of the additional terms of the transport equations of the normal Reynolds stress components.It allows us to make calculations covering the long distance of the channel length without using algebraic assumptions for various components of the Reynolds stress.
With applying a minimum number of assumptions and empiricism, the model represents a more contemporary computational approach in a turbulent particulate flow.It is also simpler and uses the stateof-the-art modelling and computational techniques and is more accurate because no approximations are employed.
integral time scales for unladen and particle-laden flows, respectively; of gas in particle-laden and unladen flows, respectively; ε and 0 ε are the dissipation rates of the turbulence kinetic energy in particle-laden and unladen flows, respectively; of slip velocity.The additional terms of Eqs (

∆
are the width and height of the control volume; the numerical constant B equals 5.2 for the smooth wall of the channel.

∆
Figures 4-9show the cross section modifications of three components of the Reynolds stress, , caused by 700-µm glass beads, calculated by the presented RSTM model at two locations of the initial period of the grid-generated turbulence decay 46 x M = and 93 as well as beyond it for 200.
∆ are uniform, which corresponds to the initial gridgenerated homogeneous isotropic turbulence that decays downstream (see Figs 2 and 3).The distributions of modification of , ∆ remain uniform downstream.At the same time, the cross section extent of the uniformity of distributions of components of the Reynolds stress and the effect of particles on turbulence decrease, since the turbulence level decreases downstream (cf.data presented for 46 x M = and 93 in Figs 4-9).The distributions of modification of , u ∆ , v ∆ and w ∆ that have taken place beyond the initial period of turbulence decay (location 200 x M ≈ in Figs 4-9) are typical of the channel turbulent particulate flow.One can see that in this case turbulence enhancement becomes slower, since here the turbulence level is

Fig. 2 .
Fig. 2. Axial turbulence decay behind the grid 4.8 M = mm: 1 and 3 are the data by Hussainov et al. (2005) got for the unladen and particle-laden flows, respectively; 2 and 4 are the numerical data obtained for the same conditions.

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 3.The calculated axial turbulence decay behind the grids: 1 and 2 are the data got for the unladen and particleladen flows, respectively, at 4.8 M = mm; 3 and 4 are the data obtained for the same conditions at 10 M = mm.x/Mx/M w ∆ that is observed towards the wall (Figs 4-9) arises from the growth of the slip velocity (see curves 1, 2, and 3 in Fig.10).The decrease in ,u ∆ , v∆ and w ∆ in the immediate vicinity of the wall is caused by the decrease in the length scale of the energy-containing vortices and thus by the increase in the dissipation of the turbulence kinetic energy.

Fig. 10 .α
Fig. 10.The cross section distributions of the axial gas and particles velocities and particles mass concentration for the grid 4.8 M = mm: 1, u U for 46; x M = 2, , u U 3, , s u U and 4,