Simulation of solid particles behaviour in a driven cavity flow

Powder Technology 191 (2009) 327–339

Contents lists available at ScienceDirect

Powder Technology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p ow t e c

Simulation of solid particles behaviour in a driven cavity flow P. Kosinski a,⁎, A. Kosinska b, A.C. Hoffmann a a b

The University of Bergen, Norway Bergen University College, Norway

a r t i c l e

i n f o

Article history: Received 2 August 2007 Received in revised form 28 May 2008 Accepted 28 October 2008 Available online 12 November 2008 Keywords: Driven cavity flow Numerical simulation Two-phase flows

a b s t r a c t The paper focuses on simulation of driven flow in a cavity containing an incompressible fluid and solid particles. The particle phase was modelled using the Eulerian–Lagrangian (E–L) approach where the solid particles are treated as points moving in the computational domain as a result of the fluid motion. In the mathematical model, the particle–particle interactions are simulated using the hard-sphere model. Different cases were considered, where the Reynolds number of the flow (defined by the velocity of the moving boundary, size of the cavity and fluid viscosity) and particle momentum response time were varied. This issue has various applications, for instance in pneumatic transport of dusts. The results are shown as snapshots of particle location at specific points in time, as well as statistics to characterize the behaviour in the dust cloud. The two-way coupling was considered making it possible to determine the influence of the particle cloud on the fluid phase in a driven cavity flow. © 2008 Elsevier B.V. All rights reserved.

1. Introduction 1.1. Driven cavity Rotating flows of viscous fluids have various industrial applications. The main focus of researchers has been lid-driven cavity flows, where the fluid is set into motion by a part of the containing boundary. This type of flows is crucial for analysing fundamental aspects of recirculating fluids: in spite of the apparently simple geometry, liddriven cavity flows may involve high degree of complexity. Shankar and Deshpande [1] give a comprehensive review on the experimental and numerical studies related to this type of flows. They emphasize the following advantages of studying driven-cavity flows: 1. results from theory, numerical simulations and experiments can be directly compared 2. the flow domain is unchanged when the Reynolds number is varied 3. this class of flows exhibit almost all phenomena that may occur in incompressible flows: eddies, secondary flows, complex 3-D patterns, instabilities and transition to turbulence Most of the papers related to the analysis of the flow in a driven cavity may be divided roughly into two groups: 1. the analysis of the physics of the flow 2. the development of new numerical schemes (where the driven cavity is used as a benchmark)

⁎ Corresponding author. E-mail address: [email protected] (P. Kosinski). 0032-5910/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2008.10.025

Taking the first of these two objectives, there are many works devoted to the issue of transition to turbulence and regarding the location of the so-called first Hopf bifurcation (e.g. Bruneau and Saad [2], Tiesinga et al. [3]). Some work has been done on rectangular rather than square cavities. Cheng and Hung [4] use a lattice Boltzmann method for the analysis of the 2D vortex structure in a rectangular lid-driven cavity for different depth-to-width ratios and for various Reynolds numbers. Also other geometries have been tested. Povitsky [5] studies the flow in a 3D cavity, where the lid is moving along the diagonal. Sheu and Tsai [6] present an investigation of 3D flow in a cubical cavity. As mentioned, many researchers use driven cavity flows for validating their numerical schemes. Some examples are: Perron et al. [7], Elman et al. [8], Zhang [9], Auteri et al. [10], Mei et al. [11], Wright and Smith [12], Liu and Leung [13], Boivin et al. [14], Neofytou [15], Tai et al. [16], Albensdoer and Kuhlmann [17]. There have been some works devoted to the issue of motion of two walls in the cavity rather than one. An example is the work of Albensdoer and Kuhlmann [18], where the linear stability of two counter-rotating vortices driven by the parallel motion of the two walls is investigated by a finite volume method. Other works are: [19–23]. Flow with temperature effects in a driven cavity has also been the subject of some interest: a temperature gradient is introduced and this leads to buoyancy driven flow. The field becomes coupled by velocity and temperature distributions. Examples are: Agrawal et al. [24], Zhou et al. [25], Alleborn et al. [26]. An interesting issue is the use of experimental techniques for analysis of driven cavity flows: Migeon et al. in [27] and in [28] investigate the time-dependent laminar flow in cavities where one wall moves after a sudden start. Ilegbusi and Mat [29] analyse mixing

328

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

of two fluids of different densities in a 2D cavity, Vogel et al. [30] study experimentally the flow in a rectangular cavity driven by the sinusoidal motion of the floor over a broad range of parameters. 1.2. Two-phase flows

as the analysis of the physical processes present there. We investigate two aspects: 1. the behaviour of the particles 2. their influence on the fluid

