Numerical simulation of circular particles migration in oscillatory Poiseuille flow

Accepted Manuscript

Numerical simulation of circular particles migration in oscillatory Poiseuille flow Wenjun Yuan , Jianqiang Deng PII: DOI: Reference:

S0045-7930(16)30313-9 10.1016/j.compfluid.2016.10.012 CAF 3300

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

29 January 2016 9 October 2016 13 October 2016

Please cite this article as: Wenjun Yuan , Jianqiang Deng , Numerical simulation of circular particles migration in oscillatory Poiseuille flow, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.10.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights The oscillation causes a strong modification of the eddy behind the particle.



The pressure distribution on the body with rotation behavior is revealed.



Distinct flow patterns during circular particles settling are summarized.



The Magnus lift influences the balance between walls and particles interplay.

AC

CE

PT

ED

M

AN US

CR IP T



ACCEPTED MANUSCRIPT

Numerical simulation of circular particles migration in oscillatory Poiseuille flow

CR IP T

Wenjun Yuan, Jianqiang Deng*

School of Chemical Engineering and Technology, Xi'an Jiaotong University,

AN US

Xi’an 710049, P. R. China

Abstract: The present work reports on the simulation of one and two circular particles migration in a two-dimensional channel subjected to oscillatory Poiseuille

M

flow using a finite element arbitrary Lagrangian–Eulerian method. The effects of

ED

oscillation frequency on terminal angular velocity, drag force, motion trajectory, vortex structure and relative displacement during sedimentation have been

PT

investigated. The results show that the oscillation can cause a strong modification of

CE

the eddy formation and shedding behind the particle, and there is a critical oscillation

AC

frequency for the flow pattern of particles motion related to the blockage ratio. When the frequency is below the critical value, the single particle increases its rotation speed with alternate oscillation and discharge of vorticity, and eventually approaches

*

Corresponding author. Tel./fax: +86 29 82663413 Email address: [email protected] (J. Deng)

ACCEPTED MANUSCRIPT

the channel wall with one side vortex shedding; otherwise the particle is rotation shift with swinging wake near the channel center. The fluctuating velocities have little impact on the turning couples on the body, while the average vertical velocity decreases with the rising oscillation frequency. Moreover, the hydrodynamic

CR IP T

interaction between the particles has a close association with the oscillatory effect. The horizontal distance reaches the maximum and the two particles separate as falling forward at the critical frequency; while the pairs in low and high frequencies

AN US

always settle in the same horizontal line despite the different rotation characteristics. The final configuration of the pairs is sensitive to the attractor between the particles, and is caused by the Magnus type of lift balancing the wall repulsion and the

M

interplay between particles.

ED

Keywords: Circular particles; Sedimentation; Arbitrary Lagrangian–Eulerian (ALE)

CE

PT

method; Direct numerical simulation; Oscillatory flow

AC

1. Introduction

Particulate flows, involving solid–fluid interactions, are of great importance to

many industrial and biomedical applications, including removal of pollutants from watercourses, paper manufacturing, pharmaceutical processing and cells clotting in blood vessels. Often, the flow pattern during particle migration is hard to predict due to the variety of the flow field. For the sudden acceleration and deceleration of

ACCEPTED MANUSCRIPT

oscillatory shear, the vortex structure around the particle is always severely affected by the unsteady flow [1]. The study of the dynamic motion of suspensions has attracted increasing attention recently because of the diverse rheological behaviors. The distinguished experimental research on particle migration in a tube with

CR IP T

Newtonian fluid has been conducted by Segré and Silberberg [2]. They remarked that the neutrally buoyant spherical particle searches for an equilibrium radial position at about 0.6 times the tube radius. Afterward, typical numerical work confirmed that due

AN US

to the balance between the wall effect and inertia effect, there is an equilibrium layer between the wall and the channel centerline for both single and multi-particle systems [3, 4]. Meanwhile, theoretical studies [5] applied the perturbation theory to analyze

M

the nonlinear effects of particle-fluid and particle-wall interactions. It has been

ED

pointed out that the mechanism controlling particle motion is the turning couples on the body, resulting from the competition of density ratio, initial state, fluid rheology,

PT

flow type as well as particle shape [6, 7]. Indeed, Shao et al. [8] concluded that the

CE

free sedimentation of a circular particle undergoes steady motion without overshoot, weak oscillatory motion, strong oscillatory motion and irregular oscillation in proper

AC

order at appropriate particle density intervals; Zimmerman [9] indicated that

the

Stokes drag for a sphere settling in a tube depends on the channel to particle blockage ratio , and the particle will either oscillate or tumble in rotation at high blockage ratios. Moreover, it has been revealed that the driving forces influence the particle lateral displacement could be could be classified into five types: an inertial shear

ACCEPTED MANUSCRIPT

induced lift [10], a lift related to rotation [11] , a wall repulsion due to lubrication [12] and in the case of shear flow, a lift caused by the velocity profile curvature [13].The stresses on particle body and the related drag-lift correlation affected by the wall induce the eddy formation, distortion and shedding behind the particle, eventually

CR IP T

giving rise to various migration trajectories and rotation behaviors [14, 15].

In addition, for the particle-particle interaction in a narrow channel, the walls tend to push the pairs together while their mutual repulsion always keeps the particles

AN US

apart. The celebrated drafting-kissing-tumbling (DKT) process for two spherical particles falling against gravity has been observed in Fortes et al.’s [16] experiment with a channel filled with Newtonian fluid. It is noted the trailing particle would be

M

sucked into the wake of the leading one with an accelerated motion to form a pair of

ED

kissing spheres, and the two particles are inclined to tumble and split apart before traveling with the same velocity [17, 18]. Nevertheless, the later numerical results

PT

pointed out that this kind of particles interplay only occurs in certain simulation

CE

configurations [19, 20], and the final states of the two particles are accustomed to be resolved with a staggered rather than cross-stream arrangement in the channel [6].

AC

