A numerical study on the migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection

International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A numerical study on the migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection Junjie Hu, Zhaoli Guo ⇑ State Key Laboratory of Coal Combustion, School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China

a r t i c l e

i n f o

Article history: Received 29 September 2016 Received in revised form 18 January 2017 Accepted 20 January 2017

Keywords: Particle migration Poiseuille flow Thermal convection Lattice Boltzmann method

a b s t r a c t The migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection is investigated with the double-population lattice Boltzmann method. Compared with the isothermal case, the migration of the particle is notably affected by thermal convection. Five competing mechanisms are responsible for the lateral migration of the particle, i.e., wall repulsion due to lubrication, inertial lift related to shear slip, a lift due to particle rotation, a lift due to the curvature of the undisturbed velocity profile and a lift induced by thermal convection. Within the parameter range examined (Gr < 2000; Rec < 200), a critical Grashof number Grc is identified, beyond which the equilibrium position is insensitive to the initial position. Furthermore, a power law dependence of the critical Grashof number Gr c on the channel Reynolds number Rec is discovered. In addition, with the channel Reynolds number Rec increasing, the forced flow relative to the particle is becoming stronger, leading to a declining effect of thermal convection. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Particle transport is of vital importance in numerous industrial and natural processes, such as fluidized beds, lubricated transport, fluvial erosion, river sediment and sand storms [1,2]. Segré and Silberberg [3] conducted a particularly important experimental study on the motion of particles in a pipe flow, and discovered that neutrally buoyant particles always migrate to a certain lateral equilibrium position, which is the so-called Segré-Silberberg effect. Karnis et al. [4] verified the Segré-Silberberg effect, and observed that the lateral equilibrium position is related to the flow rate and dimensionless particle diameter. Tachibana et al. [5] studied the lateral migration of spheres in pipe flows experimentally, and reported that the ratio of the sphere diameter to the pipe diameter is vital to reproduce the Segré-Silberberg effect. Besides experimental studies, theoretical investigations of the particle migration have been conducted. For example, Saffman et al. [6] obtained the lift on a spherical particle in a shear flow in an unbounded domain with the perturbation theory. Ho et al. [7] and Vasseur et al. [8] investigated the lateral migration of a spherical particle in a plane Poiseuille flow bounded by two infinite plane walls. Schonberg et al. [9] extended Saffman’s analysis to a small sphere in a plane Poiseuille flow for the channel Reynolds number of order one. ⇑ Corresponding author. E-mail address: [email protected] (Z. Guo). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.01.077 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

In addition, some particle-resolved numerical simulations have been carried out to verify and extend the experimental and theoretical observations. For example, Feng et al. [10] investigated the motion of a circular particle in a Poiseuille flow with the finite-element method and obtained qualitative agreement with the results of experiments and perturbation theories. Inamuro et al. [11] studied the motion of a single and two lines of neutrally buoyant circular particles in a Poiseuille flow with the lattice Boltzmann method, and the Segré-Silberberg effect is reproduced. Also with the lattice Boltzmann method, Basßag˘aog˘lu et al. investigated colloid migration in rough-walled channels [12] and particle migration in microchannels [13]. Patankar et al. [14] studied the lift force and equilibrium position of a circular particle in a Poiseuille flow by an arbitrary Lagrangian-Eulerian finite-element method. In the above studies, the underlying assumption is that the fluid and the particle are isothermal. However, thermal convection is usually involved in numerous processes, such as fluidized reactors, pneumatic transport and volcanic eruptions. It is reported that the motion of particles is notably affected by thermal convection in previous studies [15–22]. For example, Gan et al. [15] investigated the particle settling in a vertical channel with the finite-element method, and discovered that thermal convection may fundamentally change the motion of the particle. For a single particle settling in a channel, the particle may settle straight down or migrate toward a wall or oscillate laterally, and various Grashof-number regimes are identified. Instead of the well-known drafting-

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

kissing-tumbling scenario, a pair of particles tend to separate if they are colder than the fluid and aggregate if they are hotter. These effects are attributed to the competition between thermal convection and forced flow relative to the particle. Kim et al. [19] studied the motion of a circular particle in a two-dimensional vertical channel with different temperatures at the left and right walls with the fictitious domain method, and three Grashof-number regimes are identified in terms of the particle behavior, i.e., steady settling, a transient overshoot before the steady settling and thermal levitation. Yang et al. [22] investigated two cold particles settling in a vertical channel with the lattice Boltzmann method. With different initial separation distances and relative angles between the two cold particles, three regimes are identified, i.e., repulsion, attraction and transition. Compared with their isothermal counterparts, the interaction between two cold particles is significantly influenced by thermal convection. However, the previous studies mostly focus on the particle sedimentation in a channel with thermal convection. To our knowledge, studies on the migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection are rarely reported in the literature. To deepen our understanding, it is necessary to understand the effect of thermal convection on the particle migration in a Poiseuille flow, which is the aim of the present work. The rest of this paper is organized as follows. Section 2 gives the problem description, Section 3 introduces the double-population lattice Boltzmann method for flows with thermal convection, and our code is validated in Section 4. In Section 5, the effect of thermal convection on the migration of a neutrally buoyant particle in a Poiseuille flow is discussed, and some conclusions are given in Section 6.

2159

boundaries [25] is adopted to implement the velocity and pressure boundary conditions. The dimensionless parameters involved in the present study are as follows,  confinement ratio

c¼

D ¼ 0:25; W

ð1Þ

 channel Reynolds number

Rec ¼

WU m

m

;

ð2Þ

 particle Reynolds number

Rep ¼

DU s

m

;

ð3Þ

 Prandtl number

Pr ¼

m ¼ 0:7; a

ð4Þ

 Grashof number

