Particle migration in bounded shear flow of Giesekus fluids

Particle migration in bounded shear flow of Giesekus fluids

Journal Pre-proof Particle migration in bounded shear flow of Giesekus fluids Bingrui Liu , Jianzhong Lin , Xiaoke Ku , Zhaosheng Yu PII: DOI: Refere...

2MB Sizes 1 Downloads 46 Views

Journal Pre-proof

Particle migration in bounded shear flow of Giesekus fluids Bingrui Liu , Jianzhong Lin , Xiaoke Ku , Zhaosheng Yu PII: DOI: Reference:

S0377-0257(20)30001-X https://doi.org/10.1016/j.jnnfm.2020.104233 JNNFM 104233

To appear in:

Journal of Non-Newtonian Fluid Mechanics

Received date: Revised date: Accepted date:

24 July 2019 17 November 2019 29 December 2019

Please cite this article as: Bingrui Liu , Jianzhong Lin , Xiaoke Ku , Zhaosheng Yu , Particle migration in bounded shear flow of Giesekus fluids, Journal of Non-Newtonian Fluid Mechanics (2020), doi: https://doi.org/10.1016/j.jnnfm.2020.104233

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2020 Published by Elsevier B.V.

Highlights 

The pattern of particle migrations can be roughly classified into “returning” and “passing”.



The effects of some factors on the particle migrations are explored.



The separatrix between “returning” and “passing” pattern is dependent on some factors.

Particle migration in bounded shear flow of Giesekus fluids Bingrui Liu, Jianzhong Lin*, Xiaoke Ku, Zhaosheng Yu Department of Mechanics, State Key Laboratory of Fluid Power and Mechatronic Systems, Zhejiang University, 310027 Hangzhou, China

________________________________________________________________________________ ABSTRACT In this paper, the migrations of two interacting particles in a three-dimensional bounded shear flow of Giesekus fluids are numerically investigated using the direct forcing/fictitious domain method for the Weissenberg number ranging from 0.1 to 1.0, the mobility parameter α which quantifies the shear-thinning effect ranging from 0.01 to 0.7, and the ratio of the solvent viscosity to the total viscosity being 0.1. The model is first validated by comparing the numerical results with the available data in the literature. The effects of the Weissenberg number, the shear-thinning effect, and the initial vertical distance between two particles on the particle migrations are explored. Some of the results are in agreement with the experimental ones. The results show that the pattern of particle migrations can be roughly classified into “returning” and “passing”. The variations of particle velocity and pressure field in the “returning” pattern are totally different from those in the “passing” pattern. The separatrix between “returning” and “passing” pattern is dependent on the initial vertical distance between two particles, the Weissenberg number and the shear-thinning effect. With other parameters fixed, the trajectories of particle change from the “returning” pattern to the “passing” pattern as the initial vertical distance between two particles and the Weissenberg number increase, but as the shear-thinning effect decreases. Key words: particle migration; Giesekus fluid; shear flow; Weissenberg number; fictitious domain method ________________________________________________________________________________ 1. Introduction Particle migration in non-Newtonian fluids occurs widely in many natural, industrial and biomedical systems. It is necessary to understand particles’ behavior in non-Newtonian fluids for applications, e.g., in medical applications, the fluid viscoelasticity can be used to achieve isolation __________________________________ *: Corresponding author, Email: [email protected], Fax: 0086-571-87952882.

of human leukocytes [1] and focusing of red blood cells [2]. The particle migrations will show complex patterns because of the complexity of non-Newtonian fluids in rheological behavior. Therefore, the particle migration in non-Newtonian fluids is an important object that deserves more attention. In the last decades, much work has been done to investigate the influence of non-Newtonian property on the dynamics of particles. D’Avino and Maffettone [3] reviewed the existing literatures and discussed the dynamics of non-Brownian particles in viscoelastic liquids. They focused on the peculiar phenomena induced by fluid viscoelasticity and the interior mechanism. In addition, D’Avino et al. [4] reviewed the experiments and theoretical predictions and presented the current understanding of viscoelasticity-driven cross-flow migration. In Newtonian fluids, due to the linearity of the Stokes equations, the analytical solution for the problem of two particles in a creeping linear flow is obtainable [5]. There exist, however, some non-Newtonian properties such as shear-thinning and shear-thickening for particles in non-Newtonian fluids at moderate volume fractions because of the inter-particle interactions [6-7]. Compared to the particles in Newtonian fluids, there are more complex behaviors of particle migration in non-Newtonian fluids, for instance, multi-particle aggregation in sheared suspensions [8-12] and interaction patterns of two-body or three-body particles [13-16]. Much effort has been made in understanding the effect of fluid viscoelasticity on the dynamics of suspended particles. To explore the two-body interaction in viscoelastic fluid, Hwang et al. [17] presented a new finite element scheme for two-dimensional simulation of particle suspensions in simple shear flow of Oldroyd-B fluid. The results indicated that, unlike the kissing-tumbling-separation phenomenon in shear flow of Newtonian fluid, particles in Oldroyd-B fluid tend to approach and tumble (kissing-tumbling-tumbling). Yoon et al. [15] focused on the particle-particle interactions in a shear flow of Oldroyd-B fluid confined by two sliding plates in two and three dimensions. The particle behaviors were classified into three patterns, i.e., pass, return and tumble, which strongly depend on the initial separation distance in the gradient direction. If the separation distance is sufficiently large, the particles “pass” one each other, while the particles reverse the direction and “return” to the end of channel if the separation distance is sufficiently small. The transitions between “passing” and “returning” is tumbling (particles approach and rotate as an elastic dumbbell). The patterns (passing, tumbling, returning) are affected by many factors such as the confinement, elasticity and Weissenberg number. The “returning”