For instance, Aidun and Ding [21] demonstrated a cascade of period-doubling bifurcations and a chaotic state for two circular particles sedimentation at low Reynolds numbers; Mukundakrishnan et al. [22] observed that two spherical particles circle in an out-of-phase manner in a rotating flow with regard to their radial positions. Theoretical work of Hu [23] indicated that the basic mechanisms

ACCEPTED MANUSCRIPT

controlling the motion and interaction of sediments are the competition between the particle–particle and particle–wall forces, by considering the difference in sizes and initial positions of the particles [24, 25]. The different rotation behaviors, such as anomalous rotation, rotation shift and zigzag migration, with different vortex

CR IP T

structures could be observed for the various sedimentation patterns [26, 27]. Alternately, Wang et al. [28] manifested that the settling of two particles would never undergo the DKT process when the initial longitudinal distances and diameter ratios

AN US

beyond a certain threshold; Dash and Lee [29] showed that the tumbling mechanism of two spheres is dictated by the induced turning couple, which influences the lateral motions together with their repulsive hydrodynamic forces.

M

Generally, all the reports mentioned above refer to steady flow, while the

ED

wall-bounded oscillatory flow exists in many different practical situations (i.e. fluid pumping, respiratory system of living beings, etc.). The previous work [30] suggested

PT

that the fluctuating velocity field affects the wake structure behind the cylinder and

CE

would modify the stress distribution on the settling body. The rotation behavior and the drag-lift correlation in particle-fluid suspensions could be chaotic dynamics in

AC

unsteady shear flows [31, 32]. Through experimental approach, Snook et al. [33] indicated that the spherical particles preferentially migrate away from the channel wall and concentrate at the centerline when the oscillation amplitude is large enough. Recently, Sun et al. [34] and Yapici et al. [35] have investigated the circular particle dynamics in oscillatory Newtonian shear flows using numerical and theoretical

ACCEPTED MANUSCRIPT

method, respectively. Both of them observed that the oscillatory flow promotes the circular particle closer to the centerline, and the oscillation frequency affects its equilibrium position significantly. Particularly, Ingber [36] showed that an attractive static force between the two spheres is formed in oscillatory flows, and the migration

CR IP T

direction changes depending on the amplitude of the oscillation. Our exploratory work concluded that the higher amplitude leads to the smaller range scale zigzag migration for a particle dynamic in oscillatory flow, and there are different flow

AN US

patterns concerning the oscillation amplitude and frequency [37]. However, the motion mechanism for circular particle in oscillatory flows is still not clear, and little research has been devoted particles interaction.

M

The present work focuses on the settling of one and two circular particles

ED

subjected to a pressure driven flow. The impact of oscillation frequency on the flow pattern during migration has been investigated with the finite element arbitrary

PT

Lagrangian-Eulerian method. The remainder of the paper is organized as follows. In

CE

Section 2, the governing equations describing the motion of both the fluid and the particle are introduced. In Section 3, the finite element numerical method is described,

AC

and this method is calibrated and validated with the immersed boundary results. In Section 4, one and two circular particles settling in oscillatory flow are investigated, respectively. The flow microstructures, such as migration trajectory, distribution of pressure, angular velocity, drag force and terminal velocity are analyzed for the different regime behaviors. Conclusions are addressed in Section 5.

ACCEPTED MANUSCRIPT

2. Governing equations As shown in Fig.1 (a), two identical circular particles, denoted by P1 and P2 , with boundaries Pi  t  ( i  1, 2 ), of density  p , diameter D=2R, are released with

domain  with boundary Γ i ( i  1,

CR IP T

zero velocity side by side in a two-dimensional channel of width L. The flow

,4 ) is filled with Newtonian fluid of density  f

and kinematic viscosity f . The right particle is initially in channel center and the left is placed at the same distance between the wall and the right particle surface. The

AN US

gravity g and the external pressure differential P act along the x direction. For convenience, the orientation angle  and the arc length direction of the particle are defined as shown in Figs. 1(b) and 1(c), respectively. Ω

M

Γ3 D P2

PT

Γ4

θ

P (t)

g

x

ED

y

P1

Γ2

(b)

H

AC

CE

L

(c)

Γ1 (a)

(d)

Fig. 1. (a) Problem description and notations used in the calculation of the motion of two circular particles ( P1 and P2 ); (b) orientation angle θ; (c) particle arc length direction; (d) typical computational grid for numerical analysis.

ACCEPTED MANUSCRIPT

2.1 Fluid and particle motion The fluid motion satisfies the conservation of the mass and momentum:  u  0 , Du f  f f     , Dt

(1) (2)

CR IP T

where u=(u, v) is the fluid velocity vector; f is the body force; σ is and the stress tensor. The material derivative of the velocity is given by Du u   u    u . Dt t

(3)

AN US

For a Newtonian fluid, the stress tensor is defined by the simple constitutive relation

   pI   f u + uT  ,

(4)

M

where p is the pressure and I is the unit tensor.

ED

The translational motion and the rotation of the rigid particles are governed by

CE

PT

Newton's second law and Euler equations, respectively, dU M P = F+G, dt dω J = T, dt

(5) (6)

AC

where M is the particle mass; J is the inertia moment matrix; G is the body force external gravity field, which includes gravity and buoyancy; U P = dXP dt = (U x ,U y ) and ω = dΘ dt = ωk are the translational and angular velocity vector of the particle, respectively; F and T are the total force and torque acting on the particle by the surrounding viscoelastic fluid, and can be expressed as, respectively,

ACCEPTED MANUSCRIPT

F=

σ  nds ,

(7)

 x  X P   σ  n  ds ,

(8)

P(t)

T=

P(t)

where n is the unit normal vector on particle surface pointing outwards into the fluid;

gives the position of the center of the particle.

2.2 Particle collision

CR IP T

s is the surface of the particle; x is the coordinate on the surface of the particle; Xp

AN US

In order to prevent the particles from interpenetrating each other, a supplementary collision model should be adopted in present simulation. By following the spring type short-range repulsive force [38], and the particle-particle interaction

M

takes the form,

