Migration of a solid particle in the vicinity of a plane fluid–fluid interface

European Journal of Mechanics B/Fluids 30 (2011) 76–88

Contents lists available at ScienceDirect

European Journal of Mechanics B/Fluids journal homepage: www.elsevier.com/locate/ejmflu

Migration of a solid particle in the vicinity of a plane fluid–fluid interface A. Sellier a,∗ , L. Pasol b a

LadHyX, Ecole Polytechnique, 91128 Palaiseau Cédex, France

b

LPT, Espci, 10 rue Vauquelin, 75231 Paris Cédex, France

article

abstract

info

Article history: Received 26 April 2007 Received in revised form 22 September 2010 Accepted 23 September 2010 Available online 8 October 2010 Keywords: Solid particle Stokes flow Fluid–fluid interface Boundary-integral method

The slow migration of a small and solid particle in the vicinity of a gas–liquid, fluid–fluid or solid–fluid plane boundary when subject to a gravity or an external flow field is addressed. By contrast with previous works, the advocated approach holds for arbitrarily shaped particles and arbitrary external Stokes flow fields complying with the conditions on the boundary. It appeals to a few theoretically established and numerically solved boundary-integral equations on the particle’s surface. This integral formulation of the problem allows us to provide asymptotic approximations for a distant boundary and also, implementing a boundary element technique, accurate numerical results for arbitrary locations of the boundary. The results obtained for spheroids, both settling or immersed in external pure shear and straining flows, reveal that the rigid-body motion experienced by a particle deeply depends upon its shape and also upon the boundary location and properties. © 2010 Elsevier Masson SAS. All rights reserved.

1. Introduction

velocity U and the angular velocity 

For many applications such as particle separation it is still of importance to calculate the net force and torque exerted on a stationary small and solid particle or its freely suspended rigidbody when it is immersed in an unbounded Newtonian liquid and subject to a gravity field g and/or an external velocity field u∞ (with respect to a prescribed framework). For small enough particles inertial effects are negligible and such a task is then achieved within the relevant creeping flow approximation. In that direction a considerable literature exists (see, for instance, [1,2]). As soon as the particle is a sphere with radius a and center O′ simple results are available. More precisely, if the liquid has uniform viscosity µ and density ρ a sphere immersed in the external Stokes velocity field u∞ and moving with prescribed translational velocity U (the velocity of O′ ) and angular velocity  is known to experience the following net force F and torque 0 with respect to its center O′

U = u∞ (O′ ) +

[

F = 6πµa u∞ (O′ ) +

a2 6

] ∇ 2 u∞ (O′ ) − U ,

(1)

0 = 4πµa3 [∇ ∧ u∞ (O′ ) − 2]. Accordingly, if homogeneous with uniform density ρs a sphere freely suspended in the uniform gravity field g and the external Stokes flow velocity u∞ experiences the following translational

∗

Corresponding author. Fax: +33 1 69 33 30 30. E-mail address: [email protected] (A. Sellier).

0997-7546/$ – see front matter © 2010 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.euromechflu.2010.09.006

=

1 2

a2 6

2

∇ 2 u∞ (O′ ) + (ρs − ρ)a2 g, 9

(2)

∇ ∧ u∞ (O ). ′

For non-spherical particles there is no analytical equivalent to (1)–(2) but using a boundary-integral approach (see [3]) is in such circumstances quite efficient to numerically compute with good accuracy the particle’s rigid-body motion (U, ). In practice, boundaries might be encountered and likely to deeply affect the motion (U, ) predicted for the unbounded liquid. The hydrodynamic interactions occurring near a boundary furthermore depend upon its shape and nature: gas–liquid interface, fluid–fluid interface or solid–liquid boundary. Those three types of boundaries have been previously addressed for a plane boundary with the following different approaches and ranges of applications: (i) Analytical studies for a spherical particle. These works resort to the bipolar coordinates and provide the net force and torque acting on a sphere either translating or rotating at a prescribed velocity [4–9] or immersed in an external Couette [10], shear [11,12] or polynomial [13–15] flow, for a solid–liquid wall. The case of a fluid–fluid plane interface has been treated using the same method by Lee and Leal [16] including special cases such as the solid–fluid and gas–liquid boundaries and dealing only with the translation or rotation of the sphere (no ambient flow). Amazingly, the same analysis for this latter case has been recently repeated in Nguyen and Evans [17].

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

(ii) Use of singularities. The fundamental singularity that complies with the boundary conditions at a fluid–fluid interface is derived in Lee et al. [18] by extending a method proposed by Lorentz [13]. Using such singularities, which are the counterparts of the usual Stokeslets, [18,19] then provide high-order asymptotic approximations of the net force and torque applied on a sphere when it translates, rotates or is kept fixed in a pure shear or straining external flow field far enough from a gas–liquid, fluid–fluid or solid–fluid interface. Using the same approach, the case of a slender and axisymmetric body located far from a fluid–fluid interface and either translating [20] and/or rotating [21] or being embedded in ambient shear and straining flows [19] has been asymptotically solved. This work provides a new approach to systematically and accurately compute the rigid-body motion (U, ) of an arbitrarily shaped particle located close to a general plane boundary (i.e. a gas–liquid, fluid–fluid or solid–liquid boundary) in the presence of a gravity field and/or an external and arbitrary Stokes velocity field. It is organized as follows. Not only the governing problem but also six auxiliary Stokes flows induced by translations or rotations of the body and a resulting linear system governing the unknown rigid-body motion (U, ) are introduced in Section 2. Boundary-integral equations which dictate on the particle’s surface the tractions associated with these auxiliary Stokes flows are established in Section 3 together with velocity and pressure integral representations in the liquids. This boundary formulation is asymptotically solved for a distant boundary in Section 4. A numerical treatment, quite required for arbitrary particle–boundary gaps, is achieved using a boundary element technique and tested against [16] and the asymptotic analysis in Section 5. Numerical results for spheroids immersed in a uniform gravity field or pure shear or straining external velocity fields are given for different values of the liquids viscosity ratio and discussed. Finally, a few conclusions close the paper in Section 6. 2. Governing equations, relevant auxiliary flows and linear system This section presents the relevant governing low-Reynoldsnumber problem and shows how it is fruitful in determining the particle’s rigid-body motion to appeal to six auxiliary and specific Stokes flows in the liquids. 2.1. Adopted framework and governing problem

∇.σ(Vl , Pl ) = µl ∇ 2 Vl − ∇ Pl = 0 and ∇.Vl = 0 for (−1) x3 < 0 (l = 1, 2), l

V1 = V2 ,

V2 .e3 = 0,

for l = 1, 2 on Σ .

Fig. 1. A solid particle P near the x3 = 0 plane fluid–fluid interface Σ between two immiscible and viscous Newtonian liquids 1 and 2. The particle is immersed in the liquid 2 and subject to a uniform gravity g and/or external and steady Stokes flows (Vl , Pl ) in each liquid = 1, 2.

occupies the domain Dl (see Fig. 1) and the particle migrates with respect to the steady interface Σ at unknown translational velocity U (the velocity of O′ ) and angular velocity . When migrating, the particle also disturbs each flow (Vl , Pl ) so that the velocity and pressure fields at a point M in the fluid l become Vl + ul and Pl + ρl g.OM + pl , respectively. Introducing the separation distance h = O′ O.e3 > 0 between the interface and the particle’s center of mass O′ and the velocity scale V = Max(|U|, ||h) we further assume that µ2 > 0 and Re = ρ2 hV /µ2 ≪ 1. Accordingly (see [18]), the disturbances (ul , pl ), with stress tensors σ l , are quasi-steady low-Reynolds-number flows. Such flows are governed by the equations and boundary conditions

µl ∇ 2 ul = ∇ pl , ∇.ul = 0 and (ul , pl ) → (0, 0) as |x| → ∞ in Dl (l = 1, 2), u1 = u2

and

u2 .e3 = e1 .(σ 1 − σ 2 ).e3 = e2 .(σ 1 − σ 2 ).e3 = 0

(3) (4)

We impose a uniform gravity field g and wholly immerse in fluid 2 a solid particle P with smooth surface S, mass M , volume V , center of mass O′ and center of volume O′′ . The fluid l now

(5) (6) (7)

The particle has a negligible inertia and thus experiences zero net force and torque (for instance with respect to the center of mass O′ ). Exploiting (3) and designating by n the unit outward normal on S, those requirements yield the additional equalities

σ 2 .ndS = (ρ2 V − M)g,

∫S

(8) O′ M ∧ σ 2 .ndS = −M O′ O′′ ∧ g.

S

Note that we actually consider a fluid–fluid interface Σ with 0 < µl < ∞ for l = 1, 2 when requiring (4) and (6). Still for 0 < µ2 < ∞ it is however possible to also deal with either gas–liquid or solid plane surfaces Σ by adding the following cases: (i) µ1 = 0. In such circumstances Σ is a gas–liquid interface. We then discard the flows (V1 , P1 ) and (u1 , p1 ) and replace (4) and (6) with the conditions V2 .e3 = u2 .e3 = 0