behavior is associated with the wall confinement [15, 17, 18] which provides a repulsive force on the particles and causes particles to move in the opposite direction. The “tumbling” motion has never been found for Giesekus fluid with two-dimensional simulations [18]. The “tumbling” behavior is related to the Oldroyd-B viscoelastic model and resulted from elasticity instead of shear-thinning [15, 19]. Meanwhile the simulation results showed that increasing Weissenberg number leads to an increase in separatrix between “passing” and “tumbling” [15].Vazquez-Quesada et al. [20] extend the three-dimensional Smoothed Particle Hydrodynamics non-colloidal particulate model to investigate the dynamics of single and mutually interacting “noncolloidal” rigid spheres under shear flow with Oldroyd-B fluid. The effects of Reynolds number and wall confinement on the particle behavior are great, especially on the tumbling behavior. And the repulsion force between particles also affect the two-body trajectories. Chiu et al. [21] adopted distributed Lagrange multiplier/fictitious domain method to study two-particle interactions in a three-dimensional bounded shear flow of Oldroyd-B fluids, and found that, although the trajectories of “passing” and “tumbling” are similar to those in a Newtonian fluid, they lose the symmetry due to the effect of elastic force arising from viscoelastic fluids. For relatively high Weissenberg number, there is less “returning” behavior. In addition, the results indicated that the particle can tumble first and kayak later for higher Weissenberg number. Nie and Lin [16] applied the lattice Boltzmann direct forcing/fictitious domain to simulate the motions of three circular particles in a two-dimensional bounded shear flow for power-law fluid, and revealed that particles are more likely to pass for higher Reynolds numbers or lower power-law indices. In addition to the numerical simulations, Snijkers et al. [19] explored experimentally the effect of viscoelastic fluids on the interaction between two equally sized spheres in shear flow, and found that both “passing” and “returning” behaviors occurred and the “passing” trajectory was strongly asymmetric in wormlike micellar surfactant solution and shear-thinning elastic polymer solution. In addition, the start-up effect is important and only “passing” behavior is observed under the reported experimental conditions in Boger fluids. Based on the experimental results, they attributed the overall nature of interactions to the shear-thinning of viscosity, rather than the relative magnitude of normal stress differences. The other peculiar phenomenon in non-Newtonian is particle alignment (chain-like structure) with regard to which a conclusive explanation for the role of the specific rheological property of the

suspensions has not been given since this phenomenon was first visualized [8]. Loon et al. [22] carried out a series of experiments to elucidate the role of fluid properties as well as the geometrical conditions in string formation. The chain-like structure has been observed in a suspension of fd virus particles, which does exhibit shear-thinning but has no measurable normal stress difference. It revealed that shear-thinning is an essential condition for particle alignment and the elasticity is not a necessary rheological feature, although the normal stresses can promote the phenomenon. Additionally, concentration and confinement effects promote alignment probably due to increasing collision probability and inhibiting diffusion. To the contrary, Jaensson et al. [23] summarized the different conclusions about which fluid property can determine particle aggregation. They employed the finite element method to simulate the particle alignment of two and three spherical, rigid particles in a viscoelastic shear flow of Giesekus fluids. The results indicated that alignment was only observed in a constant, elastic Oldroyd-B fluid, but not observed in non-elastic shear-thinning fluids. They concluded that the normal stress differences determine the particle alignment, whereas shear-thinning can have a strong promotion. The different views need to be further clarified. As stated above, the three-dimensional numerical simulation of particle interaction and migration in a bounded shear flow of Giesekus fluids is rarely addressed in the literature. Therefore, the objective of the present paper is to simulate numerically the migration of two interacting particlesin a three-dimensional flow bounded by two sliding plates with the Giesekus fluid model. Attentions are concerned on the effect of Weissenberg number, the shear-thinning effect, the initial vertical distance between two particles on the particle migrations. The rest of this paper is organized as follows. The numerical method is introduced and validated by comparing the present results with the available data in the literature in Section 2. In Section 3, the simulation results are presented and discussed. Finally, conclusions are drawn in Section 4. 2. Numerical method 2.1. Fictitious domain method In this paper, the numerical simulations are based on the direct forcing/fictitious domain (DF-FD) method which is an improved version of the earlier distributed Lagrange multiplier/ fictitious domain (DLM/FD) method [24-26]. The DLM/FD method was originally developed by

Glowinski et al. [27]. Generally, the key point of this approach is that the interior of the particle is filled with fluid and a pseudo body-force (i.e. distributed Lagrange multiplier) is introduced over the particle inner domain to enforce the fictitious fluid to satisfy the rigid-body motion constraint. The algorithm of DF/FD will be briefly introduced thereinafter and the detailed description can be found in [28]. For the sake of simplicity, only one particle is considered in the following formulae. The particle density, volume, moment of inertia, translational velocity and angular velocity are ρs, Vp, J, U and ωp, respectively. Furthermore, the fluid density, fluid polymer viscosity and fluid solvent viscosity are ρf, ηp and ηs, respectively. Let P(t) and Ω denote the solid domain and the entire domain including both interior and exterior of the solid body, respectively. Following scales are introduced for the nondimensionalization: H for length, U0 for velocity, H/U0 for time, ρf U20 for pressure and ρf U20 ⁄H for pseudo body force. Therefore, the dimensionless FD formulations for the incompressible fluid consist of the following three parts: (1) Combined momentum equations