ED

 0,  FRp   C  1 2 2 x  x  n1 n2   x  x   ,   n2       n1  x  x   n1 n2  p

x n1  x n 2   x n1  x n 2  

,

(9)

PT

where x n1 is the mesh node on P1 and x n 2 is the mesh node on P2 , which are

CE

chosen to ensure that x n1  x n 2 is the smallest distance between the two particles;

AC

C    p  f  Ap g is the relative gravitational force scale, where Ap is the particle

area;  p and  are the stiffness parameter and the threshold distance, respectively. Therefore, a similar repulsive force to handle the collision between the particle

and wall can be written as,

ACCEPTED MANUSCRIPT

 0,  FRw   C  1 2 2 x  x  w ni   x  x   ,   ni      w  x  x   p w ni 

x w  x ni   x w  x ni  

,

(10)

where x ni is the node on particle surface closest to the wall and x w is the

CR IP T

corresponding point on the wall, i.e. x w  x ni is the smallest distance between any point on the circle and the wall .

2.3 Boundary conditions

AN US

No-slip conditions are imposed on the channel walls and on the particle boundary,

(11)

u  V     x  XP  on Pi (t )

(12)

M

u  v  0 on Γ 2 and Γ 4

ED

Furthermore, we assume the outflow section to be the traction-fee form boundary condition [39]. In addition, for the inlet boundary condition, we consider

PT

two types of flows: nonoscillatory flow and oscillatory flow. For the nonoscillatory

CE

flow, namely classical Poiseuille flow, a constant inlet pressure P0 is adopted to

AC

maintain the flow. For the oscillatory flow, the inlet pressure is assumed to be P0 (1+ Asin(2πFt )) , where A and F are the dimensionless amplitude of oscillatory

flow and the oscillation frequency measured in hertz, respectively. Note that, the developed velocity profile of Poiseuille flow in a two dimensional straight channel section is parabolic. Therefore, the velocity corresponding to the constant inlet pressure is imposed on the flow field as the initial condition. It should

ACCEPTED MANUSCRIPT

be mentioned that the real velocity profile with oscillatory pressure inlet may be far from parabolic due to unsteady flow. It requires iterative calculations to get the real information of the computational domain.

CR IP T

3. Numerical method and validation Direct simulation of particles migration in fluids has been carried out using a two-dimensional generalized Galerkin finite element (FE) method based on a

AN US

combined formulation of the fluid and particle momentum equations [23], and an arbitrary Lagrangian–Eulerian (ALE) moving mesh technique is adopted to handle the changes of the fluid domain [40]. For the discretization of the fluid–particle system, we use triangular elements with continuous quadratic interpolation (P2) and

M

linear continuous interpolation (P1) by the Delaunay–Voronoi method. The main

ED

advantage of the ALE formulation with respect to other methods is the use of a

PT

boundary-fitted mesh during the whole simulation leading to high accuracies around the particle where the largest gradients are expected. However, as the mesh becomes

CE

more and more distorted as the particle moves, it will reduce the solution accuracy of

AC

the flow field. Therefore, we should monitor the mesh quality during the particle motions as well, and regenerate a new finite-element mesh if unacceptable element distortion is detected. In this section, we briefly describe the numerical method, and validate it with the immersed boundary scheme.

ACCEPTED MANUSCRIPT

3.1 Arbitrary Lagrangian–Eulerian (ALE) mesh movement From the general kinematic theory for the ALE technique [41] , the material time derivative of the velocity at a time instant t and a given point x in the fluid domain is expressed as,

CR IP T

D u u  x, t     u  uˆ     u , Dt t 

(13)

and the referential time derivative in the domain constant keeping the coordinates  is defined as,

AN US

  . u  x, t   u  x   , t  , t   fixed t t

(14)

Note that x   ,t  is where the fluid mechanics problem generates and uˆ  x,t 

(15)

ED

M

is the domain (or mesh) velocity, it becomes, d x   , t   uˆ . dt

Therefore, when the referential domain is in accord with the spatial domain, it

PT

obtains uˆ  0 ; when the mesh velocity coincides with the particles velocity, we have

CE

uˆ  u . The referential time derivative reduces to the local Eulerian time derivative

and the Lagrangian time derivative, respectively.

AC

In general, the domain (or mesh) velocity in Eq. (15) is only constrained at the

domain boundaries, and has to follow the confining flow geometry and the particles movements. The mesh motion in the fluid satisfies an elliptic partial differential equation [40], i.e. Laplace’s equation, to guarantee its smooth variation,    k euˆ   0

in   t  ,

(16)

ACCEPTED MANUSCRIPT

where k e is a function involving the domain deformation, and it controls the region far from the particles absorbs most of the deformation while the region near the particles retains its shape. With rigid circular particles in this work, the mesh velocity is allowed to slip on

CR IP T

the surface of the particles, thus the nodes on a particle surface can move with the particle with its translational velocity but do not need to rotate with the particle. Therefore, the mesh velocity uˆ satisfies the following boundary conditions, for x Pi  t 

AN US

uˆ  Vi  i   x  Xi  ,

on  i .

uˆ  0

(17) (18)

Similarly, if the particles undergo acceleration, an acceleration field, aˆ  x,t  ,

M

of the domain satisfies,

in   t 

ED

   k eaˆ   0

(19)

and with the boundary conditions,

PT

aˆ  Vi  i   x  Xi   i  Vi

CE

a0

for x Pi  t 

on  i ,

(20) (21)

AC

where Vi  dVi dt and i  di dt . Finally, the weak formulations for the mesh velocity Eq. (16) and acceleration

Eq. (19) can be written as, Find uˆ 

mesh1

and aˆ 

mesh2

, such that

  k uˆ  uˆ d   uˆ  e

0

mesh0

(22)

ACCEPTED MANUSCRIPT

and

  k aˆ  aˆ d   aˆ  e

mesh0

,

(23)

0

where 0 is the domain occupied by the fluid at a given time instant t0;

is a

CR IP T

natural space for the velocity of the fluid–solid mixture, and the function spaces are defined as, mesh1