Wall-bounded fluid-particle flows are present in many industrially relevant systems: fluidized beds, pneumatic or slurry transport of powders, coal combustion, dust explosions, catalytic reactions, stirred tanks containing particles and many others. Of interest is the determination of particle–turbulence interaction and particle–particle interactions and both have been studied in literature. There are hundreds of works related to the interesting and still not fully resolved issue. Here we show just a few ones showing a general trend. White [31] performed analytical studies of the equilibrium of a particle lying on a bed of particles already in 1940. He analysed the problem theoretically to determine which factors are of importance. Many years later Tsuji et al. [32] performed numerical simulations to determine the motion of particles and fluid in a horizontal channel. Sommerfeld [33] demonstrated that numerical simulations of confined particulate two-phase flows require detailed modelling of particle–wall collisions. Frank et al. [34] simulated the motion of solid particles in a horizontal, twodimensional, turbulent channel flow, based on the Lagrangian approach. Kuru et al. [35] studied the mechanisms responsible for the initial growth of sand waves on the surface of a settled layer of particles experimentally and theoretically. Doron and Barnea [36] obtained data of pressure drop and flow patterns in a pipe with a flow of solid–liquid mixture experimentally. According to Sommerfeld and Huber [37] the wall roughness considerably alters the rebound behaviour of the particles and causes in average a redispersion of the particles. Hayden et al. [38] investigated particle entrainment experimentally. They measured the pick-up velocity for particles of different diameter. Portela and Oliemans [39] analysed dispersion of particles and their deposition and re-suspension at the walls. Their main interest was the particle–turbulence interaction.

2. Mathematical model

1.3. The objective of this research

In this research, the flow of the particles was modelled using the Eulerian–Lagrangian (E–L) approach where the particles are treated as points, whose motion is the result of the influence of the fluid phase. The biggest drawback of this method is the fact that in real industrial applications the number of the particles is too large to consider the behaviour of all of them. In order to overcome this problem, a number of particles are often “pooled” to form “virtual particles” with the same parameters, so that the number of the particles to be considered can be reduced. For small geometries with limited number of particles, the E–L type of model can be used accurately without introducing virtual particles. This was done in the present research. The E–L approach makes it possible to investigate the fundamental processes. The equations of motion for the jth particle may be written in the following way:

Our paper studies the flow in a square 2D cavity filled with solid particles. This is an interesting problem, which may yield much information about the interaction between fluid and particle in a wide range of practical configuration. This has not been widely studied so far. As mentioned above, the rotating flows of viscous fluids have are of interest for many researchers and engineers. This also refers to the case when there is also a particle phase present there. Here we can name many industrial applications, starting with mixing and segregation, and finally going to more scientific issues like the analysis of eddies in turbulent flow. This is exactly what happens while studying a flow in a driven cavity. An important problem while dealing with two-phase flows is the proper choice of a good numerical model. Again the driven cavity flow is a proper tool: the flow with a pure fluid has been precisely investigated by many researchers (benchmarks) and besides the physics of the flow involves problems present virtually also in many other types of flows. A numerical scheme that works correctly for simulating a flow in a driven cavity, is supposed to work also for modelling other applications. An example, shown further in the paper, are the particles that move on spiral trajectories due to the rotation of the mixture. Insufficient precision in the particle algorithm will result in migration of the particles out of their real trajectories. This also happens, for instance, in a turbulent flow in a pipeline. Thus, the objective of the paper is both to name a numerical scheme that one can use for simulating other two-phase flows, as well

2.1. The flow of the continuous phase The mathematical model for the continuous phase (fluid) is based on the assumption that the volume fraction of the solid phase (particles) is low so that the effect of the particles on the equation of continuity is negligible. Thus the flow of the continuous phase is modelled using standard Navier–Stokes equations modified to take into account momentum exchange with particles: • the equation of continuity Y

jd u = 0

ð1Þ

• the equation of momentum ρ


Y  Au Y Y Y Y + udðju Þ = −jp + μj2u−f p At

ð2Þ

Y

where: f p expresses the action of the particles on unit volume of Y fluid; u is the fluid velocity; p is the pressure; ρ is the fluid density, μ is the dynamic viscosity. The momentum equation, Eq. (2) differs from standard Navier– Y Stokes equation by one extra term, fp , which is modelled as: Y

Y

f p = ∑ fpj

ð3Þ

Y

where: fpj represents the action of each individual particle [40]. 2.2. The flow of the particles

Y

mj

dv j Y = f pj dt

ð4Þ

Y

Ij

dω j Y = T pj dt

ð5Þ Y

where: mj is the particle mass; v j is its velocity, Ij its moment of inertia, Y Y ω j its angular velocity and T pj its torque. Y The force f pj , which is also used in Eq. (2), describes all forces acting on jth particle from the surrounding fluid. This is mainly the drag force that may be expressed using the fundamental formula: Y Y  Y Y ju−v j j u−v j Y Y ð6Þ f pj = f drag = CD Ap ρ 2

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

329

modelled since the average distance between the particles was high enough to neglect this. 2.3. Numerical scheme

Fig. 1. The grid used in simulation.

where: CD is the coefficient of the drag force and Ap is the area projected on a plane normal to the flow direction. The value of the drag coefficient CD is a primarily a function of the Reynolds number, defined as: Y

Re =

Y

ρdju−v j j μ

ð7Þ

where: d is the particle diameter, μ is the viscosity One of the empirical formulas that can be used is (by Schiller and Neumann [40]): CD =

  24 1 1 + Re2=3 Re 6

ð8Þ

In this research, we only considered the drag force acting on a particle. There are also other mechanisms like buoyancy, lift and added-mass forces. They can be easily implemented into the code but their effect was negligible in this research. Y The torque Tpj can be modelled using the following formula:  5 ρ d Y Y T pj = CR jΩjΩ 2 2 Y