Gr ¼

gbD3 ðT h  T c Þ

m2

2 ½0; 2000;

ð5Þ

where U m ¼ W 2 Dp=ð8LmÞ is the undisturbed maximum velocity of the fluid at the centerline, U s is the steady velocity of the particle when a migration to the equilibrium position is achieved, m is the kinematic viscosity of the fluid, a is the thermal diffusivity of the fluid, g is the gravity acceleration, and b is the thermal expansion coefficient of the fluid.

2. Problem description As depicted in Fig. 1, a neutrally buoyant particle with diameter D is located in a channel with length L ¼ 20D and width W ¼ 4D. Initially, the particle is positioned at (x0 ; y0 ), where x0 ¼ L=2 while different values are assigned to y0 to investigate the effect of thermal convection on the particle migration. To decrease the computational cost, a moving domain technique is adopted, that is, when the particle moves a lattice unit in the longitudinal direction, the mesh is shifted a grid spacing in the same direction [23,24]. The temperature of the particle is fixed at T c ¼ 0, while the temperatures of the lower and upper walls of the channel are kept at T h ¼ 1, and the initial fluid temperature is T h as well. The no-slip boundary conditions are applied to the lower and upper walls of the channel, while the pressure boundary conditions are imposed for the inlet and outlet with a pressure difference Dp ¼ pi  po , where the non-equilibrium extrapolation method for straight

3. Numerical method 3.1. Double-population lattice Boltzmann method The double-population lattice Boltzmann method is adopted in the present study, in which the flow field and temperature field are solved respectively with two lattice Boltzmann equations (LBEs) [26],

f i ðx þ ci dt; t þ dtÞ  f i ðx; tÞ ¼ Xi ðf Þ;

ð6Þ

g i ðx þ c i dt; t þ dtÞ  g i ðx; tÞ ¼ Xi ðgÞ;

ð7Þ

where f i ðx; tÞ and g i ðx; tÞ are the distribution functions at position x and time t with discrete velocity c i for the flow field and temperature field, respectively, Xi ðf Þ and Xi ðgÞ are the corresponding discrete collision operators, and dt is the time step. In practice, the evolution of the distribution functions is decomposed into two steps, i.e., collision step þ

f i ðx; tÞ ¼ f i ðx; tÞ þ Xi ðf Þ;

ð8Þ

g þi ðx; tÞ ¼ g i ðx; tÞ þ Xi ðgÞ;

ð9Þ

and streaming step þ

f i ðx þ ci dt; t þ dtÞ ¼ f i ðx; tÞ;

ð10Þ

g i ðx þ c i dt; t þ dtÞ ¼ g þi ðx; tÞ;

ð11Þ

þ

Fig. 1. Sketch of a neutrally buoyant particle migrating in a Poiseuille flow with thermal convection.

where f i ðx; tÞ and g þ i ðx; tÞ are the post-collision distribution functions for the flow field and temperature field, respectively. The multi-relaxation-time (MRT) model is adopted for the flow field due to its superior numerical stability [23,27], while the popular single-relaxation-time (SRT) model is adopted for the temperature field since the temperature field is passively advected by the fluid flow and a simple scalar equation is obeyed if the viscous heat dissipation and compression work are negligible, i.e.,

2160

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168






1 2

 

1 Xi ðf Þ ¼ ðM1 SMÞij ðf j  f eq I  S M Fj; j Þ þ dt M

Xi ðgÞ ¼ 

1

sg

ð12Þ

ij

ðg i  g eq i Þ;

ð13Þ

eq

where f i and g eq i are the equilibrium distribution functions for the flow field and temperature field, respectively, M is the transform matrix, S is the diagonal relaxation matrix for the flow field, sg is the dimensionless relaxation time for the temperature field, I is the identity matrix, and F i is the forcing term dependent on the buoyancy force induced by thermal convection. In the present study, the D2Q9 (2-dimension and 9-velocity) model [28] is adopted, and the equilibrium distribution functions are defined as

# c i  u ðc i  uÞ2 u  u q; uÞ ¼ wi q 1 þ 2 þ  2 ; cs 2cs 2c4s
 ci  u ; g eq i ðT; uÞ ¼ wi T 1 þ c2s

The gravity force experienced by the fluid is then given by

G ¼ qg ¼ q0 ½1  bðT  T h Þg ¼ q0 g  q0 gbðT  T h Þ;

where the first term q0 g is a constant which is absorbed into the pressure term, while the second term is the buoyancy force term induced by the temperature difference, namely,

F ¼ q0 gbðT  T h Þ:

F i ¼ wi

ð15Þ

8 for i ¼ 0; > < 4=9; wi ¼ 1=9; for i ¼ 1; 2; 3; 4; > : 1=36; for i ¼ 5; 6; 7; 8;

ð16Þ

c i is the discrete velocity given by

8 for i ¼ 0; > < ð0; 0Þ; for i ¼ 1; 2; 3; 4; c i ¼ cðcos½ði  1Þp=2; sin½ði  1Þp=2Þ; > : pffiffiffi 2cðcos½ð2i  1Þp=4; sin½ð2i  1Þp=4Þ; for i ¼ 5; 6; 7; 8; ð17Þ

pffiffiffi cs ¼ c= 3 is the sound speed of the model, c ¼ dx=dt is the lattice speed, dx is the grid spacing. Generally, the transform matrix M is

1

1

1

1

1

1

1

1

1

B 4 1 1 1 1 2 2 2 2 C C B C B C B 4 2 2 2 2 1 1 1 1 C B C B 1 0 1 0 1 1 1 1 C B 0 C B 2 0 1 1 1 1 C M¼B C; B 0 2 0 C B 0 0 1 0 1 1 1 1 1 C B C B B 0 0 2 0 2 1 1 1 1 C C B C B 1 1 1 1 0 0 0 0 A @ 0 0 0 0 0 0 1 1 1 1