1 r    B  λ inΩ   2u u  u u  r  p  t Re Re Wi

(1)

u  U  ωp  r in P(t)

 r  1Vp

 r  1 J *

(2)

dU    λdx P dt

dωp dt

(3)

   r  λdx

(4)

P

(2) Continuity equation u  0

in Ω

(5)

(3) Giesekus constitutive equation B α BI T  u B  B u   u   B  0  B  I 2  t Wi Wi

in Ω

(6)

where u, p, and λ represent the fluid velocity, fluid pressure and Lagrange multiplier (pseudo body-force), respectively, r is the position vector with respect to the centre of mass of the particle,

ρr=ρs/ρf is the particle/fluid density ratio, Re=ρfU0H/η0 is the Reynolds number with η0=ηs+ηp being the total zero shear-rate viscosity, 𝑉𝑝∗ = 𝑉𝑝 ⁄𝐻 3 is the dimensionless particle volume, 𝐽∗ = 𝐽⁄𝜌𝑠 𝐻 5 is

the dimensionless moment of inertia, 𝑊𝑖 = 𝜆𝑡 𝑈0 ⁄𝐻 is the Weissenberg number with 𝜆𝑡 being the fluid relaxation time, B is the polymer configuration tensor and related to the polymer stress τ=ηp(B-I)/λt, 𝛼 is the mobility parameter which quantifies the entity of the shear-thinning (α=0 gives the Oldroyd-B model without shear-thinning). A fractional-step time scheme is employed to decouple the combined Eqs.(1)-(6) into the following three sub-problems. (1) Fluid sub-problem for 𝒖𝑓∗ and p: u  un r 2u 1   p  3G n  G t 2 Re 2

n  1



r  2u

n

2 Re

  u  0



n

(7) (8)

which is fundamentally the solution of the Navier-Stokes equation and where 𝐆 = −𝒖 · ∇𝒖 + [(1 − 𝛽)⁄𝑅𝑒𝑊𝑖 ]∇ · 𝚩 . A finite-difference-based projection method on a homogeneous halfstaggered grid can be used for spatial discretization, and all spatial derivatives are discretized with the second-order central difference scheme [28]. (2) Particle sub-problem for 𝑼𝑛+1 , 𝝎𝑛+1 , 𝜆𝑛+1 , 𝒖𝑛+1

rVp

r

 u  U n 1 Un   r  1Vp     λ n dx P t t t  

J ωnp 1 t

  r  1

J ωnp

 u    r    λ n dx P t  t 

(9)

(10)

where particle velocities 𝑼𝑛+1 and 𝝎𝑛+1 can be obtained without iteration because Eqs.(9)-(10) 𝑝 have been reformulated so that all the terms on the right-hand side are known, hence 𝝀𝑛+1 defined at the Lagrangian nodes can be determined by Eq.(11) and 𝒖𝑛+1 defined at the Eulerian nodes can be corrected from Eq.(12):

λ

n 1



U n 1  ωnp 1  r  u t



un 1  un  t λ n 1  λ n

 λn

(11)



(12)

where u and λ are defined at the Eulerian and Lagrangian nodes, respectively. Therefore, the tri-linear interpolation function is used to transfer u from Eulerian nodes to Lagrangian nodes and λ from Lagrangian nodes to the Eulerian nodes.

(3) Constitutive equation sub-problem for B Bn 1  Bn  un 1 Bn  Bn un 1  un 1 t





T

 Bn 



α Bn  I Wi



2



Bn+1 -I 0 Wi

(13)

in which 𝒖𝑛+1 has been determined from the first two sub-problems and Eq.(13) can be solved with the first-order time scheme, the central difference scheme for the velocity gradient, and the third-order up winding MUSCL scheme [29] for the convective term. The FD method has been widely used and well validated [24-26, 30-33]. In this work, we further verify the feasibility of the FD method by comparing the numerical results on particle rotation in a three-dimensional shear flow of Oldroyd-B to those obtained from the finite element method and experiments [34]. We consider the case of a particle with radius of 0.1 freely rotating in an Oldroyd-B fluid (the particle-fluid density ratio is 1) with its mass center at (0, 0, 0). The computational domain is Ω=(-1,1)×(-1,1)×(-1,1) and the grid size is h=H/256, and Re=0.1, Wi=1.0, β=ηs/η0=0.5, time step Δt=5×10-4. Angular velocities of particle are shown in Fig.1 where the simulated results show a good agreement with those of corresponding experiments.

Fig.1.Comparison of angular velocities of particle freely rotating in an Oldroyd-B fluid with its mass center at (0,0,0) for different Weissenberg number.

In order to verify the feasibility of the FD method further, the migrations and trajectories of two interacting particles in Stokes flow of Newtonian fluid are simulated and compared with the results given by Choi et al. [18] under the condition of same parameters. The trajectories are shown in Fig. 2 where we can see that the two results are very consistent.