ð9Þ

The value of the coefficient CR, used in Eq. (9), is primarily a function of the rotational Reynolds number, defined as: Y

ρd2 j Ωj μ

ð10Þ

1 Y j× u−ωj 2

ð11Þ

ReR =

When solving the Navier–Stokes equations, the computational domain was discretized using a staggered grid such that the pressure, p, was located in the centres of the computational cells, whereas the xcomponent of fluid velocity was stored in the midpoints of the vertical cell edges. Similarly, the y-component of fluid velocity was located in the midpoints of the horizontal cell edges (see Fig. 1). The general numerical scheme resembles the Chorin projection method developed by Chorin and Temam (see [42]) for a flow of a pure fluid. We modified it by use of the Adams–Bashfort time-advancement scheme, instead of the Euler scheme, in order to increase the accuracy. The reason for this modification was that we simulated transient flow rather than stationary, making a better time discretization necessary. In the addition to this we included the model for the particle phase in the algorithm. The use of explicit schemes does not impose any problems, since the particle relaxation time is greater than the time step imposed by the numerical stability of the flow equations (see also [43]). We explain the methods used for each step in the algorithm (for clarity, we show the Euler time–advancement scheme, which is easier to explain): 1. Stability condition: calculation of the time step, Δt. The current time step is found taking into account standard Courant–Friedrichs–Levy (CFL) conditions. 2. The predictor step: calculation of the tentative velocity of fluid. The equation of continuity is not enforced here yet. Gn−1 i;j = Gni;j =

1 2 Y ðn−1Þ Y ðn−1Þ  Y ðn−1Þ j u i;j −u i;j d ju i;j Re

1 2 Y ðnÞ Y ðnÞ  Y ðnÞ j u i;j −u i;j d ju i;j Re

  3 n 1 n−1 Yn ~ Gi;j − Gi;j u i;j = u i;j + Δt 2 2

Y

The value of the coefficient [40]: 8 64 > > for ReR b32 < ReR CR = 12:4 128:4 > > for 32bReR b1000 : 0:5 + ReR ReR

ð12Þ

Collisions of particles with walls were modelled using the hardsphere approach [40] and a model based on seeking so-called time of collision (see e.g. [41]). Collisions between the particles were not

ð14Þ

ð15Þ

For discretization of the convective terms, we compared different schemes: First-Order Upwind (FOU), Quadratic Upstream Interpolation for Convective Kinematics (QUICK) [44], Sharp and Monotonic Algorithm for Realistic Transport (SMART) [45] and Variable-Order Non-Oscillatory Scheme (VONOS) [46]. We chose the last one: VONOS since it led to results closest to benchmarks. Similar conclusions were drawn by Varonos and Bergeles [46]. This scheme is also the most time-consuming when simulating a

with: Ω=

ð13Þ

Fig. 2. The description of the problem.

330

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

Fig. 3. Snapshots of particle position for four points in time: (a) 501.25, (b) 502.5, (c) 503.75, (d) 505.0. St = 0.02778.

flow in a given mesh. However, its accuracy makes it possible to use coarser grids, thus leading to over-all shorter simulation time. 3. The calculation of particle motion and updating of fluid velocities due to two-way coupling. The objective is to solve the following equations (for each particle): Yn

mj

v i;j

+1

Yn

−v j

Δt

Y

Y

= f ðu; v Þ

ð16Þ

and: n+1

Y ω i;j

n

Y −ω i;j

Δt

Y

= f ðω Þ

ð17Þ

Eq. (16) was solved using the 4th order Runge–Kutta scheme (with the same time step as for the fluid flow). We also tried other

techniques, namely a 1st order semi-implicit scheme (where the drag force coefficient, CD, was treated explicitly), the Adams scheme, as well as the Runge–Kutta–Fehlberg method (4th/5th order). The first two are not accurate enough, whereas the last one did not lead to better results than the standard Runge–Kutta method. This was tested by investigating of the influence of the time step value. In order to evaluate the value of the drag force acting on the particles, the fluid velocity in the particle location has to be evaluated as well. This was done by a linear interpolation between the grid points, where the fluid velocity is known. Eq. (17) was solved using an implicit scheme. The final position of a particle was found by: Yn

ri

+1

Yn

= r i + Δt

 Yn + 1 Yn v +v 2

ð18Þ

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

4. The Poisson procedure and enforcing of the equation of continuity: calculation of pressure. At this stage the momentum equation, without the convective terms, is considered: Yn

u i;j

+1

~ − ui;j

Δt

= −jpni;j+ 1

ð19Þ

331

6. Boundary conditions for the final velocities. Since the velocities should vanish at solid boundaries, we apply the non-slip conditions for the lower, right and left boundaries. The x-component of the velocity at the upper boundary is equal to the given lid velocity. 3. Numerical experiments

with the equation of continuity jdY ui;j = 0

ð20Þ

Application to the divergence operator to Eq. (19) leads to the Poisson equation for the pressure: Δpni;j+ 1 =

1 ~ jd ui;j Δt

ð21Þ

5. The corrector procedure: calculation of the updated fluid velocity by use of Eq. (19)