and el .σ(V2 , P2 ).e3 = el .σ 2 .e3 = 0

for l = 1, 2 on Σ .

el .[σ(V1 , P1 ) − σ(V2 , P2 )].e3 = 0

on Σ ,

u2 = −V2 + U +  ∧ O′ M on S .

∫

We consider, as illustrated in Fig. 1, two immiscible Newtonian fluids 1 and 2 occupying the semi-infinite domains x3 > 0 and x3 < 0, respectively. As in Lee et al. [18], the fluid–fluid interface Σ is henceforth assumed to remain flat and motionless although each fluid l (l = 1, 2), with uniform viscosity µl and density ρl , is subject to a uniform gravity field g and/or an ambient steady Stokes flow with pressure Pl and velocity Vl with respect to a Cartesian framework (O, x1 , x2 , x3 ) attached to Σ . Each external flow (Vl , Pl ), with stress tensor σ(Vl , Pl ), is steady and free from body force. Due to the well-known conditions on a fluid–fluid interface (see [22]) these Stokes flows (V1 , P1 ) and (V2 , P2 ) satisfy the equations

77

(9)

(ii) µ1 = ∞. For this choice Σ is a solid plane wall. One then once more ignores the flows (U1 , P1 ) and (u1 , p1 ) replacing this time (4) and (6) with the relations V2 = u2 = 0 on Σ .

(10)

Henceforth we thus assume that µ1 ≥ 0 and replace (4) and (6) with (9) or (10) if µ1 = 0 or µ1 = ∞, respectively.

78

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

Within the adopted framework, one obtains the particle’s rigidbody motion (U, ) and the resulting flow disturbances (ul , pl ) for given triplet (ρ1 , ρ2 , µ2 ), viscosity ratio λ = µ1 /µ2 and gravity and/or ambient admissible (recall (3)–(4)) flows (Vl , Pl ) by solving (5)–(8). For a spherical particle with center O′ and radius a this can be analytically achieved, whatever the separation a/h, for the pure gravity-driven migration (i.e. in the absence of ambient flows (Vl , l )) owing to the solution established in Lee and Leal [16] by using bipolar coordinates. The case of a sphere located far from the interface in ambient straining or shear flows is asymptotically solved versus the small parameter a/h in Yang and Leal [19]. To the authors’ very best knowledge, the only available analysis for non-spherical particle near a flat fluid–fluid interface is given in [19–21] which provides an asymptotic solution for a straight and slender solid body translating or rotating in absence of gravity and ambient flow fields. For non-spherical particles a numerical treatment of (5)–(8) is possible but the computational works coping with one arbitrarily shaped particle by employing a collocation method [23,24] or a boundary element technique [25,26] are unfortunately restricted to the case of a solid plane wall Σ (i.e. to λ = µ1 /µ2 = ∞). A new approach is thus needed for a gas–liquid (λ = 0) or fluid–fluid (0 < λ < ∞) interface Σ . A very first idea might be to develop a Finite Element Code accurately solving Eqs. (5) in large enough but truncated fluid domains together with the mixed boundary conditions (6)–(7) for U2 prescribed on the surface S and a trial rigid-body motion X = (U, ). From the computed flow (u2 , p2 ) it is further possible to approximate (with a lack of precision inherent to the evaluation of the Cartesian derivatives of the velocity u2 ) the associated surface traction σ2 .n on S and therefore the integrals arising in (8). Those steps are likely to be iteratively worked out by changing X until conditions (8) hold. However, such a procedure not only requires a good guess value for X but would also be tremendously cpu time consuming and, above all, not accurate enough (when computing the stress tensor on S). We thus advocate and implement in this work a boundary method which is free from these disadvantages.

Accordingly, adding (12) to (13) immediately shows that

∫

(i)

u2 .σ L,2 .ndS =

S

For i = 1, 2, 3 and L = T or L = R let us introduce auxiliary (i) (i) (i) and coupled steady Stokes flows (uL,1 , pL,1 ) with stress tensor σ L,1 (i)

(i)

(i)

in fluid 1 and (uL,2 , pL,2 ) with stress tensor σ L,2 in liquid 2 which are induced by a translation or a rotation, respectively, of the solid particle along the ei direction. In other words, such flows obey (5)–(6) and the specific boundary conditions (i)

(i)

and uR,2 = ei ∧ O′ M on S .

uT ,2 = ei

(11)

(u(Li,)1 , p(Li,)1 )

satisfy (5) in D1 . The usual Both flows (u1 , p1 ) and reciprocal identity (see [2]) in fluid 1 therefore yields the relation

∫  Σ

(i)

(i)



u1 .σ L,1 − uL,1 .σ 1 .e3 dS = 0.

(12)

(i)

S

uL,2 .σ 2 .ndS .

(15)

Note that (15), obtained by applying in each liquid domain the usual reciprocal identity and exploiting the boundary conditions (6), is new and of basic importance when enforcing conditions (8). Indeed, by virtue of (11) and (15) it is straightforward to obtain

∫

ei .σ 2 .ndS =

∫S

∫ S

(i)

u2 .σ T ,2 .ndS ,

ei .[O M ∧ σ 2 .n]dS = ′

∫

S

S

(16)

(i)

u2 .σ R,2 .ndS .

Recalling boundary condition (7), (16) makes it possible to cast equalities (8) into the form of a 6-equation linear system for the unknown Cartesian velocity components Uj = U.ej and Ωj = .ej . (i)

= σ (Li,)2 .n the surface traction exerted on the (i) (i) particle’s boundary S by the specific flows (uL,l , pL,l ) and defining Designating by fL

the following coefficients i,j AL

=−

i,j

BL = −

1

µ2 1

µ2

∫ ∫S

(i)

ej .fL dS

and (17)

(ej ∧ O′ M).f(Li) dS ,

S

one indeed gets, adopting the usual tensor summation convention, the 6-equation linear system (i = 1, 2, 3) i,j

i ,j

AT Uj + BT Ωj =

  ∫ (M − ρ2 V )g.ei − V2 .f(Ti) dS µ2 ,

(18)

S i,j

2.2. Auxiliary coupled Stokes flows and resulting well-posed linear system

∫



i ,j

AR Uj + BR Ωj = − ρ2 V [O′ O′′ ∧ g].ei +

∫

(i)

V2 .fR dS

 µ2 . (19)

S

Before going on we check whether (18)–(19) is a well-posed system. For six arbitrary real numbers Ti , Ri (i = 1, 2, 3) let us (i) (i) consider in the fluid domain D2 the Stokes flow u′2 = Ti uT ,2 +Ri uR,2 (i)

with pressure p′2 = Ti pT ,2 + Ri pR,2 and stress tensor σ ′2 = Ti σ T ,2 + (i) Ri σ R,2 . Under our notations the rate E2

of dissipation of mechanical energy within fluid 2 for this flow satisfies (see [2])

∫

u2 .σ 2 .ndS − ′

E2 =

′

S

∫ Σ

u′2 .σ ′2 .e3 dS < 0.

(20)

Because (u′2 , p′2 ) obeys (6) note that u′2 .σ ′2 .e3 vanishes on the interface Σ . Recalling (11) and definition (17), (20) then becomes i,j

i,j

i,j

i ,j

Exploiting the reciprocal identity in D2 for the Stokes flows (u2 , p2 ) (i) (i) and (uL,2 , pL,2 ) also provides the additional link

E2 = −µ2 {Ti AT Tj + Ri AR Tj + Ti BT Rj + Ri BR Rj } < 0.

∫ 

(j) (j) i,j in (15) with (uL,2 , pL,2 ) and this change yields the properties AT j,i i,j j ,i i ,j j ,i

S

(i) u2 .σ L,2

−

−

∫  Σ

(i) uL,2 .σ 2 (i)



.ndS (i)



u2 .σ L,2 − uL,2 .σ 2 .e3 dS = 0.

(13) (i)

(i)

Because, as do both (u1 , p1 ) and (u2 , p2 ), the flows (uL,1 , pL,1 ) and

(u(Li,)2 , p(Li,)2 ) obey conditions (6) on the surface Σ , it is clear that     (i) (i) (i) (i) u1 .σ L,1 − u2 .σ L,2 .e3 = u1 . σ L,1 − σ L,2 .e3 = 0   = uL(i,)1 .σ 1 − uL(i,)2 .σ 2 .e3 on Σ . (14)

(j)

(21)

(j)

In addition, since (uL,2 , pL,2 ) fulfills (5)–(6) one can replace (u2 , σ 2 )

=

AT , BT = AR and BR = BR . Those relations and (21) prove that (18)–(19) has a 6 × 6 real-valued, symmetric and positivedefinite matrix and therefore a unique solution (U, ) whatever M , ρ2 , V , g and the ambient velocity V2 on S . The basic message delivered by (17) and (18)–(19) is that one can determine (U, ) by solely evaluating on the particle’s surface (i) (i) a very few quantities: the surface tractions fL = σ L,2 .n for i = 1, 2, 3 and L = T , R. We show how to compute these key vectors and the flow disturbances in liquids 1 and 2 in the next section.

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

