Numerical simulations of the separation of elastic particles in a T-shaped bifurcation

Accepted Manuscript

Numerical simulations of the separation of elastic particles in a T-shaped bifurcation Marco Trofa, Massimiliano Maria Villone, Gaetano D’Avino, Martien A. Hulsen, Paolo Antonio Netti, Pier Luca Maffettone PII: DOI: Reference:

S0377-0257(16)00028-8 10.1016/j.jnnfm.2016.01.015 JNNFM 3746

To appear in:

Journal of Non-Newtonian Fluid Mechanics

Received date: Revised date: Accepted date:

29 October 2015 16 January 2016 25 January 2016

Please cite this article as: Marco Trofa, Massimiliano Maria Villone, Gaetano D’Avino, Martien A. Hulsen, Paolo Antonio Netti, Pier Luca Maffettone, Numerical simulations of the separation of elastic particles in a T-shaped bifurcation, Journal of Non-Newtonian Fluid Mechanics (2016), doi: 10.1016/j.jnnfm.2016.01.015

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 • Elastic particle dynamics in a T-shaped bifurcation is studied by nu-

CR IP T

merical simulations.

• A separation is possible based on particle elasticity.

AC

CE

PT

ED

M

AN US

• Critical elasticity number is related to particle confinement.

1

ACCEPTED MANUSCRIPT

Marco Trofa∗

CR IP T

Numerical simulations of the separation of elastic particles in a T-shaped bifurcation Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Universit` a di Napoli Federico II, P.le Tecchio 80, 80125 Napoli, Italy Center for Advanced Biomaterials for Healthcare @CRIB, Istituto Italiano di Tecnologia, Largo Barsanti e Matteucci 53, 80125 Napoli, Italy

Massimiliano Maria Villone

AN US

Center for Advanced Biomaterials for Healthcare @CRIB, Istituto Italiano di Tecnologia, Largo Barsanti e Matteucci 53, 80125 Napoli, Italy

Gaetano D’Avino

Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Universit` a di Napoli Federico II, P.le Tecchio 80, 80125 Napoli, Italy

M

Martien A. Hulsen

Department of Mechanical Engineering, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands

ED

Paolo Antonio Netti

PT

Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Universit` a di Napoli Federico II, P.le Tecchio 80, 80125 Napoli, Italy Center for Advanced Biomaterials for Healthcare @CRIB, Istituto Italiano di Tecnologia, Largo Barsanti e Matteucci 53, 80125 Napoli, Italy

Pier Luca Maffettone

CE

Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Universit` a di Napoli Federico II, P.le Tecchio 80, 80125 Napoli, Italy

∗

AC

Corresponding author. Tel.: +39 0817682280 Email addresses: [email protected] (Marco Trofa), [email protected] (Massimiliano Maria Villone), [email protected] (Gaetano D’Avino), [email protected] (Martien A. Hulsen), [email protected] (Paolo Antonio Netti), [email protected] (Pier Luca Maffettone)

Preprint submitted to Journal of Non-Newtonian Fluid Mechanics

February 1, 2016

ACCEPTED MANUSCRIPT

Abstract

CR IP T

In this work, we study the dynamics of an elastic particle suspended in an inertialess Newtonian fluid flowing through a channel with an orthogonal

side branch (asymmetric T-shaped bifurcation) by means of 2D Arbitrary Lagrangian Eulerian (ALE) Finite Element Method (FEM) direct numerical simulations.

AN US

The simulations show that, when the fluid is equally partitioned in the

two downstream branches and the particle starts from the inflow channel centerline, a sufficiently deformable particle is deviated into the ‘side’ branch, at variance with a rigid particle, which remains in the ‘main’ outlet. The effects of the elastic capillary number and the confinement ratio on the

M

particle trajectory and deformation near the bifurcation are investigated. We

their elasticity.

ED

discuss how this device can be exploited for separating particles based on

Keywords: Suspension, Separation, Bifurcation, T-shaped channel, Elastic

PT

particle, Direct Numerical Simulations

CE

1. Introduction

Suspension flow through channel bifurcations has many important appli-

AC

cations in biological and technological processes. Several works deal with microfluidic devices consisting of a main channel with side branches, aimed at separating particles with different properties (see, for example, [1] and the references therein), either for design or analytical purposes [2]. Indeed, it is convenient to preventively separate and concentrate rare target particles 3

ACCEPTED MANUSCRIPT

before analysis. In general, suspensions made of particles with homogeneous properties are desirable. From the biological side, a clear example of suspen-

CR IP T

sion flow in bifurcations is represented by the peripheral circulatory system: here blood, which can be considered as a concentrated suspension of cells

(RBCs, WCs and platelets), flows in a network made of many capillaries. Moreover, the partitioning of blood constituents at vessel junctions plays an

important role in determining the microvessel hematocrit and has relevant

AN US

physiological consequences [3, 4].

Experimental and numerical studies have demonstrated that, when a suspension of rigid, non-Brownian particles with a diameter comparable to the vessel transversal dimension flows through an asymmetric bifurcation under inertialess conditions, the partitioning of particles between the two

M

downstream branches is different from the partitioning of the suspending fluid [5, 6, 7]. This phenomenon has been qualitatively explained by the

ED

presence of excluded volume effects, due to the existence of a layer adjacent to the channel wall (a particle radius in thickness) unaccessible to the particle

PT

centers, and particle-wall hydrodynamic interactions, which determine a difference, both in magnitude and direction, between the particle and fluid

CE

velocities [6]. Both contributions become more and more relevant as the particle confinement increases. Two relevant parameters that affect the distribution of particles between

AC

the two downstream branches are the bifurcation geometry (in particular the angle between the branches) and the volume fraction of the suspended particles [7]. For Y-shaped bifurcations, the branch with a greater volumetric flow rate receives a higher particle volume fraction as compared to the main

4

ACCEPTED MANUSCRIPT

channel, and this trend increases with particle size. For non-symmetric bifurcations, as in T-shaped or oblique channels, the particles preferentially move

CR IP T

toward the main outlet rather than the side branch. This effect is enhanced by increasing the branching angle from 90◦ to 135◦ [6]. In concentrated suspensions, the interparticle interactions become dominant and, at low confinement,