3.1. Initial data Two-dimensional flows cannot be realized experimentally, but numerical simulation 2-D of driven cavity flows offers considerable simplifications compared to simulating 3-D flows and this can lead to better understanding of some issues. This is the case in this paper, where the issue is to analyse the behaviour of particle cloud suspended in a rotating viscous flow. If B is a length scale equal to the edge length of the cavity, and U is a speed scale equal to the speed of the moving boundary, then all the lengths

Fig. 4. Snapshots of particle position for four points in time: (a) 501.25, (b) 502.5, (c) 503.75, (d) 505.0. St = 2.778.

332

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

and velocities may be normalized using a Reynolds number: Re=UB/v. In this research the value was always the same and equal to 1000. The first stage was the comparison of the code performance for simulation of pure flow (without particles) with benchmarks. Various Reynolds numbers, different grid sizes and other parameters were varied to obtain optimal parameters for, on the one hand, fit with the data and, on the other, computational effort. A grid of 64 per 64 cells was chosen as satisfactory for our applications. Subsequently, testing for two-phase flow was performed in a similar way: at first with one particle only tracked in the precalculated flow field. Here the necessity of using a lower value of the CFL coefficient (compared to pure fluid flow) was observed. Finally this coefficient was chosen as 0.125. As mentioned above, the fourth order numerical scheme was crucial here: for lower order schemes the results were not satisfactory. The residual error, ε, which is present after iteration for each time step for computation of the Poisson equation was taken as 10− 5. The maximum number of iterations per time step was limited to 700. We did not observe any crucial influence of the presence of the particles on the convergence of the numerical scheme. The simulations were performed as follows: at first the computation was started without particles until a steady flow was achieved. This was assumed to occur when the difference between the average speeds of the flow in all the computational cells in two time steps was smaller than 2 × 10− 6: Δuave =

1 Yn + 1 Yn ∑ju −u j N

ð22Þ

where: N is the total number of computational cells in the domain. According to this criterion, steady flow occurred after t = 500 and at this point in time the particles were added. At the next time step the square domain was filled with particles in the way shown in Fig. 2. The distance between the particles was high enough so that it could be assumed that the particles did not interact with each other. It must be stressed, however, that locally regions of high particle concentration may occur and this assumption may lead to less correct results in these regions. However, we do not anticipate that this will influence the final results crucially. Different initial positions of the particles were tested in order to investigate how this affects the behaviour of the particles and the fluid. This was carried out as follows: the particles were introduced in regular arrays with a subsequent small displacement in a random direction. Thus every simulation was actually run with slightly different initial conditions. Initially the particles occupied the domain x∈(0.25,0.75)∧y∈(0.25,0.75). The reason for not filling the whole domain was to avoid collisions with the walls occurring at the very beginning. Testing showed that the initial distribution of the particles did not crucially influence the statistics shown below; the actual positions of the particles were only slightly different for different simulations. All together 14 cases were considered, where the following parameters were varied:

In this research two-way coupling was used for modelling, where not only the influence of the fluid phase on the particles' motion is considered but also the influence of the particles on the velocity of the fluid phase. To compute the latter, the particle size relative to the grid cell size is important. For one-way coupling, on the other hand, where only the fluid phase influences the particle motion, the size of the particles with respect to the mesh cell size is of no importance: the particles do not alter the parameters of the fluid and thus even very coarse particles can be modelled (see e.g. [43]) Since we used the two-way coupling approach, we had to be careful with particle size relative to grid size. This was chosen as: d/Δ = 0.0064 for grid 64 per 64, where Δ is the grid cell size. The density of the particles was 1000 times higher than the density of the fluid phase (heavy particles). The restitution coefficient and the Coulomb friction coefficient, necessary for modelling of collisions with the walls was 0.9 and 0.15, respectively. 3.2. Results 3.2.1. Particle motion At first we show how the process develops in order to better illustrate and understand the issue. Fig. 3 shows snapshots of particle position for four points in time: 501.25, 502.5, 505.0 and 507.5. The Stokes number was 0.02778. This figure shows how the rotating fluid brings the particles into motion and how the particles are pushed against the upper (moving) and left-hand walls of the cavity. When the Stokes number increases to 2.778 the process is similar (see Fig. 4). However, due to higher inertia the particles are centrifuged outward earlier than for the lower value of Stokes number. Comparing the (a) plates in Figs. 3 and 4 it can be seen that the particles ever at the higher value of the Stokes number, are put into motion by the fluid, but do not react as promptly to the rotation of the fluid as at lower Stokes number. In order to better quantify the behaviour of the particles, the following statistics was performed: so-called mean square displacement (MSD) for time t of all the particles was computed: MSDðt Þ =

i 1 h ∑ ðxi ðt Þ−xoi Þ2 + ðyi ðt Þ−yoi Þ2 N

where for ith particle: xi(t), yi(t) is the position at time t and xoi, yo is the initial position.

• the Stokes number, defined as: St =

τv U ρp d2 U = B 18μ B

ð23Þ

where: τv is the particle momentum response time and ρp is particle density.The following seven values of the Stokes number were considered in simulations: 0.02788, 0.2778, 0.6944, 1.3889, 2.0835 or 2.778. Additionally a case when no particles were present at all was tested as well. • Relative distance between the particles: L/d where: L is the distance between the particles centres The following two values of the relative distance were considered: L/d = 25 or 50. This corresponded to the total number of particles: 38,809 and 9801, respectively.

