Eulerian-Eulerian modelling of particle-laden two

Powder Technology 301 (2016) 999–1007

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Eulerian-Eulerian modelling of particle-laden two-phase flow A. Kartushinsky a, S. Tisler a, J.L.G. Oliveira b,⁎, C.W.M. van der Geld c a b c

Faculty of Civil Engineering, Department of Mechanics, Tallinn University of Technology, Ehitajate tee 5, 19086 Tallinn, Estonia Technology Center, Federal University of Santa Catarina, Joinville, SC 89218-000, Brazil Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, Eindhoven, The Netherlands

a r t i c l e

i n f o

Article history: Received 3 October 2015 Received in revised form 7 July 2016 Accepted 22 July 2016 Available online 25 July 2016 Keywords: Eulerian model Particle-laden Turbulent flow Concentration profile

a b s t r a c t When typical fluid RMS velocities scale with the terminal velocities of inertial particles, turbulent particle-laden pipe flows exhibit special features. Recent experiments have shown that the taking of Eulerian averages requires special measures in such flows. An existing Eulerian numerical method has been adapted in order to accommodate for these special features. The effects of flow direction (upward or downward) and mean concentration (in the range 5.0 × 10−6 to 3.2 × 10−5) on radial particle distribution and on the mean axial velocities of both particles and fluid are investigated. Inertial particles have a Stokes number equal to 2.3, based on the particle relaxation time and viscous scales; Reynolds number is 10,300, based on the bulk velocity and the pipe diameter. To get agreement between experimental findings and model predictions, a lift force term had to be included; core peaking in up-flow and wall peaking in down-flow result from it. The difference of mean particle and fluid Eulerian velocities is found to decrease towards zero near the wall in both up-flow and down-flow. This feature required inclusion of a modified drag component in the particle axial momentum equation of the Eulerian model. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Particle-laden flows are found in many engineering fields and have therefore been investigated for several decades. The dispersion of pollutants in an urban environment, sediment transport or the fluidized catalytic cracking of carbohydrates are often studied. Predicting the behavior of laden-flows is therefore of quite some interest in engineering applications as found in blast furnaces, see Wright et al. [1] and Zhou et al. [2]. An abundant number of experimental and numerical investigations of flows loaded with particles were performed; some of them are reviewed below. Direct numerical simulation (DNS) of turbulent pipe flow with high concentration of large-size inertial particles is practically impossible. Limitations in computational speed and memory capacity cannot be overcome with DNS, at present. Prediction of particle laden flows in these conditions is possible with Eulerian-Eulerian modelling if the modelling coefficients used can be tuned to experimental results in a wide range of process conditions. This approach was followed for various turbulent two-fluid models, using experimental data of gas-particle flows or bubbly flows; see for example Peirano et al. [3], Pialat et al. [4], Konan et al. [5], Ozel et al. [6], Cheung et al. [7], Zeng and Zhou [8] and Zhou [9]. To the best of our knowledge, no such modelling attempt was performed with liquid flows with finite-sized neutrally-buoyant solid particles.

⁎ Corresponding author. E-mail address: [email protected] (J.L.G. Oliveira).

http://dx.doi.org/10.1016/j.powtec.2016.07.053 0032-5910/© 2016 Elsevier B.V. All rights reserved.

In the present study, new experimental data for both up-flow and down-flow of inertial particles in a pipe with an inner diameter of 100 mm are used to validate the Eulerian-Eulerian model developed by Kartusinsky et al. [10]. It turns out that when the ratio of a characteristic RMS fluid velocity, urms, to the terminal velocity of inertial particles, UTV, is of order 1, changes in this model are necessary to get agreement between measured and predicted mean Eulerian axial relative velocities. There is clearly a research gap in the existing Eulerian modelling; and the reason for this gap is explained as follows. There is lack of knowledge on how to predict ensemble averages of fluid properties at or near a particle location in a strongly heterogeneous turbulent flow field. Common expressions for the drag force on an inertia particle, for example, are in terms of Lagrangian relative velocities. The fluid velocity seen by an inertial particle, u(xs)|L, might on averaging be different from the Eulerian mean fluid velocity at a certain radial position in the pipe, 〈u〉; the Lagrangian average relative velocity, 〈us − u(xs)|L〉, is not equal to the difference of Eulerian averages, 〈us〉 – 〈u〉. An inertia particle might preferentially be located in fluid patches moving towards the wall, for example; “crossing trajectories” is another way of describing this difference. For this reason, some kind of weighing of fluid velocities is necessary in order to derive a proper expression for the drag force exerted on inertia particles in a Eulerian model. It is not known a priori what kind of weighing is required and this causes the research gap. The motivation for using both upflow and downflow is the striking difference in particle concentrations of these flows. The more different the trends to be modeled, the better the validation of the model is. Because of the different results of up-flow and down-flow, the only term in the Eulerian-Eulerian model that may accommodate for both trends

1000

A. Kartushinsky et al. / Powder Technology 301 (2016) 999–1007