particles are almost equally partitioned between downstream branches, regardless of an unequal division of the suspending fluid [8]. The nature of

AN US

the suspending fluid also has a relevant effect on the particle partitioning. A non-Newtonian fluid, indeed, can significantly affect the particle dynamics through the bifurcation, due to the existence of normal stresses that alter the shear rate distribution around the particles and promote a migration mechanism toward the channel centerline [9, 10].

M

Several studies have also investigated the role of particle deformability on the fluid-particle partitioning [11, 12, 13]. Indeed, the ability to separate

ED

particles based on their mechanical properties is of primary importance in medicine, since cell deformability is a clinical indicator of a wide range of

PT

diseases, e.g., malaria [14] and leukemia [15]. Blood is, then, a preferential sample, because it is composed by a heterogeneous mixture of cells that

CE

continuously respond to the physiological and pathological changes of human body. In this regard, another interesting target for this selective separation are the circulating tumor cells (CTCs), which are responsible for blood-borne

AC

metastasis and most cancer-related deaths [16]. Results for drops, red blood cells, and capsules in bifurcations are available

in the literature. The case of a drop suspended in a fluid flowing through a bifurcating channel has been numerically studied by Manga [11] by means

5

ACCEPTED MANUSCRIPT

of the boundary integral method. The author found that the drop has an increasing tendency to enter the branch at higher flow rate for low (drop to

CR IP T

fluid) viscosity ratios, high flow velocities, and high confinements. Secomb et al. [17, 12] used the finite element method to compute the motion of a

red blood cell (modeled as a viscoelastic drop) at a bifurcation. Despite the two-dimensional approximation, their calculations succeed in predicting

the experimentally observed RBC behavior, including tank-treading and cell

AN US

elongation. They found that cell trajectories deviate from the background flow streamlines due to migration toward the channel centerline upstream of the bifurcation, and due to flow perturbations caused by cell obstruction in the neighborhood of the bifurcation. In addition, flexible cells determine smaller obstruction as compared to rigid particles because they tend to flatten against

M

the flow divider. Finally, Woolfenden et al. [13] investigated numerically the behavior of a two-dimensional elastic fluid-filled capsule through the boundary

ED

integral method. They found that a capsule moving through a channel with a side branch is drawn out of the main channel when the flow rate in the side

PT

branch is sufficiently high. For equal flow rates but different branch widths, capsules tend to enter the branch with the highest centerline velocity. When

CE

the oulet channels are also of equal width, there is no significant bias toward one exit or the other. In this case, the capsule deformation may play a crucial role in determining the ultimate route taken, i.e., capsules with different

AC

elastic properties may follow different routes out of the junction in identical flow conditions. They also found that the deformation experienced by the capsules in the region of the junction strongly depends on the bifurcation angle.

6

ACCEPTED MANUSCRIPT

It so appears that the dynamics of particles in a branching channel has been already addressed for rigid particles, drops, and capsules. However, to

CR IP T

the best of our knowledge, a study on homogeneous elastic particles is missing, which is the aim of the present work. In this regard, it has to be pointed

out the importance of these systems as models for microgels, broadly used in drug delivery and cell encapsulation [18].

In this work, we carry out a 2D numerical investigation on the dynamics

AN US

of a dilute suspension of elastic particles in a Newtonian fluid flowing through a channel with an orthogonal side branch. In particular, we focus on the trajectory and deformation experienced by the particles during the whole path from the main inlet channel to one of the two outlet branches. Since the suspension is assumed to be dilute (typical in vitro conditions), single

M

particle simulations suffice to describe the suspension dynamics. Like drops and capsules, also elastic particles have been found to migrate toward the

ED

channel centerline in a pressure-driven flow [19]. Therefore, all the simulations have been performed by releasing the particle with its center lying on the

PT

centerline of the inflow channel. Results for different capillary numbers and confinement ratios are presented and discussed.

CE

2. Governing equations

AC

In this work, we consider a 2D ‘T-shaped’ channel with one inlet and two

outlets, as depicted in Fig. 1. The two outflow branches are orthogonal and all the sections of the device have height H. The outflow branch aligned with the inflow one will be referred to as the ‘main branch’ and the other outlet, orthogonal to the main branch, will be denoted as ‘side branch’. 7

ACCEPTED MANUSCRIPT

Γw ∂P(t)

𝑥

Γ0

Ω

H

P(t)

Γ1

-0.5

L0

CR IP T

𝑦 0.5

L1

-1.5

L2

R

(b) -2.5 1

2

3

4

Γ2

5

(a)

6

7

AN US

0

90°

Figure 1: Computational domain and mesh used in the simulations with β = 0.4. (a) Zoom of the smooth corners, with curvature radius R = 0.1H, between the main and the side

M

channel. (b) Refined mesh when the particle is close to the bifurcation corner.

An initially circular elastic particle with diameter dp is suspended in a

ED

Newtonian liquid entering the channel with average velocity u¯0 . We denote by Ω the union of the liquid and solid domains. As the particle moves due to the fluid flow, we denote by P (t) the time-dependent particle domain,

PT

and by ∂P (t) its boundary. Consequently, Ω − P (t) is the fluid domain. We denote by Γ0 the inflow section, by Γ1 and Γ2 the main and side outflow

CE

sections, respectively, and by Γw the channel walls. A Cartesian reference frame is set with its origin at the middle of the inlet boundary Γ0 and the

AC

x-axis orthogonal to it (see Fig. 1). As illustrated in inset (a) in Fig. 1, the corners between the main and the side channel are rounded in order to avoid numerical problems (see below). The curvature radius of such smooth corners is R = 0.1H. Due to the characteristic flow velocities and length scales occurring in 8

ACCEPTED MANUSCRIPT

microfluidic devices, the Reynolds number is generally very small [20], thus inertia is here neglected.

CR IP T

Although the use of a displacement-based formulation for solids is quite common, we apply the velocity-based approach for an elastic solid as proposed

in [21] and [22], which consists in considering the elastic solid particle as a droplet of a viscoelastic fluid with infinite relaxation time (see below). In such