ð24Þ

Fig. 5. The history of mean square displacement.

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

333

Fig. 6. History fluid velocity at point (0.25, 0.25) after the particles have been introduced for (a) L/d = 25, x-component; (b) L/d = 50, x-component; (c) L/d = 25, y-component; (d) L/d = 50, y-component.

We show the results only for the initial relative distance L/d = 25, since for the other value, 50, the mean square displacement was nearly the same. Histories of the mean square displacement are shown in Fig. 5 for different values of the Stokes number. The rate of the increase is

higher for smaller values of the Stokes number, which is caused by the fact that the particles have smaller momentum response time and react faster to the fluid flow. When the particle cloud makes one half of the full rotation, the mean square displacement obtains its maximum value and later

Fig. 7. Change in velocity at point (0.25, 0.25) after the particles has been introduced, where: (a) velocity x-component, (b) velocity y-component.

334

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

Fig. 8. Change in velocity after the particles has been introduced, where: (a) velocity x-component, (b) velocity y-component.

decreases. This is especially visible for the Stokes number equal to 0.02778, where this reaches nearly zero. This local minimum corresponds to one full rotation of the particle cloud. Afterwards, the mean square displacement increases again. However, the local minima will at higher values since the particles migrate towards the wall so that the average distance from their initial position at full turns steadily increases. At higher values of the Stokes number the mean square displacement quickly increases. As mentioned above, this is due to the higher inertia. 3.2.2. Velocity change In the model, so-called two-way coupling was implemented, which means that particles may influence the flow of fluid. This influence is shown in Fig. 6 as the history of x- and y-component of fluid velocity

monitored at coordinates (0.25, 0.25). Before the particles were introduced (time 500 on the graph) the velocity was nearly constant. After the particles have been distributed in the cavity, the velocity started to vary. This is more intense for higher Stokes numbers and for smaller distance between the particles (higher concentrations). In order to quantify this better, we show the change of the fluid velocity due to the particles at point with coordinates (0.25, 0.25). This is presented in Fig. 7, where the influence of the Stokes number and the relative distance L/d are analysed. The point in time is 505. The increase in the relative distance influences the velocity in a more intense way, as seen in Fig. 7. This discussion refers to only one monitor point. In Fig. 8 we show two contour plots with the distribution of the change of the velocity after the particles have been introduced. Fig. 8a corresponds to the

Table 1 Vorticity after t = 502.5, L/d = 50, x = 0.5 y

St = 2.778

2.0835

1.3889

0.6944

0.2778

0.02788

No particles

0.125 0.250 0.375 0.500 0.625 0.750 0.875

−1.2005869 2.4231493 2.0556857 2.0066930 2.0044907 2.1099805 1.5522362

−1.2004307 2.4232284 2.0559721 2.0070184 2.0044989 2.1101444 1.5523280

−1.2002590 2.4232982 2.0563327 2.0074059 2.0049126 2.1104372 1.5526645

−1.2008743 2.4232256 2.0567243 2.0075265 2.0053163 2.1108572 1.5520807

−1.2010152 2.4230531 2.0571006 2.0082365 2.0059570 2.1114345 1.5514413

−1.2011228 2.4230630 2.0571340 2.0087603 2.0064175 2.1118911 1.5512882

−1.2011603 2.4230530 2.0571154 2.0088180 2.0064642 2.1120477 1.5512440

Table 2 Vorticity after t = 502.5, L/d = 25, x = 0.5 y

St = 2.778

2.0835

1.3889

0.6944

0.2778

0.02788

No particles

0.125 0.250 0.375 0.500 0.625 0.750 0.875

− 1.1978681 2.4239834 2.0520025 2.0017470 1.9989109 2.1049760 1.5551696

−1.1979032 2.4238944 2.0527709 2.0018957 1.9995036 2.1056488 1.5554178

− 1.1976680 2.4241781 2.0536898 2.0031996 2.0003613 2.1065125 1.5561843

−1.2001474 2.4238016 2.0554129 2.0052342 2.0023682 2.1081790 1.5538114

− 1.2005945 2.4231185 2.0569564 2.0068139 2.0043849 2.1100537 1.5527497

− 1.2010815 2.4230689 2.0571313 2.0086273 2.0062497 2.1117712 1.5514082

−1.2011603 2.4230530 2.0571154 2.0088180 2.0064642 2.1120477 1.5512440

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

335

Fig. 9. Histories of enstrophy for (a) L/d = 25 and (b) L/d = 50.

Stokes number 1.3889, and Fig. 8b to the Stokes number 2.778. Both figures are for the value of the L/d = 25. 3.2.3. Vorticity profiles Vorticity can be defined as the curl of the fluid velocity: Y

Y

ω = j×u

ð25Þ

and for 2D situation has only one component (in the z-direction): Auy Auy − ω= Ax Ay

3.2.4. Enstrophy and total kinetic energy Bruneau and Saad [2] suggest using global quantities for comparison for analysis of flow processes in a driven cavity. Among other things they use the enstrophy and the total kinetic energy of the system defined as: Enstrophy: Z=

ð26Þ