at the same time is the lift force. An adaptation of the Saffman lift force to non-Stokesian regimes will be proposed and fitted to the measurement data. The outline of the paper is as follows. A detailed description of the Eulerian method is presented in Section 2, including particle and fluid equations and boundary conditions. Computational method is given in Section 3, while the reference experimental conditions are shown in Section 4. Results are presented and analyzed in Section 5 and concluding remarks are finally given in Section 6. 1.1. Experimental work on particle-laden flows Most of the experimental works in particle-laden flows deal with two types of flow: one in which turbulence is insignificant, characterized by | urms/UTV | ≪ 1, and one in which particles act as flow tracers (|urms/UTV | ≫ 1). In the first case, little interaction of particles with turbulent fluctuations occurs and in the second case particle turbulence is the same as fluid turbulence. The last case is the basis of PIV measurements techniques and is not considered any further here. In both cases, knowledge of Eulerian averaged quantities is sufficient to provide reasonable fluid and particle statistics. However, this does not hold for particles or bubbles which possess | urms/UTV | ≈ O (1). This case will be considered in detail in the present study. Most of the experimental studies of particle-laden flows concerned air with solid particles, small as compared to the Kolmogorov length scale and with a high mass density. These are turbulent wall-bounded flows with | urms/UTV | ≪ 1. Good examples are Tsuji & Morikawa [11], Tsuji et al. [12], Kulick et al. [13], Kussin & Sommerfeld [14], Caraman et al. [15], Yang & Shy [16], Benson et al. [17] and Borée and Caraman [18]. Experimental data for inertial particles with a size exceeding the Kolmogorov length scale are scarce. Turbulent transport of near-neutrally buoyant particles in water was only recently experimentally examined, see Volk et al. [19], but only for homogenous and isotropic turbulence. Only very few particle laden liquid flow experiments were performed in turbulent flow in a channel; see Sato & Hishida [20] and Suzuki et al. [21]. In this cases |urms/UTV | ≪ 1 and measurements are in a square channel, with only a restricted class of inertial particles. Turbulent liquid flows in pipes with almost neutrally buoyant solid particles were to the best of our knowledge not performed, and the experiment of Oliveira [22] is therefore of particular interest. 1.2. Numerical work on particle-laden flow Also most of the numerical works in particle-laden flows can be characterized by |urms/UTV | ≪ 1. Models to predict turbulence augmentation or attenuation by particles in two-phase flows are usually simplified, see Elgobashi [23]. Numerical work of particle-laden flows is mostly limited to point-force predictions; see Poelma et al. [24]. Fully numerical computations of finite particles require significant memory and processing speed and only simple geometries with limited numbers of particles were studied, see Tryggvason et al. [25] and Dijkhuizen et al. [26]. Numerical works in particle-laden flows were therefore focused on small solid particles, compared to the Kolmogorov length scale. The flow satisfies the so-called two-way coupling approximation, see Ferrante & Elghobashi [27], Balachandar & Eaton [28] and Soldati & Marchioli [29]. Simulations involving finite-size particles are not usual. Pan & Banerjee [30] investigated the effect of finite-size particles in a turbulent channel flow. They showed that turbulent stresses become larger when particle size is larger than the dissipative scale. Some studies with |urms/UTV | ≈ O (1) have recently become available with the aid of direct numerical simulations (DNS) analysis. Investigation of heavy finite-size particles in open-channel has been performed in the dilute regime by Kidanemariam et al. [31] and Kidanemariam & Uhlmann [32]. Even more recently, Picano et al. [33]

modeled turbulent channel flow with dense suspensions of neutrally buoyant spheres. To the best of our knowledge, no Eulerian modelling attempts were made for |urms/UTV | ≈ O (1). 1.3. Objectives The aim of the present work is to develop and validate a numerical method able of reproducing turbulent particle-laden pipe flow statistics. Particularly, the model must deal with laden turbulent flows with the ratio of the RMS fluid velocity to the terminal velocity of inertial particles of order 1. The prediction is focused on the radial distribution of inertial particles in both upflow and downflow, and on the mean axial velocity of particles and fluid. The experimental data of Oliveira [22] will assist the method validation. 2. Numerical method 2.1. Description of CFD model The present CFD model utilises the Reynolds-Averaged NavierStokes (RANS) equations, along with the eddy viscosity hypothesis, to model both the liquid and particulate phases. For the liquid phase, the equations are closed using the standard two-equation k–ε model, corrected for coupling between the two phases, i.e. two-way coupling. In the present model cylindrical coordinates are utilized, where u and v are the mean velocity in the axial (x) and radial (r) directions, respectively. Subscripts p and s refer to the conditions of a single particle and particulate phase, respectively. The model assumes the flow through the pipe to be axisymmetric, steady and incompressible. For the particulate phase, a statistical PDF model by Zaichik et al. [34] is used to close the equations; see the appendix. 2.2. Transport equations of the continuous phase The first transport equation used in the current CFD model is the continuity equation of the liquid: ∂rρu ∂rρv þ ¼0 ∂x ∂r

ð1Þ

where ρ is the mass density of the liquid. The momentum equation of the liquid phase in the axial direction is as follows: ∂rρu2 ∂rρuv ∂ ∂u ∂ ∂u ∂v ∂rp − ¼ 2r ðμ t þ μ Þ þ r ðμ t þ μ Þ þ − þ ∂r ∂x ∂x ∂r ∂r ∂x ∂x ∂x2 sffiffiffiffiffiffi3 ! ! ð2Þ u−us ρ ðv−vs Þ ∂u5 −g 1− þ f mei ραr4 0 τp ρp τsaf ∂r where p is pressure, μ is the dynamic laminar viscosity of the liquid, and μt is the dynamic turbulent viscosity. α is the particle mass concentration calculated as: α = β (ρp/ρ); where β is the particle volume fraction, ρp is the mass density of solid particles and g is the acceleration of gravity. The variable us stands for the axial velocity of particulate phase. Let δ be the particle diameter and the particle Reynolds number, Res, be defined by:

Res ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρδ ðu−us Þ2 þ ðv−vs Þ2 μ

ð3Þ

where vs represents the radial velocity of the particulate phase. The parameter τ'p is the particle response time that specifies the drag. τ'p is defined as: τ'p = τp/CD; with τp, the particle relaxation time represented by Eq. (4), Albrecht et al. [35], and CD, the drag coefficient for the non-Stokesian regime as given by Schiller and Naumann [36], see Eq. (5), where f (δ/2h) represents a correction term for the presence


of the wall. h

τp ¼

i ρp δ2 1 þ 0:5ρ=ρp

ð4Þ

18μ

CD ¼

24 Res

1 2=3 1 þ Res f ðδ=2 hÞ; Res b1000 6

ð5Þ

1001

1.44 and cε2 = 1.92. The expression for the production of turbulent energy is provided by (13). The dynamic turbulent viscosity, μt, is computed by (14) with the constant cμ = 0.09. " 2 # 2 v2 ∂u ∂v2 ∂u ∂v Pk ¼ μ t 2 þ2 þ2 þ þ r ∂x ∂r ∂r ∂x