2

(24)



 aˆ     ; aˆ  Vi  i   x  Xi   i  Vi on Pi  t  ; aˆ  0 on  i , (25) 2

mesh0





 uˆ   0  ; uˆ  0 on  . 2

AN US

mesh2



 uˆ     ; uˆ  Vi  i   x  Xi  on Pi  t  ; uˆ  0 on  i ,

(26)

where    corresponds to the space for the two-dimensional velocity field in the 2

M

fluid.

ED

3.2 Mesh generation and time discretization In this work, we measure the mesh quality by checking the change in the

PT

element volume and aspect ratio in comparison with its value in the initial

f1  max  f1e  and f 2  max  f 2e  , 1e N el

1e N el

(27)

AC

CE

undeformed mesh. The two monitoring parameters are defined as,

where N el is the total number of mesh elements. f1e and f 2e are the changes of the element volume and aspect ratio, respectively, and can be expressed as, f1e  log V e V0e  and f 2e  log  S e S0e 

(28)

with V e and V0e are the volume of the eth element and the initial undeformed mesh,

ACCEPTED MANUSCRIPT

respectively; S e and S0e are the aspect ratio of the eth element and the initial undeformed mesh, respectively. The aspect ratio is given by, (l e )3 S  e , V e

(29)

CR IP T

where l e the maximum length of the sides of the eth element. Once the element volume or the aspect ratio is greater than 1.39 (i.e. f1e or f 2e is larger than four times or smaller than 1 4 of its original value), the mesh is considered too distorted and needs to be regenerated.

AN US

Afterwards, the intermediate velocity field U * , linear momentum P * and angular momentum L* must be projected to satisfy the incompressibility condition and the non-penetration condition,

M

 U n1  0

(31)

ED

U n1  n  U solid, n1  n on P(t ) .

(30)

The detailed decomposition rule for the intermediate velocity and the projection

PT

method could be found in the reference [42] and [43], respectively.

CE

Regarding the time discretization, the momentum balance and continuity equations are decoupled by a finite-difference scheme. As the forces in the

AC

fluid-particle system are internal, the hydrodynamic forces and moments acted on the particle could be eliminated through simplifying the formularies. Therefore, there is no need to compute them explicitly. The velocity field of the system is determined implicitly, while the position of the particle and the grid nodes are updated explicitly [40]. This explicit–implicit scheme is numerically stable and a second-order order

ACCEPTED MANUSCRIPT

backward difference formula (BDF) is used for time integration. The second-order backward difference scheme reads, 3 n1 1 t  2t n  t n1  g1  t n1 , Un1   g w  t n1 , Un1  , n  1 , 2 2

(32)

CR IP T

where g1  t, U  t   accounts for the diffusive and convective terms, g w  t, U  t   is the term that depends on the domain velocity. This scheme can be initialized by a one-step second-order Crank–Nicolson method [44]. However, the formula does not satisfy the geometric conservation law (GCL) [45, 46] in general, when applied to the

AN US

conservative formulation. Therefore, when g=g1  g w is substituted into Eq. (32), the second order BDF scheme is achieved by,

3 n1 1 n1 3 t n1 1 tn n t  2t  t   n g  t , U  t  dt   n1 g  t , U  t  dt , 2 2 2 t 2 t

(33)

M

Moreover, we can modify Eq. (32) by using a more accurate quadrature rule to

ED

integrate the two terms in the right hand side. Unfortunately, if this procedure is

PT

carried out on the whole g, the coerciveness of the form g1 will be destroyed since the second integral appears with a negative sign. Thus, to impose the GCL

AC

g w only.

CE

compliance to the second order BDF scheme we are obliged to operate on the term

3.3 Validation

In this work, the channel to the particle blockage ratio is defined as β = L D . The major axis is D=0.1cm; the density and the kinematic viscosity of fluid are

f  1.0g cm3 and f  1 10-3 Pa  s , respectively; the particle density is  p  1.5f ;

ACCEPTED MANUSCRIPT

the initial pressure gradient is P0  1800Pa . Moreover, for the dimensionless parameters, the amplitude is A  0.5 ; the stiffness parameter is  p  0.25 . Mesh and time convergences are checked in all the simulations. The computational domain height H is also validated to make sure that the particle motion

CR IP T

is not affected by the inflow and outflow boundaries. Furthermore, when the particle reaches quite close to the wall, an extra refinement between the channel wall and the particle is added for the spatial convergence. In Table 1, some parameters of the

AN US

meshes used are presented by referring to a single circular particle with blockage ratio   6 and initial position d0  5 24 . Table 1

M

Mesh parameters for a single circular sedimentation,   6, d0  5 24. M1

M2

M3

On the particle boundary

24

28

32

Total mesh

6809

8765

10272

PT

ED

Mesh label

CE

Meanwhile, the sensibility of the mesh and time step size to the migration trajectory are reported in Figs. 2(a) and (b), respectively. There are two different

AC

oscillation frequencies ( F  10 Hz and 100 Hz ) are presented using the mesh parameters in Table 1. It can be observed that, in both frequencies, a mesh with 28 elements (M2) on the particle boundary is sufficient to achieve convergence, although the comparison with M1 shows very small deviations. Nevertheless, regarding the time convergence, we found that a smaller step size should be chosen as the high

ACCEPTED MANUSCRIPT

oscillation frequency of the flow field and the fast dynamics involving in the small particle-wall distance. Hence, for the numerical simulations presented in the next sections, the time step size is chosen to be Δt  0.001s .

0.8

CR IP T

0.8

F=10 Hz

F=10 Hz

y/L

0.6

y/L

0.6

0.4

0.4 M1, Δt=0.001s M2, Δt=0.001s M3, Δt=0.001s

0.2 0

1

2

3

t (s)

M2, Δt=0.005s M2, Δt=0.001s M2, Δt=0.0005s

0.2

0

1

2

3

t (s)

(b)

M

(a)

F=100 Hz

AN US

F=100 Hz