a way, we maintain a conforming mesh across the interface and we can apply

AN US

a mesh updating scheme (ALE) that filters out the tank treading motion, thus avoiding big mesh distortion within a short time.

Assuming also incompressibility of both the suspending and the suspended phase, the continuity and momentum balance equations governing the dynamics of both the liquid and the solid read

(1)

∇·σ =0

(2)

ED

M

∇·u=0

with u the velocity vector and σ the total stress tensor. For σ we can write

CE

PT

σ = −pI + 2ηD σ = −pI + τ

on Ω − P (t) on P (t)

(3) (4)

where p is a pressure, I is the identity tensor, η is the viscosity of the

AC

suspending liquid, D = (∇u + (∇u)T )/2 is the symmetric part of the velocity gradient tensor, and τ is the extra-stress in the elastic solid, for which we use the upper-convected Maxwell constitutive equation ∇

τ +

τ = 2GD λ 9

(5)

ACCEPTED MANUSCRIPT

In Eq. (5), G is the shear modulus of the elastic material and λ is the relaxation time, but, since this tends to infinity, the term τ /λ is negligible

CR IP T

and the constitutive equation for the solid reduces to ∇

τ = 2GD

(6)

which is the neo-Hookean constitutive equation in terms of the velocity

gradient. The symbol (∇ ) denotes the upper-convected time derivative,

AN US

defined, as usual, as

∂τ + u · ∇τ − (∇u)T · τ − τ · ∇u ∂t

∇

τ≡

(7)

For computational reasons we prefer to express the extra-stress in terms of the conformation tensor c defined through the relation

M

τ = G(c − I)

(8)

ED

from which the neo-Hookean constitutive equation becomes ∇

c= 0

(9)

PT

The continuity and momentum balance equations for the suspending medium and the suspended particle are solved with the following boundary

AC

CE

conditions on the channel walls, inlet and outlets

u=

"

3 u¯0 2

u = 0 on Γw #  2 ! 2y 1− , 0, 0 H

(10) on Γ0

(11)

σ · m = −p1 m

on Γ1

(12)

σ · m = −p2 m

on Γ2

(13)

10

ACCEPTED MANUSCRIPT

Equation (10) is the no-slip condition for the suspending fluid on the channel walls. Equation (11) imposes the well-known parabolic velocity profile for a

CR IP T

Newtonian fluid in a straight channel given by Poiseuille law [23]. Finally, Eqs. (12) and (13) are the outflow conditions for the fluid in the main and

the side branch, respectively, with m the outwardly directed unit vector normal to the outflow sections, and p1 , p2 the ambient pressures that the

fluid ‘feels’ when crossing Γ1 and Γ2 , respectively. Actually, more relevant

AN US

than the pressure values p1 and p2 is the pressure difference ∆p = p1 − p2 , which sets the relative weight of the two outlet flow rates. Fixed ∆p, different p1 - and p2 -values only shift the pressure field of the fluid in the channel. The boundary conditions on the particle-matrix interface ∂P (t) read (14)

(σ · n)|m = (σ · n)|p

(15)

M

u|m = u|p

ED

where n is the unit vector normal to the particle-matrix interface and directed toward the suspending liquid. Equations (14) and (15) express the continuity

PT

of velocity and traction across the above mentioned interface. Since both the particle and the suspending medium are inertialess, no

CE

initial conditions on the velocities are required. Only the initial condition for the conformation tensor in the solid needs to be specified, and we assume

AC

that the particle is initially stress-free, or, in symbols, c|t=0 = I

(16)

The governing equations are made dimensionless by choosing H as the

characteristic length, u¯0 as the characteristic velocity, H/¯ u0 as the characteristic time, and η u¯0 /H as the characteristic stress. Then, the following 11

ACCEPTED MANUSCRIPT

dimensionless parameters will appear in the equations: the confinement ratio β = dp /H, the elastic capillary number Ca = η u¯0 /(HG), and the dimension-

CR IP T

less pressure difference ∆P = ∆p/(η u¯0 /H). All the results presented in the following are given in terms of dimensionless quantities, the usual superscript ∗

being omitted for the sake of simplicity.

3. Numerical method and code validation

AN US

3.1. Numerical method

The governing equations presented in section 2 are solved through the Arbitrary Lagrangian Eulerian (ALE) Finite Element Method (FEM). We adopt the log-conformation representation (LCR) [24], which is a transformation of the original equation for the conformation tensor c given by Eq. (9)

M

to an equivalent equation for s = log c

∂s + u · ∇s + h(∇u, s) = 0 ∂t

ED

(17)

where the symbol (., .) denotes the inner product. Details on the expression for

PT

h(∇u, s) can be found in [24]. Solving the equation for s instead of the one for c leads to major stability improvements. Furthermore, we use the Streamline-

CE

Upwind/Petrov-Galerkin (SUPG) technique [25] for the constitutive equation. Hence, the weak form of the balance equations and the constitutive equation

AC

is

(q, ∇ · u) = 0 ((∇v)T , τ ) − (∇ · v, p) = 0

 w + τ (u − um ) · ∇w,

δs + (u − um ) · ∇s + h(∇u, s) δt 12

(18) 

(19) =0

(20)

ACCEPTED MANUSCRIPT

with q, v, and w the test functions for the pressure p, velocity u and the log-conformation tensor s, respectively. In Eq. (20), τ is the SUPG parameter,

CR IP T

for which we take τ = h/2U , with h a characteristic element length and U the magnitude of the velocity vector u, evaluated in each integration point separately. In the same equation, um is the velocity of the mesh nodes

and δ/δt the grid derivative, both coming from the ALE formulation [26]. Furthermore, proper inner products (., .) are defined on the full domain Ω.

AN US

The continuity and momentum balance equations are decoupled from the

constitutive equation through the method of D’Avino and Hulsen [27], where an implicit stress formulation is developed. The discretized continuity and momentum balance equations read

M

∇ · un+1 = 0

∇ · τ n+1 − ∇pn+1 = 0

(21) (22)

ED

In Eq. (22), the extra-stress tensor term is replaced by the space-continuous but time-discretized form of the constitutive equation through an Euler scheme