∫ 

3. Advocated boundary formulation

S

This section recalls in Section 3.1 the relevant Green’s tensor for the addressed problem (previously employed by Lee et al. [18]) and also derives in Section 3.2 single-layer integral representations of the fluid velocity and pressure in each liquid domain. Those basic representations yield in Section 3.3 boundary-integral equations which permit us to obtain the particle’s rigid-body motion and, if needed, the disturbed velocity and pressure fields in each liquid.

=

(j)

(j)

79



u2 .σ 2,1 − u2,1 .σ 2 .ndS (x)

∫ 

(j)

(j)



u2 .σ 2,1 − u2,1 .σ 2 .e3 dS (x).

Σ

(28)

Since y lies in the liquid domain D1 the reciprocal identity for the (j) (j) (j) flows (u1 , p1 ) and (u1,1 , p1,1 ) with stress tensors σ 1 and σ 1,1 only holds in D1 \ {y} and a standard limit procedure this time yields

∫ 

(j)

(j)



3.1. The Green tensor

[u1 .ej ](y) = −

For two given integers m in {1, 2} and j in {1, 2, 3} we introduce (j) (j) at the field point x(x1 , x2 , x3 ) the coupled Stokes flows (u1,m , p1,m )

Combining (28) with (29) and exploiting (see (14)) the conditions (6) and (23) (for m = 1), one ends up with the boundary-integral representation

(j)

(j)

if x3 > 0 and (u2,m , p2,m ) if x3 < 0 induced by a concentrated force point of strength ej located at the pole y(y1 , y2 , y3 ) such that (−1)m y3 < 0 and subject to all the boundary conditions at the (j) interface Σ . These two flows, with stress tensors σ l,m , thus obey

µl ∇ 2 u(l,jm) 

(j)

(j)

= 

ul,m , pl,m

(j) u1,m

=

(j) u2,m .e3

∇ p(l,jm)

− δ3d (x − y)ej ,

∇.ul(,jm)

= 0, (22)

→ (0, 0) at ∞ in Dl ,

[u1 .ej ](y) = −

∫

with δ3d the three-dimensional delta function. The related Green (j) velocity tensor Gkj (x, y)ek ⊗ ej and pressure vector pj (x, y) = pl,m (j)

Gkj (x, y) = Gjk (y, x).

8πµ2 Gkj (x, y) =

R1kj

(x, y),

8π pj (x, y) =

Pj1

(x, y)

if x3 > 0 and y3 < 0,

(25)

and, designating by δkj the kronecker delta, one arrives at 8πµ2 Gkj (x, y) =

δkj [(x − y).ek ][(x − y).ej ] + |x − y| |x − y|3 2 + Rkj (x, y), x3 < 0, y3 < 0,

S

u2 .σ 2,1 .ndS (x)

∫ 

∫S 

+

Σ

0=

∫  Σ

(j)

(j)



u2,2 .σ 2 − u2 .σ 2,2 .ndS (x) (j)

(j)



u2,2 .σ 2 − u2 .σ 2,2 .e3 dS (x),

(j)

(j)



u1,2 .σ 1 − u1 .σ 1,2 .e3 dS (x)

(31) (32)

and it follows that (30) is replaced with the representation

[u2 .ej ](y) = −

∫

(j)

S

∫

u2,2 .σ 2 .ndS (x) (j)

+ S

u2 .σ 2,2 .ndS (x) for y in D2 .

(33)

As shown in Appendix B, boundary condition (7) together with the definition of Gkj (x, y) and property (24) permit us to rewrite (30) and (33) as

[ul .ek ](x) = −

∫

[f(y).ej ]Gkj (x, y)dS (y) S

for x in Dl and l = 1, 2 (26)

[(x − y).ej ] + Pj2 (x, y) if x3 < 0 and y3 < 0 (27) |x − y|3

with δkj the kronecker delta and functions R1kj , R2kj , Pj1 and Pj2 displayed in Appendix A. Those functions dependent upon the viscosity ratio λ = µ1 /µ2 and turn out to be regular as the observation point x tends to the pole y. When the pole y lies far from Σ (i.e. when y3 → −∞) these functions moreover vanish so that (26)–(27) recover the usual free-space Green velocity tensor and pressure vector Cartesian components (see [1]). 3.2. Velocity and pressure integral representations We now obtain integral representations of the disturbances

(u1 , p1 ) and (u2 , p2 ) governed by (3)–(5) by extending the material presented in Pozrikidis [3] for a solid boundary Σ . Let us first select in (22)–(23) a pole y located in the liquid domain D1 , i.e. such that y3 > 0 (case m = 1). In the liquid domain D2 both (u2 , p2 ) (j) (j) (j) and (u2,1 , p2,1 ) are Stokes flows with stress tensors σ 2 and σ 2,1 , respectively. The reciprocal identity then implies that

(j)

+

Now if y lies in the liquid domain D2 the reader may easily check that the counterparts of (28)–(29) read

(24)

When the pole y lies in liquid 2 (i.e. when y3 < 0) one actually arrives at

(29)

(30)

(j)

such that Gkj (x, y) = ul,m .ek and pj (x, y) = pl,m if (−1)l x3 < 0 and (−1)m y3 < 0 have been analytically obtained in Lee et al. [18]. Using those results shows that

∫

for y in D1 .

    = e1 . σ (1j,)m − σ 2(j,)m .e3 = e2 . σ (1j,)m − σ 2(j,)m .e3 = 0 (23)

on Σ ,

(j) u2,1 .σ 2 .ndS (x)

S

[u2 .ej ](y) = −

(j) u2,m ,

8π pj (x, y) = 2

Σ

u1,1 .σ 1 − u1 .σ 1,1 .e3 dS (x).

(34)

where, recalling that in D2 the ambient Stokes flow (V2 , P2 ) has stress tensor σ(V2 , P2 ), we introduced the surface traction f on S as f = [σ(V2 , P2 ) + σ 2 ].n. Hence, the flow disturbances u1 and u2 are produced by a source distribution on the particle’s surface whose density is −f = −[σ(V2 , P2 ) + σ 2 ].n. Accordingly, the pressure disturbances p1 and p2 read pl (x) = −

∫

[f(y).ej ]pj (x, y)dS (y) S

for x in Dl and l = 1, 2.

(35)

As revealed by (34), the Cartesian components of the velocity ul admit a single-layer integral representation in the entire liquid domain Dl . According to (34)–(35) the evaluation of the flows (u1 , p1 ) and (u2 , p2 ) solely require to compute the previously defined surface traction f on the particle’s surface S . How such a key task is achieved is explained in the next subsection. 3.3. Relevant boundary-integral equations and proposed procedure It is easy to check that (34) (but not (35)) also holds for y located on the surface S. By virtue of (7), one thus immediately arrives at

80

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

the Fredholm boundary-integral equation of the first kind

where G0kj is the usual Oseen–Green tensor prevailing in an

[U +  ∧ (x − OO′ ) − V2 ].ek ∫ = − [f(y).ej ]Gkj (x, y)dS (y) for x on S

unbounded liquid 2 (here obtained by setting the term R2kj to zero (i),n

in (26)) and the velocities uL (36)

S

which provides f once the rigid-body motion (U, ) is known. The (i) (i) basic point is that our auxiliary Stokes flows (uL,l , pL,l ) also satisfy (5)–(6) with a boundary condition (11) analogous to (7) (set V2 = 0 and adequately select U and ). Thus (34)–(36) hold for those (i) (i) flows and the needed surface traction fL = σ L,2 .n on S obeys the boundary-integral equation

[u(Li,)2 .ek ](x) = −

∫ 

(i)



fL .ej (y)Gkj (x, y)dS (y)

for x on S .

(37)

S

In view of the previous material, the advocated strategy to solve the governing problem (5)–(8) thus consists in the following steps: (i)

(i) One first obtains the required surface tractions fL on S by inverting six boundary-integral equations (37) with conditions (11). The solution is obtained up to a constant multiple of n (as obtained in Pozrikidis [3] for the unbounded flow). The unknown rigid-body motion (U, ) is subsequently determined by solving the well-posed linear system (18)–(19), under definition (17). (ii) Injecting in (36) the previously computed value of (U, ) one further gains the surface traction f occurring in (34)–(35) by inverting the boundary-integral equation (36). If needed, the disturbances (u1 , p1 ) and (u2 , p2 ) can be subsequently evaluated at any field point x in the liquids by using the integral representations (34) and (35).

(i),0