ð13Þ

2

In Eq. (5), the variable h stands for the distance of the particle center to the pipe wall. Faxén [37] derived a wall correction f for creeping flow in the following form: " f ðδ=2 hÞ ¼

1−C 1

3

9 δ 1 δ þ 16 2h 8 2h

−

4

45 δ 256 2h

−

5 #!−1

1 δ 16 2h

ð6Þ with C1 equal to 1. The parameter C1 was introduced in (6) in order to be able to vary the effect of the wall in subsequent computations. The bigger C1, the bigger the influence of the wall is. The momentum equation of the liquid phase in the transverse direction is: ∂rρuv ∂rρv2 ∂ ∂v ∂u ∂ ∂v ∂p þ 2r ðμ t þ μ Þ −r − þ þ ¼ r ðμ t þ μ Þ ∂x2 ∂x ∂x ∂r ∂r ∂r ∂r ∂r sffiffiffiffiffiffi3 ! ð7Þ v−v ð u−u Þ ∂u s s 5 −f ραr 4 mei τ0p τsaf ∂r The expression in brackets represents the effects of drag and lift owing to the presence of the particulate on the momentum equation of the liquid phase. The variable τSaf represents the time scale of the Saffman [38] lift force, Eq. (8). The term fmei allows the application of the Saffman lift force to non-Stokesian regimes, Mei [39]; see Eq. (9). Alternatively, the lift force as derived by Auton [40] could replace the modified Saffman lift expression in Eq. (7).

τSaf

h i ρp δ 1 þ 0:5ρ=ρp ¼ pffiffiffiffiffiffi ρμ

f mei ¼

ð8Þ