Fig.2. Comparison of trajectories of two particles in Newtonian fluid for different vertical separations.

2.2. Collision model In order to prevent the mutual penetration of particles when two particles collide, the following soft-sphere collision model is adopted:





Fij  F0 1  dij dc nij

(14)

where Fij, dij and nij represent the repulsive force acting on particle j from particle i, the distance between two particles and the unit normal vector pointing from the center of particle i to that of particle j, respectively, dc indicates a cutoff distance and the repulsive force is activated when dij
and 2h in the streamwise and longitudinal directions, respectively. The channel length, width and height are denoted by L, W and H, and the computational domain spans over [-W/2, W/2], [-H/2, H/2], [-L/2, L/2] in x, y, z directions, respectively. A Cartesian reference is set with its origin at the channel center. The periodic boundary conditions are imposed in both streamwise and crosswise directions, and the no-slip boundary condition is applied at channel walls and particle surface.

Fig.3. Schematic diagram of two interacting particles in a bounded shear flow.

Fig.4. Schematic diagram of particles in y-z plane.

Moreover, the constant velocity of both walls in the opposite direction is U0/2 and thus the shear rate is 𝛾̇ = 𝑈0 ⁄𝐻 . The block ratio ε is defined as ε=dp/H where dp is the particle diameter. In order to balance the computational efficiency and accuracy, we select L=3H and W=2H. In the simulation, the mesh size is 1/128, time step is 5×10-4, the Reynolds number is 0.1 and block ratio ε is 0.2. 3. Results and discussion 3.1. Trajectories of particle In this section, the trajectories of two particles under the conditions of Wi=1.0, α=0.01, Re=0.1

and β=ηs/η0=0.1 are simulated to investigate the particle migration induced by fluid viscoelasticity. The centers of mass of particles are initially located at (0, -h, l) and (0, h, -l), respectively. The particles are initially at rest with a fixed horizontal (streamwise) distance l= 0.5 but a varied vertical (longitudinal) separation h and then released. In order to explore the effect of h, the simulations with different h have been performed. Figure 5 shows the trajectories of particle projected on the y-z plane of the flow. As shown in Fig.5(a), when initial h is small, two particles initially move toward each other and then move in the opposite direction. We call this behavior as “returning” pattern. When initial h is large, the different behavior is shown in Fig.5(b) where two particles initially move toward each other along the horizontal direction with a fixed h and then suddenly move in the opposite direction toward the walls under the action of repulsive force, finally move in the opposite direction along the horizontal direction after moving away from the walls slightly under the repulsive action of walls. We call this behavior as “passing” pattern.

(a) “returning” pattern (b) “passing” pattern Fig.5. Schematic representations of two-particle interactions on the shear plane.

Figure 6 shows the trajectories of particle for different initial vertical separations D (D is defined as D=h/a with a being the particle radius). The solid and dash lines denote the trajectories of particle 1 and 2, and dots and diamonds represent the initial and final position of particles, respectively. We can see that there exist two patterns, i.e., “returning” and “passing”. The trajectories of particles are in agreement with those of experiments [19]. Comparing the results in present simulation with those in Oldroyd-B fluid [15, 17, 21], the trajectories of particles do not show the “tumbling” in Giesekus fluid under the given conditions, which is in accordance with the experimental results [19]. Two particles show the “returning” pattern when D is small, and the trajectories of two particles are anti-symmetric to the axes of y=0 and z=0. In the “passing” pattern,

all particles with different initial positions have same trajectories when they move toward the walls and then move along the horizontal direction. This means that the interaction between particles is negligible when two particles move far away from each other. It is worth noting that the value of D corresponding to the separatrix between “returning” and “passing” is 0.17 under these specific parameters.

Fig. 6. Trajectories of particle for different D (l=0.5, Wi=1.0, α=0.01, β=0.1).

Variations of velocity for particle 1 with time for different D are shown in Fig.7. We can see that the effect of D on the particle velocity is significant for the “passing” pattern. The larger the D is, the faster the particle velocity reaches its maximum and the sooner the particle reaches the final position. With regard to the “passing” pattern, although the time required for a particle to reach its maximum and minimum velocity is different for different D, the values of maximum, minimum and terminal vertical velocity as well as the values of maximum and terminal horizontal velocity are almost identical. The vertical velocity hardly changes with time for D=0.1 as shown in Fig.7(a). The horizontal velocity decreases when the two particles move toward each other as shown in Fig.7(b), and the larger D is, the greater the reduction is.

(a) vertical velocity (b) horizontal velocity Fig.7. Variation of velocity for particle 1 with time for different D (l=0.5, Wi=1.0, α=0.01, β=0.1).

3.2. Influence of fluid rheological property 3.2.1. Effect of Weissenberg number The trajectories of particle projected on the y-z plane for different initial vertical separations D and Weissenberg numbers (Wi) are shown in Fig.8. Weissenberg number compares the elastic forces to the viscous forces. We can see that whether the particle trajectory is “returning” pattern or “passing” pattern depends on D, however, the values of D corresponding to the separatrix between “returning” and “passing” pattern are different for different Wi. This can be more clearly illustrated in Fig.9 where the phase diagram for Wi and D has been produced, and circles and triangles represent “passing” and “returning” pattern, respectively. From this figure we can obtain the separatrix between “returning” and “passing” pattern. Increasing Wi causes the value of critical D to decrease (critical D is the value of D corresponding to the separatrix between “returning” and “passing” pattern), i.e., the fluid elasticity causes the particle trajectory changes to the “passing” pattern at a smaller D.