This parameter has been widely used for the flow analysis within a driven cavity. The presence of particles will influence this value. Since the process is not steady the vorticity will also be a function of time. Profiles of vorticity along a vertical line going through the cavity centre are presented in Tables 1 and 2. They correspond to various values of the Stokes number and the relative distance L/d.

1 ∫Ω jjωjj2 dΩ 2

ð27Þ

Total kinetic energy: E=

1 ∫Ω jjujj2 dΩ 2

ð28Þ

The values are stationary before the particles are introduced into the cavity. After this, they begin to vary in time. In Figs. 9 and 10 the histories of these parameters are shown for different cases.

Fig. 10. Histories of total kinetic energy for (a) L/d = 25 and (b) L/d = 50.

336

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

Fig. 11. The change of the enstrophy (a) and the total kinetic energy (b) due to particles.

The particles influence these global quantities. How the Stokes number and the relative distance between the particles influence the change of enstrophy and the total kinetic energy is shown in Fig. 11. The first conclusion is that, as expected, an increase of the Stokes number and decrease in the relative initial distance between the particles will increase the change of both enstrophy and the total kinetic energy of the system. In the simulations, when the relative distance, L/d, decreased by a factor of 2, the total number of particles in the system increased 3.9 times. It can be noticed in Fig. 11 that when L/d increases by a factor of two, the change in enstrophy increases about 3–4 times and the change in total kinetic energy about 4 times. This means that there is a nearly linear relation, especially for the total kinetic energy. The influence of the Stokes number was, however, far from linear. When the Stokes number increased 10 times, the change of enstrophy increased 6.89 times for L/d = 25, and 6.25 for L/d = 50.

Similarly, when the Stokes number increased 10 times, the change of the total kinetic energy increased 5.4 times for L/d = 25, and 5.0 for L/d = 50. 4. Comparison with experiments We now compare the performance our code with experiments shown in the paper by Tsorng et al. [47]. This paper discusses to a 3D measurement of a particle suspended in a fluid enclosed in a transparent cubic cavity. All the initial data, presented below, are dimensional as shown in their paper. The experimental setup consisted of a cavity with side 10 cm, filled with a fluid of viscosity 37.2 mm2/s and specific gravity 1.21. At the top, a conveyer belt system was installed moving with a speed 17.5 cm/s. This data corresponded to the Reynolds number 470. The density of the particle was nearly equal to the density of the fluid so that the

Fig. 12. Trajectory of one particle in the driven cavity. Comparison of numerical simulation (a) with experiments (b) (from [47], Fig. 11a, with kind permission of Springer Science and Business Media).

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

neutral buoyancy was obtained. The same parameters were included into our code. Before the experiments, the particle was plunged into the cavity so that it touched the lid and its initial location was approximately at the middle of the lid. Fig. 12 shows the comparison of numerical simulation and the experiments as a trajectory of the particle. The total simulation (experiment) time was 25 s, whereas the trajectory sketched as a bold line corresponds to total time 8 s.

337

Fig. 13 shows how the process develops with emphasis on the fluid flow as well as the particle behaviour caused by this flow. Snapshots of the flow field are shown for six different points in time: 0.75, 1.00, 1.25, 2.00, 3.00 and 3.75 s. This illustrates the behaviour of the particle: it can be observed that the transient structure of the flow at the beginning of the simulation is of importance. The formation of the vortex clearly seen in Fig. 13 is responsible for the drift of the particle. This explains the shape of the particle trajectory shown in Fig. 12, a shape which is obtained both numerically and experimentally.

Fig. 13. Location of one particle in the field flow for various points in time: a) 0.75 s; b) 1.00 s; c) 1.25 s; d) 2.00 s; e) 3.00 s; f) 3.75 s.

338

P. Kosinski et al. / Powder Technology 191 (2009) 327–339

Until time reaches 8 s, both numerical results and the experiments match very well. Nevertheless, after this time, some differences appear that may be a result of the following: • The actual flow is 3D and these three-dimensional effects do not make it possible to simulate the real flow for a long period of time (where 2D modelling is assumed) • It is only possible to estimate the approximate initial location of the particle in the experiments • There is some lubrication flow that occurs through the small gaps between the conveyor belt and the top of the cavity. This may slightly influence the flow pattern in the cavity. We suppose that the first issue is of main importance, since Tsorng et al. report a motion of the particle in the lateral direction after a few seconds after the beginning of the process. This cannot be simulated using 2D assumption. A similar conclusion can be drawn while analysing for instance reference [48], where studies of rotating filtration are reported. The authors investigate Taylor–Couette flow involving a rotation inner cylinder and a stationary outer cylinder. This configuration resembles our research involving, as it does, a rotational flow with particles. The main difference between this reference and our research is the threedimensional flow leading to very interesting conclusions and results: the flow in the third direction is responsible for a special segregation of the particles. Tsorng et al. also present results from a case in which more than one particle were distributed in the system. These results are, however, difficult to compare directly with numerical observations. 5. Concluding remarks The paper focuses on simulation of two-phase flow in a driven cavity for the Reynolds number equal to 1000. Different cases were considered, where the Stokes number and the initial distance between the particles varied. It was shown how the particle cloud behaves in the system and simultaneously how the particles influence the flow of the fluid. The simulations show that the particles tend to migrate towards the walls and this process is more intense for particles with a higher value of the Stokes number. This was shown by presenting snapshots of particle position, as well as the mean square displacement of the particles. The particle cloud influences the flow. In this paper we showed that when the Stokes number increases and distance decreases, the change in fluid velocity is more clear. This was quantitatively described, but we must emphasize that more tests are needed to describe the process more thoroughly: for instance, it is interesting to compare different Reynolds numbers and to investigate the influence of particle–particle interactions. Also global parameters like enstrophy and the total kinetic energy were compared during various tests. The trends for the change of the parameters shown in Figs. 6–11 are small and it might be thought that they are within the numerical accuracy of the solver. However, we emphasize that the effects of the parameters tested, like, for example, the Stokes number are qualitatively correct. Thus we expect that numerical errors do not obscure the physics of the process. An interesting issue is the comparison with experiments. The model and the code are able to simulate a real flow, but at the same time the assumption that the flow is planar does not lead to accurate results for longer simulation times. 3D simulations are of importance for achieving better agreement for longer simulation times: this is, however, time-consuming because of the amount of particles that are present in the system. Acknowledgement This research was financed through the Hyperion project by the Research Council of Norway and Norsk Hydro.