PT

implicit in the velocity [28]. For this, we evaluate the weak form at the new

CE

time step tn+1 and replace the extra-stress τ in Eq. (22) with

AC

τ n+1 = G(cn − I)+   G∆t − (un+1 − um,n+1 ) · ∇cn + (∇un+1 )T · cn + cn · ∇un+1 (23)

As a result, Eqs. (21) and (22) are a Stokes-like system in the velocity field un+1 and pressure field pn+1 , both at the new time step tn+1 . After computing un+1 and pn+1 , the log-conformation field sn+1 is obtained by Eq. (20) discretized by a first or second-order scheme based on backwards 13

ACCEPTED MANUSCRIPT

differencing (Gear) 1st order (only for the first step)

w + τ (un+1 − um,n+1 ) · ∇w,

sn+1 − sn ∆t



ˆn+1 ) = 0 (25) + (un+1 − um,n+1 ) · ∇sn+1 + h(∇un+1 , s

2nd order

AN US



(24)

CR IP T

ˆn+1 = sn s

ˆn+1 = 2sn − sn−1 s w + τ (un+1 − um,n+1 ) · ∇w,

3 s 2 n+1

− 2sn + 12 sn−1 ∆t

 ˆ + (un+1 − um,n+1 ) · ∇sn+1 + h(∇un+1 , sn+1 ) = 0 (27)

M



(26)

ˆn+1 denotes the prediction of s at the new timestep. where s

ED

Finally the conformation tensor cn+1 is obtained by projection

PT

with e a test function for c.


e, cn+1 − exp(sn+1 )

(28)

The adopted implicit-stress formulation is second-order accurate in time

CE

and gives convergent solutions without using other stabilization techniques [29], e.g., the DEVSS method [30], thus leading to a significant reduction of

AC

the computational time. Since both the suspending liquid and the suspended particle are numeri-

cally treated as fluids, the interface between them needs to be tracked: to do this, a FEM with second-order time discretization is defined on the interface. Here, the normal component of the mesh velocity equals the normal 14

ACCEPTED MANUSCRIPT

component of the physical velocity, whereas the tangential velocity of the mesh nodes is such that the distribution of the elements on the interface is

CR IP T

optimized. This approach lets the mesh get rid of the tank-treading motion of the particle, thus greatly reducing the distortion of the ALE mesh as compared to a Lagrangian description of the interface. A detailed description

of the approach used to track the interface, with several test applications, can be found in [22].

AN US

Both the fluid and the solid domains are discretized through a mesh of quadratic triangles, an example being shown in Fig. 1. Since the particle is

transported along the device by the suspending liquid, the mesh elements progressively warp. Hence, every time it overcomes a distortion threshold, the

the new mesh [26, 31].

M

mesh needs to be regenerated, and the physical fields have to be projected on

As hinted above, the corners between the orthogonal outlets are smooth,

ED

with curvature radius R = 0.1H. This avoids the occurrence of flow singularities due to sharp corners. Moreover, since great deformations are expected

PT

when the particle approaches the bifurcation corner, refinement steps are introduced. At every refinement, each of the elements on the particle-matrix

CE

interface and on the bifurcation corner is divided in two in order to ensure the presence of a minimum number of mesh elements in the gap between the

AC

particle and the corner (see inset (b) in Fig. 1). 3.2. Code validation In the computational domain, the length of the inflow branch has been

chosen sufficiently larger than the particle diameter to allow the particle assume its equilibrium deformed shape before ‘feeling’ the effect of the bi15

ACCEPTED MANUSCRIPT

furcation. In particular, we have verified that a length of the inflow branch L0 = 4H satisfies the above mentioned condition for every dp considered. In

CR IP T

addition, the length of the outflow branches is chosen sufficiently larger than the particle diameter to avoid any spurious effects on the particle dynamics

in the bifurcation region due to the imposed outflow boundary conditions. A length of the outflow branches L1 = L2 = 2H fulfills this requirement.

Mesh and time convergence are checked for all the simulation results

AN US

reported in this work. The results of a typical convergence test are reported in Fig. 2. Figure 2a shows the trajectories of a particle (with reference to its

center of volume, dashed lines) and some snapshots of its shape (solid lines) for β = 0.6, Ca = 0.07, ∆P = −0.4, and for three different combinations of mesh resolution and time step. In particular, the blue curves refer to

M

mesh M1, with 40 elements initially present on the particle boundary, and ∆t = 0.002, the red curves refer to mesh M1 and ∆t = 0.001, and the green

ED

curves refer to mesh M2, with 100 elements initially on ∂P , and ∆t = 0.002. In all the three cases, the mesh is refined twice as the particle approaches the

PT

bifurcation corner. Figure 2b shows the temporal trends of the variation of particle perimeter ∆Λ with respect to the initial one for the same values of

CE

the parameters as in panel (a). (For the definition of ∆Λ, see section 4.) From both Figs. 2a and b, a fair superposition of the three indicators considered (i.e., trajectory, shape, and perimeter variation) makes us conclude that, for

AC

β = 0.6, Ca = 0.07, and ∆P = −0.4, a mesh with 40 elements on the particle boundary and a time step ∆t = 0.002 is sufficient to reach mesh and time step convergence of the numerical results. For the other β-, Ca-, and ∆P -values considered in this work, analogous tests have been performed.

16

ACCEPTED MANUSCRIPT

0.5

0.0

CR IP T

◆ ● ■ ■ ● ◆

◆ ● ■

-0.5

-1.0

●

Mesh M1, Δt = 0.002

■

Mesh M1, Δt = 0.001

◆

Mesh M2, Δt = 0.002

3.0

3.5

AN US

◆ ● ■

4.0

4.5

5.0

(a)

25

ED

CE

PT

ΔΛ%

M

○○○○○○○○○○ ○○◆ ◆ ◆ ◆ ◆ ○○ ○◆ ◆ ○ ○ ○◆ ◆ ○ ○ Mesh M1, Δt = 0.002 ○◆ 20 ○ ◆ ○ ○ ◆ ○ Mesh M1, Δt = 0.001 ○ ○ ◆ ○ ◆ ◆ Mesh M2, Δt = 0.002 15 ○ ○ ◆ ○ 10 ○ ◆ ○ ○ ◆ ○ ○ 5 ○○◆ ◆ ○ ○ ○◆ ○ ◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○◆ ○○ ○○◆ ○○◆ ◆ 0◆ ○○

AC

0

1

2

3

4

5

6

t (b)

Figure 2: (a) Trajectories (dashed lines) and shapes (solid lines) of a particle for β = 0.6, Ca = 0.07, ∆P = −0.4, and three different combinations of mesh resolution and time step.

The trajectories refer to the particle center of volume. (b) Relative variation of particle perimeter as a function of time for the same β, Ca, ∆P , mesh resolutions and time steps as in panel (a).

17

ACCEPTED MANUSCRIPT

About mesh convergence, it is worth remarking that, for very confined geometries, i.e., for β > 0.6, higher curvatures than the ones shown in Fig.

CR IP T

2a are detected, which makes it necessary to employ finer initial meshes, with about 100 elements on ∂P . However, when this is the case, only one refinement step in the vicinity of the corner turns out to be sufficient.

For what matters time convergence, it has been found that, for low confinements, namely, for β ≤ 0.2, a ‘deflation’ of the suspended particle

AN US

happens for ∆t = 0.002. In other words, under these conditions the particle

area decreases in time. Hence, for low β-values, ∆t has to be reduced down to 5 · 10−4 to achieve convergent numerical results. 4. Results

M

In this work, we investigate the effect of particle elasticity on the dynamics (deformation and trajectory) of a particle suspended in an inertialess

ED

Newtonian fluid flowing through a channel with an orthogonal side branch. In all the simulations, a single circular particle is initially placed with its

PT

center of volume on the centerline of the inlet channel. We select this configuration because it is known from the literature that elastic particles in

CE

pressure-driven flow naturally migrate toward the channel centerline due to their deformability [19]. Furthermore, as remarked in the previous section,

AC

the length L0 of the main inlet channel is selected sufficiently great in order to assure that the particle achieves a steady-state shape before entering the bifurcation. Therefore, under these assumptions, the inlet channel can be seen as the final part of a very long straight channel. All the simulations are carried out by fixing the outlet pressure difference 18

ACCEPTED MANUSCRIPT

∆P such that flow rates through the outlet sections are the same. This implies that, in a fluid without particles, the streamline corresponding to the inlet

CR IP T

channel centerline ends up on the bifurcation corner. Such a streamline divides the fluid in two streams equally partitioned between the two branches. Under these conditions, a rigid particle is driven into the main branch, regardless of the confinement [10]. One aim of this work is to investigate whether particle

deformability may alter this scenario and deviate the particle toward the side

AN US

branch.

We present simulation results by varying the capillary number Ca between 0.01 and 0.1 (order of magnitude of the maximum value experimentally achievable in microfluidics) and the confinement ratio β in the range 0.1 - 0.8. Such a β-range has been selected because at very low confinement ratios

M

(β < 0.1), i.e., when the particle is very small as compared to the channel, the particle simply follows the fluid streamline, and the effect of the deformation

ED

(on migration and separation at the corner) is negligible. On the other hand, at very high confinement (β > 0.8), the simulations fail to converge due to the

PT

extreme curvature at particle edges close to the channel walls. In particular, in the following we show in detail particle dynamics, in terms of deformation

CE

and trajectories, for β = 0.2 and β = 0.6. Let us start with the low confinement ratio β = 0.2. Figure 3a shows the

trajectories and the shapes of a particle for capillary numbers ranging from

AC

0.01 to 0.05. The dashed lines indicate the position of the particle center of volume. We observe that the trajectories superimpose until the particle is very close to the corner. In the vicinity of the corner, depending on the capillary number, the particle can enter either the main or the side branch.

19

ACCEPTED MANUSCRIPT

Specifically, for low capillary numbers (i.e., low deformability and/or low flow rate), the particle is driven to the main branch, in agreement with the rigid

CR IP T

particle case [10]. On the other hand, for a sufficiently high capillary number (i.e., high deformability and/or high flow rate), the particle deviates to the

side branch. Hence, particle deformability promotes the deviation to the side branch, and a critical capillary number exists such that the particle switches from one outlet to the other. Similarly to the trajectories, the shapes also are

AN US

undistinguishable until the particle gets very close to the corner. Notice also

that the shape in the inflow channel remains essentially circular due to quite low confinement ratio. However, as the particle almost touches the corner, different shapes are clearly observed as the capillary number is varied. This is evident by looking at Fig. 4 where the shapes of the particle at different

M

times are shown. The selected times correspond to the positions indicated by the symbols in Fig. 3a. In each panel of Fig. 4, the shapes are superimposed

ED

such that the centers of volume coincide. While very slight deviations are noted in Fig. 4a (particle far from the corner), the shapes gradually depart

PT

from the circle as the particle approaches the corner, showing elliptical or bean-like shapes. As expected, the most deformed shape is observed for the

CE

highest capillary number and for the minimum particle-corner distance (black line in Fig. 4d). To quantify the deformation experienced by the particle, we evaluate the

AC

time evolution of the variation of its perimeter with respect to the initial circular shape. The relative perimeter variation is evaluated as ∆Λ =

Λ(t) − πdp πdp

(29)

where Λ(t) is the (time-dependent) perimeter of the deformed shape. Recalling 20

ACCEPTED MANUSCRIPT

0.2

0.0

-0.2 ● ■ ◆ ▲ ▼ -0.4

CR IP T

▼ ▲ ◆ ● ■

● ■ ● ◆ ▲ ▼ ■ ◆ ▼▲

3.6

3.8

AN US

-0.6 4.0

4.2

4.4

4.6

4.8

5.0

(a)

▼

6

ΔΛ%

3

■

Ca = 0.02

◆

Ca = 0.03

▲

Ca = 0.04

▼

Ca = 0.05

PT

2

Ca = 0.01

M

4

●

ED

5

1

CE

0

AC

0

1

▼

▲

■ ■

◆ ■ ●

2

◆

◆

▼ ▲

▼ ▲ ◆ ■ ●

▲

●

● 3

4

5

6

t (b)

Figure 3: (a) Trajectories (dashed lines) and shapes (solid lines) of a particle with β = 0.2 and Ca between 0.01 and 0.05. The dashed lines refer to the particle center of volume. (b) Relative variation of particle perimeter as a function of time.

21

0.10

0.10

0.05

0.05

a) t = 2.0