(a) Wi=0.1

(b) Wi=0.3

(c) Wi=0.5

(d) Wi=0.7

(e) Wi=1.0 Fig.8. Trajectories of particle for different D and Wi (l=0.5,α=0.01, β=0.1).

Fig. 9. Phase diagram for Wi and D and Wi (l=0.5, α=0.01, β=0.1).

Fig. 10.Trajectories of particle for different D (l=0.5, α=0.01, β=0.1).

For further illustrate the effect of Wi on the particle migration, the trajectories of particle for different D and Wi are shown in Fig.10 where only the trajectories of particle 1 are given because the motion of two particles is centrally symmetric. It can be seen that the trajectories of particle show the “returning” pattern and “passing” pattern at D=0.1 and 0.3, respectively, for all values of Wi. However, the trajectories of particle show the “returning” pattern and “passing” pattern when Wi≤0.5 and Wi≥0.7, respectively, at D=0.2, which illustrates that the trajectories of particle tend to the “passing” pattern as the fluid elasticity is enhanced. Actually the value of Wi corresponding to the separatrix between “returning” and “passing” is 0.6 under these specific parameters. With regards to the effect of Wi on the “returning” pattern, the larger Wi is, the smaller the minimum distance between particles 1 and 2 is, and the larger the terminal vertical shift of particle (the

vertical distance between the final and initial position) is. For the “passing” pattern as shown in Fig.10, the terminal vertical shift of particle increases with decreasing Wi, which is identical to the experimental result [19]. These indicate that the particles have a large migration range under large elastic forces. Figure 11 shows the velocity variation of particle 1 with time in the “returning” pattern, and the corresponding trajectory of particle is shown in the case of D=0.1 in Fig.10. The positive and negative sign on ordinate represents the direction of particle motion. We can see that, as shown in Fig.11(a), the vertical velocity increases when two particles move toward each other, and then decreases when two particles move in the opposite direction. The larger Wi is, the smaller the velocity variation of particle is. However, for the horizontal velocity as shown in Fig.11(b), the opposite is true. These indicate that the elastic force of fluid enhances the particle horizontal velocity, but reduces the vertical velocity under the combined effect of wall. The velocity variation of particle 1 with time in the “passing” pattern is shown in Fig.12 and the corresponding trajectory of particle is shown in the case of D=0.3 in Fig.10, which shows a more complex pattern than that in Fig.11. The vertical velocity of particle is zero in the initial and final stages as shown in Fig.12(a). In the period of 10
(a) vertical velocity (b) horizontal velocity Fig. 11.Velocity variation of particle 1 with time in the “returning” pattern for different W i(l=0.5,α=0.01, β=0.1).

(a) vertical velocity (b) horizontal velocity Fig. 12.Velocity variation of particle 1 with time in the “passing” pattern for different Wi (l=0.5,α=0.01, β=0.1).

For the case of D=0.2 in Fig.9, the trajectories of particle show the “returning” pattern and “passing” pattern for Wi=0.1 and 1.0, respectively, and the corresponding pressure contours at different times are shown in Fig.13 and 14, respectively. It can be seen that the pressure values are similar in the initial stage, but the difference is obvious in the follow-up stage. From the figures we can see the relationship between the particle trajectory and pressure distribution of the flow.

(a)

(b)

(c) Fig. 13.Pressure contours in the “returning” pattern for Wi=0.1(l=0.5, α=0.01, β=0.1, D=0.2).

(a)

(c)

(b)

(d)

(e) Fig. 14. Pressure contours in the “passing” pattern for Wi=1.0 (l=0.5, α=0.01, β=0.1, D=0.2).

In order to understand the physical behavior of particles in confined shear flow of viscoelastic fluid, we also examine the polymer extensions for the “returning” and “passing”. Figure 15 shows the instantaneous distribution of the polymer extension when particles behave the “passing”. Firstly, the strong elongational flow generated throughout the process of “passing”. Secondly, the polymer extension becomes larger as two particles approach each other. Additionally, there seem to be a typical orientation for the highly stretched regions. While as we can see in Fig.16, there is no obvious elongational flow and the polymer extension is very small and almost unchanged in the process of “returning”. What is more, a non-uniform and anisotropic microstructure becomes more pronounced in the flow with higher Weissenberg number by comparing Fig.15 and Fig.16.

(a)

(b)

(c)

(d)

(e) Fig. 15. Polymer extension tr(B) in the “passing” pattern for Wi=1.0 (l=0.5,α=0.01, β=0.1, D=0.2).

(a)

(b)

(c) Fig. 16. Polymer extension tr(B) in the “returning” pattern for Wi=0.1(l=0.5,α=0.01, β=0.1, D=0.2).

3.2.2. Shear-thinning effect As shown in Giesekus constitutive equation (6), the mobility parameter 𝛼 quantifies the shearthinning property of fluid. The shear-thinning effect becomes stronger with the increase of α. The trajectories of particle for different D and α are shown in Fig.17 where only the trajectories of particle 1 are given. We can see that the trajectories of particle show the “returning” pattern at D=0.1 for all values of α, and at D=0.2 for all values of α except α=0.01. However, the trajectories