( 1−0:3314ϕ0:5 expð−0:1 Res Þ þ 0:3314ϕ0:5 0:0524ðϕRes Þ0:5

Res ≤40 Res N40

ð14Þ

2.3. Boundary conditions As the flow is assumed to be asymmetrical, the following boundary conditions are employed at the pipe axis (r = 0): ∂u ∂k ∂ε ∂us ∂α ∂ks ¼ ¼ ¼ ¼ ¼ v ¼ vs ¼ ¼0 ∂r ∂r ∂r ∂r ∂r ∂r

ð15Þ

For the liquid, based on the control volume method - see Fertziger and Perić [42] and Perić and Scheuerer [43], the liquid shear stress at the pipe wall is: 8 > >
> ≅ : ρν t ∂r ln ðΣ yþ Þ ρν

ð16Þ

where y+ is the distance from the pipe boundary in wall units, Δr is the computational grid spacing, Ψ = 0.41 is the von Karman constant, Σ = eΨB = 8.4317 (B is an empirical constant related to the thickness of visk0.5, is calculated cous sublayer, B ~ 5.2). The friction velocity, v+ = c0.25 μ using the wall function method, with the assumption of local equilibrium between production and dissipation of the turbulent energy. At the pipe boundary, the production of the turbulent energy is:

ð9Þ

ð10Þ

k ε

The effect of solid particles in k–ε model is given by the last term on the RHS of Eqs. (11) and (12), see Schwarzkopf et al. [41]. The transport equations of the dispersed phase are presented in the appendix.

Pk ¼ τw

where the parameter ϕ is given by Eq. (10). ϕ varies in the range: 0.005 b ϕ b 0.4. ∂u δ ϕ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ∂r 2 ðu−us Þ þ ðv−vs Þ

μ t ¼ ρcμ

pffiffiffi ∂u k=ðψΔr Þ ≅2τw c0:25 μ ∂r

ð17Þ

In Eq. (17), τw is the quantitative parameter implemented into the wall function definition and determined by the flow regime in the vicinity of a wall, either laminar or turbulent regime. At the pipe boundary, the dissipation rate of the turbulence energy is: k ε ¼ 2c0:75 μ

1:5

=ðψΔrÞ

ð18Þ

Equations for the turbulent kinetic energy, k, and for the dissipation of the turbulence energy of the liquid phase, ε, are presented in Eq. (11) and (12).

For the particulate phase at the pipe boundary, the following boundary conditions are employed (Zaichik et al. [44]; Zaichik, [45]; Zaichik et al. [34]):

∂rρuk ∂rρvk ∂ μ ∂k ∂ μ ∂k þ ¼ r t þμ þ r t þμ þ Pk ∂x ∂r ∂x σ"k ∂x ∂r σ k# ∂r 2 2 ðu−us Þ þ ðv−vs Þ þ ραr τ0p

Ds

ð11Þ

∂rρuε ∂rρvε ∂ μ ∂ε ∂ μ ∂ε ε þ ¼ r t þμ þ r t þμ þ ðcε1 P k −cε2 ε Þþ ∂x ∂r ∂x " σ ε ∂x ∂r #σ ε ∂r k 2 2 1:416 μ Res ðu−us Þ þ ðv−vs Þ 0:058cε2 αr τ0p 3πδ2 ð12Þ where the k–ε model constants are as follows: σk = 1, σε = 1.3, cε1 =

∂α ¼ 0; ∂r

τ 0p

vs ¼ 0;

!0:5 ∂us 1−χex 1−χ 2 − us ð19Þ ¼2 1 þ χex 1 þ χ π‹v0 2s › ∂r

where χ is the reflection coefficient, which is the probability of particles recoiling off the boundaries and back to the flow. In the current model χ = 1/3, see Binder & Hanratty [46]. ex is the coefficient of restitution in the axial direction, modeled as:

ex ¼

8 > < 1−ξx > :

5 7

0≤ξx ≤ ξx N

2 7

2 7

ð20Þ

1002


This coefficient ex depends on the attack angle of a particle with a wall. Here, the parameter ξx = ηx (1 + ey) tan(φx), where ey is the coefficient of restitution normal to the wall; ηx is the coefficient of friction between the particles and the wall; φx = tan− 1 (vs/us) is the attack angle between the trajectory of a particle and the wall. The coefficient of restitution reflects the loss of a particle momentum as a particle hits the wall. Here, ey = 1 and ηx = 0.39. At the pipe inlet (x = 0), the following conditions are imposed:

The fluid time-scale τf in the Stokes number, St, is based on viscous scales; see info below Table 1. The evaluation of the particles timescale, τp, is given by Eq. (4). A relaxation time of τp ≈ 4 ms is obtained for the tracers. The terminal velocity specified in Table 1 is attained in quiescent fluid when gravitational and drag forces are in equilibrium:

us ð0; r Þ ¼ kslip uð0; r Þ;

where CD is the drag coefficient, which is a function of the particle Reynolds number, Res. In Table 1, Res is based on the particle diameter and the terminal velocity; Res = δ|UTV |/ν. In the Stokes regime, CD is given by Eq. (24). The drag coefficient for the non-Stokesian regime is calculated by Eq. (5).

vs ð0; r Þ ¼ kslip vð0; r Þ;

α ð0; r Þ ¼ α 0

ð21Þ

where kslip is the initial slip coefficient, and α0 is the initial flow mass loading (set to match the experimental conditions). Finally, at both the pipe inlet and outlet, the following boundary conditions are applied: ∂u ∂k ∂ε ∂us ∂vs ∂α ∂ks ¼ ¼ ¼ ¼ ¼ ¼v¼ ¼ 0: ∂x ∂x ∂x ∂x ∂x ∂x ∂x

ð22Þ

3. Computational method The numerical procedure was performed in two steps. In the first step, the single-phase flow field (subscript “0”) was obtained numerically by solving the equations of the liquid in the absence of the particles. In the second step, the two-phase flow was computed using the results obtained for the single-phase flow field as the initial conditions. The inlet boundary conditions for the particulate phase were set by taking into account some velocity slip, that was carried out by introducing the initial slip coefficient kslip, which is usually kslip ~ 0.8–1. The control volume method was applied to solve the 2D partial differential equations written for both the liquid (Eqs. 1, 2, 6, 10 and 11) and the particulate phase (Eqs. A.1 – A.4), respectively, along with the boundary conditions listed above. These governing equations were solved using the implicit lower and upper (ILU) matrix decomposition method with the flux-blending differed-correction and upwinddifferencing schemes by Perić and Scheuerer [43]. This method was previously utilized in the calculations of the particulate turbulent flows in pipes and channels of the round, square or rectangular cross-sections. The calculations were performed in the dimensional form for all the flow conditions. The number of the control volumes was 1,154,304, which encompassed 150 pipe diameters. 4. Experimental data for numerical code validation Particle and fluid statistics of particle-laden pipe flows obtained by the present numerical method are compared with the experimental datasets of Oliveira [22]. Vertical flows were investigated by 3D particle tracking velocimetry (3D-PTV) in a dedicated test rig in Eindhoven University of Technology, The Netherlands. All tests were performed at the same bulk flow Reynolds number (Reb), 10,300, based on the pipe diameter and bulk velocity. Information about the experimental setup, the particle tracking algorithm used to obtain trajectories of tracers and inertial particles and the trajectory analysis method of individual particle trajectories is given by Oliveira et al. [47] and Oliveira et al. [48]. Properties of polystyrene particles applied in the particle-laden experiments are summarized in Table 1.

UTV ¼ f4 ðρp−ρf Þ δ g=3 CD ρ f g

CD ¼

0:5

ð23Þ

24 ; Res b1 Res

ð24Þ

The value for UTV was obtained by an iterative computation using Eq. (23) and Eqs. (24) or (5). The Kolmogorov scale for fully developed single-phase pipe flow at Reb = 10,300 was computed from DNS results by Veenman [50]. The Kolmogorov length is about 0.6 mm in the pipe core and 0.2 mm at the wall region. The bulk flow velocity, Ub, was approximately 100 mm/s, and the ratio Ub/UTV was of the order of 102 for seeding particles, see Table 1. Since Ub ≫ UTV, τp b τf and dp b η, the employed seeding particles worked well as flow tracers. For inertial particles, the ratio Ub/UTV is of the order of 10, τp N τf and dp N η for the entire cross-section. Therefore, inertial particles have significant inertial characteristics, as was also exhibited by the differences in the Lagrangian statistical properties, auto- and cross-correlations of velocities and accelerations for tracers and inertial particles; see Oliveira et al. [48]. 4.1. Experimental conditions Downward and upward vertical flows were measured at Reb = 10,300 by Oliveira [22]. The bulk velocity of each flow, Ub, was adapted to temperature changes to keep Reb ≈ 10,300. Particle-laden flows with mean concentrations of inertial particles, ‹β›, ranging from 5.0 × 10−6 to 3.2 × 10−5, have been tested. Properties of inertial particles were selected with the aim of testing particle-laden flows which have a characteristic root-mean-square velocity representative of the turbulent carrier phase, urms, and the terminal velocity of the dispersed phase, UTV, of same order of magnitude: urms/UTV ≈ O (1). Further explanations will be provided in Section 6. A particle-laden experimental case is represented here by a number (1 or 2), which indicates the flow direction and a letter (A, B or C), indicating the applied mean concentration of inertial particles. A summary of all tested particle-laden cases is shown in Table 2. 5. Results 5.1. Concentration profiles of inertial particles In Fig. 1 the volumetric concentration profiles of inertial particles, β, obtained by the present numerical method (solid lines) are compared with experimental results (squares) obtained by Oliveira [22]. The

Table 1 Properties of inertial particles applied in the particle-laden experiment of Oliveira [22]. Particles

Mass density [kg/m3]

Diameter dp [mm]

Terminal Velocity, UTVa [mm/s]

Resa

St = τp/τfb

Flow tracers Inertial particles

1050 1050

0.2 0.8

1.0 10.2

0.18 7.76

0.14 2.31

a

Settling velocity of a particle in an infinite, stagnant pool of water. Fluid time-scale is based on viscous scales as given by: τf = ν/u2τ. For Reb b 105, the wall shear velocity can be estimated as uτ = (U2bf/8)1/2 with f = a Re-m b , m = 0.25 and a = 0.316; see Hinze [49]. τf is roughly 28 ms. b

A. Kartushinsky et al. / Powder Technology 301 (2016) 999–1007 Table 2 Summary of particle-laden flow experiments of Oliveira [22]. All experimental runs have been performed at Reb = 10,300. Case

Flow Direction

〈β〉 × 10−5

1A 1B 1C 2A 2B

Upward Upward Upward Downward Downward

0.5 1.4 3.2 1.8 2.8

volume-averaged concentration, 〈β〉, ranged from 0.5 × 10− 5 to 3.2 × 10−5; see Table 2. The effect of the flow orientation with respect to gravity on the volumetric concentration profiles of inertial particles is clear in Fig. 1. Close to the pipe walls in upward flows, for r/R N 0.85, a reduction in β occurs; see cases 1 A/1B/1C in Fig. 1a, b and c, respectively. For downward flows, there is a peak at r/R ≈ 0.97; see cases 2 A/2B in Fig. 1d and e. Therefore, the direction of the vertical flow, upward or downward, strongly affects the concentration profiles in the near-wall zone. The numerical method is able to reproduce, to a certain extent, the wall-peaking tendency of downflows and the core-peaking trend of upflows; see solid lines in Fig. 1. The agreement is generally good; only for cases 2 A and 2B the concentration in the range 0.84 b r/ R b 0.96 is somewhat overpredicted. Only with a lift force of a kind the different trends for up- and downflow can be reproduced. Here, the adaptation of the Saffman lift force to non-Stokesian regimes, as provided by Mei [39], was adopted. The lift force depends crucially on the gradient of mean axial velocity in radial direction, ∂ u/∂r; see the last term on the RHS of Eqs. (2) and (A.3). It is concluded that the lift force is essential for the prediction of particulate two-phase flows.

1003

a

b

c

5.2. Mean axial velocity profiles The effects of flow direction and mean concentration on the mean axial velocity profiles of fluid, u, and that of inertial particles, us, are presented in Fig. 2. The bulk flow velocity, ub, was adjusted to keep the Reynolds number, Reb, for each experiment the same, 10,300. Inertial particle and fluid velocity profiles are normalized in these figures by the corresponding ub. In Fig. 2a–e, the fluid velocity profiles obtained by the numerical method are represented by dashed lines and the particle velocity profiles by solid lines. Experimental results attained by Oliveira [22] are represented by symbols: circles for fluid velocity profiles and squares for particle mean axial velocities. In Fig. 2a–c, the mean axial velocity profiles of fluid and inertial particles for cases 1A/1B/1C are presented. No significant changes in the experimental axial velocity profiles are noticed when the mean volumetric concentration, 〈β〉, is varied from 0.5 × 10−5 to 3.2 × 10−5 in upward flows. Numerical results follow the experimental results closely, except for the velocity profiles of particles close to the wall, i.e. for r/R N 0.8. The numerical results nevertheless reproduce the reduction in the mean relative axial velocity, u–us, in the near-wall region; however, this reduction occurs in the experiments at r/R values of about 0.8 already. The profiles of fluid and inertial particles, u and us, for downflow cases with 〈β〉 = 1.8 × 10−5 and 2.8 × 10−5 (cases 2A/2B) are presented in Fig. 2d and e, respectively. Of course, fluid particles now lag behind inertial particles. Similar results to cases 1 A/1B/1C regarding the reduction in the mean relative axial velocity were obtained. Once again, the reduction in u–us occurs at r/R values of about 0.8 in the experiments. These trends are further examined below. In order to elucidate the trends observed in Fig. 2, the profiles of u–us as a function of the radial coordinate are presented in Fig. 3. Numerical results for particle-laden flows summarized in Table 2 coincide and therefore only one curve is given in Fig. 3. This is probably related to the low concentration of inertial particles applied in the simulations. The numerical results with C1 = 1 (Faxén correction) were found to underestimate the effect of the wall; see solid lines in Fig. 3. Only very

d

e

Fig. 1. Effect of (upward or downward) flow direction on the volumetric concentration profiles of inertial particles, β (r/R). Fig. 1a, b and c represent upward cases and Fig. 1d and e, downward cases.

1004


a

b

c

close to the wall the relative velocity drops down. Fig. 3 also shows that the drag force coefficient at the pipe center is overestimated by about 10%. With a value of 10 for C1, better agreement between experiments and predications is found, see the dashed lines in Fig. 3. This shows that the wall effect in turbulent pipe flow cannot be estimated by the Faxén correction term. The reason for this is a complex interaction between eddies and particle motion in the entire region r/R N 0.7. The ensemble average of the Lagrangian drag coefficient of Schiller & Naumann that must be taken in the Eulerian model can only be taken with sufficient predictive capacity of crossing trajectories and preferential locations of inertia particles with respect to the fluid. The normal averaging does not apply. A drag force solely based on the difference of these two velocities would be zero and could not compensate the average body force on the particle, which is nonzero. Another way of averaging is therefore required. The other way of averaging is accounted for by the tuning coefficient C1. Several studies revealed that the rise velocity of particles is reduced by turbulence due to the increased residence time of particles in the upward side of large vortices, see Aliseda and Lasheras [51], for example. In the Eulerian model, the decrease of the axial component of the velocity slip is caused by three reasons: 1) addition of the lift force on the RHS of Eqs. (2) and (A.3) (the last term); 2) direct contact of the dispersed phase with the wall; namely, the effect of the friction coefficient ηx, Eq. (20); 3) increase of the drag, in the example of Eq. (5) by the increase of coefficient C1. The sample computations of Fig. 3 show the necessity of a modified drag in Eulerian modelling of particle-laden pipe flow. Alternatively, a modified drag or lift component relying on particle inertia and turbulence features, e.g. urms/UTV, could be used as well. The tuning coefficient is expected to approach 1 if the ratio |urms/UTV | is far bigger or much b 1 and to go to the value found in our study for values of this ratio around 1. 5.3. Second order velocity statistics