0.00

-0.05

-0.05

-0.10

-0.10

-0.10 -0.05

0.00

0.05

0.10

-0.10 -0.05

ED

0.05

PT

-0.05

CE

-0.10 -0.05

0.00

0.05

0.10

0.05

0.10

0.05

c) t = 3.6

0.00

0.00

0.10

M

0.10

-0.10

b) t = 2.8

AN US

0.00

CR IP T

ACCEPTED MANUSCRIPT

d) t = 4.4

0.00 -0.05 -0.10

0.05

0.10

-0.10 -0.05

0.00

Figure 4: Particle shapes for β = 0.2 and different time steps (the same shown in Fig 3a).

AC

For the legend see Fig 3b.

22

ACCEPTED MANUSCRIPT

that the area of the particle does not change, high ∆Λ-values imply large deviations from the circular shape. In Fig. 3b, the trends of the relative

CR IP T

perimeter variation for the just examined capillary numbers are shown. The symbols on each curve denote the same time-values as in Fig. 3a, i.e., the ones appearing in Fig. 4. First of all, we notice that no perimeter variation occurs as the particle moves through the inlet, confirming that the shape

remains essentially circular. As the particle approaches the corner, the curves

AN US

increase regardless of the capillary number. Of course, higher Ca-values promote stronger perimeter variations. All the curves achieve a maximum and then decrease. This is due to the fact that, after attaining the maximum deformation, the particle enters one of the two outlet branches and the