by

q¼

X f i;

qu ¼

X

i

ci f i þ

i

dt F; 2

T¼

X

gi :

ð24Þ

i

Through the Chapman-Enskog expansion, the macroscopic equations can be derived from the two LBEs, i.e.,

@q þ r  ðquÞ ¼ 0; @t

ð25Þ

@u 1 þ ðu  rÞu ¼  rp þ mr2 u þ F; @t q

ð26Þ

@T þ r  ðTuÞ ¼ r  ðarTÞ; @t

ð27Þ

where the fluid pressure is

p ¼ qc2s ;

ð28Þ

the kinematic viscosity of the fluid is related to the relaxation time

ss as




m ¼ c2s ss 




ð18Þ

 1 dt; 2

ð29Þ

a ¼ c2s sg 

 1 dt: 2

ð30Þ

3.2. Unified iterative scheme for moving boundaries

ð19Þ

where sq and sj are related to the conserved moments and their values do not influence the computation, which are set to be sq ¼ sj ¼ ss , se and s can be determined by the linear stability analysis [23,27], which are set to be se ¼ 0:9 and s ¼ 0:8, sq is determined by the magic relation, i.e., sq ¼ ð8ss  1Þ=ð16ss  8Þ, to enforce the no-slip boundary condition for straight boundaries [29,30], and ss is related to the kinematic viscosity of the fluid and is set to be ss ¼ 0:75 in the present study. To couple the flow field and temperature field, the Boussinesq approximation is adopted, in which the fluid properties are assumed to be constant except that the density is linearly dependent on the temperature, i.e.,

q ¼ q0 ½1  bðT  T h Þ;

ð23Þ

and the thermal diffusivity of the fluid is related to the relaxation time sg as

and the corresponding diagonal relaxation matrix S is

S ¼ diagðsq ; se ; s ; sj ; sq ; sj ; sq ; ss ; ss Þ1 ;

  c i  u ðc i  uÞ  F: þ c i c2s c4s

The fluid density q, velocity u and temperature T are obtained

ð14Þ

where q; u and T are the fluid density, velocity and temperature, respectively, wi is the weight coefficient of the model which is given by

1

ð22Þ

To implement the buoyancy force, the forcing scheme proposed by Guo et al. [31] is adopted, i.e.,

"

eq fi ð

0

ð21Þ

ð20Þ

where q0 is the fluid density at the intimal fluid temperature T h .

In the implementation of the lattice Boltzmann method to moving boundaries, two key issues should be addressed, i.e., how to treat the curved surface of a solid particle on a uniform Cartesian grid and how to initialize a fresh node coming from the moving particle. In previous studies, these two issues are usually considered separately [32,33]. Recently, Hu et al. [34] developed a unified iterative scheme (UIS) to treat these two issues consistently, which is demonstrated to be a superior treatment for moving boundaries. As depicted in Fig. 2, after the collision step and streaming step, the distribution functions at fluid nodes and boundary nodes with links to other fluid nodes are updated. For a boundary node or a fresh node, the unknown distribution functions are determined following two steps, i.e., prediction step and correction step. In the prediction step, the unknown distribution functions are approximated. Here, for simplicity, we extrapolate the unknown distribution functions along the direction of a discrete velocity c i which maximizes the quantity ci  n, where n is the outward unit normal vector at point xw . As depicted in Fig. 2, the unknown distribution functions at node xf (it can be a boundary node or a fresh node) are predicted with f i ðxf Þ ¼ f i ðxff Þ.

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

2161

Fig. 3. Sketch of the Galilean invariant momentum exchange method. The momentum carried by the fluid particle towards the solid surface C before the þ collision is ðc i  uw Þf i ðxf ; tÞ relative to the solid surface C, while the relative momentum is ðci  uw Þf i ðxf ; t þ dtÞ after the collision with the solid surface C.

Fig. 2. Sketch of the unified iterative scheme for moving boundaries. C is the solid surface, xff is the fluid node, xf is the boundary node, xw is the intersection point between the lattice link and the solid surface, and xs is the solid node.

To reflect the boundary effect, a correction step is followed, where the distribution functions are decomposed into their equilibrium parts and non-equilibrium parts. After the prediction step, the unknown distribution functions at node xf are approximated, therefore, the density qf and velocity uf at node xf can be obtained,

field, the unknown distribution functions for the temperature field at node xf (it can be a boundary node or a fresh fluid node) are firstly predicted with g i ðxf Þ ¼ g i ðxff Þ. After the prediction step, the unknown distribution functions for the temperature field at node xf are approximated, as a result, the temperature T f at node xf can be obtained, as well as the equilibrium parts g eq i ðxf Þ and nonequilibrium parts g ne i ðxf Þ. Since the temperature T w and velocity uw at point xw are known, the equilibrium parts g eq i ðxw Þ can be determined. The non-equilibrium parts can be achieved by a first-order extrapolation with those at node xf , which is similar to those for the flow field. After the fictitious distribution functions at point xw are constructed, a first-order interpolation is adopted to correct the unknown distribution functions for the temperature field at node xf , just like Eq. (31). Similarly, to improve the numerical accuracy and decrease the inconsistency, an enforced iteration is adopted.

eq