(i),0 = ei , uR = ei ∧ O′ M, [ ]   c δkj + dδk3 δj3 i,j,0 (i),1 uL .ek (x) = − AL , a [ ]   c δkj + dδk3 δj3 i,j,1 (i),2 uL .ek (x) = − AL a ∫ 1 ′ ′ ′ ′ e[δk3 δj3 yj − δj3 δk3 yk ] + 2 a S   + [δkj (f + v) + (g + v)δk3 δj3 ]y′3 f(Li),0 .ej (y)dS (y)  e   i,j,0 ′ ′ ′ ′ − 2 AiL,3,0 δk3 xk − δk3 AL δj3 xj a [ ] f + v + (g + v)δk3 i,k,0 − AL x′3 2

uT

(40) (41)

(42)

a

′ with δk3 = 1 − δk3 and, for λ = µ1 /µ2 , the following coefficients

c (λ) = e(λ) = f (λ) =

v(λ) =

4. Asymptotic analysis for a distant interface

satisfy

2 − 3λ 32π (1 + λ) 2 + 3λ 64π (1 + λ) 2 − 5λ

,

λ

32π (1 + λ)

3(λ + 2) 32π (1 + λ)

, (43)

, ,

64π (1 + λ)

d(λ) = −

g (λ) = −

6 + 5λ 64π (1 + λ)

, (44)

. (i),n

As soon as the particle, with typical length scale a and center of volume O′ , lies far from the interface one may asymptotically solve the boundary-integral equation (37) and the system (18)–(19) versus the small parameter ϵ = a/h ≪ 1 with h = O′ O.e3 . By contrast with the analysis carried out for a sphere in Lee et al. [18], the present section establishes asymptotic results for arbitrarily shaped particles and also provides analytical results for ellipsoids admitting three planes of symmetries normal to e1 , e2 and e3 .

According to (40)–(42), fL is the surface traction exerted on the particle if immersed in an unbounded liquid 2 and subject on its (i),n boundary to the velocity uL . Such a velocity is divergence free (for n = 2 this is due to the equality f + g + 2(e + v) = 0 obtained from (43)–(44)) and therefore the boundary-integral equation (39) (i),n admits a solution fL (up to a constant multiple of n). 4.2. First-order solution for arbitrarily shaped particles (i),1

Because uL We use Cartesian coordinates (O′ , x′1 , x′2 , x′3 ) centered at the particle’s center of volume O′ with x1 = x′1 , x2 = x′2 and x3 = x′3 − h and for y(y1 , y2 , y3 ) on the particle’s surface S similarly set y′1 = y1 , y′2 = y2 and y′3 = y3 + h. If both x and y lie on S one has |x′k | = O(a), |y′k | = O(a) and for a particle distant from Σ also |x′3 | ≪ h, |y′3 | ≪ h. Exploiting those properties, (26) and Appendix A makes it possible to expand (37) versus the small parameter ϵ = a/h therefore arriving at the asymptotic approximation (i)

(i),0

fL = fL

(i),1

+ ϵ fL

2 (i),2

+ ϵ fL

+ O(ϵ 3 ).

(38)

Curtailing the details and introducing for n = 0, 1, 2 the i,j,n i,j,n (i) coefficients AL and BL by replacing in (17) surface traction fL (i),n

with fL , one ends up with the pyramidal system of boundaryintegral equations

−

∫ 

(i),n

fL

   .ej (y)G0kj (x, y)dS (y) = u(Li),n .ek (x)

S

for x on S and n = 0, 1, 2

turns out to be uniform it is straightforward to

(i),1

4.1. A pyramidal system

(39)

(i),0

express fL in terms of the tractions fT whatever the geometry of the particle. More precisely, (41) yields (up to a constant multiple of n) the solution (i),1

fL

=−

1



a

i,1,0 (1),0 fT

c (λ)[AL

+ [c (λ) + d(λ)] (i),1

Hence, fL (i),0

fL

+ AiL,2,0 f(T2),0 ]

i,3,0 (3),0 AL fT



.

(45)

only requires the knowledge of the surface tractions

exerted on the particle when it translates or rotates in an (i)

unbounded liquid 2. Injecting in (17)–(19) the approximation fL = (i),0

(i),1

fL + ϵ fL yields (U, ) = (U ,  ) + ϵ(U ,  ) with (U , 0 ) the particle’s rigid-body motion in an unbounded liquid 2 and ϵ(U1 , 1 ) the first-order correction. By superposition, two cases arise: (i) Pure gravity-driven migration. In this case there is no external flow and (U1 , 1 ) ̸= (0, 0) so that the freely suspended particle experiences O(a/h) and thus long range interactions with the distant interface. (ii) Pure external flow. In the absence of gravity the very specific solution (45) when put into (17)–(19) easily shows that 0

0

1

1

0

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

(U1 , 1 ) = (0, 0). Thus, a particle freely suspended in an external flow experiences short range particle–interface interactions which scale as (a/h)2 far from the interface. Finally, note that the net force and net torque exerted on a particle moving under the action of a gravity and/or an external flow field are affected (see (17)) up to order ϵ = a/h by the distant interface. 4.3. Second-order solution for ellipsoids

(i),0

fT

= −µ2 CT(i) s(M )ei ,

(i),0

fR

(i)

(46)

(i)

s(M ) = {(x′1 /a21 )2 + (x′2 /a22 )2 + (x′3 /a23 )2 }−1/2 .

(47)

Denoting by Ve = 4π a1 a2 a3 /3 the volume of the ellipsoid, one thus arrives at i,j,0

= BiT,j,0 = 0,

i,j,0

= 3Ve CT(i) δij ,

i,j,0

= Ve CR(i) (a21 + a22 + a23 − a2i )δij .

BR

AR

(i),1

Exploiting (45) then provides fT (i),1

fR

i,j,1

AT

i,j,1

AR

(48)

= −[c (λ)+δi3 d(λ)]AiT,i,0 f(Ti),0 /a,

= 0 and =−

δij a

[c (λ) + δi3 d(λ)](AiT,j,0 )2 ,

(49)

= BTi,j,1 = BiR,j,1 = 0. (i),2

It is too heavy to get the second-order approximation fL

but

i ,j ,2 easy to calculate the related approximation AL

by resorting to the usual reciprocal identity for Stokes flows in an unbounded liquid 2. For instance, we get

∫

(i),2

∫

(j),0

.uT(i),2 dS , ∫ ∫ (ej ∧ O′ M).fR(i),2 dS = fR(j),0 .uR(i),2 dS . ej .fT

fT

dS =

S

S

(50)

S

S

From (42), (46) and (49) it thus easily follows that i,j,2

AT

δij

=

a2

[c (λ) + δi3 d(λ)]2 (AiT,j,0 )3 ,

i,j,2

BR

= 0.

(51)

As detailed in Appendix C, replacing the subscript R with T in (50) i,j,2 j,i,2 permits one to evaluate each quantity BT = AR . The only nonzero coefficients read 2,1,2 AR

1,2,2

AR

=

1,2,2 BT

(2)

= −CR

1,1,0 Ve AT

= B2T ,1,2 = CR(1) Ve A2T ,2,0





a21 e(λ) + a23 [f (λ) + v(λ)] a2

a22 e

(λ) +

a23

[f (λ) + v(λ)] a2

i,j





, (52)

. (53)

i ,j

In summary, each coefficient AL and BL , introduced in (17), is of order ϵ 3 except the following ones (without summing over indices i in (54)) i,i

i,i,0

AT = AT

{1 − Ai ϵ + (Ai ϵ)2 } + O(ϵ 3 ), i,i,0

Ai = [c (λ) + δi3 d(λ)]AT

/a,

2,1

= BT1,2 = ϵ 2 BT1,2,2 + O(ϵ 3 ),

1,2

= BT2,1 = ϵ 2 BT2,1,2 + O(ϵ 3 ).

AR AR

Ui = Ui0 {1 + Ai ϵ} + O(ϵ 3 ), i,i,0

Ui0 = (ρs − ρ2 )Ve g.ei /[µ2 AT U20 a

{r1 ϵ 2 + O(ϵ 3 )},

2,2,0

r1 = −

Ω2 =

AT

U10 a

(56)

],

(57)

{a22 e(λ) + a23 [f (λ) + v(λ)]} , a(a22 + a23 )

{r2 ϵ 2 + O(ϵ 3 )},

1,1,0

= −µ2 CR(i) s(M )ei ∧ O′ M

with coefficients CT and CR given in Appendix C and

AT

If homogeneous with uniform density ρs the ellipsoid then experiences far from the interface a pure gravity-driven rigid-body (U, ) obtained by solving (18)–(19) for O′ = O′′ and V2 = 0 with (54)–(55) and given by

Ω1 =

Result (42) prevents one from gaining in closed form the second-order approximation when the particle’s geometry is arbitrary. As detailed in this section, analytical results however exist for ellipsoidal particles. If the particle’s domain has inequality x′12 /a21 + x′22 /a22 + x′32 /a23 ≤ 1, with x′i = O′ M.ei , it is known (see [27]) that (without summing over indices i in (46))

81

(54)

(55)

{a21 e(λ) + a23 [f (λ) + v(λ)]} , a(a21 + a23 )  0 3 |U |ϵ . Ω3 = O

r2 =

AT

(58)

a