d

e

Numerical and experimental RMS results of the turbulent kinetic energy, k1/2, for the liquid and the dispersed phase are compared in Fig. 4. Results are scaled with the bulk flow velocity. In all experimental test conditions, the loading of particles is too small to affect turbulence statistics really. For this reason, only one specific case is presented. The largest turbulent kinetic energy is observed at r/R ~ 0.96 and the trends of simulations and measurements are the same. Neutrally-buoyant particles, indicated by triangles in Fig. 4, and fluid particles possess similar values of velocity fluctuations. Differences in turbulent behavior occur near the wall, where the particle diameter exceeds the Kolmogorov length scale: dp/η is roughly 1.3 in the pipe centerline and increases to 4 in the wall region. The agreement between numerical and experimental data is mainly satisfactory in the range r/R N 0.5 (or about 75% of the cross-section area). Although differences occur elsewhere, the numerical model is able to reproduce the approximately homogeneous turbulence behavior at pipe core, as well as the turbulence inhomogeneity at the wall region. 6. Conclusions A Eulerian numerical method has been developed to simulate turbulent particle-laden pipe flows for flow conditions when typical fluid Fig. 2. Effect of upward or downward flow direction and mean particle concentration on the mean velocity profiles of fluid, u, and inertial particles, us. The velocities are normalized by the bulk velocity of each flow, ub. Fig. 2a, b and c represent upward cases and Fig. 2d and e, downward cases.