as well as the equilibrium parts f i ðxf Þ and non-equilibrium parts ne f i ðxf Þ. Since the density at point xw can be well approximated with that at node xf for a nearly incompressible flow, i.e., qw ¼ qf , and the velocity at point xw can be obtained, i.e., eq uw ¼ up þ xp  ðxw  xc Þ, the equilibrium parts f i ðxw Þ can be determined, where up and xp are the translational velocity and angular velocity of the solid body, respectively, and xc is the center of the solid body. The non-equilibrium parts can be achieved by a ne ne first-order extrapolation with those at node xf , i.e., f i ðxw Þ ¼ f i ðxf Þ. So far, the fictitious distribution functions at point xw are coneq ne structed, i.e., f i ðxw Þ ¼ f i ðxw Þ þ f i ðxw Þ. To correct the unknown distribution functions at node xf , a first-order interpolation is adopted,

f i ðxf Þ ¼ ðqfi ðxff Þ þ f i ðxw ÞÞ=ðq þ 1Þ;

ð31Þ

where q ¼ jxf  xw j=jxf  xs j is the fraction of the intersected link in the fluid region. Different from previous boundary conditions and refilling algorithms [32,33], UIS is a multi-step method, where an enforced iteration is adopted to enforce the no-slip boundary condition and decrease the inconsistency between the constructed distribution functions and those evolutionary ones. In each iteration step, the correction step is executed, and the unknown distribution functions, density and velocity at node xf are updated. The enforced iteration is executed until the error of the velocity at node xf is less

3.3. Hydrodynamic force on a solid particle To update the velocity and position of the particle, it is crucial to evaluate the force exerted on the particle by the fluid accurately. Here, due to its efficiency and robustness [35], the momentum exchange method (MEM) is adopted to evaluate the hydrodynamic force. The original MEM was proposed by Ladd [36,37]. However, this scheme is not Galilean invariant, and a Galilean invariant momentum exchange method (GIMEM) was proposed [38]. Recently, Peng et al. [35] demonstrated that the Galilean invariance is crucial for correct moving fluid-solid interaction. In GIMEM, the momentum transfer is correlated to the relative velocity of the fluid particle to solid surface instead of the absolute velocity of the fluid particle, which is more physically reasonable [35]. As depicted in Fig. 3, the momentum carried by the fluid particle towards the solid surface C before the collision is þ ðc i  uw Þf i ðxf ; tÞ relative to the solid surface C, while the relative momentum is ðc i  uw Þf i ðxf ; t þ dtÞ after the collision with the þ

solid surface C, where f i ðxf ; tÞ is the post-collision distribution functions at time t; f i ðxf ; t þ dtÞ is the updated distribution functions at time t þ dt which is constructed with UIS for curved boundaries (Section 3.2). Therefore, the hydrodynamic force exerted on the particle is

Fh ¼

xf

6

than a criterion (e.g., 10 ) or the iteration step number is greater than a limit (e.g., 100). Although UIS is presented to treat moving boundaries consistently without thermal convection, i.e., the fluid and the particle are isothermal, it is straightforward to be extended to treat moving boundaries with thermal convection. Similar to those for the flow

XX

 þ ðc i  uw Þf i ðxf ; tÞ  ðc i  uw Þf i ðxf ; t þ dtÞ ;

ð32Þ

i

and the torque exerted on the particle is

Th ¼

XX  

þ ðxw  xc Þ  ðc i  uw Þf i ðxf ;tÞ  ðci  uw Þf i ðxf ;t þ dtÞ : xf

i

ð33Þ

2162

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

Fig. 4. Time history of the lateral position of the particle with different initial positions.

3.4. Particle movement The particle movement is governed by the Newtonian motion equations, i.e.,

dup ¼ F h þ ðqp  qf ÞgAp ; dt dxp Ip ¼ Th; dt dxp ¼ up ; dt dhp ¼ xp ; dt

mp

ð34Þ

Fig. 6. Time history of the lateral position of the particle.

Table 1 Equilibrium positions of the particle at Gr ¼ 1000 and 2000. Gr

Present

Yu et al. [16]

Feng et al. [17]

Yang et al. [22]

1000 2000

2.91 2.76

2.89 2.74

2.90 2.73

2.89 2.79

ð35Þ i dt h tþdt F þ ðqp  qf ÞgAp ; mp h dt ¼ xtp þ T htþdt ; Ip

ð36Þ

uptþdt ¼ utp þ

ð38Þ

ð37Þ

xtþdt p

ð39Þ

where mp ¼ qp Ap and Ip ¼ mp D2 =8 are the mass and moment inertia of the particle, respectively, xp and hp are the position and angular position of the particle, respectively, qp and qf are the particle den-

xtþdt ¼ xtp þ dtuptþdt ; p

ð40Þ

htþdt p

ð41Þ

sity and fluid density, respectively, Ap ¼ pD2 =4 is the area of the particle. To solve the above equations, the first-order forward Euler scheme is adopted, i.e.,

where the superscript t and t þ dt denote the corresponding physical quantity at time t and t þ dt, respectively.

¼

htp

tþdt ; p

þ dt x

4. Code validation To validate our code, simulations of a neutrally buoyant isothermal particle moving in a Poiseuille flow and a cold particle settling in a vertical channel are conducted.

4.1. A neutrally buoyant isothermal particle moving in a Poiseuille flow

Fig. 5. Sketch of a cold particle settling in a vertical channel.

Segré and Silberberg [3] were first to study the problem experimentally, and discovered that a neutrally buoyant particle always migrates to a lateral equilibrium position between the centerline and the wall. As depicted in Fig. 1, the fluid and the particle are isothermal, i.e., T h ¼ T c ¼ 0. For convenience of comparison, we should take the same channel Reynolds number Rec , as a result, the pressure difference is set to be Dp ¼ 0:00267 [39], corresponding to a channel Reynolds number Rec ¼ 96:12. Fig. 4 shows the time history of the lateral position of the particle with different initial positions. It is observed that the particle always migrates to a certain lateral equilibrium position regardless of its initial positions, as discovered by Segré and Silberberg [3]. At the same time, the grid independence is validated by varying the grid resolution, i.e., D ¼ 25dx and D ¼ 35dx, while the channel Reynolds number Rec is kept fixed, and almost the same results are obtained with these two grid resolutions. Given the computational cost, the coarse grid resolution D ¼ 25dx is adopted in the following simulations unless otherwise specified. Li et al. [39] reported

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