References [1] P. Shankar, M. Deshpande, Fluid mechanics in the driven cavity, Annual Review of Fluid Mechanics 32 (2000) 93–136. [2] C.H. Bruneau, M. Saad, The 2D lid-driven cavity problem revisited, Computers & Fluids 35 (2006) 326–348. [3] G. Tiesinga, F.W. Wubs, A.E.P. Veldman, Bifurcation analysis of incompressible flow in a driven cavity by the Newton–Picard method, Journal of Computational and Applied Mathematics 140 (2002) 751–772. [4] M. Cheng, K.C. Hung, Vortex structure of steady flow in a rectangular cavity, Computers & Fluids 35 (2006) 1046–1062. [5] A. Povitsky, Three-dimensional flow in a cavity at yaw, Nonlinear Analysis 63 (2005) 1573–1584. [6] T.W.H. Sheu, S.F. Tsai, Flow topology in a steady three-dimensional lid-driven cavity, Computers & Fluids 31 (2002) 911–934. [7] S. Perron, S. Boivin, J.M. Herard, A finite volume method to solve the 3D Navier– Stokes equations on unstructured collocated meshes, Computers & Fluids 33 (2004) 1305–1333. [8] H.C. Elman, V. Howle, J.N. Shadid, R.S. Tuminaro, A parallel block multi-level preconditioner for the 3D incompressible Navier–Stokes equations, Journal of Computational Physics 187 (2003) 504–523. [9] J. Zhang, Numerical simulation of 2D square driven cavity using fourth-order compact finite difference schemes, Computers and Mathematics with Applications 45 (2003) 43–52. [10] F. Auteri, N. Parolini, L. Quarapelle, Numerical investigation on the stability of singular driven cavity flow, Journal of Computational Physics 183 (2002) 1–25. [11] R. Mei, W. Shyy, D. Yu, L.S. Luo, Lattice Boltzmann method for 3-D flows with curved boundary, Journal of Computational Physics 161 (2000) 680–699. [12] J.A. Wright, R.W. Smith, An edge-based method for incompressible Navier–Stokes equations on polygonal meshes, Journal of Computational Physics 169 (2001) 24–43. [13] C.H. Liu, D.Y.C. Leung, Development of a finite element solution for the unsteady Navier–Stokes equations using projection method and fractional-θ-scheme, Computer Methods in Applied Mechanics and Engineering 190 (2001) 4301–4317. [14] S. Boivin, F. Cayre, J.M. Herard, A finite volume method to solve the Navier–Stokes equations for incompressible flows on unstructured meshes, International Journal of Thermal Science 39 (2000) 806–825. [15] P. Neofytou, A 3rd order upwind finite volume method for generalised Newtonian fluid flows, Advances in Engineering Software 36 (2005) 664–680. [16] C.H. Tai, Y. Zhao, K.M. Liew, Parallel-multigrid computation of unsteady incompressible viscous flows using a matrix-free implicit method and high-resolution characteristics-based scheme, Computer Methods in Applied Mechanics and Engineering 194 (2005) 3949–3983. [17] S. Albensdoer, H.C. Kuhlmann, Accurate three-dimensional lid-driven cavity flow, Journal of Computational Physics 206 (2005) 536–558. [18] S. Albensdoer, H.C. Kuhlmann, Three-dimensional instability of two counterrotating vortices in a rectangular cavity driven by parallel wall motion, European Journal of Mechanics B/Fluids 21 (2002) 307–316. [19] H.C. Kuhlmann, M. Wanschura, H.J. Rath, Flow in two-sided lid-driven cavities: non-uniqueness, instabilities, and cellular structures, Journal Of Fluid Mechanics 336 (1997) 267–299. [20] H.C. Kuhlmann, M. Wanschura, H.J. Rath, Elliptic instability in two-sided lid-driven cavity flow, European Journal of Mechanics B/Fluids 17 (1998) 561–569. [21] C.H. Blohm, H.C. Kuhlmann, The two-sided lid-driven cavity: experiments on stationary and time-dependent flows, Journal Of Fluid Mechanics 450 (2002) 67–95. [22] S. Albensoeder, H.C. Kuhlmann, Linear stability of rectangular cavity flows driven by anti-parallel motion of two facing walls, Journal Of Fluid Mechanics 458 (2002) 153–180. [23] S. Albensoeder, H.C. Kuhlmann, H.J. Rath, Multiplicity of steady two-dimensional plows in two-sided lid-driven cavities, Theoretical And Computational Fluid Dynamics 14 (2001) 223–241. [24] L. Agrawal, J.C. Mandal, A.G. Marathe, Computations of laminar and turbulent mixed convection in a driven cavity using pseudo-compressibility approach, Computers & Fluids 30 (2001) 607–620. [25] Y. Zhou, R. Zhang, H. Staroselsky, H. Chen, Numerical simulation of laminar and turbulent buoyancy-driven flows using a lattice Boltzmann based algorithm, International Journal of Heat and Mass Transfer 47 (2004) 4869–4879. [26] N. Alleborn, H. Raszillier, F. Durst, Lid-driven cavity with heat and mass transport, International Journal of Heat and Mass Transfer 42 (1999) 833–853. [27] C. Migeon, A. Texier, G. Pineau, Effects of lid-driven cavity shape on the flow establishment phase, Journal of Fluids and Structures 14 (2000) 469–488. [28] C. Migeon, G. Pineau, A. Texier, Three-dimensionality development inside standard parallelepipedic lid-driven cavities at Re = 1000, Journal of Fluids and Structures 17 (2003) 717–738. [29] O.J. Ilegbusi, M.D. Mat, A comparison of predictions and measurements of kinematic mixing of two fluids in a 2D enclosure, Applied Mathematics Modelling 24 (2000) 199–213. [30] M.J. Vogel, A.H. Hirsa, J.M. Lopez, Spatio-temporal dynamics of a periodically driven cavity flow, Journal of Fluid Mechanics 478 (2003) 197–226. [31] C.M. White, The equlibrium of grains on the bed of a stream, Proceedings of the Royal Society of London. Mathematical physical sciences (1940) 332–338. [32] Y. Tsuji, Y. Morikawa, T. Tanaka, N. Nakatsukasa, M. Nakatani, Numerical simulation of gas–solid two-phase flow in a two-dimensional horizontal channel, International Journal of Multiphase Flow 13 (5) (1987) 671–684. [33] M. Sommerfeld, Modelling of particle–wall collisions in confined gas-particle flows, International Journal of Multiphase Flow 18 (1992) 905–926.