of particle show the “passing” pattern at D=0.3 for all values of α except α=0.7. These illustrates that the trajectories of particle tend to the “returning” pattern as the shear-thinning effect is enhanced. We have calculated that the values of α corresponding to the separatrix are 0.05 at D=0.2, and 0.6 at D=0.3, respectively. For the “returning” pattern at D=0.2, the smaller α is, the shorter the minimum distance between the particle 1 and 2 is. For the “passing” pattern at D=0.3, the larger α is, the closer the particle is to the wall, i.e., the shear-thinning effect tends to make the particles move towards the wall, which is consistent with experimental results [35].

Fig.17. Trajectories of particle for different D and α (l=0.5, β=0.1, Wi=1.0).

Fig.18. Phase diagram for α and D (l=0.5, β=0.1, Wi=1.0).

The phase diagram for α and D is shown in Fig.18 where we can obtain the separatrix between “returning” and “passing” pattern. It can be seen that increasing α causes the value of critical D to increase, i.e., the trajectories of particle tend to the “returning” pattern as the shear-thinning effect is enhanced. Figure 19 shows the velocity variation of particle 1 with time in the “returning” pattern, and the corresponding trajectory of particle is shown in the case of D=0.1 in Fig.17. The positive and negative sign on ordinate represents the direction of particle motion. The vertical velocity first increases and then decreases as shown in Fig.19(a), the smaller α is, the smaller the velocity variation of particle is. However, for the horizontal velocity as shown in Fig.19(b), the opposite is true. The smaller α is, the larger the velocity variation of particle is. These indicate that the shearthinning effect reduces the particle horizontal velocity, but enhances the vertical velocity.

(a) vertical velocity (b) horizontal velocity Fig. 19. Velocity variation of particle 1 with time in the “returning” pattern for different α(l=0.5,Wi=1.0, β=0.1).

The velocity variation of particle 1 with time in the “passing” pattern is shown in Fig.20 and the corresponding trajectory of particle is shown in the case of D=0.3 in Fig.17. The vertical velocity of particle is zero in the initial and final stages as shown in Fig.20(a). In the intermediate stage, the vertical velocity increases to maximum and then decreases to zero, and then increases to maximum in the opposite direction and finally decreases to zero. The smaller α is, the earlier this process begins and ends, and smaller the maximum velocity is. The maximum vertical velocities are almost equal for different α. For the horizontal velocity as shown in Fig.20(b), the velocity also changes in a process of decrease-increase-decrease-increase. The smaller α is, the earlier this process begins and ends, and the larger the maximum velocity is. The shear-thinning effect has a larger impact on the “passing” pattern than the “returning” pattern by comparing Fig.19and 20.

(a) vertical velocity (b) horizontal velocity Fig. 20. Velocity variation of particle 1 with time in the “passing” pattern for different α (l=0.5,Wi=1.0, β=0.1).

Figure 21 shows the vertical velocity variation of particle 1 with time in the “returning” pattern

(α=0.3, D=0.2) and around the separatrix between “returning” and “passing” pattern (α=0.1,D=0.2 and α=0.3, D=0.25). It can be seen that the velocities change in a process of increase- decreaseincrease-decrease for the cases around the separatrix between “returning” and “passing” pattern, which is qualitatively different from that for the “returning” pattern. Therefore, the vertical velocity variation of particle can be used to predict the transition from the “returning” pattern to the “passing” pattern.

Fig. 21. Vertical velocity variation of particle 1 with time for different α and D(l=0.5,Wi=1.0, β=0.1).

3.3. Influence of other factors Different from the particle migration in Oldroyd-B fluid, the “tumbling” pattern is not observed in Giesekus fluid in this work. Actually, no matter in Newtonian [36] fluid or Oldroyd-B fluid [15, 21], the particles migrations are mainly in the form of “passing” pattern, while the “returning” and “tumbling” patterns appear seldom. As we can see from the work of Kulkarni and Morris [36], the Reynolds number has a great impact on the particle trajectories, especially for the “passing” pattern and “tumbling” pattern. The in-plane spiraling trajectories (“tumbling” pattern) have not been observed in [36] for Re≥0.1. It is worth noting that the Reynolds number in [36] is defined as 𝑅𝑒 = 𝜌𝛾̇ 𝑎2⁄𝜇 . According to this definition, the Reynolds number in the present work is 0.001. Similar phenomenon has been observed in the work of Bachelor and Green [37]. Specifically, the closed trajectory of particles can be observed at Re=0, and vanish at finite inertia. Some researchers [36, 38, 39, 40] attribute the vanishing of the closed trajectory to the loss of the fore-aft symmetry as Reynolds number increases. Meanwhile, the fore-aft asymmetric interactions appear for “passing” pattern at small Re. The wall confinement (H/a) also has a great impact on the particle trajectory of two particles. As