At this stage it is worth pointing out the following properties of the derived asymptotic expansions (54)–(58): i,j,0 (i) For a sphere with radius a one gets AT = 6π aδij and using (54)–(55) immediately recovers the results obtained in Lee et al. [18]. (ii) The proposed approximations only require to evaluate the (i) (i) coefficients CT and CR with the formulae displayed in Appendix C. Note that the case of either prolate or oblate spheroids then yields simple results also available in Appendix C. (i) (iii) For i = 1, 2 the coefficient Ai = 3c (λ)Ve CT /a vanishes for (i)

the critical viscosity ratio λc = 2/3. Since CT > 0 (see [1]), the ellipsoid subject to a gravity parallel to a distant interface settles faster or slower than in an unbounded fluid for λ < 2/3 and λ > 2/3, respectively. Note that the ellipsoid then rotates parallel to the interface with a weak angular velocity of order ϵ 2 . Of course, if the homogeneous ellipsoid does not admit a plane of symmetry parallel to the distant interface results (54)–(58) do not hold. Symmetry actually shows that we can in practice asymptotically track as time evolves the gravity-driven motion of the homogeneous ellipsoid considered in this subsection by applying (54)–(58) if g is normal to the interface or a1 = a3 and g ∧ e1 = 0 or a2 = a3 and g ∧ e2 = 0. 5. Numerical results for a close interface When the particle lies close to the interface we determine its rigid-body motion by implementing the numerical procedure advocated in Section 3.3. This section also benchmarks the adopted numerical strategy and discusses our numerical results. 5.1. Numerical strategy and comparisons Each integral equation (37) is numerically inverted by using a standard boundary element method (see, among others, Beskos [28] and Bonnet [29]). For the sake of accuracy we use on the particle’s surface an N-node mesh consisting of isoparametric and curvilinear triangular elements. The weak singularity produced by the first two terms on the right-hand side of (26) is adequately handled by using polar coordinates, as explained in Bonnet [29]. The discretized counterpart of (37), a 3N-equation linear system with non-symmetric and fully populated square matrix, is finally solved by Gaussian elimination. In order to ascertain the accuracy of the implemented procedure, comparisons have been made with the analytical results obtained in Lee and Leal [16] for a sphere

82

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

Table 1 Computed normalized coefficients t11 and t33 versus the number of collocation points on the sphere’s surface for five values of the viscosity ratio λ and two separation parameters h = 1.1a and h = 2a. For comparisons the analytical results provided in Lee and Leal [16] are also displayed. tii t11 t11 t11 t11 t33 t33 t33 t33 t11 t11 t11 t11 t33 t33 t33 t33

h/a

N 74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16]

λ=0

1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 2 2 2 2 2 2 2 2

0.7380 0.7412 0.7413 0.7413 4.3429 4.6054 4.6275 4.6255 0.8372 0.8387 0.8388 0.8389 1.5913 1.5960 1.5966 1.5967

λ= 0.1

λ=1

λ = 10

0.7792 0.7804 0.7806 0.7806 4.5386 4.8279 4.8525 4.8503 0.8673 0.8689 0.8691 0.8691 1.6275 1.6325 1.6332 1.6331

1.0496 1.0517 1.0519 1.0519 5.7942 6.3080 6.3539 6.3494 1.0378 1.0401 1.0404 1.0404 1.8147 1.8208 1.8216 1.8216

1.8118 1.8215 1.8206 1.8207 8.3871 9.7636 9.8966 9.8812 1.2996 1.3033 1.3037 1.3037 2.0538 2.0616 2.0625 2.0626

λ=∞ 2.2475 2.2674 2.2643 2.2643 9.3805 11.288 11.482 11.459 1.3781 1.3823 1.3827 1.3828 2.1162 2.1244 2.1254 2.1255

and with the asymptotic results established in Section 4.2 for spheroidal particles distant from the interface. For a sphere with radius a and center O′ symmetries show that 1,2 2,1 2,1 1,2 i,i i,i only AT , BR and BT = AR = −BT = −AR are non-zero with 2 ,2

1 ,1

2,2

1 ,1

Table 2 Computed normalized coefficients c11 , c12 and c33 versus the number of collocation points on the sphere’s surface for five values of the viscosity ratio λ and two separation parameters h = 1.1a and h = 2a. For comparisons the analytical results provided in Lee and Leal [16] are also displayed. cij

N

h/a

λ=0

λ = 0.1

λ=1

λ = 10

λ=∞

c11 c11 c11 c11 c12 c12 c12 c12 c33 c33 c33 c33 c11 c11 c11 c11 c12 c12 c12 c12 c33 c33 c33 c33

74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16] 74 242 1058 [16]

1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 2 2 2 2 2 2 2 2 2 2 2 2

1.0631 1.0637 1.0654 1.0655 0.1038 0.1044 0.1044 0.1044 0.9250 0.9208 0.9205 0.9205 0.9968 0.9960 0.9958 0.9958 0.0386 0.0389 0.0389 0.0389 0.9885 0.9848 0.9846 0.9847

1.0833 1.0840 1.0858 1.0859 0.1018 0.1025 0.1025 0.1025 0.9375 0.9334 0.9331 0.9331 1.0007 0.9998 0.9996 0.9996 0.0363 0.0365 0.0366 0.0366 0.9912 0.9875 0.9874 0.9874

1.1967 1.1986 1.2005 1.2008 0.0811 0.0814 0.0812 0.0812 1.0038 1.0001 1.0000 1.0000 1.0185 1.0178 1.0176 1.0176 0.0229 0.0230 0.0229 0.0230 1.0038 1.0001 0.999 1.0000

1.3801 1.3865 1.3887 1.3889 −0.0143 −0.0169 −0.0192 −0.0192 1.0941 1.0910 1.0915 1.0915 1.0380 1.0373 1.0371 1.0371 0.0016 0.0015 0.0015 0.0015 1.0167 1.0131 1.0130 1.0130

1.4433 1.4522 1.4546 1.4549 −0.0792 −0.0849 −0.0887 −0.0876 1.1192 1.1164 1.1170 1.1171 1.0426 1.0420 1.0418 1.0418 −0.0048 −0.0051 −0.0050 −0.0050 1.0196 1.0160 1.0159 1.0159

AT = AT and BR = BR . It is thus sufficient to compute the normalized coefficients 1 ,1

t11 =

AT

6π a

3 ,3

,

t33 =

3 ,3

c33 =

BR

8π a

, 3

AT

6π a

1,1

,

c11 =

1,2

c12 =

BR

8π a3

, (59)

BT

8 π a2

versus the distance h = O′ O.e3 and the viscosity ratio λ = µ1 /µ2 . This is numerically achieved for λ = 0, 0.1, 1, 10, ∞ when the sphere is either close to (h = 1.1a) or in the vicinity (h = 2a) of the interface by using three increasingly refined N-node meshes on its surface. As illustrated in Tables 1 and 2, the results nicely agree with Lee and Leal [16] and suggest that using 242 collocation points on the sphere’s boundary provides a sufficient relative accuracy of order 10−3 in the range h ≥ 1.1a. It is also worth testing the numerics against formulae (56)–(58). When the ellipsoid happens to be a spheroid with a1 = a2 and a3 = sa1 we select a = a3 and thus ϵ = a3 /O′ O.e3 ≫ 1. In this case the theoretical values of A1 = A2 , A3 and r2 = −r1 solely depend upon the ‘‘slenderness ratio’’ s and λ and read A1 = r2 =

2 − 3λ 2(1 + λ)(χ ′ + α1 )

,

A3 = −

2 + 3λ + s2 (2 − 3λ) 4(1 + s2 )(1 + λ)(χ ′ + α1 )

2 + 3λ

(1 + λ)(χ ′ + s2 α3 )

,

(60)

(61)

with χ ′ , α1 and α3 given in Appendix C versus s. Calculating (U10 , U30 ) for ϵ = 10−10 and (U1 , a3 Ω2 , U3 ) for ϵ = 10−2 provides numerical approximations of Ai and r2 as (Ui /Ui0 − 1)ϵ −1 and a3 Ω2 ϵ −2 /U10 , respectively. As shown in Table 3, these computed values perfectly match with (60)–(61).

Fig. 2. Three addressed spheroids (a1 = a2 ) near the interface Σ . (a) Oblate spheroid with a1 = 2a3 . (b) Sphere with radius as . (c) Prolate spheroid with a3 = 3a1 /2. The ratio a1 /as is given in cases (a) and (c) as explained in Section 5.2.1.

5.2.1. Gravity-driven motion If only the gravity g acts such a spheroid, with uniform density ρs , translates without rotating at the velocity Ui0 ei given by (56) when it lies far from the interface, i.e. for h/a3 → ∞. For symmetry reasons it appears that (U, ) = (U3 e3 , 0) if g ∧ e3 = 0 whereas U and  are parallel with g and g ∧ e3 , respectively if g.e3 = 0. Accordingly, it is sufficient to compute versus the separation ratio h/a3 the normalized velocity components u1 = u3 =