2163

Fig. 7. Lateral equilibrium position of the particle yeq against Gr at Rec ¼ 96:12 : Gr 2 ½0; 2000 (a), Gr 2 ½0; 100 (b), Gr 2 ½60; 70 (c).

Fig. 8. Particle Reynolds number Rep against Gr at Rec ¼ 96:12 : Gr 2 ½0; 2000 (a), Gr 2 ½0; 100 (b), Gr 2 ½60; 70 (c).

Fig. 9. Time evolution of the lateral position of the particle at Gr ¼ 64 and 65.

that the lateral equilibrium position is 0.2874, and our result is 0.2872, which is in good agreement. 4.2. A cold particle settling in a vertical channel No thermal effect is involved in the above case. To validate our code for flows with thermal convection, a cold particle settling in a vertical channel is simulated. This problem was firstly investigated by Gan et al. [15] with the finite-element method, and was regarded as a benchmark by many researchers to validate their numerical methods [16–22]. As depicted in Fig. 5, a particle with diameter D ¼ 25dx is located in a vertical channel with width W ¼ 4D and height H ¼ 320D. Initially, the particle is positioned

Fig. 10. Time evolution of the shear stress on the lower wall of the channel at Gr ¼ 64 and 65.

off the centerline by half the diameter. The density ratio between the particle and the fluid is qr ¼ qp =qf ¼ 1:00232. The reference Reynolds number Rer ¼ DU r =m is set to be 40.5, where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U r ¼ pðD=2Þðqr  1Þg is the reference velocity. The temperature of the particle is fixed at T c ¼ 0, while temperatures of the side walls of the channel are kept at T h ¼ 1, and the initial fluid temperature is T h as well. Fig. 6 shows the time history of the lateral position of the particle. Gan et al. [15] reported that the centerline of the channel is no longer the equilibrium position of the particle for 810 < Gr < 2150, instead, a steady settling of the particle close to one of the side walls is observed, which is in agreement with our result qualitatively. Furthermore, the equilibrium positions of the

2164

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

Fig. 11. Vortex structure around the particle at Gr ¼ 64 and 65.

Fig. 12. Time evolution of the slip velocity at Gr ¼ 64 and 65.

Fig. 13. Time evolution of the angular velocity of the particle at Gr ¼ 64 and 65.

particle at Gr ¼ 1000 and 2000 are listed in Table 1, together with those in the literature. It is observed that the agreement between our results and those in the literature is satisfactory.

5. Results and discussions Now, we come to the problem described in Section 2. The pressure difference is set to be Dp ¼ 0:0027, and the corresponding

Fig. 14. Time evolution of the lateral position of the particle at Gr ¼ 64 and 65 with and without rotating.

channel Reynolds number is Rec ¼ 96:12. To investigate the effect of thermal convection on the particle migration, different initial positions are studied, i.e., y0 =W ¼ 0:25, 0.5 and 0.75. The lateral equilibrium position of the particle yeq against Gr is depicted in Fig. 7. At the very start, a rapid decreasing of yeq with Gr increasing is observed, following by a slow decreasing until a steady lateral equilibrium position around 0.2146 is achieved. In addition, it is observed that the equilibrium position is dependent on the initial position for small Gr, while it is becoming independent for large one, indicating that there is a critical Grashof number Gr c , beyond which the equilibrium position is insensitive to the initial position. If the initial position of the particle is above the centerline, i.e., y0 =W ¼ 0:75, the equilibrium position is approaching the centerline for Gr < 60. However, it is interesting that a sudden drop of the equilibrium position at Gr ¼ 70 is observed. Furthermore, the sudden drop is observed between Gr ¼ 64 and 65 (Fig. 7(c)). To rule out the possibility of numerical artifacts, simulations with a finer grid D ¼ 35dx are conducted, and almost the same results are obtained (Fig. 7(c)). Recently, Peng et al. [35] pointed out that the results may depend on the specific implementation for moving fluid-solid interaction with the lattice Boltzmann method, such as the implementation of the boundary condition and refilling algorithm. To validate our results are robust, simulations with an alternative implementation of the boundary

2165

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

Fig. 15. Temperature field around the particle at Gr ¼ 64 and 65.

Table 2 Parameters.

Case Case Case Case Case

Fig. 16. Five competing mechanisms are responsible for the lateral migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection.

condition and refilling algorithm are conducted, where the boundary nodes are treated with the non-equilibrium extrapolation (NEE) method for curved boundaries [32] and fresh fluid nodes are treated with the averaged extrapolation (AE) refilling algorithm [33]. Once again, we obtain almost the same results (Fig. 7(c)). Fig. 8 shows the particle Reynolds number Rep against Gr, and the similar instability phenomenon is observed (Fig. 8(c)). Generally, a slight decreasing of Rep with Gr increasing is observed, reflecting the fact that the steady velocity of the particle is decreasing with Gr increasing. Since the undisturbed velocity profile of the

1 2 3 4 5

Dp

Um

Rec

0.0007 0.0013 0.0027 0.0040 0.0053

0.02 0.04 0.08 0.12 0.16

24.03 48.06 96.12 144.1 192.2