P. Kosinski et al. / Powder Technology 191 (2009) 327–339 [34] T. Frank, K.P. Schade, D. Petrak, Numerical simulation and experimental investigation of a gas-solid two-phase flow in a horizontal channel, International Journal of Multiphase Flow 19 (1993) 187–198. [35] W.C. Kuru, D.T. Leighton, M.J. McCready, Formation of waves on a horizontal erodible bed of particles, Journal of Multiphase Flow 21 (6) (1995) 1123–1140. [36] P. Doron, D. Barnea, Pressure drop and limit deposit velocity for solid–liquid flow in pipes, Chemical Engineering Science 50 (10) (1995) 1595–1604. [37] M. Sommerfeld, N. Huber, Experimental analysis and modelling of particle–wall collisions, International Journal of Multiphase Flow 25 (1999) 1457–1489. [38] K.S. Hayden, K. Park, J.S. Curtis, Effect of particle characteristics on particle pickup velocity, Powder Technology 131 (1) (2003) 7–14. [39] L.M. Portela, R.V.A. Oliemans, Eulerian–Lagrangian DNS/LES of particle–turbulence interactions in wall-bounded flows, International Journal for Numerical Methods in Fluids 26 (2003) 719–727. [40] C. Crowe, M. Sommerfeld, Y. Tsuji, Multiphase Flows with Droplets and Particles, CRC Press LLC, 1998. [41] S. Sundaram, L.R. Collins, Numerical considerations in simulating a turbulent suspension of finite-volume particles, Journal of Computational Physics 124 (1995) 337–350. [42] M. Griebel, T. Dornseifer, T. Neunhoeffer, Numerical Simulation in Fluid Dynamics. A Practical Introduction, Society for Industrial and Applied Mathematics, 1997.

339

[43] L.M. Portela, R.V.A. Oliemans, Eulerian–Lagrangian DNS/LES of particle–turbulence interactions in wall-bounded flows, International Journal for Numerical Methods in Fluids 43 (2003) 1045–1065. [44] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Computational Methods and Applications in Mechanical Engineering 19 (1979) 59–98. [45] P.H. Gaskell, A.K.C. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, International Journal for Numerical Methods in Fluids 8 (1988) 617–641. [46] A. Varonos, G. Bergeles, Development and assessment of a Variable-Order Nonoscillatory Scheme for convection term discretization, International Journal for Numerical Methods in Fluids 26 (1998) 1–16. [47] S.J. Tsorng, H. Capart, J.S. Lai, D.L. Young, Three-dimensional tracking of the long time trajectories of suspended particles in a lid-driven cavity flow, Experiments in Fluids 40 (2006) 314–328. [48] S.T. Wereley, R.M. Lueptow, Inertial particle motion in a Taylor Couette rotating filter, Physics of Fluids 11 (1999) 325–333.