deformation is progressively reduced. Finally, notice that the green curve,

M

corresponding to Ca = 0.03, shows a smoother variation near the maximum because, for such a capillary number, the particle spends more time near the

ED

corner as compared to lower or higher Ca-values. Indeed, Ca = 0.03 is close to the above mentioned critical capillary number that ‘discriminates’ particles

PT

going to the main or the side branch. We now move to a more confined case. Figure 5a reports the trajectories

CE

described by the particle center of volume for β = 0.6 and Ca ranging from 0.01 to 0.09, along with the shapes assumed by the particle. These shapes are also shown in Fig. 6 for the times corresponding to the symbols in Fig. 5a.

AC

Similarly to the previous case, the shapes in this figure are centered with respect to the particle center of volume. The results reported in Figs. 5 and 6 show that: (i) the trajectories start to separate as soon as the particle deviates from the rectilinear path in the main flow channel; (ii) even when

23

ACCEPTED MANUSCRIPT

0.5

0.0

CR IP T

■▲ ◆ ● ▼ ● ■ ◆ ▲ ▼

●

■ ● ■ ◆ ▼▲ ◆ ▲ ▼

AN US

-0.5

-1.0 3.0

3.5

4.0

4.5

5.0

(a)

▼

15

Ca = 0.01

■

Ca = 0.03

◆

Ca = 0.05

▲

Ca = 0.07

▼

Ca = 0.09

PT

ΔΛ%

20

●

10

CE AC

1

▼ ◆

▲ ◆

5

0 0

▲

ED

25

M

30

▼ ▲ ◆ ■ ● 2

▼ ▲ ◆ ■ ●

■

■ ● 3

● 4

5

6

t (b)

Figure 5: (a) Trajectories (dashed lines) and shapes (solid lines) of a particle with β = 0.6 and Ca between 0.01 and 0.09. The dashed lines refer to the particle center of volume. (b) Relative variation of particle perimeter as a function of time.

24

0.4

0.4

0.2

0.2

a) t = 2.0

-0.2

-0.2

-0.4

-0.4

-0.4

-0.2

0.0

0.4

-0.4

-0.2

0.0

0.2

0.4

0.2

0.4

d) t = 4.4

0.0

PT

CE -0.4

0.0

0.2

c) t = 3.6

0.0

-0.2

0.4

ED

0.2

-0.4

0.2

M

0.4

-0.2

b) t = 2.8

0.0

AN US

0.0

CR IP T

ACCEPTED MANUSCRIPT

-0.2

-0.4 0.2

0.4

-0.4

-0.2

0.0

AC

Figure 6: Particle shapes for β = 0.6 and different time steps (the same shown in Fig 5a). For the legend see Fig 5b.

25

ACCEPTED MANUSCRIPT

the particle travels in the main flow channel, a non-circular shape can be detected for any investigated capillary number; for instance, a bullet-like

CR IP T

shape is visible for the highest Ca-value; (iii) near the corner, a remarkable deformation with a concavity toward the corner is observed due to the flow

elongational components, even for the lowest capillary number (blue lines); (iv) as for the case with β = 0.2, at high capillary numbers (orange and black

curves), the particle ends up in the side branch, whereas, at low Ca-values, it

AN US

moves to the main outflow branch; hence, also in this case a critical capillary number (between 0.05 and 0.07) exists. Of course, the stronger deformations observed for β = 0.6 as compared to the case at β = 0.2 are due to the more significant hydrodynamic interactions between the particle and the channel walls, resulting in higher shear rate gradients around the particle. All these

M

observations are reflected in the trends of the relative perimeter variation reported in Fig. 5b. First of all, notice that, during the initial times, the

ED

curves for Ca > 0.03 increase (≈ 1 − 2%) due to the deformation experienced by the particle even when traveling within the inlet far from the bifurcation.

PT

However, after this initial deformation, a constant trend can be identified for any Ca-value, i.e., the particle has reached a steady-state shape, confirming

CE

that the main channel length has been correctly selected. As soon as the particle approaches the bifurcation, the perimeter slightly goes down, achieves a minimum, then it remarkably increases. As for the case with β = 0.2, the

AC

deviation from the circular shape is larger for higher capillary numbers, and the maximum deformation occurs when the particle is almost in contact with the corner. After that, the center of volume moves in one of the two outlets and the deformation is reduced.