Fig. 2. Particle trajectories for (a) the meshes labeled as ‘M1’, ‘M2’ and ‘M3’ with time step size

ED

Δt=0.001s and (b) different time step size Δt with ‘M2’. Two oscillation frequencies are presented.

PT

Finally, for validation, a single circular particle free settling under gravity is simulated under the same conditions as used by Luo et al. [26] who solved the

CE

problem with an immersed boundary technique, as shown in Fig. 3. It is clear that

AC

both the lateral displacement and the angular velocity of the particle in our results accord quite well with the data from immersed boundary method. Therefore, the ALE method to investigate the particle migration is accurate and acceptable.

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 3. Free sedimentation of a circular particle in a channel for  p f  1.5 . (a) Lateral displacement and (b) angular velocity of the particle.

M

4. Results and discussion

Fig. 4 shows that a single particle settling seeks an equilibrium position at the

ED

channel side in the nonoscillatory flow independent of its initial position. It can be

PT

seen that all the migrations appear iconic zigzag motion in microscopic observation

CE

due to the tentative interaction of the wall and the local shear rate. When the particle is initially released at the opposite side of its final site, it directly migrates inward

AC

toward the final equilibrium position; whereas for the one near the equilibrium position, it moves to the channel centerline first, and would be repulsed back to the final position after a period of time. This behavior accords quite well with the previous numerical results [7, 20], that the final equilibrium position of a circular particle in the narrow channel has nothing to do with the particle initial position. Thus,

ACCEPTED MANUSCRIPT

in the following investigations, we perform a single circular particle ( P1 ) sedimentation in oscillatory flows with its initial position the same distance between the wall and the right particle surface, i.e. d0  5 24 . This flow configuration is the same as where P1 located in Fig. 1(a), so that we can have a contrast with the two

ED

M

AN US

CR IP T

particles settling process to clarify the particle-wall effect.

CE

initial positions.

PT

Fig. 4. Lateral migration of a circular particle (   6 ) in the nonoscillatory flow with different

AC

4.1 A single particle settling in oscillatory flow In order to have a systematic understanding of the flow pattern during

sedimentation in the oscillatory flow, we have performed a series of simulations to examine the effect of oscillation frequency. However, we choose the typical migration behaviors to discuss in detail, and most of the parameters analyzed are

ACCEPTED MANUSCRIPT

AN US

CR IP T

from a steady state or a time periodic solution that the particle reaches.

Fig. 5. Terminal state of (a) angular velocity and (b) drag force range of particles for different

M

oscillation frequencies at  =6 .

As shown in Fig. 5, the terminal state of angular velocity and drag force range

ED

for different oscillation frequencies are presented. Clearly, the angular velocity shoots

PT

up in the low frequencies and peaks at a critical value, and then the particle tends to swing or nearly to be irrotational at high oscillation frequencies (Fig. 5(a)). In

CE

addition, the drag force range rises at a near-liner trend and the average value remains

AC

at the same level as the frequency increases (Fig. 5(b)). Therefore, judging from the feature of particle rotation in oscillatory flows, there may be different migration patterns at certain oscillation frequency intervals. Owing to the increasing drag force, we can gather that the lateral displacement of the particle, the stresses on the particle surface as well as the particle migration velocities present different characteristics

ACCEPTED MANUSCRIPT

with the unique flow regions, which will be further studied in the following discussions. The typical transversal displacements of a single circular particle migration in oscillatory flows are shown in Fig. 6. It can be found that the particle is eventually at

CR IP T

the channel side at the low oscillation frequency (F=10 Hz), which is identical to the stable rotation in the nonoscillatory flow. However, in the high frequency oscillatory flows (F=80 Hz and 100Hz), the circular particle acquires an equilibrium position

AN US

close to the centerline at last. Combined with the lateral migration trajectories and the low rotation speeds in these two situations, the particle is nearly in a tumbling state at low oscillation frequencies and in a rotation shift state in the high frequency

M

oscillatory flows. In particular, the particle rapidly approaches close to the wall with

ED

the large range zigzag migration at an intermediate frequency (F=30 Hz). Both Xia et al. [47] and Luo et al. [26] indicated that the ‘rotation shift’ is closely associated with

PT

vortex shedding, while the ‘zigzag migration’ is caused by a combination of the wall

CE

effect and the vortex shedding. Therefore, the sedimentation flow in the channel is dominant in this flow configuration, and the wall effect plays a principal role when

AC

the particle is close to the wall. Obviously, the micro phenomena, such as ‘rotation shift’ and ‘zigzag migration’ could be observed for one particle sedimentation, depending on the mean terminal settling Reynolds number and the initial position. It is clear that different oscillation frequencies of flow field correspond to different types of particle motion, and would have the same effects with a certain particle for

ACCEPTED MANUSCRIPT

AN US

CR IP T

different initial stations.

Fig. 6. Typical lateral migration trajectories of a single circular particle settling in the oscillatory

M

flow with different oscillation frequencies at  =6 .

ED

To evaluate the effect of blockage ratio, the lateral migration trajectories of a single circular particle settling at different oscillation frequencies are presented for

PT

different blockage ratios (  =5 and  =8 ), as shown in Figs. 7(a) and (b), respectively.

CE

It can be noted that the finial equilibrium position of the particle is more near to the

AC

channel with the bigger blockage ratio, and the particle will also reach the final position faster in a more narrow channel. These observations are consistent well with previous study by Feng et al. [7] for particle migration in Poiseuille flow. In addition, when the channel is wide enough, i.e.  =8 (Fig. 7(b)), the particle eventually moves close to the left wall at low frequencies, which is quite different from the regular pattern with small blockage ratios (  =5 and 6 ). It seems that the particle

ACCEPTED MANUSCRIPT

rotates in an opposite direction in the wide channel due to the relative smaller velocity profile curvature. Further, the large range zigzag migration at an intermediate frequency (F=30Hz) could also be observed in our simulations for different blockage ratios. Unfortunately, the special zigzag migration for  =8 occurs at the frequency

CR IP T