a

Appendix A. Transport equations of the dispersed phase

(u - us ) / | ub |

0.15

Exp Num Num*

0.05

sffiffiffiffiffiffi3 ! ! u−u ρ ð v−v Þ ∂u5 s s ραr 4 −g 1− þ f mei τ 0p ρp τ saf ∂r

ðA:1Þ

2

upward case 1B 0.5 r/R

1

b 0.15

Exp Num Num*

0.1

0.05 downward case 2A 0 0

The particulate phase mass conservation, axial and radial momentum equations are presented in Eqs. (A.1), (A.2) and (A.3), respectively. ∂rραus ∂rραvs ∂ ∂α ∂ ∂α þ rρDs þ ¼ rρDs ∂x ∂x ∂r ∂r ∂x ∂r ∂rραu2s ∂rραus vs ∂ ∂us ∂ ∂us ∂vs þ þ ¼ 2rμ s α þ rμ s α þ ∂x ∂x ∂r ∂x ∂r ∂r ∂x

0.1

0 0

(u - us ) / | ub |

1005

0.5 r/R

1

Fig. 3. Mean relative axial velocity, u–us, as a function of the radial coordinate, for upflows (a) and downflows (b). The solid line (Num) stands for numerical results with C1 = 1 (Faxén correction), whereas the dashed line (Num*) represents numerical results with C1 = 10.

RMS velocities scale with the terminal velocities of inertial particles, i.e. when |urms/UTV | ≈ O (1), and for more usual conditions when this ratio is much b1. The radial distributions of inertial particles with Stokes number equal to 2.3 (based on the ratio of particle relaxation time to viscous scales) in up-flow and down-flow as well as the mean axial velocities of both particles and liquid could be closely reproduced. Comparison was made with new experimental data at bulk Reynolds number of 10,300. The main trends of inertial particle concentrations were captured by the present code. The core peaking in upflow and the wall peaking in downflow were reproduced with the application of the Saffman lift force, following Mei [39]. The actual Eulerian relative velocity in entire region for r/R N 0.7 is decreased by the interaction of large eddies and motion of inertial particles. This paper has shown a way to account for this interaction in Eulerian modelling of this type of two-phase flows.

∂rραus vs ∂rραv2s ∂ ∂vs ∂us ∂ ∂vs þ 2rμ s þ ¼ rμ s þ þ ∂x ∂r ∂x ∂r ∂x ∂r ∂r 2 sffiffiffiffiffiffi3 ! v−vs ðu−us Þ ∂u5 −f fmei ραr4 τ0p τ saf ∂r

ðA:2Þ

ðA:3Þ

where Ds is the coefficient of turbulent diffusion of particles and μs is the turbulent viscosity of the particulate phase. The transport equation for the particle turbulence kinetic energy, ks, is shown in (A.4) for two-dimensional non-stationary flow conditions. ∂ks ∂ks 5 ∂ ∂ks ∂ ∂ks ¼ þ − αDs rαDs ρα us þ vs 9 ∂x ∂x r∂r ∂r ∂x ∂r 2 ∂us ∂vs vs ∂us ∂vs vs þ ραks þ μ s α þ þ þ þ r r 3 " ∂x ∂r ∂x ∂r # 2 2 2 v 2 ∂us ∂vs ∂us ∂vs 2ð f u k−ks Þ s þ þ2 þ2 þ μsα 2 þ r τ0p ∂x ∂r ∂r ∂x ðA:4Þ Eq. (A.4) describes the balance between convective, turbulent diffusive, production and dissipation terms of the particle turbulence kinetic energy due to the particles interaction with the fluidic eddies. Here, the velocity correlations of the particulate phase are calculated with the help of the particles turbulent energy: ‹u'2s › = (2/3) ks and ‹v'2s › = (2/ 3) ks. These terms in x- and r- directions of the momentum equations of the particulate phase express the effect of turbophoresis force.

Acknowledgements The authors gratefully acknowledge support of this work by Norway-Estonian project EMP230 and by Brazilian National Council of Research (CNPq) through the project call “Science without borders”, protocol number: Proc. 405700/2013-0.

Fig. 4. Square root of the turbulent kinetic energy, k1/2. Results are scaled with the bulk flow velocity, Ub. Squares and triangles represent experimental data for fluid and particles, whereas solid and dashed lines represent numerical data for fluid and particles. Results are presented with data of case 2A.

1006