5.2. Numerical results and discussion We consider, as illustrated in Fig. 2, homogeneous and spheroidal solid particles in the vicinity of the interface Σ with a1 = a2 and h = O′ O.e3 = O(a3 ) so that the asymptotic analysis of Section 4 breaks down and numerical treatment is required.

U.e1 U10 U.e1 U30

,

w2 =

[

.e2 U30

] as

if g = ge1 ; (62)

if g = ge3 .

Note that in scaling .e2 in (62) we used U30 (not U10 ) and also introduced a length as . This latter is such that the sphere with radius as and density ρs admits far from the interface and for g = ge3 the same velocity as the spheroid, i.e. U30 = 2(ρ −

ρs )a2s g /[9µ2 ]. Using (56), (48) yields CT(3) = 3/[2a2s ]. This link

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

83

Table 3 Comparisons between the theoretical triplet (A1 , r2 , A3 ) and the computed one (A′1 , r2′ , A′3 ) for λ = 0.3, 2/3, 2 using 1058 collocation points on a sphere (s = 1) or an oblate spheroid (s = 1/2) and 1634 collocation points on a prolate spheroid (s = 3/2).

(s, λ)

A1

A′1

r2

r2′

A3

A′3

(0.5, 0.3) (0.5, 2/3) (0.5, 2.0) (1, 0.3) (1, 2/3) (1, 2.0) (1.5, 0.3) (1.5, 2/3) (1.5, 2.0)

0.25153 0 −0.39635 0.15865 0 −0.25000 0.12631 0 −0.19904

0.25158 7.16E−5 −0.39625 0.15866 1.25E−5 −0.24998 0.12630 −8.87E−6 −0.19904

0.58081 0.57075 0.55490 0.14423 0.11250 0.06250 0.06330 0.03675 −0.00510

0.58022 0.57013 0.55422 0.14380 0.11207 0.06207 0.06323 0.03668 −0.00517

−1.51465 −1.62954 −1.81061 −0.83654 −0.90000 −1.00000 −0.55083 −0.66100 −0.73444

−1.51438 −1.62923 −1.81022 −0.83651 −0.89995 −0.99995 −0.55083 −0.66100 −0.73443

Fig. 3. Velocity component u3 for λ = 0 (—), λ = 0.5 (◦), λ = 1 (∗), λ = 5 (•) and λ = ∞ (solid line). The smallest or largest value of u3 for a given pair (d, λ) is obtained for the oblate or prolate spheroid, respectively.

and (77), (81)–(83) provide a1 (and then a3 ) when (as , a3 /a1 ) is prescribed. Three particles, depicted in Fig. 2(a–c) and having the same value for U30 , are considered here: an oblate spheroid with a1 = 2a3 ∼ 1.3456as , the sphere with radius as and a prolate spheroid with 2a3 /3 = a1 ∼ 0.8570as . We plot in Fig. 3 for those spheroids and λ = 0, 0.5, 1, 5, ∞ the quantity u3 versus the separation parameter d = h/a3 in the range [1/0.9, 10]. Each spheroid is seen to behave as the sphere: it is slowed down by the interface for any prescribed viscosity ratio λ (u3 < 1 decreases as d drops) and the velocity decreases as λ increases. As predicted by (56), the velocity slowly asymptotes to unity as d becomes large. By contrast, the interface–particle interactions become significant in the range d ≤ 5. These interactions are found for a given value of d = O(1) to deeply depend upon the particle shape and the viscosity ratio. For example at d = 2 one obtains u3 ∼ 0.2 (i.e. very strong interactions) for the oblate spheroid and λ = ∞ (solid wall case) whereas u3 ∼ 0.7 (i.e. reasonable interactions) for the prolate spheroid and λ = 0 (free surface case). For a given pair (d, λ) the prolate spheroid or the sphere moves faster than the sphere or the oblate spheroid, respectively. This is because a1 and thus the fluid resistance to the particle’s translation decreases as the slenderness ratio s = a3 /a1 increases. Plotting in Fig. 4 the velocity components u1 and w2 shows that the prolate spheroid exhibits the same behaviour as the sphere if g = ge1 . Again (remind (56)) u1 slowly tends to its limit u1∞ far from the interface. If u1∞ = 1 for the sphere one obtains u1∞ = 0.9225 for the prolate spheroid with U10 ̸= U30 . The magnitude of interface–particle interactions monotonically decreases as the particle–interface gap increases but, by contrast with u3 , those interactions may either speed up (u1 > u1∞ ) or slow down (u1 < u1∞ ) the solid. For a prescribed location d there exists a critical

Fig. 4. Velocity component u1 (a) and w2 (b) for a sphere (letter s next to each plot) or a prolate spheroid if λ = 0 (—), λ = 0.5 (◦), λ = 1 (∗), λ = 5 (•) and λ = ∞ (solid line). The horizontal solid lines in (a) indicate the asymptotic value of u1 for d = h/a3 large.

viscosity ratio λc = λc (d) at which the particle ignores (as regards its translational velocity) the interface (u1 = u1∞ ) with u1 smaller or larger than u1∞ for λc > λc or λc < λc , respectively. The plots in Fig. 4(a) show that for λ large or small the interactions become strong close to the interface whereas u1 is close to u1∞ whatever d as soon as λ = O(1). This behaviour, well illustrated by the curves associated with λ = 0.5, 1, indicates that λc = O(1) (as seen in Section 4.2, λc ∼ 2/3 for d roughly exceeding 10). The angular velocity w2 in seen in Fig. 4(b) to exhibit the same trends with this time another critical viscosity ratio at which no

84

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

5.2.2. Case of pure shear and straining external flows We now immersed, in the absence of gravity, the previous spheroids in pure ambient shear flow V2sh and straining flow V2st such that KE ̸= 0 and V2sh = Kx3 e1 ,

V2st = E (x1 e1 + x2 e2 − 2x3 e3 ).

(63)

The above straining flow is only applied for a non-solid boundary Σ (due to unsatisfied no-slip conditions) and if O′ O is not parallel with e3 a simple shift permits one to reduce the problem to the selected straining flow because a particle immersed in a uniform external flow experiences the velocity of that flow (appeal to (18)–(19)). For 0 < λ < ∞ the associated external flows in liquid 1 are readily V1sh = Kx3 e1 /λ and V1st = E (x1 e1 + x2 e2 − 2x3 e3 ). If at rest or freely suspended the spheroid then experiences net force F and torque 0 (with respect to its center O′ ) or a rigid-body migration (U, ), respectively. For symmetry reasons it appears that (U, , F, 0) = (U1 e1 , Ω2 e2 , F1 e1 , Γ2 e2 ) for V2 = V2sh whilst (U, , F, 0) = (U3 e3 , 0, F3 e3 , 0) for V2 = V2st . Using the key reciprocal identity (15) one obtains

∫

(1)

V2sh .fT dS ,

F1 = −

Γ2 = −

(2)

V2sh .fR dS , S

S

∫

∫

(64)

(3)

V2st .fT dS .

F3 = − S

The limit values U10 , U30 , F10 , F30 and Γ20 for a spheroid with s = a3 /a1 located at h = O′ O.e3 in the absence of the interface are calculated by appealing to (18), (19), (46), (47), (64) and (77). The results read U10 = −Kh, U30 = 2hE and

Ω20 = Γ20 Fig. 5. Velocity component u1 (a) and w2 (b) for the oblate spheroid if λ = 0 (—), λ = 0.3 (), λ = 0.5 (◦), λ = 0.7 (△), λ = 1 (∗), λ = 2 (N), λ = 5 (•) and λ = ∞ (solid line). The horizontal solid line in (a) indicates the asymptotic value of u1 for d = h/a3 large.

rotation occurs for the given separation d. As predicted by (58), w2 quickly decays as d increases. By virtue of (61), r2 vanishes for λw = 2(1 + s2 )/[3(s2 − 1)] provided that s = a3 /a1 > 1 (prolate spheroid). For the selected prolate spheroid s = 3/2 and thus λw = 26/15. Finally, observe that except when close enough to a very viscous interface the prolate spheroid translates and rotates slower than the sphere for the same value of d. Since its rotation is found to exhibit a specific behaviour, the case of the oblate spheroid is reported in another figure and for additional values of λ. The velocity u1 admits the same behaviour as the sphere and the prolate spheroid (with this time u1∞ = 1.1421) and the results displayed in Fig. 5(a) are not discussed. By contrast, the oblate spheroid does not at all rotate as a sphere close to the interface and for large values of λ. Indeed, w2 this time tends to zero for a vanishing spheroid–interface gap. This is understood if one thinks about the limit case of a disk located parallel to a solid wall: the no-slip and therefore zero normal velocity condition at a solid wall clearly prevent the disk from rotating. The velocity w2 is always positive even for large values of λ and such a property, already observed close to a solid wall in [24], is predicted far from the interface by (61) which shows that r2 is strictly positive whatever λ for the oblate spheroid. The angular velocity is positive, vanishes close to and far from the interface and admits for any given viscosity ratio λ a maximum value at h/a3 ∼ 1.33.