F=60Hz for a smaller range, and the critical frequency for this phenomenon appears to be rose, regardless of the increasing or decreasing of the blockage ratio. Consequently, the sizes of the particle and the channel width affect the interaction

AN US

between the particle and the wall, and have a significant influence on the particle dynamics in the oscillatory flow.

0.8

M

0.8

PT

0.4

0.2

AC

0

CE

F=0 Hz F=60 Hz

0.4

0.2

F=20 Hz F=100 Hz

F=30 Hz

2

3

1

0.6

y/L

y/L

ED

0.6

t (s) (a)

F=0 Hz F=60 Hz

0

F=20 Hz F=100 Hz

1

2

F=30 Hz

3

t (s) (b)

Fig. 7. Lateral migration trajectories of a single circular particle settling in different channels with (a)  =5 and (b)  =8 at different oscillation frequencies.

Further studies of pressure difference on the particle surface and vortex structure behind the particle are shown in Fig. 8. As previously stated, the results reported here

ACCEPTED MANUSCRIPT

are for the situation when after the initial transient dies out, and the particle reaches the steady state or the time periodic solution. Figs. 8(a) and (b) manifest that the oscillatory flow has little impact on the periodic oscillation in the wake and discharge of vorticity from behind the particle at low frequencies. The steady pressure

CR IP T

difference on particle surface is balanced against the local shear rate to sustain a constant rotation for this simulation. Because of the torque exerted by flow field, the maximum differential pressure decreases spatially from vertical direction to

AN US

horizontal direction. For the micro transversal displacement in this case, the pressure distribution is centrally symmetric. Parallel results, at the intermediate frequency, due to the wall effect the periodic vortex shedding only appears at the left side of the

M

channel, causing the particle to get close to the right wall periodically. The pressure

ED

difference decreases to balance with the local shear rate when the particle is moving away from the wall (Fig. 8(c)). The lift force related to particle rotation is responsible

PT

for this oscillatory motion. In contrast, the vortex shedding has been accelerated and

CE

is almost in the same rhythm of the flow field at high oscillation frequency. As great lateral migration could not occur in such a small time period, the particle has no

AC

evident shift relative to the transient variations of pressure difference on its surface (Fig. 8(d)).

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 8. Pressure difference on particle surface and vortex structure with time steps from a motion

CE

period for different oscillation frequencies at  =6 .

AC

For corresponding to the drag force fluctuating, the terminal velocities in a U x versus U y plot are shown in Fig.9. Due to the different rhythm of particle migration and flow field oscillation, the fluctuating of velocities inevitably occur in different motion cycle. The horizontal velocity range is essentially unchanged at F=10 Hz and 80 Hz on account of the asymmetric vortex shedding; while at the intermediate frequency, the particle moves faster in the horizontal direction for the one side vortex

ACCEPTED MANUSCRIPT

shedding. It can also be noted that the average vertical velocities are almost the same when the equilibrium positions are close to the centerline, and the particle attached to the channel is always accompanied by a slower settling velocity. In general, for the channel to particle blockage ratio in our simulations, the conditions for Hopf

CR IP T

bifurcation are met. The particle is in a tumbling state at the channel side for the nonoscillatory flow, which is shown as the closed circular orbits in Fig.9 (a). As the alternate eddy formation, swaying, shedding and moving in the instantaneous flow

AN US

field, it can be inferred that there’s not enough time for the particle to feed back to the lateral motion at the high frequency. The velocity trajectory of the rotating zigzag motion is presented as the overlapping structure at low frequency (Fig.9 (b)) and

M

turns into the multi-layer reticular structure as the oscillation frequency increases

ED

(Fig.9 (c)). However, in the high frequency flow, the typical offset oscillation is obtained, which appears as the up and down orbits (Fig.9 (d)). The multitudinous

PT

velocities plots are obtained from the unsteady rotation behavior and the increasing

CE

drag force range at different oscillation frequencies, which is significant influenced

AC

by the response time of the fluid-structure interaction.

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Fig. 9. Trajectories of terminal velocities with a few cycles of motion at different oscillation

CE

frequencies for  =6 .

AC

4.2 Two particles settling in oscillatory flow For a pair of particles settling from a horizontal line, the walls tend to push them

together and their mutual repulsion keeps them apart, the typical DKT process is reproduced with free settling in our simulations. However, as shown in Fig. 10(a), the Karman vortex street of single particle settling gives way to a wavy but still attached wake, and the particle-particle interplay is finally resolved with a cross-stream

ACCEPTED MANUSCRIPT

arrangement in the channel flow. According to Figs. 10(b), (c) and (d), the particles retain their flow details at the low oscillation frequency, and suffer the oscillatory effect at relative high frequencies. The weak oscillatory motion is reproduced by the confinement of the wall at F=20 Hz, and the accelerated vortex shedding promotes a

CR IP T

longer attached wake at F=80 Hz. As the rotation gives rise to a smaller Magnus type of lift at this situation, the particles tend to move away from the wall, and eventually keep at the same level. In addition, although the particles are at the same weight in

AN US

our simulations, the in-depth study shows that the average forces acting on the particle clusters are bigger than that on single settlement for the larger blockage ratio. It means that the two-particles will obtain an equilibrium position closer to the wall

AC

CE

PT

ED

M

with lower sedimentation velocity than the separated single particle.

Fig. 10. The effect of oscillation frequency on the terminal vortex structures of the two circular particles settling at  =6 .

ACCEPTED MANUSCRIPT

The relative distances of the conjugate diameter endpoints changes over time under different oscillation frequencies are then illustrated in Fig. 11. It is noted that, during the sedimentation, there is a transition from primary base state to periodic spinning or swaying of the circle through the so called ‘hopf bifurcation’. The

CR IP T

rotation or sway of the particle would remain constant as time goes on. This periodic motion agrees qualitatively with the previous theoretical results [48]. Moreover, at relative low frequencies, the particle could rotate in a circle, and the orientation is

AN US

prone to be continual altering when the oscillation of the flow field goes up. The angular velocity is expected to increase with the oscillation frequency within certain realms. However, as revealed in Fig. 10, in high frequency (F=80 Hz), the body could