The coefficient Ds and the viscosity μs are given by Eqs. (A.5) and (A.6), respectively. Ds ¼

2 0 τ ðks þ g u kÞ 3 p

τ0p ks μ s ¼ ρ f u νt þ 3

ðA:5Þ

ðA:6Þ

where gu and fu are coefficients that characterize the particles involvement into the turbulent motion according to the PDF approach by Zaichik and Alipchenkov [52], and vt is the coefficient of turbulent viscosity. gu and fu are obtained with Eqs. (A.7) and (A.8). 1

fu ¼

1þ

gu ¼

ðA:7Þ

τ0p T Lp

τ0p T Lp = 1þ 0 τp T Lp

ðA:8Þ

TLp is the particle Lagrange integral time scale of the fluid along the particle trajectory when the moving particles interact with various eddies. TLp is calculated with the aid of a corrected Stokes number function as in accordance to Wang & Stock [53]: TLp = f(StE) TL. The particle Stokes number, StE, is calculated by Eq. (A.9), and the corrected Stokes function, by Eq. (A.10). St E ¼

τp C 0D T E

ðA:9Þ 2

TE 6 1− f ðSt E Þ ¼ 6 TL 4

3 TE −1 7 TL 7 0:4ð1þ0:01St E Þ 5 ð1 þ St E Þ

ðA:10Þ

where TL is the Lagrangian integral time scale of the liquid, TL = 0.3 (k/ ε), Sommerfeld [54]; and TE is the Eulerian integral time-scale of the liquid, TE = 2.81TL; see Zaichik et al. [34]. The particles turbulent energy ks is calculated by the differential equation Eq. (A.4) of the PDF model by Zaichik et al. [34]. The diffusion and turbulent viscosity coefficients, Eqs. (A.5) and (A.6), are calculated by the algebraic expressions derived by the PDF model by Zaichik et al. [34]. References [1] B. Wright, P. Zulli, Z.Y. Zhou, A.B. Yu, Gas-solid flow in an ironmaking blast furnace I: physical modelling, Powder Technol. 208 (2011) 86–97. [2] Z.Y. Zhou, H.P. Zhu, B. Wright, A.B. Yu, P. Zulli, Gas-solid flow in an ironmaking blast furnace II: discrete particle simulation, Powder Technol. 208 (2011) 72–85. [3] E. Peirano, V. Delloume, F. Johnsson, B. Leckner, O. Simonin, Numerical simulation of the fluid dynamics of a freely bubbling fluidized bed: influence of the air supply system, Powder Technol. 122 (2002) 69–82. [4] X. Pialat, O. Simonin, P. Villedieu, A hybrid Eulerian–Lagrangian method to simulate the dispersed phase in turbulent gas-particle flows, Int. J. Multiphase Flow 33 (2007) 766–788. [5] N.A. Konan, O. Kannengieser, O. Simonin, Stochastic modeling of the multiple rebound effects for particle–rough wall collisions, Int. J. Multiphase Flow 35 (2009) 933–945. [6] A. Ozel, P. Fede, O. Simonin, Development of filtered Euler–Euler two-phase model for circulating fluidised bed: high resolution simulation, formulation and a priori analyses, Int. J. Multiphase Flow 55 (2013) 43–63. [7] S.C.P. Cheung, G.H. Yeoh, F.S. Qi, J.Y. Tu, Classification of bubbles in vertical gas– liquid flow: part 2 – a model evaluation, Int. J. Multiphase Flow 39 (2012) 135–147. [8] Z.X. Zeng, L.X. Zhou, A two-scale second-order moment particle turbulence model and simulation of dense gas–particle flows in a riser, Powder Technol. 162 (2006) 27–32. [9] L.X. Zhou, Advances in studies on two-phase turbulence in dispersed multiphase flows, Int. J. Multiphase Flow 36 (2010) 100–108.