fluid is parabolic with its maximum at the centerline, while the equilibrium position of the particle is approaching the wall with Gr increasing, as a result, the velocity of the particle closer to the wall is smaller, leading to a smaller Rep . For Gr < 60, it is observed that Rep is increasing with Gr increasing for the initial position y0 =W ¼ 0:75, since the equilibrium position of the particle is approaching the centerline (Fig. 7(b)), resulting in a larger particle velocity. To explore the interesting instability phenomenon, the time evolution of the lateral position of the particle at Gr ¼ 64 and 65 is depicted in Fig. 9. For the initial position y0 =W ¼ 0:75, a bifurcation phenomenon around tm=D2 ¼ 1 is observed. Before t m=D2 ¼ 1, the trajectory of the particle is almost the same at Gr ¼ 64 and 65, however, quite different behavior is observed after the critical time. Specially, a continuing migration of the particle toward the lower equilibrium position at Gr ¼ 65 is observed, while a turning around of the particle at Gr ¼ 64 is observed, following by a migration toward the upper equilibrium position. In addition, another initial position y0 =W ¼ 0:6 is simulated, and the similar bifurcation phenomenon around tm=D2 ¼ 1 is observed. For the isothermal case, Feng et al. [10] concluded that the lateral migration of a neutrally buoyant particle in a Poiseuille flow is determined by four competing mechanisms, i.e., wall repulsion due to lubrication, inertial lift related to shear slip, a lift due to particle rotation and a lift due to the curvature of the undisturbed velocity profile. With thermal convection involved, an additional mechanism should be included, i.e., a lift induced by thermal convection. To dig out the reason of the interesting bifurcation phenomenon, the above mechanisms are analyzed one by one in the following

Fig. 17. Lateral equilibrium positions of the particle yeq against Gr at Rec ¼ 96:12: a cold particle migrating in the hot fluid and a hot particle migrating in the cold fluid.

part. Around tm=D2 ¼ 1, the position of the particle is below the centerline, since the wall repulsion is always inclined to push the particle toward the centerline, as a result, the direction of the wall repulsion is upward. Since it is difficult to quantify the wall

2166

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

where l ¼ mqf is the dynamic viscosity of the fluid. As depicted in Fig. 10, compared with that at Gr ¼ 64, a sudden drop of the shear

Fig. 18. Critical Grashof number Gr c against Rec . The symbols are the numerical results, and the dashed line is the fitting curve.

repulsion due to lubrication directly, an alternative method is to obtain the shear stress on the lower wall of the channel, which is evaluated as




sab ¼ lðru þ ruT Þ ¼  1 

 1 X eq ðf  f i Þc ai c bi ; 2ss i i

ð42Þ

stress around t m=D2 ¼ 4:5 at Gr ¼ 65 is observed. To explain the sudden drop, the quite different vortex structure around the particle at Gr ¼ 64 and 65 is shown in Fig. 11. At Gr ¼ 65, vortexes are formed between the lower side of the particle and the lower wall, while vortexes are not notable at Gr ¼ 64. The stronger interaction between the vortex and the wall may be responsible for the drop of the shear stress. As a result, we may speculate that a much weaker wall repulsion is exerted on the particle at Gr ¼ 65 than that at Gr ¼ 64. To analyze the contribution of the inertial lift related to shear slip, a slip velocity uslip is defined, which is the difference between the particle velocity and the undisturbed velocity of the fluid at the center of the particle. As depicted in Fig. 12, the slip velocity is negative, that is, the particle lags the fluid. Feng et al. [10] pointed out that the particle is close to the centerline if it lags the fluid while close to a wall if it leads the fluid. Therefore, the direction of the inertial lift related to shear slip is upward as well. In addition, since the slip velocity at Gr ¼ 64 and 65 is approximately equal, thus, the contribution of the inertial lift to the lateral migration of the particle is comparable. Fig. 13 shows the time evolution of the angular velocity of the particle at Gr ¼ 64 and 65. Around t m=D2 ¼ 1, a shift of the

Fig. 19. Pressure distribution around the particle at the equilibrium state.

2167

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168 Table 3 Steady equilibrium positions of the particle yeq;s at different Rec . Rec

24.03

48.06

96.12

144.1

192.2

yiso eq;s

0.2883

0.2879

0.2872

0.2867

0.2866

yeq;s

0.1943 0.0940

0.2145 0.0734

0.2146 0.0726

0.2146 0.0721

0.2148 0.0718

yiso eq;s  yeq;s