26

ACCEPTED MANUSCRIPT

0.07 0.06

Side

● ●

Cac

● 0.04

● ●

0.03 ●

Main

Main

0.02

0.00 0.0

Side

●

AN US

0.01

●

CR IP T

0.05

● ●

0.2

0.4

0.6

0.8

1.0

β

M

Figure 7: Dependence of the critical capillary number Cac on the confinement β.

We just mentioned that, for both the confinement ratios investigated,

ED

a critical capillary number Cac exists such that for Ca > Cac the particle deviates to the side branch, otherwise it ends in the main outflow branch. This critical value is reported in Fig. 7 for several confinement ratios from

PT

β = 0.1 to β = 0.8. At low confinement ratios a relatively small capillary number is sufficient to drive the particle to the side branch. The value of

CE

Cac increases of almost an order of magnitude with growing confinement. A non-monotonic dependence of Cac on the confinement β is found, with a

AC

maximum achieved at β ≈ 0.7. This is probably due to strong interactions of bigger particles with the channel walls. To summarize, for a given confinement, all the particles with Ca < Cac

remain in the main branch. For a fixed fluid, geometry and flow rate, such a condition is equivalent to say that particles with elastic modulus G higher 27

ACCEPTED MANUSCRIPT

than a critical value Gc end up in the main outlet. In contrast, particles with Ca > Cac (i.e., with elastic modulus G < Gc ) deviate to the side branch.

CR IP T

This phenomenon can be exploited to separate softer from harder particles, or even to select a group of particles with an elastic modulus in a certain

range [G1 ,G2 ]. Indeed, in a first step, the flow rate could be selected such

that particles with G < G2 deviate to the side branch. In a second step, these particles could be reflowed with a lower flow rate (thus, lowering Gc ),

AN US

so making the particles with G < G1 < G2 move to the side branch. Hence, particles with elasticity in the range of interest could be collected from the main branch. Of course, the inverse procedure would be also possible, first ‘cutting’ the softer particles and then the harder ones. The optimal procedure depends on the elastic modulus distribution of the suspension to treat.

M

In conclusion, the simulation results show that a simple device geometry, i.e., a straight channel with an orthogonal side branch, can be used to separate

ED

elastic particles based on their deformability. In addition, the separation can be tuned through variation of the flow rate. The only relevant technical

PT

‘drawback’ is connected to the fine control of the outlet flow rates, that need to be as close as possible. Indeed, an unbalance of the flow rates may affect

CE

the particle trajectories driving all of them toward one outlet regardless of their deformability. Finally, we remark that the 2D nature of our simulations can only give

AC

qualitative results on the dynamics of elastic particles in channel bifurcations. However, we are confident that, although more complex shapes could be expected in 3D near the corner bifurcation, the main phenomenology remains the same.

28

ACCEPTED MANUSCRIPT

5. Conclusions The dynamics of a single two-dimensional elastic particle immersed in

CR IP T

a Newtonian fluid flowing in a channel with an orthogonal side branch is

investigated through direct numerical simulations. The model equations are

solved through the Arbitrary Lagrangian Eulerian (ALE) Finite Element Method (FEM). The pressure difference between the two outlets is selected

such that the flow rates through the exit branches are equal. Simulations are

AN US

carried out by releasing a circular particle with its center lying on the main

channel centerline and by varying the capillary number and the confinement ratio.

As expected, high capillary numbers and confinement ratios promote large deformations, especially when the particle approaches the corner. The results

M

show that elasticity of the particle can deviate it toward the side branch, at variance with rigid particles that end up in the main exit channel. Hence,

ED

for a fixed confinement ratio, we can identify a critical capillary number that divides particles moving to the main branch (low elasticity and/or flow rate)

PT

and particles going to the side outlet (high elasticity and/or flow rate). Such a critical value increases with the confinement ratio, achieving a maximum at

CE

a strong confinement, then decreasing at even stronger confinement ratios. Although the simulations are limited to 2D, the results presented in this

AC

work indicate that a device with a simple geometry can be used to separate elastic particles based on their deformability. This separation technique is highly flexible as the critical elastic modulus can be tuned by manipulating the total flow rate.

29

ACCEPTED MANUSCRIPT

References [1] N. Pamme, Continuous flow separations in microfluidic devices, Lab

CR IP T

Chip 7 (12) (2007) 1644–1659. doi:10.1039/b712784g.

[2] S. Zheng, J.-Q. Liu, Y.-C. Tai, Streamline-based microfluidic devices

for erythrocytes and leukocytes separation, J. Microelectromech. Syst. 17 (4) (2008) 1029–1038. doi:10.1109/jmems.2008.924274.

AN US

[3] A. R. Pries, K. Ley, M. Claassen, P. Gaehtgens, Red cell distribution

at microvascular bifurcations, Microvasc. Res. 38 (1) (1989) 81–101. doi:10.1016/0026-2862(89)90018-6.

[4] C. Pozrikidis, J. Davis, Blood flow through capillary networks, in: S. M.

M

Becker, A. V. Kuznetsov (Eds.), Transport in Biological Media, Academic Press, New York, 2013, Ch. 6, pp. 213–252.

ED

[5] D. M. Audet, W. L. Olbricht, The motion of model cells at capillary bifurcations, Microvasc. Res. 33 (3) (1987) 377–396. doi:10.1016/

PT

0026-2862(87)90029-x.

[6] B. W. Roberts, W. L. Olbricht, Flow-induced particulate separations,

CE

AIChE J. 49 (11) (2003) 2842–2849. doi:10.1002/aic.690491116.

AC

[7] B. W. Roberts, W. L. Olbricht, The distribution of freely suspended particles at microfluidic bifurcations, AIChE J. 52 (1) (2006) 199–206. doi:10.1002/aic.10613.

[8] C. Xi, N. C. Shapley, Flows of concentrated suspensions through an

30

ACCEPTED MANUSCRIPT

asymmetric bifurcation, J. Rheol. 52 (2) (2008) 625–647. doi:10.1122/ 1.2833469.

CR IP T

[9] M.-G. Koh, C. Kim, Separation of particles of different sizes from non-