[10] A.I. Kartusinsky, E.E. Michaelides, M.T. Hussainov, Y. Rudi, Effects of the variation of mass loading and particle density in gas–solid particle flow in pipes, Powder Technol. 193 (2009) 176–181. [11] Y. Tsuji, Y. Morikawa, LDV measurements of an air–solid two-phase flow in a horizontal pipe, J. Fluid Mech. 120 (1982) 385–409. [12] Y. Tsuji, Y. Morikawa, H.H. Shiomi, LDV measurements of an air–solid two-phase flow in a vertical pipe, J. Fluid Mech. 139 (1984) 417–434. [13] J.D. Kulick, J.R. Fessler, J.K. Eaton, Particle response and turbulence modification in fully developed channel flow, J. Fluid Mech. 277 (1994) 109–134. [14] J. Kussin, M. Sommerfeld, Experimental studies on particle behaviour and turbulence modification in horizontal channel flow with different wall roughness, Exp. Fluids 33 (2002) 143–159. [15] N. Caraman, J. Borée, O. Simonin, “Effect of collisions on the dispersed phase fluctuation in a dilute tube flow: experimental and theoretical analysis, Phys. Fluids 15 (12) (2003). [16] T.S. Yang, S.S. Shy, Two-way interaction between solid particles and homogeneous air turbulence: particle settling rate and turbulence modification measurements, J. Fluid Mech. 526 (2005) 171–216. [17] M. Benson, T. Tanaka, J.K. Eaton, The effects of wall roughness on particle velocities in a turbulent channel flow, ASME J. Fluids Eng. 127 (2005) 250–256. [18] J. Borée, N. Caraman, Dilute bidispersed tube flow: role of interclass collisions at increased loadings, Phys. Fluids 17 (2005) 055108. [19] R. Volk, E. Calzavarini, E. Lévêque, J.-F. Pinton, Dynamics of inertial particles in a turbulent von Kármán flow, J. Fluid Mech. 668 (2011) 223–235. [20] Y. Sato, K. Hishida, Transport process of turbulence energy in particle-laden turbulent flow, Int. J. Heat Fluid Flow 17 (1996) 202–210. [21] Y. Suzuki, M. Ikenoya, N. Kasagi, Simultaneous measurements of fluid and dispersed phases in a particle-laden turbulent channel flow with the aid of 3-D PTV, Exp. Fluids (2000) 185–193. [22] J.L. Oliveira, G. 3D-PTV of particle-laden turbulent pipe flowspH.D. thesis Technische Universiteit Eindhoven, 2012. [23] S. Elgobashi, On predicting particle-laden turbulent flows, Appl. Sci. Res. 52 (1994) 309–329. [24] C. Poelma, J. Westerweel, Ooms, Turbulence statistics from optical whole-field measurements in particle-laden turbulence, Exp. Fluids 40 (3) (2006) 347–363. [25] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759. [26] W. Dijkhuizen, M. van Sint Annaland, J.A.M. Kuipers, Numerical and experimental investigation of the lift force on single bubbles, Chem. Eng. Sci. 65 (2010) 1274–1287. [27] A. Ferrante, S. Elgobashi, On the physical mechanisms of two-way coupling in particle-laden isotropic turbulence, Phys. Fluids 15 (2) (2003) 315–329. [28] S. Balachandar, J.K. Eaton, Turbulent dispersed multiphase flow, Annu. Rev. Fluid Mech. 42 (2010) (2010) 111–133. [29] A. Soldati, C. Marchioli, Physics and modelling of turbulent particle deposition and entrainment: review of a systematic study, Int. J. Multiphase Flow 35 (9) (2009) 827–839. [30] Y. Pan, S. Banerjee, Numerical simulation of particle interactions with wall turbulence, Phys. Fluids 8 (10) (1996) 2733–2755. [31] A.G. Kidanemariam, C. Chan-Braun, T. Doychev, M. Uhlmann, Direct numerical simulation of horizontal open channel flow with finite-size, heavy particles at low solid volume fraction, New J. Phys. 15 (2) (2013) 025031. [32] A.G. Kidanemariam, M. Uhlmann, Direct numerical simulation of pattern formation in subaqueous sediment, J. Fluid Mech. 750 (2014). [33] F. Picano, W.-P. Breugem, L. Brandt, Turbulent channel flow of dense suspensions of neutrally buoyant spheres, J. Fluid Mech. 764 (2015) 463–487. [34] L.I. Zaichik, V.M. Alipchenkov, E.G. Sinaiski, Particles in Turbulent Flows, Wiley-VCH, 2008. [35] H.-E. Albrecht, M. Borys, N. Damaschke, C. Tropea, Laser Doppler and Phase Doppler Measurement Techniques, Springer-Verlag, 2003. [36] L. Schiller, Z. Naumann, A drag coefficient correlation, Z. Ver. Deutsch. Ing. 77 (1935). [37] H. Faxén, Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist, Ann. Phys. 373 (10) (1922) 89–119. [38] P.G. Saffman, The lift on a small sphere in a slow shear flow, J. Fluid Mech. 22 (1965) 385–400. [39] R. Mei, An approximation expression for the shear lift on a spherical particle at finite Reynolds number, Int. J. Multiphase Flow 18 (1) (1992) 145–147. [40] T.R. Auton, The lift force on a spherical rotational flow, J. Fluid Mech. 183 (1987) 199–218. [41] J.D. Schwarzkopf, T.C. Crowe, P. Dutta, A model for particle laden turbulent flows, Proc. ASME FEDSM-2009-78385 2009, pp. 1–4 (August 2–6, Vail, Colorado, USA). [42] J.H. Fertziger, M. Perić, Computational Methods for Fluid Dynamics, Springer, 1996 356 pp. [43] M. Perić, G. Scheuerer, CAST – A finite volume method for predicting two-dimensional flow and heat transfer phenomena, GRS - Technische Notiz SRR– 89–01, 1989. [44] L.I. Zaichik, I.N. Gusev, E.I. Guseva, Turbulent flow and precipitation of particles in channels, Engineering Turbulence Modelling and Experiments, Elsevier, Amsterdam 1990, pp. 907–916. [45] L.I. Zaichik, Transport and deposition of inertial particles in turbulent pipe flow, ERCOFTAC Bull. 336 (1998) 16–19. [46] J.L. Binder, T.J. Hanratty, A diffusion model for droplet deposition in gas/liquid annular flow, Int. J. Multiphase Flow 17 (1) (1991) 1–11.

A. Kartushinsky et al. / Powder Technology 301 (2016) 999–1007 [47] J.L.G. Oliveira, C.W.M. van der Geld, J.G.M. Kuerten, Lagrangian and Eulerian statistics of pipe flows measured with 3D-PTV at moderate and high Reynolds numbers flow, Turbul. Combust. 91 (1) (2013) 105–137. [48] J.L.G. Oliveira, C.W.M. van der Geld, J.G.M. Kuerten, “Lagrangian velocity and acceleration statistics of fluid and inertial particles measured in pipe flow with 3D particle tracking velocimetry, Int. J. Multiphase Flow 73 (2015) 97–107. [49] J.O. Hinze, Turbulence, McGraw-Hill, 1975. [50] M.P.B. Veenman, Statistical Analysis of Turbulent Pipe Flow: A Numerical ApproachPh.D. thesis Technische Universiteit Eindhoven, 2004. [51] A. Aliseda, J.C. Lasheras, Preferential concentration and rise velocity reduction of bubbles immersed in a homogeneous and isotropic turbulent flow, Phys. Fluids 23 (2011) 093301.

1007

[52] L.I. Zaichik, V.M. Alipchenkov, Transport and algebraic models of particles kinetic stresses in turbulent flows, Proc. 2nd Int. Symposium on Turbulence and Shear Flow Phenomena Stockholm 2001, pp. 193–198. [53] L.-P. Wang, D.E. Stock, Dispersion of heavy particles in turbulent motion, J. Atmos. Sci. 50 (13) (1993) 1897–1913. [54] M. Sommerfeld, Modelling of particle-wall collisions in confined gas-particle flows, Int. J. Multiphase Flow 18 (1992) 905–926.