direction of the angular velocity is observed, i.e., from positive to negative. Due to the Magnus effect, a lift is exerted on the particle, whose direction is downward if the angular velocity is positive while upward if the angular velocity is negative. To demonstrate the effect of the lift due to particle rotation, a comparative study is conducted, where the particle is not allowed to rotate. As depicted in Fig. 14, it is interesting that a continuing migration of the particle toward the lower equilibrium position at Gr ¼ 64 and 65 is observed if the particle is without rotating, instead of the bifurcation phenomenon around tm=D2 ¼ 1. Thus, the lift due to particle rotation is a crucial contribution to the interesting bifurcation phenomenon. At the same time, it is observed that the lateral equilibrium position is closer to the centerline without rotating, which is in agreement with the report by Feng et al. [10]. Then, we analyze the effect of the lift due to the curvature of the undisturbed velocity profile. Around t m=D2 ¼ 1, the velocity of the particle in the longitudinal direction is 0:6759U m , and the undisturbed fluid velocity at the center of the particle is 0:9471U m ; 0:7696U m at its lower tip and 0:9996U m at its upper tip. It is clear that a higher velocity of the fluid relative to the particle at the upper side is created by the curvature of the undisturbed velocity profile, as a result, the particle is sucked toward the centerline by the lower pressure at the upper side induced by the stronger flow. As the particle is cold while the fluid is hot, as a result, the fluid around the particle will be cooled by the particle. According to the Boussinesq assumption, the fluid density will increase with its temperature decreasing, and a downward thermal convection will impinge on the upper side of the particle. Fig. 15 shows the temperature field around the particle at Gr ¼ 64 and 65. At Gr ¼ 64, the cold fluid layer is mostly distributed at the left side of the particle, resulting in a weaker lift induced by thermal convection, while the cold fluid layer is concentrated at the upper side of the particle at Gr ¼ 65, leading to a more effective lift. From the above analysis, five competing mechanisms are responsible for the lateral migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection, as depicted in Fig. 16. Among these five competing mechanisms, three of them are directed upward, i.e., wall repulsion due to lubrication, inertial lift related to shear slip and a lift due to the curvature of the undisturbed velocity profile, while two of them are directed downward, i.e., a lift due to particle rotation and a lift induced by thermal convection. At Gr ¼ 65, the upward wall repulsion is reduced significantly while the downward lift due to thermal convection is enforced effectively, as a result, the downward lift is just able to conquer the upward lift, and a continuing migration of the particle toward the lower equilibrium position is observed. However, the downward lift is unable to conquer the upward lift eventually at Gr ¼ 64, and a turning around of the particle is observed, leading to the interesting bifurcation phenomenon. In the above discussion, the particle is cold while the fluid is hot. For the opposite case, i.e., the hot particle migrating in the cold fluid, the same conclusion should be expected, since the direction of thermal convection is just opposite. As depicted in Fig. 17, it is observed that the lateral equilibrium positions of the particle are symmetric about the centerline in these two cases. In addition, a bifurcation phenomenon

between Gr ¼ 64 and 65 for the hot particle migrating in the cold fluid is observed as well, which is just the same to that for the cold particle migrating in the hot fluid. From the above discussion, a critical Grashof number Gr c is observed, beyond which the equilibrium position is insensitive to the initial position. In addition, it is reported that the lateral migration of a neutrally buoyant particle is affected significantly by the flow rate [4,11], which is represented by the channel Reynolds number Rec . Next, we will investigate the effect of Rec on the critical Grashof number Gr c and steady equilibrium position yeq;s , and the parameters are listed in Table 2. In simulations, we vary the pressure difference Dp to obtain different channel Reynolds numbers while maintaining a small Mach number. The critical Grashof number Gr c against Rec is depicted in Fig. 18. A monotonic increasing of Gr c with Rec increasing is observed. Furthermore, a least-square fitting is conducted, and an approximate power law dependence of Gr c on Rec is identified, i.e., Gr c ¼ 0:01385Re1:841 , where the relative error is within 10%. c Fig. 19 shows the pressure distribution around the particle at the equilibrium state. It is observed that the pressure at the lower side of the particle is increasing with Rec increasing, leading to a stronger wall repulsion due to lubrication. To conquer the increasing wall repulsion, the stronger thermal convection is desired, corresponding to an increasing Gr c . The steady equilibrium positions of the particle at different Rec are listed in Table 3, together with those without thermal convection, i.e., the particle and the fluid are isothermal. For the isothermal case, it is observed that the equilibrium position is closer to the wall with Rec increasing, which is in agreement with the report by Feng et al. [10]. Compared with that without thermal convection, the equilibrium position with thermal convection is always closer to the wall, which is related to the downward lift induced by thermal convection. In addition, it is observed that the difference of the equilibrium position without and with thermal convection yiso eq;s  yeq;s is declining with Rec increasing. Gan et al. [15] attributed the special phenomenon with thermal convection to the competition between the forced flow relative to the particle and thermal convection, since the forced flow relative to the particle is becoming stronger with Rec increasing, as a result, a declining effect of thermal convection is observed. 6. Conclusions The migration of a neutrally buoyant particle in a Poiseuille flow with thermal convection is investigated with the doublepopulation lattice Boltzmann method. Compared with the isothermal case, the migration of the particle is notably affected by thermal convection. Five competing mechanisms are analyzed to illuminate the lateral migration of the particle with thermal convection, i.e., wall repulsion due to lubrication, inertial lift related to shear slip, a lift due to particle rotation, a lift due to the curvature of the undisturbed velocity profile and a lift induced by thermal convection. Within the parameter range examined (Gr < 2000; Rec < 200), the following are what we conclude:  due to thermal convection, the equilibrium position of the particle is asymmetric about the centerline and closer to the wall;

2168