M

not turn in a circle due to the rapid change of the flow field. The rotating speeds of

AC

CE

PT

ED

particles are quite low for the small swing angle at this situation.

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 11. The relative distance of the conjugate diameter endpoints changes with time at different

CE

oscillation frequencies at  =6 .

AC

In addition, as the distance between the particle-particle surfaces is smaller than the distance between particle surface and nearest wall surface for both particles, the interparticle repulsions initially dominate the particles migration. Therefore, once released from rest, the right particle is rapidly repulsed towards the right wall and the left particle gradually moves to the left wall. Due to the curvature effect, the two particles rotate at the same sense in the opposite direction and will never appear in the

ACCEPTED MANUSCRIPT

AN US

CR IP T

same side any more.

Fig. 12. Relative (a) horizontal distance and (b) vertical distance of the two circular particles

M

settling at different oscillation frequencies at  =6 .

ED

Finally, the relative horizontal distance and vertical distance of the two particles settling at different oscillation frequencies are presented in Figs. 12(a) and (b),

PT

respectively. It can be observed that, the horizontal distance peaks at the critical value,

CE

and declines at a high oscillation frequency. As demonstrated previously, at both the low and high frequencies, the two particles momentarily form an attractor in the

AC

vertical direction, and have traded the lead from time to time, eventually damped by viscosity with the same level. While at the critical frequency, the vertical distance of the pairs increases with the time. The leading particle ( P2 ) accelerates to extricate the trailing one ( P1 ) as the oscillatory effect pushed it to the low shear rate region. The attractor is destroyed by the bigger horizontal distance. Consequently, the particles

ACCEPTED MANUSCRIPT

will not align in a horizontal line any more and separate while falling in the channel. The dynamic response of the streamwise oscillation magnifies the fluid damping effect originating from the relative velocity changing between the particle and the fluid. The oscillation frequency has an impact on the particles settling, and these

CR IP T

complex behaviors are related to the particle–particle and particle–wall forces.

5. Conclusion

AN US

The present study investigates the dynamic behaviors of one and two circular particles sedimentation in a narrow channel subjected to the oscillatory pressure driven flow with a finite element arbitrary Lagrangian-Eulerian method. It is found

M

that there is a critical oscillation frequency for the particles settling. When the frequency is smaller than the critical value, the single particle increases its rotation

ED

speed, getting close to the wall with one side vortex shedding; while at high

PT

frequency, the particle is almost irrotational with the accelerated alternate shedding

CE

near the channel center. According to the analysis of the terminal angular velocity, drag force and settling velocity for the different regime behaviors, the Magnus type of

AC

lift and the wall effect related to the blockage ratio are revealed to be responsible for the complex behaviors of particles in oscillatory flow. Moreover, the quintessential attractor of two particles settling disappears at the critical frequency, the horizontal distance reaches the maximum and the two particles will separate. Currently, the vortex shedding resonance of a particle settling is suggested to explain this particular

ACCEPTED MANUSCRIPT

phenomenon. It has certain significance on the transportation and separation of particles, and deserves further investigations.

Acknowledgements

CR IP T

The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China under Grant No 21376187.

AN US

Reference

[1] Ern P. Wake-Induced Oscillatory Paths of Bodies Freely Rising or Falling in Fluids. Annual Review of Fluid Mechanics. 2012;44:97-121.

M

[2] Segre G, Silberberg A. Radial particle displacements in Poiseuille flow of suspensions. Nature. 1961;189:209-&.

ED

[3] Matas J-P, Morris JF, Guazzelli É. Inertial migration of rigid spherical particles in

PT

Poiseuille flow. J Fluid Mech. 2004;515:171-95. [4] Jebakumar AS, Premnath KN, Abraham J. Lattice Boltzmann method simulations

CE

of Stokes number effects on particle trajectories in a wall-bounded flow. Comput

AC

Fluids. 2016;124:208-19. [5] Yang B, Wang J, Joseph D, Hu HH, Pan T-W, Glowinski R. Migration of a sphere in tube flow. J Fluid Mech. 2005;540:109-31.

[6] Feng J, Hu HH, Joseph DD. Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid Part 1. Sedimentation. J Fluid Mech.

ACCEPTED MANUSCRIPT

1994;261:95-134. [7] Feng J, Hu HH, Joseph DD. 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. 1994;277:271-301.

walls. J Zhejiang Univ Sci. 2004;5:111-6.

CR IP T

[8] Shao XM, Lin JZ, Yu ZS. Sedimentation of a single particle between two parallel

[9] Zimmerman WB. On the resistance of a spherical particle settling in a tube of

AN US

viscous fluid. Int J Eng Sci. 2004;42:1753-78.

[10] Holzer A, Sornmerfeld M. Lattice Boltzmann simulations to determine drag, lift and torque acting on non-spherical particles. Comput Fluids. 2009;38:572-89.

M

[11] Bagchi P, Balachandar S. Effect of free rotation on the motion of a solid sphere

ED

in linear shear flow at moderate Re. Phys Fluids. 2002;14:2719-37. [12] Hu HH. Motion of a circular-cylinder in a viscous-liquid between parallel plates.

PT

Theor Comp Fluid Dyn. 1995;7:441-55.

CE

[13] Pan TW, Glowinski R. Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow. J Comput Phys. 2002;181:260-79.

AC

[14] Huang PY, Feng J, Joseph DD. The turning couples on an elliptic particle settling in a vertical channel. J Fluid Mech. 1994;271:1-16.

[15] Derksen J, Larsen R. Drag and lift forces on random assemblies of wall-attached spheres in low-Reynolds-number shear flow. J Fluid Mech. 2011;673:548-73. [16] Fortes AF, Joseph DD, Lundgren TS. Nonlinear mechanics of fluidization of

ACCEPTED MANUSCRIPT

beds of spherical particles. J Fluid Mech. 1987;177:467-83. [17] Zhang H, Tan YQ, Shu S, Niu XD, Trias FX, Yang DM, et al. Numerical investigation