=

Ks2 1 + s2

,

16π K µ2 a31 s3 3(χ ′ − α1 )

16π K µ2 a1 sh

, χ ′ + α1 32π E µ2 a1 sh F30 = . χ ′ + s2 α3

F10 = −

,

(65)

The influence of the interface is thus quantified by introducing the normalized velocities u1 , w2 , u3 ; forces f1 , f3 and torque γ2 u1 =

γ2 =

U1

,

w2 =

, 0

u3 =

U10

Γ2 Γ2

Ω2 Ω20 U3

,

, 0

U3

f1 = f3 =

F1 F10

,

F3

(66)

. 0

F3

As soon as the spheroid becomes distant from the interface the asymptotics established in Section 4 and (64) easily provides as leading approximations u1 = 1 + O(ϵ 2 ), w2 = 1 + O(ϵ 2 ), u3 = 1 + O(ϵ 2 ) and also f1 = 1 − A1 ϵ + O(ϵ 2 ), f3 = 1 − A3 ϵ + O(ϵ 2 )

γ2 = 1 + (1 + t 2 )r2 ϵ + O(ϵ 2 ),

(67)

with ϵ = a3 /h and A1 , A2 , r2 given by (60)–(61). For a sphere, (67) agrees with Yang and Leal [19] where refined expansions of f1 , γ2 and f3 are given up to order O(ϵ 4 ). Such expansions are found to be in perfect agreement with our numerical results for ϵ ≤ 0.5 (i.e. for d = h/a3 ≥ 2) but of course break down close to the interface. The largest differences with the computed values when d < 2 occur for small or large values of the viscosity ratio λ for γ2 or (f1 , f3 ), respectively. This feature is illustrated in Fig. 6(a) for λ = 0.1 and λ = 10. Note that Yang and Leal [19] provide additional plots for the (approximated) functions f1 , f3 and γ2 but not for the normalized velocities u1 , u3 and w2 . We therefore also draw in Fig. 6(b) those normalized velocities versus the separation parameter d = h/a3 for a few relevant values of λ. As asymptotically predicted, those velocities quickly asymptote to unity for d ≫ 1 (actually, as soon as d ≥ 5). Both u3 and

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

Fig. 6. Numerical results for a sphere. (a) Computed values of f1 (◦) and f3 (•) for λ = 10 and γ2 (∗) for λ = 0.1. Each dashed line reports the result predicted by the high-order asymptotic analysis of Lee et al. [18]. (b) Velocity u1 for λ = 0.3 (), λ = 1 (∗), λ = 5 (•) and velocities u1 (no symbols) u3 (∇) and w2 (△) for λ = 0 (dashed lines) and λ = ∞ (solid lines).

w2 are less than unity and decrease as the sphere approaches the interface. For a given sphere’s location u3 and w2 weakly depend upon λ and decrease or increase with λ, respectively. For a prescribed value of d the velocity u1 is by contrast deeply sensitive to the viscosity ratio with u1 = 1 for a critical value λc and u1 > 1 or u1 < 1 if λ < λc or λ > λc , respectively. Each spheroid is found to behave as a sphere when immersed in the pure straining flow: as depicted in Fig. 7, the quantities f3 − 1 and 1/[u3 − 1] are indeed strictly positive and increase with either λ or 1/d. The interactions with a close interface strongly slow down a spheroid and deeply depend upon its nature (oblate or prolate one): note that for d ≤ 2 the coefficients f3 (u3 ) are much larger (smaller) for the oblate spheroid than for the prolate one. Our numerical results for the ambient shear flow are given in Figs. 8 and 9. The coefficients f1 (see Fig. 8(a)) and u1 (see Fig. 9) exhibit similar trends for the oblate and prolate spheroids. Both quickly asymptote the unity far from the interface (here as soon as d ≥ 5), are positive and either larger or smaller than unity depending on the value of (d, λ). More precisely, a spheroid close to the interface may either translate faster (weak enough viscosity ratio) or slower (strong enough viscosity ratio) in the absence of interface. The normalized torque γ2 (see Fig. 8(b)) and angular velocity w2 (see Fig. 9) adopt, by contrast, different behaviors for the oblate and prolate spheroids. For a given value of d the

85

Fig. 7. Computed coefficients f3 and u3 for the prolate (dashed lines) and oblate (solid lines) spheroids for λ = 0 (), λ = 0 (∗) and λ = ∞ (). (a) Functions (f3 − 1)/5 for the oblate spheroid and f3 − 1 for the prolate spheroid. (b) Normalized velocity u3 .

normalized torque γ2 decreases or increases as λ increases for the prolate or oblate spheroid, respectively. In addition, the angular velocity w2 is strictly positive for the prolate spheroid but close enough to the interface vanishes and even may be negative for the oblate spheroid. There exists for d close to 1.1 a critical viscosity ratio at which the oblate spheroid translates without rotating parallel to the interface and faster in the absence of interface. 6. Conclusions A new procedure has been both proposed and tested to accurately compute at a reasonable cpu time cost the rigid-body motion of a solid and small particle located in the vicinity of a fluid–fluid interface (the cases of a gas–liquid and a solid–fluid boundary being obtained as special cases). The advocated method makes use of a few boundary-integral equations on the particle’s surface and permits one to cope in essence with arbitrarily shaped particles and arbitrary gravity field and/or external Stokes velocity fields. The proposed treatment allowed us to provide far from the boundary asymptotic approximations of the net force and torque exerted on a translating or rotating particle and of its incurred rigid-body migration when subject to a gravity or external flow. Different asymptotics then emerge. For instance, extending the results obtained by Lee et al. [18] for a sphere, a non-spherical

86

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

Fig. 8. Coefficients f1 and γ2 for the prolate (dashed lines) or oblate (solid lines) spheroid and λ = 0 (), λ = 0.5 (◦), λ = 0 (∗), λ = 5 (•) and λ = ∞ (). (a) f1 . (b) γ2 .

particle with typical size a subject to an external Stokes velocity field in the absence of gravity has translational velocity [1 + O(ϵ 2 )]U0 and angular velocity [1 + O(ϵ 2 )]0 with a/ϵ ≫ 1 its center-to-interface distance. The pure sedimentation of this particle however deeply depends on both its shape and the gravity orientation. For instance, homogeneous ellipsoids with a plane of symmetry parallel to a distant interface Σ have angular velocity of order O(ϵ 3 ) and translational velocity [1 + A3 ϵ]U03 if g is normal to Σ . For g parallel to Σ and aligned with one semi-axis and ei these spheroids have angular velocity of order O(ϵ 2 ) and a translational velocity [1 + Ai ϵ]U0i . The coefficients A1 and A2 have signs of 2 − 3λ and therefore vanish for the critical viscosity ratio λc = 2/3 for which the hydrodynamic interactions with the distance interface vanish even more quickly. For a close interface a numerical implementation has been achieved by using a boundary element code and tested against the analytical results of [18] for a sphere and asymptotic results for a distant spheroid. The computations for a sphere and either oblate or prolate spheroids in the vicinity of the interface show that the rigid-body motion of a particle deeply depends upon its shape, its position and the gravity and/or external velocity fields. Finally, one should note that the challenging case of two weakly non-Newtonian viscoelastic fluids might be treated in a similar fashion by resorting to the procedure advocated in Ho and Leal [30] and Chan and Leal [31].

Fig. 9. Coefficients u1 (solid lines) and w2 (dashed lines) for λ = 0(), λ = 0.5(◦), λ = 0(∗), λ = 5(•) and λ = ∞(). (a) Prolate spheroid. (b) Oblate spheroid.

Appendix A This appendix provides the regular functions R1kj , R2kj , Pj1 and Pj2 appearing in (25)–(27) for a source point y(y1 , y2 , y3 ) located in liquid 2, i.e. such that y3 < 0. For a field point x(x1 , x2 , x3 ) we denote by x′ (x1 , x2 , −x3 ) its image with respect to Σ and introduce the quantities ri = (x − y).ei , ri′ = (x′ − y).ei , r = |x − y| and r ′ = |x′ − y|. Recalling that λ = µ1 /µ2 , setting Jkj = δkj − 2δk3 δj3 and using the usual tensor summation convention over repeated indices, it is found, see [18], that

δkj + Jkj [(δki + Jki )ri ]rj (1 + λ)R1kj (x, y) = + r r3   x3 2 − 2 3 δj3 rk − δk3 rj + y3 − 3y3 rk rj /r for x3 > 0, r   4λ rj − δj3 x3 3y3 r3 rj Pj1 (x, y) = + 1+λ r3 r5 for x3 > 0,

(68)

(69) Jkj − λδkj

[(Jki − λδki )ri′ ]rj′

(1 + λ)R2kj (x, y) = + r′ r3   x3 + 2λ ′3 Jk3 rj′ − Jkj y3 − [Jki ri′ ]δj3 + 3y3 [Jki ri′ ]rj′ /r ′2 r