Newtonian suspension by using branched capillaries, J. Chem. Eng. Jpn. 40 (11) (2007) 964–972. doi:10.1252/jcej.07we006.

[10] G. D’Avino, M. A. Hulsen, P. L. Maffettone, Separation of particles in

non-Newtonian fluids flowing in T-shaped microchannels, Adv. Model.

AN US

Simul. Eng. Sci. 2 (9) (2015) 1–23. doi:10.1186/s40323-015-0033-9. [11] M. Manga, Dynamics of drops in branched tubes, J. Fluid Mech. 315 (1996) 105–117. doi:10.1017/s0022112096002352.

[12] J. O. Barber, J. P. Alberding, J. M. Restrepo, T. W. Secomb, Simulated

M

two-dimensional red blood cell motion, deformation, and partitioning in

ED

microvessel bifurcations, Ann. Biomed. Eng. 36 (10) (2008) 1690–1698. doi:10.1007/s10439-008-9546-4.

PT

[13] H. C. Woolfenden, M. G. Blyth, Motion of a two-dimensional elastic capsule in a branching channel flow, J. Fluid Mech. 669 (2011) 3–31.

CE

doi:10.1017/s0022112010004829. [14] M. Diez-Silva, M. Dao, J. Han, C.-T. Lim, S. Suresh, Shape and biome-

AC

chanical characteristics of human red blood cells in health and disease, MRS Bull. 35 (5) (2010) 382–388. doi:10.1557/mrs2010.571.

[15] M. J. Rosenbluth, W. A. Lam, D. A. Fletcher, Force microscopy of nonadherent cells: A comparison of leukemia cell deformability, Biophys. J. 90 (8) (2006) 2994–3003. doi:10.1529/biophysj.105.067496. 31

ACCEPTED MANUSCRIPT

[16] M. Yu, S. Stott, M. Toner, S. Maheswaran, D. A. Haber, Circulating tumor cells: approaches to isolation and characterization, J. Cell Biol

CR IP T

192 (3) (2011) 373–382. doi:10.1083/jcb.201010021. [17] T. W. Secomb, B. Styp-Rekowska, A. R. Pries, Two-dimensional simulation of red blood cell deformation and lateral migration in microvessels, Ann. Biomed. Eng. 35 (5) (2007) 755–765. doi:10.1007/

AN US

s10439-007-9275-0.

[18] I. Galaev, B. Mattiasson (Eds.), Smart polymers: applications in biotechnology and biomedicine, CRC Press, Boca Raton, 2007. [19] C. K. W. Tam, W. A. Hyman, Transverse motion of an elastic sphere

s0022112073001497.

M

in a shear field, J. Fluid Mech. 59 (1) (1973) 177–185. doi:10.1017/

ED

[20] T. M. Squires, S. R. Quake, Microfluidics: Fluid physics at the nanoliter scale, Rev. Mod. Phys. 77 (3) (2005) 977–1026. doi:10.1103/

PT

revmodphys.77.977.

[21] T. Gao, H. H. Hu, Deformation of elastic particles in viscous shear flow,

CE

J. Comput. Phys. 228 (6) (2009) 2132–2151. doi:10.1016/j.jcp.2008. 11.029.

AC

[22] M. M. Villone, M. A. Hulsen, P. D. Anderson, P. L. Maffettone, Simulations of deformable systems in fluids under shear flow using an arbitrary Lagrangian Eulerian technique, Comput. Fluids 90 (2014) 88–100. doi:10.1016/j.compfluid.2013.11.016.

32

ACCEPTED MANUSCRIPT

[23] R. B. Bird, E. S. Warren, N. L. Edwin, Transport phenomena 2nd Ed., John Wiley & Sons, New York, 2002.

CR IP T

[24] M. A. Hulsen, R. Fattal, R. Kupferman, Flow of viscoelastic fluids past a cylinder at high Weissenberg number: Stabilized simulations using

matrix logarithms, J. Non-Newtonian Fluid Mech. 127 (1) (2005) 27–39. doi:10.1016/j.jnnfm.2005.01.002.

AN US

[25] A. N. Brooks, T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis

on the incompressible Navier-Stokes equations, Comput. Methods in Appl. Mech. Eng. 32 (1-3) (1982) 199–259. doi:10.1016/0045-7825(82) 90071-8.

M

[26] H. H. Hu, N. A. Patankar, M. Y. Zhu, Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian-Eulerian technique, J.

ED

Comput. Phys. 169 (2) (2001) 427–462. doi:10.1006/jcph.2000.6592.

PT

[27] G. D’Avino, M. A. Hulsen, Decoupled second-order transient schemes for the flow of viscoelastic fluids without a viscous solvent contribution, J.

CE

Non-Newtonian Fluid Mech. 165 (23) (2010) 1602–1612. doi:10.1016/ j.jnnfm.2010.08.007.

AC

[28] A. C. B. Bogaerds, M. A. Hulsen, G. W. M. Peters, F. P. T. Baaijens, Time dependent finite element analysis of the linear stability of viscoelastic flows with interfaces, J. Non-Newtonian Fluid Mech. 116 (1) (2003)

33–54. doi:10.1016/s0377-0257(03)00099-5.

33

ACCEPTED MANUSCRIPT

[29] M. G. H. M. Baltussen, M. A. Hulsen, G. W. M. Peters, Numerical simulation of the fountain flow instability in injection molding, J. Non-

jnnfm.2010.03.001.

CR IP T

Newtonian Fluid Mech. 165 (11-12) (2010) 631–640. doi:10.1016/j.

[30] R. Gu´enette, M. Fortin, A new mixed finite element method for computing viscoelastic flows, J. Non-Newtonian Fluid Mech. 60 (1) (1995)

AN US

27–52. doi:10.1016/0377-0257(95)01372-3.

[31] N. O. Jaensson, M. A. Hulsen, P. D. Anderson, Stokes-Cahn-Hilliard formulations and simulations of two-phase flows with suspended rigid particles, Comput. Fluids 111 (2015) 1–17. doi:10.1016/j.compfluid.

AC

CE

PT

ED

M

2014.12.023.

34