on

the role of discrete element

method

in

combined

LBM-IBM-DEM modeling. Comput Fluids. 2014;94:37-48.

CR IP T

[18] Dash SM, Lee TS, Lim TT, Huang H. A flexible forcing three dimension IB-LBM scheme for flow past stationary and moving spheres. Comput Fluids. 2014;95:159-70.

AN US

[19] Shao XM, Liu Y, Yu ZS. Interactions between two sedimenting particles with different sizes. Applied Mathematics & Mechanics. 2005;26:407-14. [20] Liu ML. Numerical simulation of particle sedimentation in 3D rectangular

M

channel. Appl Math Mech-Engl. 2011;32:1147-58.

ED

[21] Aidun CK, Ding EJ. Dynamics of particle sedimentation in a vertical channel: Period-doubling bifurcation and chaotic state. Phys Fluids. 2003;15:1612-21.

PT

[22] Mukundakrishnan K, Hu HH, Ayyaswamy PS. The dynamics of two spherical

CE

particles in a confined rotating flow: pedalling motion. J Fluid Mech. 2008;599:169-204.

AC

[23] Hu HH. Direct simulation of flows of solid-liquid mixtures. Int J Multiphase Flow. 1996;22:335-52.

[24] Sun R, Chwang AT. Interactions between two touching spherical particles in sedimentation. Phys Rev E. 2007;76:046316. [25] Wang YL. Simulation of sedimentation of two circular particles with collision

ACCEPTED MANUSCRIPT

considered in vertical channel. Appl Math Mech-Engl. 2006;27:983-91. [26] Luo K, Wei A, Wang Z, Fan J. Fully-resolved DNS study of rotation behaviors of one and two particles settling near a vertical wall. Powder Technol. 2013;245:115-25.

CR IP T

[27] Tan JH, Luo K, Fan JR. Three-dimensional particle-resolved direct numerical simulation study of behaviors of a particle settling near a vertical wall. Powder Technol. 2015;275:290-303.

AN US

[28] Wang L, Guo ZL, Mi JC. Drafting, kissing and tumbling process of two particles with different sizes. Comput Fluids. 2014;96:20-34.

[29] Dash S, Lee T. Two spheres sedimentation dynamics in a viscous liquid column.

M

Comput Fluids. 2015;123:218-34.

ED

[30] Konstantinidis E, Balabani S, Yianneskis M. The effect of flow perturbations on the near wake characteristics of a circular cylinder. J Fluid Struct.

PT

2003;18:367-86.

CE

[31] Kim S, Lee KB, Lee CG. Theoretical approach on the turbulence intensity of the carrier fluid in dilute two-phase flows. Int Commun Heat Mass. 2005;32:435-44.

AC

[32] Nilsen C, Andersson HI. Chaotic rotation of inertial spheroids in oscillating shear flow. Phys Fluids. 2013;25:013303.

[33] Snook B, Butler JE, Guazzelli É. Dynamics of shear-induced migration of spherical particles in oscillatory pipe flow. J Fluid Mech. 2015;786:128-53. [34] Sun B, Yu ZS, Shao XM. Inertial migration of a circular particle in

ACCEPTED MANUSCRIPT

nonoscillatory and oscillatory pressure-driven flows at moderately high Reynolds numbers. Fluid Dyn Res. 2009;41:055501. [35] Yapici K, Powell RL, Phillips RJ. Particle migration and suspension structure in steady and oscillatory plane Poiseuille flow. Phys Fluids. 2009;21:053302.

CR IP T

[36] Ingber MS. Combined static and hydrodynamic interactions of two rough spheres in nonlinear shear flow. Journal of Rheology (1978-present). 2010;54:707-18.

AN US

[37] Yuan WJ, Deng JQ, Cao Z, Mei M. Sedimentation of an elliptical particle in periodic oscillatory pressure driven flow. Fluid Dyn Res. 2015;47:065504. [38] Glowinski R, Pan TW, Hesla TI, Joseph DD, Periaux J. A fictitious domain

M

approach to the direct numerical simulation of incompressible viscous flow past

2001;169:363-426.

ED

moving rigid bodies: Application to particulate flow. J Comput Phys.

PT

[39] Gresho PM. Incompressible fluid dynamics: Some fundamental formulation

CE

issues. Annual Review of Fluid Mechanics. 1991;23:413-53. [40] Hu HH, Patankar NA, Zhu MY. Direct numerical simulations of fluid–solid

AC

systems using the arbitrary Lagrangian–Eulerian technique. J Comput Phys. 2001;169:427-62.

[41] Hughes TJ, Liu WK, Zimmermann TK. Lagrangian-Eulerian finite element formulation for incompressible viscous flows. Comput Method Appl M. 1981;29:329-49.

ACCEPTED MANUSCRIPT

[42] Bell JB, Colella P, Glaz HM. A second-order projection method for the incompressible Navier-Stokes equations. J Comput Phys. 1989;85:257-83. [43] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys. 1997;135:118-25.

CR IP T

[44] Formaggia L, Nobile F. Stability analysis of second-order time accurate schemes for ALE–FEM. Comput Method Appl M. 2004;193:4097-116.

[45] Mackenzie J, Mekwi W. An unconditionally stable second-order accurate

AN US

ALE–FEM scheme for two-dimensional convection–diffusion problems. IMA Journal of Numerical Analysis. 2012;32:888-905.

[46] Boffi D, Gastaldi L. Stability and geometric conservation laws for ALE

M

formulations. Comput Method Appl M. 2004;193:4717-39.

ED

[47] Xia ZH, Connington KW, Rapaka S, Yue PT, Feng JJ, Chen SY. Flow patterns in the sedimentation of an elliptical particle. J Fluid Mech. 2009;625:249-72.

PT

[48] Szeri AJ, Milliken WJ, Leal LG. Rigid particles suspended in time-dependent

CE

flows: Irregular versus vegular motion, disorder versus order. J Fluid Mech.

AC

1992;237:33-56.