shown in the work of Yoon et al. [15], the vertical distance between two particles corresponding to the separatrix between “returning” and “passing” pattern increases when H/a is increased to 16, while the “returning” pattern is never observed when H/a is further increased to 32. Additionally, Vazquez-Quesada et al. [20] speculated that the weak confinement could promote “tumbling” pattern. Moreover, Vazquez-Quesada et al. [20] investigated the effect of collision model on the particle motion, and found that different values for a parameter in same collision model could produce different “passing” pattern. Nie et al. [16] studied the impact of different collision models on the particle trajectory, and observed the small differences in the particle trajectory when different collision strategies were used. Vazquez-Quesada et al. and Nie et al. concluded that the essential patterns of the particle trajectory are not altered, regardless the different parameters in one collision model or different collision models are used. Here, we further study the effect of contact force F0 on the particle trajectory. The values of F0 are selected as 10, 100 and 1000, respectively. In order to explore whether the “tumbling” pattern appears under the different values of F0, we choose other parameters (i.e., l=0.5, α=0.01, β=0.1, D=0.2) which corresponds to the state close to the transition between “returning” and “passing” pattern. As can be seen in Fig.22, the values of F0 have effect on the particle trajectory, but no effect on the pattern of particle trajectory.

Fig. 22.Trajectories of two particles at different F0 (l=0.5, Wi=1.0, β=0.1, α=0.01, D=0.2).

4. Conclusions The migrations of two particles in a bounded shear flow of Giesekus fluids are numerically investigated using the direct forcing/fictitious domain method. Attentions are mainly paid to the effect of Weissenberg number, the shear-thinning effect, the initial vertical distance between two particles on the particle migrations. The model is validated by comparing the numerical results with the available data in the literature. Some of the results are in agreement with the experimental ones. Conclusions are summarized as follow: When two particles are released in the flow with a certain initial horizontal and vertical distance, the trajectories of particle will show “returning” pattern and “passing” pattern, but do not show the “tumbling” pattern which appeared in Oldroyd-B fluid. The variations of particle velocity and pressure field in the “returning” pattern are totally different from those in the “passing” pattern. The separatrix between “returning” and “passing” pattern is dependent on the initial horizontal and vertical distance between two particles, the Weissenberg number and the shear-thinning effect. With other parameters fixed, the trajectories of particle change from the “returning” pattern to the “passing” pattern as the initial vertical distance (D) between two particles and the Weissenberg number (Wi) increase, but as the shear-thinning effect (α) decreases. The value of D corresponding to the separatrix between “returning” and “passing” is 0.17 under l=0.5, Wi=1.0, α=0.01, β=0.1, the value of Wi corresponding to the separatrix is 0.6 under l=0.5, D=0.2, α=0.01, β=0.1, and the values of α corresponding to the separatrix are 0.05 at D=0.2, and 0.6 at D=0.3, respectively, under l=0.5, Wi =1.0, β=0.1. Acknowledgements This work was supported by the Major Program of the National Natural Science Foundation of China with Grant Nos. 11632016 and 91634103.

Conflict of Interest None

References [1] J.Sroka, A. Kordecka, P. Włosiak, Z. Madeja,W. Korohoda, Separation methods for isolation of human polymer phonuclear leukocytes affect their motile activity, Eur. J. Cell Biol. 88(2009) 531-539. [2] Y.W.Kim, J.Y.Yoo, Three-dimensional focusing of red blood cells in microchannel flows for bio-sensing applications, Biosensors and Bioelectronics 24 (2009) 3677-3682. [3] G. D’Avino, P.L.Maffettone, Particle dynamics in viscoelastic liquids, J. Non-Newtonian Fluid Mech. 215(2015) 80-104. [4] G. D’Avino, F.Greco, P.L.Maffettone, Particle migration due to viscoelasticity of the suspending liquid and its relevance in microfluidic devices, Annu. Rev. Fluid Mech. 49 (2017) 341-360. [5] G.K.Batchelor, J.T.Green, The hydrodynamic interaction of two small freely-moving spheres in a linear flow field, J. Fluid Mech. 56 (1972) 375. [6] R.L.Hoffman, Discontinuous and dilatant viscosity behavior in concentrated suspensions. I. Observation of a flow instability, Trans. Soc. Rheol. 16 (1972)155-173. [7] J.J. Stickel, R.L.Powell, Fluid mechanics and rheology of dense suspensions, Annu. Rev. Fluid Mech. 37 (2005) 129-149. [8] J. Michele, R.Patzold, R.Donis, Alignment and aggregation effects in suspensions of spheres in non-Newtonian media, Rheol.Acta16(1977)317-321. [9] L.Petit, B. Noetinger, Shear-induced structures in macroscopic dispersions, Rheol. Acta 27 (1988) 437-441. [10] M.K.Lyon, D.W.Mead, R.E.Elliott, L.G.Leal, Structure formation in moderately concentrated viscoelastic suspensions in simple shear flow, J. Rheol. 45 (2001) 881-890. [11] D. Won, C. Kim, Alignment and aggregation of spherical particles in viscoelastic fluid under shear flow, J. Non-Newtonian Fluid Mech. 117 (2004) 141-146. [12] R. Scirocco, J. Vermant, J. Mewis, Effect of viscoelasticity of suspending fluid on structure formation in suspensions, J. Non-Newtonian Fluid Mech. 117 (2004) 183-192. [13] F.Gauthier, H.L.Goldsmith, S.G.Mason, Particle motions in non-Newtonian media. I: Couette flow, Rheol. Acta10 (1971) 344-364. [14] M. R. Hashemi, R. Fatehi, M.T. Manzari, SPH simulation of interacting solid bodies suspended in a shear flow of an Oldroyd-B fluid, J. Non-Newtonian Fluid Mech. 166 (2011) 1239-1252.