J. Hu, Z. Guo / International Journal of Heat and Mass Transfer 108 (2017) 2158–2168

 the equilibrium position of the particle is dependent on the initial position for small Gr, however, the equilibrium position is insensitive to the initial position once beyond a critical Grashof number Gr c ;  a power law dependence of the critical Grashof number Gr c on the channel Reynolds number Rec is discovered;  with the channel Reynolds number Rec increasing, the forced flow relative to the particle is becoming stronger, resulting in a declining effect of thermal convection, i.e., the equilibrium position of the particle with thermal convection is closer to that without thermal convection. In the present study, the two-dimensional case is analyzed and discussed. And a similar three-dimensional case can be modeled as well, which we will consider carefully in the near future. Acknowledgments This work is subsidized by the National Natural Science Foundation of China (51125024 and 51390494). References [1] B.K. Cook, D.R. Noble, J.R. Williams, A direct simulation method for particlefluid systems, Eng. Comput. 21 (2004) 151–168. [2] Z.G. Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulations of particulate flows, J. Comput. Phys. 202 (2005) 20–51. [3] G. Segré, A. Silberberg, Radial particle displacements in Poiseuille flow of suspensions, Nature 189 (1961) 209–210. [4] A. Karnis, H.L. Goldsmith, S.G. Mason, The flow of suspensions through tubes: V. Inertial effects, Can. J. Chem. Eng. 44 (1966) 181–193. [5] M. Tachibana, On the behaviour of a sphere in the laminar tube flows, Rheol. Acta 12 (1973) 58–69. [6] P.G. Saffman, The lift on a small sphere in a slow shear flow, J. Fluid Mech. 22 (1965) 385–400. [7] B.P. Ho, L.G. Leal, Inertial migration of rigid spheres in two-dimensional unidirectional flows, J. Fluid Mech. 65 (1974) 365–400. [8] P. Vasseur, R.G. Cox, The lateral migration of a spherical particle in twodimensional shear flows, J. Fluid Mech. 78 (1976) 385–413. [9] J.A. Schonberg, E.J. Hinch, Inertial migration of a sphere in Poiseuille flow, J. Fluid Mech. 203 (1989) 517–524. [10] J. Feng, H.H. Hu, D.D. Joseph, Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2: Couette and Poiseuille flows, J. Fluid Mech. 277 (1994) 271–301. [11] T. Inamuro, K. Maeba, F. Ogino, Flow between parallel walls containing the lines of neutrally buoyant circular cylinders, Int. J. Multiphase Flow 26 (2000) 1981–2004. [12] H. Basßag˘aog˘lu, P. Meakin, S. Succi, G.R. Redden, T.R. Ginn, Two-dimensional lattice Boltzmann simulation of colloid migration in rough-walled narrow flow channels, Phys. Rev. E 77 (2008) 031405. [13] H. Basßag˘aog˘lu, S. Allwein, S. Succi, H. Dixon, J.T. Carrola, S. Stothoff, Two-and three-dimensional lattice Boltzmann simulations of particle migration in microchannels, Microfluid. Nanofluid. 15 (2013) 785–796. [14] N.A. Patankar, P.Y. Huang, T. Ko, D.D. Joseph, Lift-off of a single particle in Newtonian and viscoelastic fluids by direct numerical simulation, J. Fluid Mech. 438 (2001) 67–100.

[15] H. Gan, J. Chang, J.J. Feng, H.H. Hu, Direct numerical simulation of the sedimentation of solid particles with thermal convection, J. Fluid Mech. 481 (2003) 385–411. [16] Z. Yu, X. Shao, A. Wachs, A fictitious domain method for particulate flows with heat transfer, J. Comput. Phys. 217 (2006) 424–452. [17] Z.G. Feng, E.E. Michaelides, Heat transfer in particulate flows with direct numerical simulation (DNS), Int. J. Heat Mass Transfer 52 (2009) 777–786. [18] S.K. Kang, Y.A. Hassan, A direct-forcing immersed boundary method for the thermal lattice Boltzmann method, Comput. Fluids 49 (2011) 36–45. [19] D.Y. Kim, H.S. Yoon, M.Y. Ha, Particle behavior in a vertical channel with thermal convection in the low Grashof number regime, Comput. Fluids 48 (2011) 183–191. [20] Z.G. Feng, S.G. Musong, Direct numerical simulation of heat and mass transfer of spheres in a fluidized bed, Powder Technol. 262 (2014) 62–70. [21] A. Eshghinejadfard, D. Thévenin, Numerical simulation of heat transfer in particulate flows using a thermal immersed boundary lattice Boltzmann method, Int. J. Heat Fluid Flow 60 (2016) 31–46. [22] B. Yang, S. Chen, C. Cao, Z. Liu, C. Zheng, Lattice Boltzmann simulation of two cold particles settling in Newtonian fluid with thermal convection, Int. J. Heat Mass Transfer 93 (2016) 477–490. [23] L. Wang, Z. Guo, B. Shi, C. Zheng, Evaluation of three lattice Boltzmann models for particulate flows, Commun. Comput. Phys. 13 (2013) 1151–1172. [24] S. Tao, J. Hu, Z. Guo, An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows, Comput. Fluids 133 (2016) 1–14. [25] Z. Guo, C. Zheng, B. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Phys. 11 (2002) 366–374. [26] Z. Guo, B. Shi, C. Zheng, A coupled lattice BGK model for the Boussinesq equations, Int. J. Numer. Methods Fluids 39 (2002) 325–342. [27] P. Lallemand, L. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546. [28] H. Qian Y, D. D’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479–484. [29] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmannlattice fluid with boundaries, Phys. Fluids 13 (2001) 3452–3459. [30] C. Pan, L. Luo, C. Miller, An evaluation of lattice Boltzmann schemes for porous medium flow simulation, Comput. Fluids 35 (2006) 898–909. [31] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308. [32] Z. Guo, C. Zheng, B. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (2002) 2007–2010. [33] H. Li, H. Fang, Z. Lin, S. Xu, S. Chen, Lattice Boltzmann simulation on particle suspensions in a two-dimensional symmetric stenotic artery, Phys. Rev. E 69 (2004) 031919. [34] J. Hu, S. Tao, Z. Guo, An efficient unified iterative scheme for moving boundaries in lattice Boltzmann method, Comput. Fluids 144 (2017) 34–43. [35] C. Peng, Y. Teng, B. Hwang, Z. Guo, L. Wang, Implementation issues and benchmarking of lattice Boltzmann method for moving rigid particle simulations in a viscous flow, Comput. Math. Appl. 72 (2016) 349–374. [36] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1: Theoretical foundation, J. Fluid Mech. 271 (1994) 285–309. [37] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2: Numerical results, J. Fluid Mech. 271 (1994) 311– 339. [38] B. Wen, C. Zhang, Y. Tu, C. Wang, H. Fang, Galilean invariant fluidsolid interfacial dynamics in lattice Boltzmann simulations, J. Comput. Phys. 266 (2014) 161–170. [39] H. Li, X. Lu, H. Fang, Y. Qian, Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimensions, Phys. Rev. E 70 (2004) 026701.