for x3 < 0,

(70)

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

2rj′

Pj2 (x, y) =

4λ

−



δj3 x3 + rj′

1+λ for x3 < 0. r ′3

r ′3

′

3x3 r3′ rj

+

r ′5 (71)

setting χ ′ = χ /a21 one arrives by elementary algebra at χ ′ = 2 and αi = 2/3 for a sphere (t = 1) and for oblate (t < 1) or prolate (t > 1) spheroids at the formulae

χ′ = √

t

χ = √

2t

Inside the particle P the flows and the rigid-body motion U +  ∧ O′ M are Stokes flows with stress (j) tensors σ(V2 , P2 ), σ 2,m and −P δij ei ⊗ ej where the pressure P is uniform. Applying in P the reciprocal theorem to these flows thus yields

∫ S

∫

(U +  ∧ O′ M).σ (2j,)m .ndS (x) = −P (j)

S

V2 .σ 2,m .ndS (x) =

∫

∫

(j)

S

u2,m .ndS (x),

(j)

S

u2,m .σ(V2 , P2 ).ndS (x).

(72) (73)

(j) u2,m

Since is divergence free in P the integral on the right-hand side of (72) vanishes and therefore one obtains, by virtue of (7) and (73),

∫ S

(j) u2 .σ 2,m .ndS (x)

∫

(j)

=− S

u2,m .σ(V2 , P2 ).ndS (x).

(74)

Setting f = [σ(V2 , P2 ) + σ 2 ].n on S and using the definition of Gkj (x, y) and (74), (30) and (33) then become

∫

[ul .ej ](y) = −

[f(x).ek ]Gkj (x, y)dS (x).

(75)

S

Finally, switching (x, k) and (y, j) in (75) and exploiting property (24) provide the announced result (34). Appendix C Upon introducing the function ∆(t ) = {(a21 + t )(a22 + t )(a23 + t )}2 and the quantities

χ = a1 a2 a3

∞

∫ 0

dt

∆(t )

,

α i = a1 a2 a3

∞

∫ 0

dt

(

a2i

+ t )∆(t )

(76)

(i)

(i)

the coefficients CT and CR read (see [27]), without summation over indices i in (77), (i)

(i)

CT = 4/[χ + a2i αi ],

CR = 4/[χ − a2i αi ].

(77)

i ,j ,2 BT

Each term is evaluated by appealing to the reciprocal identity and (42). More precisely, one gets i,j,2

BT

∫

(j),0

fR

=−

.uT(i),2 dS

S

(j)

∫

(ej ∧ O′ M).ek D(ki),l x′l s(M )dS

= CR

(78)

S

with the definition e i,i,0 ′ (i),l ′ Dk = 2 [AT δi3 δk3 δil − A3T ,3,0 δi3 δk3 δkl ] a

[ −

f + v + δi3 (g + v) a2

]

i,i,0

AT

δki δl3 .

(79)

Appealing to the identity

∫



t2 − 1

Appendix B

(V2 , P2 ), (u(2j,)m , p2(j,)m )

(ej ∧ O′ M).ek x′l s(M )dS = ϵlkj a2l Ve ,

(80)

S

one finally obtains (52)–(53). We conclude this appendix by paying special attention to spheroids such that a1 = a2 and a3 = ta1 . Upon

87

′

√

1 − t2

α1 = α2 =

log 2t 2 + 2t

arctan

t2 t2

−1

[ 1−



t2 − 1 − 1

1 − t2 t

χ′ 2t 2

]

,



if t > 1,

(81)

 if t > 1,

α3 =

1 t2

−1

(82)

[χ ′ − 2].

(83)

References [1] J. Happel, H. Brenner, Low Reynolds Number Hydrodynamics, Martinus Nijhoff, 1973. [2] S. Kim, S.J. Karrila, Microhydrodynamics: Principles and Selected Applications, Butterworth, 1991. [3] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge University Press, 1992. [4] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, Chem. Engng Sci. 16 (1961) 242–251. [5] A.J. Goldman, R.G. Cox, H. Brenner, Slow viscous motion of a sphere parallel to a plane wall. I. Motion through a quiescent fluid, Chem. Engng Sci. 22 (1967) 637–651. [6] W.R. Dean, M.E. O’Neill, A slow motion of viscous liquid caused by the rotation of a solid sphere, Mathematika 10 (1963) 13–24. [7] M.E. O’Neill, A slow motion of viscous liquid caused by a slowly moving solid sphere, Mathematika 11 (1964) 67–74. [8] M.E. O’Neill, A slow motion of viscous liquid caused by a slowly moving solid sphere: an addendum, Mathematika 14 (1967) 170–172. [9] M.E. O’Neill, K. Stewartson, On the slow motion of a sphere parallel to a nearby plane wall, J. Fluid Mech. 27 (1967) 705–724. [10] A.J. Goldman, R.G. Cox, H. Brenner, Slow viscous motion of a sphere parallel to a plane wall. II, Couette flow, Chem. Engng Sci. 22 (1967) 653–660. [11] H. Tozeren, R. Skalak, Stress in a suspension near rigid boundaries, J. Fluid Mech. 82 (1977) 289–307. [12] M. Chaoui, F. Feuillebois, F, Creeping flow around a sphere in shear flow close to a wall, Q. J. Mech. Appl. Math. 56 (2003) 381–410. [13] S.L. Goren, M.E. O’Neill, On the hydrodynamic resistance to a particle of a dilute suspension when in the neighborhood of a large obstacle, Chem. Engng Sci. 26 (1971) 325–338. [14] L. Pasol, A. Sellier, F. Feuillebois, A sphere in a second degree polynomial creeping flow parallel to a wall, Q. J. Mech. Appl. Math. 59 (2006) 587–614. [15] L. Pasol, M. Chaoui, S. Yahiaoui, F. Feuillebois, F, Analytical solutions for a spherical particle near a wall in axisymmetrical polynomial creeping flows, Phys. Fluids 17 (2005) 073602.1–073602.13. [16] S.H. Lee, L.G. Leal, Motion of a sphere in the presence of a plane interface. Part 2. An exact solution in bipolar co-ordinates, J. Fluid Mech. 98 (1980) 193–224. [17] A.V. Nguyen, G.M. Evans, Exact and global rational approximate expressions for resistance coefficients for a colloidal solid sphere moving in a quiescent liquid parallel to a slip gas–liquid interface, J. Colloid Interface Sci. 273 (2004) 262–270. [18] S.H. Lee, R.S. Chadwick, L.G. Leal, Motion of a sphere in the presence of a plane interface. Part 1. An approximate solution by generalization of the method of Lorentz, J. Fluid Mech. 98 (1979) 705–726. [19] S. Yang, L.G. Leal, Particle motion in Stokes flow near a plane fluid–fluid interface. Part 2. Linear shear and axisymmetric straining flows, J. Fluid Mech. 149 (1984) 275–304. [20] G.R. Fulford, J.R. Blake, On the motion of a slender body near an interface between two immiscible liquids at very low Reynolds number, J. Fluid Mech. 127 (1983) 203–217. [21] S. Yang, L.G. Leal, Particle motion in Stokes flow near a plane fluid–fluid interface. Part 1. Slender body in a quiescent fluid, J. Fluid Mech. 136 (1983) 393–421. [22] D.D. Joseph, Y.Y. Renardy, Fundamentals of Two-fluid Dynamics. Part I: Mathematical Theory and Applications, Springer-Verlag, 1991. [23] R. Hsu, P. Ganatos, The motion of a rigid body in viscous fluid bounded by a plane wall, J. Fluid Mech. 207 (1989) 29–72. [24] R. Hsu, P. Ganatos, Gravitational and zero-drag motion of a spheroid adjacent to an inclined plane at low Reynolds number, J. Fluid Mech. 268 (1994) 267–292. [25] L. Elasmi, M. Berzig, F. Feuillebois, Stokes flow for the axisymmetric motion of several spherical particles perpendicular to a plane wall, Z. Angew. Math. Phys. 54 (2003) 1–24. [26] A. Sellier, Settling motion of interacting solid particles in the vicinity of a plane solid boundary, C. R. Mécanique 333 (2) (2005) 111–116. [27] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. A 102 (1922) 161–179.

88

A. Sellier, L. Pasol / European Journal of Mechanics B/Fluids 30 (2011) 76–88

[28] D.E. Beskos, Introduction to Boundary Element Methods, Elsevier Science Publishers, 1998. [29] M. Bonnet, Boundary Integral Equation Methods for Solids and Fluids, John Wiley & Sons Ltd., 1999.

[30] B.P. Ho, L.G. Leal, Migration of rigid spheres in a two-dimensional unidirectional shear flow of a second-order fluid, J. Fluid Mech. 76 (1976) 783–799. [31] P.C.-H. Chan, L.G. Leal, A note on the motion of a spherical particle in a general quadratic flow of a second-order fluid, J. Fluid Mech. 82 (1977) 549–559.