[15] S.Yoon, M.A.Walkley, O.G.Harlen, Two particle interactions in a confined viscoelastic fluid under shear, J. Non-Newtonian Fluid Mech. 185-186 (2012) 39-48. [16] D. Nie, J. Z. Lin, Behavior of three circular particles in a confined power-law fluid under shear, J.Non-Newtonian Fluid Mech. 221 (2015) 76-94. [17] W.R.Hwang, M.A.Hulsen, H.E.H. Meijer, Direct simulations of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames, J. Non-Newtonian Fluid Mech. 121 (2004) 15-33. [18] Y.J.Choi, M.A.Hulsen, H.E.H.Meijer, An extended finite element method for the simulation of particulate viscoelastic flows, J. Non-Newtonian Fluid Mech. 165 (2010) 607-624. [19] F. Snijkers, R. Pasquino, J. Vermant, Hydrodynamic interactions between two equally sized spheres in viscoelastic fluids in shear flow, Langmuir 29 (2013) 5701-5713. [20] A. Vazquez-Quesada and M. Ellero, SPH modeling and simulation of spherical particles interacting in a viscoelastic matrix, Phys. Fluids 29 (2017) 121609. [21] S.H.Chiu, T.W.Pan, R.Glowinski, A 3D DLM/FD method for simulating the motion of spheres in a bounded shear flow of Oldroyd-B fluids, Comput. Fluids172 (2018) 661-673. [22] S.Van Loon, J.Fransaer, C.Clasen, J.Vermant, String formation in sheared suspensions in rheologically complex media: the essential role of shear-thinning, J.Rheol. 58 (2014) 237-254. [23] N.O.Jaensson, M.A.Hulsen, P.D.Anderson, Direct numerical simulation of particle alignment in viscoelastic fluids, J. Non-Newtonian Fluid Mech. 235 (2016) 125-142. [24] Z.Yu, N.Phan-Thien, Y.R. Fan, R.I.Tanner, Viscoelastic mobility problem of a system of particles, J. Non-Newtonian Fluid Mech. 104 (2002) 87-124. [25] Z.Yu, A.Wachs, Y.Peysson, Numerical simulation of particle sedimentation in shear-thinning fluids with a fictitious domain method, J. Non-Newtonian Fluid Mech.136 (2006) 126-139. [26] Z.Yu, A.Wachs, A fictitious domain method for dynamic simulation of particle sedimentation in Bingham fluids, J. Non-Newtonian Fluid Mech. 145 (2007) 78-91. [27] R.Glowinski, T.W.Pan, T.I.Hesla, D.D.Joseph, A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiphase Flow 25 (1999) 755-794. [28] Z.Yu, X.Shao, A direct-forcing fictitious domain method for particulate flows, J. Comput. Phys. 227 (2007) 292-314. [29] B.V.Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to

Godunov's method, J. Comput. Phys. 32 (1979) 101-136. [30] Z.Yu, N.Phan-Thien, R.I.Tanner, Rotation of a spheroid in a Couette flow at moderate Reynolds numbers, Physical review E 76,026310(2007). [31] X.Shao, Z.Yu, B.Sun, Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers, Phys. Fluids 20,103307(2008). [32] P.Wang, Z.Yu, J.Lin, Numerical simulations of particle migration in rectangular channel flow of Giesekus viscoelastic fluids, J. Non-Newtonian Fluid Mech. 262 (2018) 142-148. [33] Z.Yu, P.Wang, J.Z. Lin, H.H.Hu, Equilibrium positions of the elasto-inertial particle migration in rectangular channel flow of Oldroyd-B viscoelastic fluids, J. Fluid Mech. 868 (2019) 316-340. [34] F.Snijkers, G.D’Avino, P.L.Maffettone, F.Greco, M.A.Hulsen, J.Vermant, Effect of viscoelasticity on the rotation of a sphere in shear flow, J. Non-Newtonian Fluid Mech. 166 (2011) 363-372. [35] F.DelGiudice, G.D’avino, F.Greco, P. A.Netti, P.Maffettone, Effect of fluid rheology on particle migration in a square-shaped microchannel, Microfluid. Nanofluid. 19 (2015) 1-10. [36] P. M. Kulkarni, J. F. Morris, Pair-sphere trajectories in finite-Reynolds-number shear flow, J. Fluid Mech. 596 (2008) 413-435. [37] G. K. Batchelor, J. T. Green, The hydrodynamic interaction of two small freely-moving spheres in a linear flow field, J. Fluid. Mech. 56 (1972) 375-400. [38] C. J. Lin, J. H. Peery, W. R. Schowalter, Simple shear flow round a rigid sphere: Inertial effects and suspension rheology, J. Fluid Mech. 44 (1970) 1-17. [39] G. Subramanian, D. L. Koch, Centrifugal forces alter streamline topology and greatly enhance the rate of heat and mass transfer from neutrally buoyant particles to a shear flow, Phys. Rev. Lett. 96(2006) 134503. [40] G. Subramanian, D. L. Koch, Inertial effects on the transfer of heat or mass from neutrally buoyant spheres in a steady linear velocity field, Phys. Fluids 18 (2006) 073302.