Properties and solution techniques for a mixed type boundary integral equation arising in creeping flow problems

Computers & Fluids 64 (2012) 141–156

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Properties and solution techniques for a mixed type boundary integral equation arising in creeping flow problems A. Ramachandran ⇑, K. Tsiglifis 1, L.G. Leal Department of Chemical Engineering, University of California at Santa Barbara, Santa Barbara, CA 93106, USA

a r t i c l e

i n f o

Article history: Received 3 June 2011 Received in revised form 12 February 2012 Accepted 8 April 2012 Available online 17 May 2012 Keywords: Creeping flow Boundary integral method Slip Ill-posedness Hebeker operator Richardson iteration Condition number

a b s t r a c t We consider the properties of a creeping flow, boundary integral operator which is the sum of a single layer potential arising from a surface density, and two double layer potentials, one each for the tangential and normal components of the same density with different coefficient weights. Specific cases of this operator are encountered quite frequently in the creeping flow literature, especially in the solution of problems with slip at an interface. We examine, in particular, the eigenspectrum of this mixed-type operator. As expected, when the coefficient of either the tangential or normal component of the double layer term is either much smaller in magnitude than unity or identically equal to zero, the solution of this operator for the unknown density is rendered ill-posed, and thus, solution with the Richardson iteration technique requires a prohibitively large number of iterations. However, a surprising result is that the solution via inversion (e.g. Gaussian elimination, LU decomposition, etc.), for the component of the density whose coefficient in the double layer term of the operator is non-zero, can actually be a well-posed problem. We verify our theoretical developments by considering a spherical viscous drop with a slipping interface suspended in a viscous liquid undergoing a uniaxial extensional flow, the interfacial slip being proportional to the local tangential stress. In this case, the operator is equivalent to the mixed type operator with the normal coefficient set to zero. We demonstrate that the solution of this operator for the tangential stress is a well-posed problem, and we delineate the parametric space for a stable solution scheme. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Physical phenomena often involve the motion of solid objects through fluids, or the motion of deformable interfaces separating two liquids, such as the lipid bilayer or the liquid–liquid interface (whether clean or decorated with surfactants). In order to study the dynamics of such phenomena, we require knowledge of the velocity field at the interface. In the creeping flow limit, given an ambient velocity field u1(x0), the disturbance velocity and stress fields in the bulk fluids are determined solely by the instantaneous interfacial position. A powerful technique for determination of the interfacial velocities for a given shape of the interface is the Boundary Integral (BI) method [1–4]. In this technique, the creeping flow equations are reduced to an integral form giving the velocity at any point x [1–4] in the suspending fluid in terms of the stress and velocity at the boundaries of the flow domain. When evaluated in the limit as x ? x0 on the

⇑ Corresponding author. Present address: Department of Chemical Engineering and Applied Chemistry, University of Toronto, Toronto, Ontario, Canada M5S 3E5. E-mail address: [email protected] (A. Ramachandran). 1 Present address: Department of Mechanical & Industrial Engineering, University of Thessaly, Leoforos Athinon, Pedion Areos, 38834 Volos, Greece. 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.04.008

boundary, this reduces to an equation for the surface velocity field u(x0) exterior to the interface,

1 1 uðx0 Þ ¼  2 8p

Z

Gðx; x0 Þ  fðxÞdSðxÞ þ

S

1 8p

Z

PV

uðxÞ S

 Tðx; x0 Þ  nðxÞdSðxÞ þ u1 ðx0 Þ;

ð1:1Þ

and if applicable, to similar equation for the surface velocity field ^ ðx0 Þ interior to the interface. u

1 1 ^ ðx0 Þ ¼ u 2 8p k

Z

Gðx; x0 Þ  ^fðxÞdSðxÞ 

S

 Tðx; x0 Þ  nðxÞdSðxÞ:

1 8p

Z

PV

^ ðxÞ u

S

ð1:2Þ

The first and second integrals appearing in both equations above are known as the single and double layer potentials respectively, by analogy with electrostatics. The surface S represents the combined surfaces of all the closed interfaces in the suspending ^  n are denoted as f(x) and fluid. The traction vectors r  n and r ^fðxÞ, respectively, where r and r ^ are the stress tensors corresponding to the Newtonian suspending and interior fluids, and n is the outward pointing normal. The tensorial functions G(x, x0) and T(x, x0) are fundamental solutions to the creeping flow equation corresponding to the application of a point force at a position x0 in an unbounded medium:

142

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

I rr Gðx; x0 Þ ¼ þ 3 ; r r rrr Tðx; x0 Þ ¼ 6 5 ; r

ð1:3Þ ð1:4Þ

Here I is the unit tensor, r = x  x0 and r = jrj. The superscript PV symbolizes the principal value of the double layer potential, implying that the integral is evaluated with the point x0 located on the surface S. In subsequent text, we will drop the PV symbol, and it will be assumed that x0 is on the surface S. The solution of Eq. (1.1) and/or (1.2) requires boundary conditions for either f or u on the interface S. For flows involving solid bodies, the velocity at the solid–liquid interface is restricted to rigid body motion, i.e. a uniform translational velocity plus a rigidbody rotation. Therefore, the BI equation (1.1) reduces to the determination of the stress f at the solid–liquid interface via the solution of an integral equation of the first kind.

1 8p

Z

Gðx; x0 Þ  fðxÞdSðxÞ ¼ bðx0 Þ;

ð1:5Þ

S

The Eq. (1.5) has been successfully employed in the past in several different creeping flow situations [5–13]. Although the simplest implementation of the BI method for the simplest class of problems involving rigid, no-slip boundaries leads to the first-kind integral equation presented in 1.5, for almost all other circumstances (which generally involve deformable interfaces), the BI method yields either a second-kind integral equation, or mixed type integral equation that involves both single and double layer contributions. For example, the formulation of the BI equation for a fluid–fluid interface was shown by Rallison and Acrivos [15] to lead to an integral equation of the second kind,

uðx0 Þ 

b 4p

Z

PV

uðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ ¼ bðx0 Þ:

ð1:6Þ

S

The eigenspectrum of this equation has been studied extensively, and the stability of solving the equation for a given b via Picard iterations has been examined in the literature [3,4,16]. The mixed-type integral equation has been encountered in the literature in two important contexts. The first is the resolution of the ill-posedness of the Fredholm integral equation of the first kind in Eq. (1.5). To circumvent this problem of ill-posedness, alternative, well-behaved formulations of the BI equations for creeping flow have been examined in the literature. Generally these formulations are based upon the idea that writing the BI equations in terms of fictitious variables, related to the ‘‘natural’’ variables u and f, one can get an integral equation of the second kind for a solid–liquid interface. This approach is called the completed double layer boundary element method (CDLBEM) [3]. One example of this formulation is the method of Hebeker, for the no-slip, rigid boundaries [4,14,17–20]. The Hebeker approach requires the solution of an equation involving the fictitious variable q (instead of the natural variable f):

qðx0 Þ 1 þ 2 8p

Z S

PV

qðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ þ

bH 8p

Z

Gðx; x0 Þ

S

 qðxÞdSðxÞ ¼ bðx0 Þ:

ð1:7Þ

The scalar bH is a positive constant. The above hybrid operator, which involves both single and double layer terms, effectively allays the twin issues of the ill-conditioned nature of the first kind, single layer operator and the range deficiency of the second kind, double layer operator. It has been successfully used [17–20] to calculate stress on solid surfaces and to study drop-solid interactions even for very close separations between drops and solid surfaces, which requires fine mesh resolution.

A second class of problems which yields mixed-type BI equations is that involving slip of the classic, Navier type at either an interface or a rigid boundary [21–23]. The problem of slip at a solid–liquid interface [22,23] requires the solution of the following integral equation for the surface stress f:

Z bðtÞ HðxÞ  fðxÞ bðtÞ PV  ½HðxÞ  fðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ 2 8p S Z 1 Gðx; x0 Þ  fðxÞdSðxÞ ¼ bðx0 Þ; þ 8p S

ð1:8Þ

where b is a known quantity, and the tensor H is defined as

HðxÞ ¼ I  nðxÞnðxÞ:

ð1:9Þ

Slip is also known to occur at the interface between two high molecular weight polymeric liquids due to the low viscosity of the diffuse, low-entanglement-density interface between the two liquids [24–30]. The solution procedure for this problem also requires calculation of the tangential component of the interfacial stress, and the form of the governing equation is identical to (1.8) except for a positive sign in front of the second term. The integral Eq. (1.8) is of mixed type - first kind for normal component, second kind for tangential component. Also, while the constant bH in the Hebeker operator can be chosen arbitrarily, in the mixed-type operator, the constant b(t) cannot be set arbitrarily; it is fixed by the choice of the viscosity ratio of the fluids separated by the interface, and the slip coefficient [21–23]. A mixed-type integral equation would also arise in creeping flows involving bilayer interfaces, when intermonolayer friction [31] is included in the description of the physics of bilayer deformation. Motivated by the occurrence of different mixed-type BI operators in the literature, we examine the properties of a general mixed-type integral operator that may encompass all of these cases. This operator J± is defined as

n J q ¼

bðnÞ ½qðx0 Þ  nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ  qðx0 Þ

o

2 n o ðnÞ b ½qðxÞ  nðxÞnðxÞ þ bðtÞ HðxÞ  qðxÞ S Z 1  Tðx; x0 Þ  nðxÞdSðxÞ þ Gðx; x0 Þ  qðxÞdSðxÞ ¼ bðx0 Þ: 8p S ð1:10Þ

1  8p

Z

PV

The superscripts (n) and (t) of the constants b(n) and b(t) denote normal and tangential components. Depending on these coefficients, the above equation can be an integral equation of the second kind for either the normal component of q, or the tangential component of q, or both. It can also be an integral equation of the first kind if both b(n) and b(t) are set to zero. We note that the operator J+ reduces to the Hebeker equation for the variable q when b(n) = b(t), with the only difference that the constant now precedes the double layer terms instead of the single layer ones. It also reduces, when b(n) is set to zero, to the equation for the (unknown) tangential stress that appears in the case of boundary or interface slip [c.o. Eq. (1.7)]. We will henceforth refer to the operator J± with the normal coefficientb(n) set to zero as the slip operator JS±. The objective of this paper is to evaluate some commonly-used procedures for solving the problem Eq. (1.10) involving the operator J±. Because of the generality of the operator J±, we believe that the study will help researchers design/select stable BI routines involving any mixed type operator for a wide variety of physically-occurring creeping flows. The operator J± is linear; therefore, its behavior and solvability can be easily understood by scrutinizing its eigenvalues and eigenvectors. This relationship is explained in Section 2. We examine whether its eigenvalues are real or complex, and whether these eigenvalues are bounded. We also

143

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

examine the effect of the shape of the integration surface on the eigenvalues. Using the eigenspectrum elucidated in Section 2, we examine, in Section 3, the stability of a successive substitution method, specifically the Richardson iteration method [32], in solving the equation for q. In Section 4, we investigate the well-posedness of the operator J± in determining q for a given b, and delineate the ranges of the various parameters such as b(n), b(t), mesh size and type, and the shape of the integration surface, for the problem to be well-posed. We also explain why the use of the slip operator JS± has been successful in all prior investigations mentioned above, despite the fact that the operator has a large condition number and appears to be ill-posed. Finally, in Section 5, we consider the problem of a spherical drop with interfacial slip placed in a uniaxial extensional flow. Using the results developed in Sections 2–4, we predict the parameter range over which numerical schemes based on the first kind and mixed kind formulation for the surface stress are stable, and confirm these predictions with numerical computations. We summarize our findings in the Section 6. 2. Properties of the operator J± Since one of the major objectives of this manuscript is to understand the solvability of Eq. (1.10), we begin this section by briefly reviewing the relationship between the ill-posedness of an operator and its eigenvalues and eigenfunctions [see [3] for a more detailed discussion of operator ill-posedness specifically in the context of linear integral operators]. Let us first define the eigenfunctions and eigenvalues of the linear operator J±. To get the eigenspectrum, we need to solve the problem

J  /ðx0 Þ ¼ K/ðx0 Þ;

ð2:1Þ

where K and /(x0) represent the eigenvalue and eigenvector of the operator respectively. Note that in addition to the double layer parameters b(n) and b(t), the eigenspectrum is also a function of the surface S over which the integration is carried out. To determine the well-posed or ill-posed nature of the problem in Eq. (1.10), we need to understand the effect of introducing a ~ to the right hand side of this equation on small perturbation b ~ that satisfies the the computed q, i.e. we need to determine the q equation

~ ~ Þ ¼ b þ b: J  ðq þ q

ð2:2Þ

We can then use the definition of the eigenvalue and eigenfunction ~ as of J+ in Eq. (2.1) to write q

~¼ q

1 ~ X bj j¼1

Kj

/j ðx0 Þ

ð2:6Þ

We can see from Eq. (2.6), that if a small error e/j (x0) is introduced on the right hand side of Eq. (1.10), the solution q(x0) is augmented by the quantity e/j(x0)/Kj. Thus, the inverse of any eigenvalue of J+ may be viewed as the amplification factor for errors of the form of the corresponding eigenfunction of J±. The smaller the eigenvalue, the larger the amplification of the error of its mode. This approach also helps us understand separately the amplification of the normal and tangential components of the error. Quite clearly, the ill-posedness of a linear operator, and therefore its solvability depends crucially on its eigenspectrum. Before we proceed to determine the eigenspectrum of the operator J±, we should mention that, in the course of performing the calculations that follow, we have discovered that almost all properties of the operator J are similar in trend, and on occasions, even quantitatively the same as those of J+. A notable exception is their null space and low frequency behavior, and this is described in detail in Appendix A along with details of how to deal with the zero eigenvalues by deflating the operator. For the remainder of the paper, however, we have restricted our discussions to the operator J+, only pointing out those properties of J that are different from J+. Let us now consider the eigenvalues and eigenfunctions of J+. Recall that the eigenspectrum actually depends on the shape of the integration surface S in J+. For simplicity, we will begin with a single spherical surface of unit radius, which allows some analytical progress. 2.1. Eigenspectrum of the operator J+ on a sphere The co-ordinate system used for describing the spherical surface is shown in Fig. 1. In Appendix B, we have delineated the procedure for determining the eigenvalues and eigenfunctions of the eigenvalue problem (2.1) for the operator J±. In short, the procedure is to establish a relationship between the velocity and stress functions from Lamb’s general solution of the creeping flow equations expressed in terms of decaying and growing harmonics (outside and inside the sphere, respectively) consistent with the eigenfunction problem (2.1) [3]. We have performed this calculation only for

Owing to the linearity of the operator J+, the above equation simplifies to

~ ~ ¼ b: J q

ð2:3Þ

A simple metric for the degree of amplification of errors by an operator is to consider the ratio of the relative error in q to the rel~ ~ k=kqkÞ=ðkbk=kbkÞ, ative error in b, i.e. ðkq where k  k is a suitable norm. The maximum value of this ratio is called the condition number [33], which, with the choice of norm as the L2 norm, is simply the ratio of the magnitudes of the largest eigenvalue of J± to its smallest one.

C ¼

jKmax j jKmin j

1 X j¼1

~ / ðx0 Þ: b j j

1

x

r = ρ 2 + x2 ð2:4Þ

The larger the value of the condition number, the greater the sensitivity of the solution to errors in b. A more fundamental way of examining the amplification of er~ into a linear combination of the rors is to decompose the error b eigenfunctions of J±.

~¼ b

s = l, θ = π

Axis of symmetry

ð2:5Þ

t

n

ρ L

θ

C s = 0, θ = 0

s

Fig. 1. Schematic of the spheroidal geometry used for the calculations in this paper. The aspect ratio of the spheroid is L. The cylindrical co-ordinates x(s) and q(s) are functions of the curvilinear co-ordinate s measured from one end of the object along the curve C. The total length of the curve C is l. The unit vectors n and t are the local normal and tangent vectors to the curve C. For a sphere, C is simply a circle, and the aspect ratio a is unity. In this case, we prefer to use the spherical co-ordinates r and h (we ignore the azimuthal co-ordinate due to the assumed axisymmetry).

144

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

operator for axisymmetric modes on a sphere have a rather simple analytical form:

axisymmetric modes due to the simplification involved; however the conclusions and general trends from this analysis can ordinarily be carried over to the 3-D case [3]. The procedure outlined in Appendix B produces two sets of eigenvalues and eigenfunctions: KAm ; /Am ðx0 Þ and KBm ; /Bm ðx0 Þ, where m is the frequency in the meridional direction of the sphere. Analytical solutions are possible for the eigenspectrum, and have been computed using the symbolic toolbox in MATLAB, but these are too long to be reproduced here. Instead, we calculate the eigenvalues from these analytical expressions and show them in Fig. 2. The asymptotic forms of the eigenvalues for small and large m for the operator J+ are presented below

6bðnÞ þ 3bðtÞ þ 1 ;...; 15  bðtÞ bðnÞ jbðtÞ  bðnÞ j 1 1 þ KAm m1   þ þ 4 4m 8m2 4 4 maxðbðnÞ ; bðtÞ Þ 1 1 ¼  ; þ 2 4m 8m2 ðtÞ ðnÞ  2 b b jbðtÞ  bðnÞ j 1 1 þ KB1 ¼ ; . . . ; KBm m1   þ  3 4 4m 8m2 4 4 minðbðnÞ ; bðtÞ Þ 1 1 ¼  ; þ 2 4m 8m2

ðm þ 1Þ½2bH ðm  1Þ þ 1 ; m ¼ 1; 2; 3; . . . ð2m þ 1Þð2m  1Þ b ð2m2 þ 4m þ 3Þ þ m ¼ H ; m ¼ 0; 1; 2; . . . ð2m þ 1Þð2m þ 3Þ

KH m A ¼ KH m B

The eigenvalues presented in Fig. 2 and Eqs. (2.7) and (2.8) provide some key pieces of information. The eigenvalues are all real, and this can be proven rigorously for a spherical surface using an energy dissipation argument (see Appendix C). Both sets of eigenvalues decrease smoothly with m [except possibly for the first few eigenvalues (small m)], and saturate to constant values of max (b(n)/ 2,b(t)/2) and min (b(n)/2,b(t)/2) for large m, provided these constants are not zero. If either b(n) or b(t) is zero, (e.g. the slip operator J Sþ , which has b(n) = 0) then one eigenvalue continues to decrease monotonically with m for large m, following a 1/4m behavior. Let us now examine the eigenfunctions of the operator J+ on a sphere. These eigenfunctions are of the form (see Appendix B)

KA0 ¼ bðnÞ ; KA1 ¼

ð2:7Þ

/m ¼ C rm Pm ðcos hÞir þ C hm

ð2:10Þ

β ( n ) = 0.1, β (t ) = 0.1 (blue)

β ( n ) = 0,

β ( n ) = 0.1, β (t ) = 0 (red)

β ( n ) = 0,

β (t ) = 0.1 by red curve)

β ( n ) = 0, β (t ) = 0 (black)

β ( n ) = 0, β (t ) = 0 (black)

ν

ν

(a) β ( n ) = 0,

β ( n ) = 0.1, β (t ) = 0 (red) β ( n ) = 0.1,

β ( n ) = 0.1, β (t ) = 0.1 (blue) β ( n ) = 0,

β ( n ) = 0,

β (t ) = 0 (black)

β (t ) = 0 (black)

β ( ) = 0.1 (green) t

β (t ) = 0.1 (blue) β ( n ) = 0, β (t ) = 0.1 (green)

ν

Dˆν A

ΛνB

β (t ) = 0.1 (green)

(green, hidden

Cˆν A

β ( n ) = 0.1,

ð2:9Þ

^ m ¼ /m / /0 m

β ( n ) = 0.1, β (t ) = 0.1 (blue) β ( n ) = 0.1, β (t ) = 0 (red)

Λν

dPm ðcos hÞ ih dh

The eigenfunction /m derived above can be normalized as

For the special case of the Hebeker operator, where the tangential and normal coefficients of the double layer term are the same and equal to bH, the solutions for the eigenvalues of the Hebeker

A

CˆνB

β ( n ) = 0.1, β (t ) = 0 (red)

ν

(b)

β ( n ) = 0, β ( n ) = 0, β (t ) = 0 (black) β ( t ) = 0.1 (green)

β ( n ) = 0.1,

β (t ) = 0.1 (blue)

β (t ) = 0.1 (blue)

β ( n ) = 0.1, β (t ) = 0 (red)

ν

ð2:8Þ

β ( n ) = 0, β ( n ) = 0.1, (t ) (t ) β = 0 (red) β = 0 (black) Dˆν B

β ( n ) = 0, β (t ) = 0.1 (green)

(c)

ν

b m and (c) coefficients of tangential component D b m for the operator J+ on a sphere as functions of the Fig. 2. The (a) eigenvalues Km, (b) coefficients of the normal component C spatial frequency in the meridional direction.

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

The normalization factor /0 is chosen such that neither the tangen^m tial nor the radial component of the normalized eigenfunction / can exceed a magnitude of 1 for h 2 [0, p].

n o /0m ¼ max /r0m ; /h0m

ð2:11Þ

where

       A;B  /r0m ¼ maxC A;B r m P m ðcos hÞ ¼ C r m  h      dPm ðcos hÞ dPm ðcos hÞ  A;B    ¼ C /h0m ¼ maxC A;B  max hm hm   h h  dh dh

ð2:12Þ

^ is The normalized eigenvector /

b m Pm ðcos hÞir þ D bm ^m ¼ C /

dPm ðcos hÞ=dh ih max½dPm ðcos hÞ=dh

ð2:13Þ

h

The reason for the choice of the above normalization will become clear in Section 4, where we analyze the well-posedness of the operator J±. We will now make some observations of the behaviors of the b m and the tangential coefficient D b m of the nornormal coefficient C ^ m , which will be useful in Section 4. For the malized eigenvector / b m and D b m are shown in Fig. 2b and c operator J+, the coefficients C as functions of the spatial frequency m for different combinations of the parameters b(n) and b(t) of the double layer term. The figure shows that when b(n) and b(t) are both non-zero (see the b(n) = 0.1 and b(t) = 0.1 case, blue curves), the normal and tangential coeffib m and D b m respectively) for both sets of eigenvectors apcients ( C proach a non-zero value asymptotically as the frequency m becomes high.

145

In the case of the slip operator J Sþ , (normal coefficient, b(n), is zero, but the tangential coefficient b(t) is non-zero), the normal b A corresponding to the eigenvalue KA decreases monocoefficient C m m tonically with m as m2, but this monotonic decrease only begins after a certain value of m (see the green curves in Fig. 2). On the b A assumes a constant value other hand, the tangential coefficient D m as m increases. For the second set of eigenvalues KBm which de1 creases monotonically as m for large m, the normal component b B is invariant with m, while the tangential component coefficient C m b B decreases asymptotically as m2 for large m, but again, coefficient D m only after a cutoff value of m. This cutoff value of n increases as b(t) becomes smaller. If the tangential coefficient b(t) is zero, but b(n) is not (the opposite of the case above, and represented by the red curves in Fig. 2), then the behaviors of the normal and tangential components in the b A and D b B reach constant values, limit of large m are reversed: C m m 2 A B b b while D m and C m decay as m after a cutoff value. If b(n) and b(t) are both set to zero (black curves in Fig. 2), then all four coefficients assume constant values with m. 2.2. Eigenspectrum of the operator J+ on other axisymmetric surfaces In this section, we examine the eigenspectrum of J± for nonspherical, axisymmetric surfaces. In particular, we consider spheroids of aspect ratio L, as shown in Fig. 1. The eigenspectrum was determined numerically in MATLAB by (a) discretizing the curve C describing the axisymmetric surface into nodes, (b) interpolating the shapes and the eigenfunctions with cubic splines,

Fig. 3. The eigenvalues of the operator J+ as a function of the number of zeros, n0, in the normal component of the corresponding eigenvectors for different aspect ratios [circles L = 1 (spheres), squares L = 5 (oblate spheroid), triangles L = 0.2 (prolate spheroid)]. Subfigures (a)–(d) show these trends for (b(n), b(t)) combinations of (0, 0), (0, 0.1), (0.1, 0) and (0.1, 0.1) respectively. The surface in each case was discretized with 601 equispaced nodes.

146

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

(c) computing the integrals over each element using Gaussian quadrature (after appropriate desingularization of leading order integrand singularities, e.g. see [21]) and (d) using the eig command in MATLAB to determine non-trivial solutions of the homogeneous system of linear equations resulting from the above described discretization of Eq. (1.9). In Fig. 3, we have shown the eigenvalues of the operator J+ for different aspect ratios, and different combinations of the double layer parameters b(n) and b(t). Since there is no obvious choice for the frequency, the abscissae for these plots is taken to be the number of zeros, m0, in the radial component of the eigenvector corresponding to an eigenvalue. We observe trends similar to those already noted for a sphere. When both b(n) and b(t) are zero, all the eigenvalues decrease monotonically as the frequency is increased [subfigure (a)]. When one of b(n) or b(t) is zero [subfigures (b) and (c)], one set of eigenvalues decreases monotonically, while the other asymptotes to a constant value. The eigenvalues asymptote to constant values for non-zero b(n) and b(t) [subfigure (d)]. Although we have not separately shown the eigenfunctions, they also show the same qualitative behavior as the eigenfunctions on a sphere. 3. Convergence in a successive substitution numerical scheme 3.1. The condition number For the operator J+, we will not consider the m = 0 eigenvalue for defining the condition number, because we recognize that this eigenvalue can be adjusted by Wielandt deflation (see Appendix A). The condition number for the operator J+ on a sphere, assuming that the monotonically increasing or decreasing behavior begins with the m = 1 eigenvalue, may be approximated using Eq. (2.7) as

  ðnÞ ðtÞ max jbðnÞ j; j maxðb2 ;b Þj ; 23   Cþ  : ðnÞ ðtÞ min j minðb2 ;b Þj ; 23

ð3:1Þ

3.2. Stability of the method of successive substitutions Let us now use the condition number to predict whether it is reasonable to solve Eq. (1.10) via a Richardson iterative technique (method of successive substitutions). Consider the solution of Eq. (1.10) via the following Richardson iteration [32]

(i)

ð3:2Þ

(i+1)

Here, q and q are the solutions at the ith and (i + 1)th iteration, and k is a scalar whose choice will be discussed shortly. If we expand the vector b in terms of the eigenfunctions of J+ as

b¼

X bj /j j

ðiÞ

q

" # X 1  ð1  kKj Þi ¼ bj /j Kj j

ð3:3Þ

ð3:4Þ

The above iterative technique will converge to the solution if the condition j1  kKjj < 1 is satisfied for all eigenvalues, or equivalently, if

0 < kKj < 2 for all j:

ð3:5Þ

An upshot of the above condition is that if a Richardson iteration of the form (3.2) is to converge, all the eigenvalues Kj must have the same sign. If the eigenvalues are both positive and negative, then no single value of k can satisfy condition (3.5). (This would be true, for example, for the operator J+ when either b(n) or b(t) is negative). Eqs. (3.5) and (3.4) also indicate that if the rate of convergence is to be high, then kKj must be as close to unity as possible for all eigenvalues. This is where the condition number C [Eq. (2.4)] becomes important. If C is close to one, then all the eigenvalues lie within a narrow band, and k can be tuned so that both kmax (Kj) and kmin (Kj) are close to unity. The Hebeker operator is thus a great candidate for solution via Richardson iterations; because bH is an adjustable parameter, its eigenvalues [see Eq. (2.8)] can be made to lie within a narrow range. On the other hand, in the slip operator J Sþ , the constant b(t) is fixed by the physics of the problem, and we have no freedom to tune this constant to satisfy kKj  1 for all eigenvalues. This can lead to a poor rate of convergence while solving the slip operator at high condition numbers using the Richardson iteration method. For example, for the operator J±, small values of b(n), or b(t), or both, or highly disparate values of b(n) and b(t) would lead to slow convergence rates, especially when the forcing function b in Eq. (1.10) comprises eigenmodes of high frequencies. The Richardson iterative technique should, therefore, be avoided for such operators. As an example, let us apply the Richardson iteration technique to solve the problem J+q = b on a sphere, with b taken as

b ¼ cosðmhÞir þ sinðmhÞih ;

One can see that if either b(n) or b(t), the coefficients of the normal and tangential components of the double layer term in the our operator, become small (much less than 4/3), C þ becomes large. Accurate solutions for such linear systems require resolution of details corresponding to both the largest and smallest eigenvalues, and are, therefore, expensive to solve numerically. Also, if the values of b(n) and b(t) are very different [e.g. b(n) = 0.1 and b(t) = 10], the condition number is high, as indicated by the large v asymptotes in Eq. (2.7). Small values of b(n) or b(t), or highly disparate values of b(n) or b(t) lead to large condition numbers, which render the inverse problem for this operator numerically expensive. Similar trends can be expected for the operators J+ on non-spherical surfaces.

qðiþ1Þ ¼ qðiÞ  k½J þ qðiÞ  b

and assume q(0) = 0 (a Neumann series), then we can write q(i) as

ð3:6Þ

(0)

and the initial guess q chosen to be 0. Since we realize that convergence issues arise only at high frequencies, we choose a relatively high frequency for the forcing function  m = 40. The discretization of the integrals is performed with the same procedure as described in section 2C. The number of nodes used to represent the axisymmetric curve C of the sphere (see Fig. 1) is 401, which allows each wavelength of the forcing function to be described by about 20 points. In this exercise, we will examine the variation of the error d with the iteration number, where d is defined as sum over all node points of the L2 norm of the difference between the solution q(i) at the ith iteration, and the actual solution q⁄, computed by Gaussian elimination.

d¼

 401  X  ðiÞ  qp  qp  p¼1

2

ð3:7Þ

For each combination of the operator parameters b(n) and b(t), we will employ a value of k of 0.1, which is approximately the value beyond which convergence could not be observed in the simulations in Fig. 4. The reason that this value is nearly constant for all the simulations is that the maximum eigenvalue, which determines the value of the largest allowable k, varies weakly with b(n) when b(n) is less than 0.5. In our simulations, b(n) is always 0.5 or smaller. In Fig. 4, we show d as a function of the Richardson iteration number for different values of the operator parameters b(n) and b(t). First, we note that if one of b(n) and b(t) is negative [e.g. b(n) =  0.5 and b(t) = 0.5], d increases monotonically with iterations,

147

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156 β ( n ) = 0,

β ( n ) = −0.5,

β (t ) = 0.5,

β (t ) = 0.5,

β (t ) = 0.5,

k = 0.05

k = 0.05

k = 0.05

0

δ

10

β ( n ) = −0.5,

β ( n ) = 0.05, β ( n ) = 0.5,

β (t ) = 0.5,

β (t ) = 0.5,

k = 0.05

k = 0.10

10

-5

β ( n ) = 0.5, β (t ) = 0.5, k = 0.05

0

100

β ( n ) = 0,

β ( n ) = 0.05,

β (t ) = 0.5,

(t )

β = 0.5, k = 0.11

200

k = 0.05

300

400

500

Iteration # Fig. 4. The variation of the error d defined in Eq. (3.7) with the Richardson iteration number as a function of the parameters b(n) and b(t) of the operator J+ for the forcing function b defined in Eq. (3.6) with m = 40. For each combination of b(n) and b(t), the Richardson iteration parameter k is taken to be the maximum value beyond which the method fails to converge (except for the combination b(n) = 0.5 and b(t) = 0.5, where convergence is not observed for any value of k).

and as can be seen in Fig. 4, this result is independent of how small the value of k is. This is consistent with our earlier assertion that convergence with Richardson iterations cannot be achieved when both b(n) or b(t) are not positive. Next, if we fix b(t) at 0.5, and vary b(n) from 0.5 to 0, we observe that more and more iterations are required to reach the same error; the rate of convergence decreases progressively. Again, this agrees with our expectation that for small b(n), the eigenvalues are smaller, leading to slower convergence. 4. Well-posedness of J+ In this section, we take a more detailed look at the amplification ~ in the forcing function into a of errors by decomposing the error b linear combination of eigenfunctions, and by analyzing the result~ in the desired function q. As per our discussion in Secing error q tion 2, the inverse of the operator eigenvalues represent the error amplification factor for the corresponding eigenmode present in ~ We will refine the definition of error amplification factor here. b. 4.1. Well-posedness of the operator J+ on a sphere For the single layer operator [b(n) = b(t) = 0], we have seen that the inverses of both sets of eigenvalues grow without bound as v increases [for e.g., see Eq. (2.8) with bH set to zero], thus leading to stronger amplification for errors of higher frequencies. As described in the introduction, it was this ill-posedness of the single layer operator that caused researchers to explore operators with fictitious variables for solid particle flow problems. The asymptotic forms of the inverse eigenvalues of the Hebeker operator [b(n) = b(t) = bH] in Eq. (2.8), for example, suggest that as long as bH is an O(1) scalar, the amplification factor is always an O(1) quantity. This indicates that even very high frequency errors will not be amplified upon inversion of the Hebeker operator. This makes the Hebeker operator an excellent choice for solving creeping flow problems involving solid surfaces (e.g. see [17–20]). Let us turn our attention to the slip operator JS+ [b(n) = 0, b(t) non-zero]. As mentioned before, this is a mixed type operator; only the tangential component of the operand is described by an integral equation of the second kind; the normal component still satisfies an integral equation of the first kind. It may also be viewed as being intermediate between the single layer operator and the Hebeker operator. Quite understandably then, and also seen in Section 2, one set of inverse eigenvalues assumes a constant value (2/ b(t)) for large v like the Hebeker operator, while the other set behaves like the single layer operator eigenvalues and grows monotonically as 4m for large m. However, because of the presence of inverse eigenvalues that grow monotonically with m, the amplifica-

tion factor is expected to grow continuously as finer and finer meshes are employed to access higher frequency modes (large m). Ultimately, we expect that the solution should be completely swamped by the amplification of even relatively small, high frequency errors, suggesting that Eq. (1.10) for the operator JS± does not represent a well-posed problem. However, this is opposite to what is observed in the quite successful numerical implementations of the operator JS+ (and JS) in the literature [21–23]! The resolution of this contradiction lies in the realization that, in the simulations involving inversion of the slip operator JS+, we only require the tangential component  1 of the stress. Thus, even though the inverse eigenvalues KAm grow monotonically with m for a given error, the product of the tangential component of the eigenvectors with the inverse eigenvalue may still reach a stable value or even decrease for asymptotically large m. The above paradox and remarks also apply for the operator J+ with a non-zero b(n) and with b(t) set to zero. In this case, even though we again have eigenvalues that decrease monotonically with frequency and are expected to corrupt the tangential component of the inverted result, the normal component is well-behaved. To understand this better, we consider the slip operator JS+ on a ~ on the RHS of Eq. (1.10) can be decomsphere. Since any error b posed into a linear combination of the eigenfunctions of the operator JS+, it is instructive to examine the relative magnitudes of the tangential and radial components of the eigenfunctions as a function of h, the meridional co-ordinate on the surface of the axisymmetric spherical surface under consideration (see Fig. 1):

~ ¼ e/ ^ A: b m

ð4:1Þ

^ A [see Because of the normalization chosen for the eigenfunctions / m Eq. (2.11)], the characteristic magnitude of either the normal or tan~ for h 2 [0, p] is less than or gential component of the error vector b equal to e. The error in the solution upon inverting the operator JS+ is

~ J 1 Sþ b ¼ eir

bA C m A

Km

Pm ðcos hÞ þ eih

bA D m

dP m ðcos hÞ=dh ; KAm maxh ½dPm ðcos hÞ=dh

ð4:2Þ

Since we are interested in the tangential component, the amplification factor we should be examining is

AAhm ¼

bA D m

KAm

ð4:3Þ

These amplification factors are shown in Fig. 5. AAhm asymptotes to a constant value for large m and non-zero b(t). This is expected, beb A and KA reach cause we have already seen in Section 2 that both D m m constant values for high frequencies. However it should be noted that this constant value increases with a decrease in b(t) as 2/b(t). For example, for b(t) = 0.01, the asymptotic value of AAhm for large m is as high as 20. It is also interesting to examine the normal component AAr of the error:

AArm ¼

bA C m

KAm

ð4:4Þ

b A decreases as m2 for large v, The results in Section 2 show that C m A while Km reaches a constant value of 2/b(t). Thus AArm decays as m2 for large m, but not before it exhibits a maximum for intermediate m, as may be seen in Fig. 5. The smaller the value of b(t), the larger the maximum in AArm and the greater the value of m at which this maximum occurs. By an analogous procedure, we can examine errors of the form

~ ¼ e/ ^ B; b m

ð4:5Þ

148

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

β

A rν

and

(n)

A

= 0, β

θν

(t )

, =0 A

θν

,

β ( n ) = 0, β (t ) = 0.01

A

ν

A

A rν ,

,

θν

β ( n ) = 0, β (t ) = 0.1

β ( n ) = 0, β (t ) = 0.01

A rν ,

β ( n ) = 0, β (t ) = 0.1

ν Fig. 5. The amplification factors for the radial and tangential components,

AArm

and AAhm respectively, as defined in Eq. (4.6), corresponding to the error defined in Eq. (4.1).

^ B . The amplification factors involving the other set of eigenvectors, / m for this set of eigenvectors can be defined in an analogous fashion:

ABhm ¼

bB D m B

Km

;

ABrm ¼

bB C m

KBm

ð4:6Þ

:

The amplification factor for the tangential component ABhn is shown in Fig. 6 for different values of the operator’s tangential coefficient b(t). It may be seen that for non-zero b(t), the tangential amplification factor ABhm increases monotonically for small m, but reaches a maximum and decreases with further increase in m, with the crossover value of m decreasing with increasing b(t). As we saw in Section 2, for high frequencies, the tangential component of the b B , varies as m2, while the corresponding eigenvalue eigenvector, D m KBm decreases as m1. Therefore, the amplification factor for the tangential component ABhm varies asymptotically as m1 for large m, i.e. it decreases as the frequency of the errors increases. On the other hand, the normal amplification factor ABrm increases with m (see b B , asympFig. 6), since the normal component of the eigenvector, C m totically reaches a constant at large m, while KBn decreases as m1. Combining the observations for the behaviors of all the amplification factors AAhm ; ABhm ; AArm , and ABrm , we see that the problem of determining the tangential component by inverting the slip operator JS+ is well-posed only when the tangential coefficient of the operator, b(t), is O(1). For small values of b(t), the problem becomes increasingly ill-posed, because the maximum values of the tangential amplification factors corresponding to both sets of eigenvalues increase with decreasing b(t). It is notable, however, that the maximum value of the tangential amplification factor ABhm is smaller than the maximum value attained by AAhm for the same value of b(t). Also, the maximum occurs for a finite value of the frequency for ABhm , while for AAhm , the maximum occurs at asymptotically large ~ we expect that frequencies. Therefore, given an arbitrary error b,

the amplification of the tangential error will be dominated by the eigenspectrum of type A rather than type B. On the other hand, the problem of determining the normal component by inverting JS+ is always ill-posed, because ABrm grows continuously with the frequency n. If the constant b(t) is set to zero, but b(n) is non-zero, then the above comments are reversed for the normal and tangential components; the problem of determining the normal component by inverting J+ is well-posed provided b(n) is O(1) but the problem of determining the tangential component is always ill-posed. 4.2. Well-posedness of the operator J+ on non-spherical, axisymmetric surfaces In this section, we examine the effect of introducing an error on the right hand side of the operator J+ for more general, axisymmetric, integration surfaces. Since we do not have analytical expressions for eigenvalues and eigenvectors, the results presented in this section are from numerical simulations. To study the amplification of errors, we choose the following ~ for the mth node. form of the error b

~ ¼ e cosðmpÞ½nðsm Þ þ tðsm Þ b

Here, m varies from 1 to N + 1, while the vectors n and t are the local normal and tangent unit vectors along the curve C. As may be seen, ~ oscillates in magnitude between ±e over each successive the error b node describing the curve. This form ensures the imposition of an error of the maximum possible frequency, and is particularly useful when we consider non-uniform meshes. The amplification factors for the tangential and normal components of the inverse of the er~ can be defined as ror b

AðtÞ ¼ B rν ,

β ( n ) = 0, β (t ) = 0.1

B rν

β

(n)

( blue ) , = 0,

β (t ) = 0.01

B rν

( red ) ,

β ( n ) = 0, B

β (t ) = 0

θν

,

β ( n ) = 0, β (t ) = 0

B

ν

B

θν

,

β = 0, β (t ) = 0.01 (n)

B

θν

,

β ( n ) = 0, β (t ) = 0.1

ν Fig. 6. The amplification factors for the radial and tangential components, ABrm and ABhm respectively, as defined in Eqs. (4.3) and (4.4), corresponding to the error defined in Eq. (4.1).

ð4:7Þ

1

e

~ maxðjt  J 1 bjÞ;

AðnÞ ¼

1

e

~ maxðjn  J 1 bjÞ

ð4:8Þ

First, we check whether the numerical solutions for amplification factors match the predictions of the analytical theory developed for a sphere in the previous subsection. For the numerical implementation, we divide the curve C into N equal parts of elemental length Ds = l/N. We expect that the behaviors of AðnÞ and AðtÞ should be dominated by those of AArm and AAhm in Fig. 5 respectively, and this is precisely what we observe in Fig. 7. For non-zero tangential coefficients (b(t) > 0), the amplification factor for the normal component grows continuously with the number of points N used to discretize the spherical surface, while that of the tangential component attains an asymptotic value for large N. When bðtÞ ¼ 0, both AðnÞ and AðtÞ increase linearly with the number of points for large N. We can also examine the amplification factors as a function of b(t) for a fixed discretization, and this is done in Fig. 8. It may be seen that for a fixed number of N uniformly-spaced points, the

149

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

10

4

N = 400

β (t ) = 0

10

N = 40

~ 1/ β ( ) t

3

2

β (t ) = 0.01

(t )

10

N = 100

10 (t )

2

β (t ) = 0.1

10

1

10

20

40

80

160

320

640

1

10 -2 10

Number of points, N

(a)

-1

β

10

(t )

(a) 3000

10

3

N = 400

2500

(n)

2000 ( n)

10

2

10

1500 1000

20

40

80

160

320

640

Number of points, N

(b) Fig. 7. The amplification factors for the (a) tangential and (b) normal components of the error given in Eq. (4.7) for the operator J when the surface of integration is a sphere. In subfigure (b), the curves for the different values of b(t)[b(t) = 0 (red curve), b(t) = 0.01 (blue curve), b(t) = 0.1 (black curve)] are almost coincident. For these simulations, b(n) was set to 0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tangential amplification factor is a decreasing function of b(t), and for large b(t) and large m, the variation is 1/b(t). The normal amplification factors are virtually independent of b(t) for fixed N. Let us now examine the effect of using a non-uniform mesh to describe the surface S, a situation that arises frequently in problems such as the stretching or relaxation of single drops. The surface S is again a sphere, but the discretization used is a combination of two uniformly meshed lengths, l1 ¼ lf and l2 ¼ lð1  f Þ, where l is the total curvilinear length of the curve describing the axisymmetric surface S. The element lengths are Ds1 = 1/400 and D s2 = 1/40 in the two regions respectively. When the parameter f is zero, the mesh is uniform but coarse; when f is 1, the mesh is uniform again, but finer. For intermediate values of f , the mesh is non-uniform. Shown in Fig. 9 are the amplification factors for the tangential and normal components of the error in Eq. (4.7) as a function of f . It can be seen in Fig. 9 that amplification factors for the fine, uniform mesh, f ¼ 1, are greater than those for the coarse, uniform mesh, f ¼ 0, as expected from the results in Fig. 7. However, the transition from the amplification factors for f ¼ 0 to those for f ¼ 1 occurs for a relatively small value of f . This implies that introduction of even a small region with a large number of mesh points on the surface S will yield error amplification factors that correspond to the mesh spacing in this highly discretized region. ~ in In Fig. 10, we have shown amplification factors for the error b Eq. (4.7) with b(t) = 0.1 when the axisymmetric surface of integration is a prolate or oblate spheroid. For the calculations in this figure, we have used a uniform mesh to describe the surface. We see that even for high aspect ratio prolate or oblate spheroids, the trends are basically the same as those seen for a sphere in the previous subsection. The error in the tangential component of

N = 100

500

N = 40 0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

β (t )

(b)

Fig. 8. The amplification factors for the (a) tangential and (b) normal components of the error given in Eq. (4.7) for the operator J over a sphere, as functions of the parameter b(t). b(n) was zero for these simulations.

the inverted result reaches an asymptotic value as the discretization of the surface is improved, while the error in the normal component continues to increase without bound. 4.3. Well-posedness of the operator J on two spheres Let us now examine the well-posedness of the operator J when the surface S consists of two spheres of unit radius separated by a center to center distance of 2(1 + D), or an end-end separation of 2D. In Fig. 11, we have shown the amplification factors for the error ~ in Eq. (4.7) as a function of the distance D. In all the calculations b shown, the non-uniform mesh described in the previous section with f = 0.075,Ds1 = 1/5120 and D s2 = 1/80 was used, and with the highly discretized regions appearing in the portions of the surfaces that face each other. This choice was necessary to maintain a node-spacing that was smaller than the distance 2D between the spheres, which ensures accurate calculation of the boundary integrals (e.g. see [34] for BI calculations of two drops undergoing a head-on collision). Fig. 11 shows that with the above choice of discretization, the amplification factors for the error are unaffected by the distance D, except when the spheres are extremely close. This behavior is similar to that of the double layer operator, as elucidated by Kim and Karrila [3]. The eigenvalues of the double layer operator for each of the two spheres overlap each other, their difference showing a D7 dependence. It is only when the two spheres are extremely close that the eigenvalues of the two operators show splitting. However, the magnitude of the tangential and normal amplification factors are determined, as expected from the results in Fig. 8, to be of the order of the error based on the smaller element size Ds1 = 1/5120 on the surface.

150

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

10

4

10

3

20 L = 0.5

L =1

18

β (t ) = 0

16

(t )

10

2

L

L = 0.1

β (t ) = 0.1 10

L = 10

L=2

14

10 40

1

0

0.2

0.4

0.6

1

12

0.8

1

80

160

320

640

320

640

Number of points, N

(a)

(a) 3000

β (t ) = 0 (red)

2500

10

β (t ) = 0.1 (black)

4

L = 0.1

2000

L = 0.5

( n ) 1500

10

3

10

2

L=2

1000 500 0

0

0.2

0.4

0.6

0.8

1

L = 10

40

80

A key conclusion from Fig. 11 is that although the amplification factors for the errors are almost independent of the separation D between the spheres even for small separations, this was achieved by maintaining a sufficiently fine mesh spacing in the contact region between the spheres. If the first kind operator is employed (b(t) = b(n) = 0), then the additional discretization required in the contact region as the spheres approach each other will lead to increasingly stronger amplification of any high frequency errors. For the mixed type operator, the tangential amplification factors are less sensitive to mesh refinement, and therefore improving the discretization increases the amplification of errors only weakly. We should mention here that, as before, the trends for the case when tangential coefficient, b(t), of the double layer is zero and normal component, b(n), is non-zero, are the same as the trends above, with the only change that results for the tangential and normal components are reversed.

5. Comparison of numerical schemes for a drop with interfacial slip placed in a uniaxial extensional flow In this section, we investigate the solution procedure for the problem of a spherical drop placed in a uniaxial extensional flow in the presence of interfacial slip. We will use the theoretical developments in the previous sections to predict when an algorithm based on a first kind or mixed type equation will converge, and confirm these predictions with numerical experiments. The

160

Number of points, N

(b) Fig. 9. The amplification factors for the (a) tangential and (b) normal components of the error given in Eq. (4.7) for the operator J over a sphere, when the discretization is improved over a fraction f of the total curvilinear length l of the axisymmetric curve C representing the sphere (see Fig. 6). When f ¼ 0, the curve C is discretized uniformly using 40 points, while for f ¼ 1, C is discretized uniformly using 400 points. For intermediate values of f, C is divided into two parts of length lf and lð1  f Þ. The portion of length lf is discretized uniformly using an element length of l/400, while the one of length lð1  f Þ is discretized uniformly using an element length of l/40. b(n) was set to 0 for these calculations.

L =1

(b) Fig. 10. The amplification factors for the (a) tangential and (b) normal components of the error given in Eq. (4.7) for the operator J when the surface of integration is either a prolate (L = 0.1) or oblate (L = 10) spheroid [see inset of subfigure (b)]. The constants b(n) and b(t) in the operator J were set to 0 and 0.1 respectively in all these calculations.

development of the BI formulation for this problem is detailed elsewhere [21], and we only summarize the equations here. The problem under consideration is depicted in Fig. 12. A drop of dimensionless radius 1 is placed at the stagnation point of a uniaxial extensional flow, u1(x), given in dimensionless form as, with by

u1 ðxÞ ¼ zez 

q 2

eq ;

ð5:1Þ

where ez is the unit vector along the axis of symmetry, and eq is the unit vector in the radial direction. The co-ordinate variables q and z are scaled with the undeformed drop radius a and the velocity is scaled with the characteristic velocity Ga, G being the strain rate of the ambient flow. The drop is initially spherical in shape, and has a viscosity that is k times that of the suspending fluid. The interface is characterized by a capillary number Ca = lGa/c (where l is the suspending fluid viscosity and c is the interfacial tension) and a dimensionless slip coefficient a = a0 l/a, where a0 is the dimensional slip coefficient [21]. The velocity field u(x0) on the surface of the drop, which is used in the kinematic condition to advance the interface shape with time, is determined by the following equation:

uðx0 Þ þ ¼

1 ð1  kÞ 4p ð1 þ kÞ

Z

PV

uðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ S

Z 2 1 1 Gðx; x0 Þ  nðxÞDf ðxÞdSðxÞ u1 ðx0 Þ  1þk 4p ð1 þ kÞ S   Z ka 1 þ ½HðxÞ  fðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ : HðxÞ  fðxÞ þ 4p S ð1 þ kÞ ð5:2Þ

151

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

10

5

Df ðxÞ ¼ β

10

(t )

=0

4

β (t ) = 0.01 2

β (t ) = 0.1 1

10 -3 10

10

-2

10

ð5:3Þ

We have assumed that there no gradients in surface tension along the drop surface. As may be seen in Eq. (5.2), the interfacial velocity can be determined if we have an expression for the tangential stress field at the drop surface. This may be obtained in several different ways, one of which is to rearrange the Ladyzhenskaya Eq. (1.1) for the velocity in the external fluid, evaluated at the interface, in the form

(t ) 103

10

1 rn Ca

-1

10

0

10

1

Z

Gðx; x0 Þ  fðxÞdSðxÞ ¼ 8pu1 ðx0 Þ  4puðx0 Þ þ

(a) 10

akHðxÞ  fðxÞ 1

2D

1

4

10

-3

10

-2

10

-1

10

0

10

1

D

(b)

ak 8p

Z

PV

½HðxÞ  fðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ 2 S Z 1 Gðx; x0 Þ  fðxÞdSðxÞ þ 8p S Z uðx0 Þ 1 k þ Gðx; x0 Þ  nðxÞDf ðxÞdSðxÞ þ ¼k 2 8p S 8p Z PV uðxÞ  Tðx; x0 Þ  nðxÞdSðxÞ: þ

ð5:5Þ

S

Fig. 11. The amplification factors for the (a) tangential and (b) normal components of the error given in Eq. (4.7) for the operator J when the surface of integration is a pair of identical spheres of unit radius separated by a center to center distance of 2(1 + D) [see inset of subfigure (b)]. The red, blue and green curves are the relationships for b(t) = 0, b(t) = 0.01 and b(t) = 0.1, with b(n) taken as 0 for all curves. A non-uniform grid on each sphere was employed in the calculations in this figure. The mesh spacing was 2p/5120 (fine) for 7.5% of the total length of the axisymmetric curve describing each sphere, beginning from the regions of the spheres facing each other. The remainder of each spherical surface was discretized with a mesh spacing of 2p/80 (coarse). This choice of discretization ensured that the mesh spacing was sufficiently smaller than the separation between the spheres. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ez n

t

λ

ð5:4Þ

With the velocity and the interface shape specified, this is an integral equation of the first kind for the stress f. The operator on the left hand side is simply the J+ operator with b(t) and b(n) set to zero. An alternative method of getting f (see [21] for details) is via the mixed-type equation

β (t ) = 0, β (t ) = 0.01, β (t ) = 0.1

10

uðxÞ

 Tðx; x0 Þ  nðxÞdSðxÞ

5

( n)

PV

S

S

D

Z

u∞ = −

ρ 2

e ρ + ze z

z

ρ

Capillary number Ca, Slip coefficient α

eρ

1

Uniaxial extensional flow C Fig. 12. The axisymmetric stretching of an initially spherical drop of unit radius by a uniaxial extensional flow u1. The viscosity of the drop is k times the viscosity of the suspending fluid. The interface is surfactant free, and is characterized by a capillary number Ca and a slip coefficient a. The vector ez is the axis of symmetry, while eq is radial unit vector in this axisymmetric cylindrical co-ordinate system. The normal and tangent vectors to the axisymmetric curve C representing the sphere are n and t respectively.

Here, n(x) is the normal vector to the surface of the drop. Df(x) is the force jump across the interface, given by

Comparison of the operator on the left hand side of the above equation with our general operator J+ in (1.10) indicates that b(n) = 0 and b(t) = ak. If we solve either (5.4) or (5.5) to determine f,we can extract the tangential component H(x)  f(x). For the two methods being compared, the numerical scheme is the same, and is described below: (1) Guess the velocity field u. (2) Given the shape of the drop, compute the magnitude of the normal stress jump Df from Eq. (5.3). (3) Compute the stress field f either using Eq. (5.4) or (5.5). (4) Knowing the tangential component of the stress field, and Df (which depends only on the shape of the drop), compute a new guess u(new) for the velocity field by solving Eq. (5.2) using the Richardson iterative technique. (5) Compute the error d = ku(new)  uk. If the error is smaller than a specified tolerance d0 (taken as 108 in our simulations), go to step 6. Otherwise go to step 3, using u(new) as the guess for the velocity field. (6) Update the shape of the drop by substituting the converged interfacial velocity field u in the kinematic condition (see [20]), update the time step, and jump to step 2. We perform our analysis for a spherical drop in an axisymmetric framework. The curve representing the sphere is discretized into N segments using N + 1 equally spaced points, with each segment described by a cubic splines. The calculation of the stress field f either using Eq. (5.4) or Eq. (5.5) was performed using LU factorization followed by Gaussian elimination for the L and U factors, but this decomposition was required to be done only once every time step. Note that we could have used the Richardson iterative method to calculate the stress field. However, for the values of the slip coefficient and viscosity ratio used here, the Richardson iteration requires small values of the constant k [see Eq. (3.2)], which leads to a large number of iterations, and ultimately consumes more CPU time than the LU decomposition method.

152

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

Therefore, the results presented below are computed only using the LU method for the stress field f. We monitor the convergence of the iterative scheme in steps 2– 5 above by comparing the tangential component of the interfacial stress vector f(x)  t(x) with that of the analytical solution f⁄(x)  t(x) for a spherical drop. The error is defined as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z  2 e¼ ðf  t  f  tÞ dA;

α = 0.03

ε

ð5:6Þ

S

where the integration is carried out over the surface of the sphere. The vector drop f⁄(x)  t(x) can be derived analytically from the results of Ramachandran and Leal [35] to be

15kð19k þ 16Þ



f ðxÞ  tðxÞ ¼ 

α = 0.02 ( Squares )

½38k2 ð5a þ 1Þ þ kð240a þ 89Þ þ 48

sin h cos h

ð5:7Þ

We can derive an expression for the net amplification factor per iteration for the two formulations using the theory presented in Sections 2 and 4. Examination of Eq. (5.2) for u and Eq. (5.4) for f suggests that the amplification factor per iteration, A, is

A¼

ak ðtÞ A Au ; 1þk

ð5:8Þ

where AðtÞ is the tangential amplification factor for the first kind Eq. (5.4), and Au is the amplification factor corresponding to the second kind Eq. (5.2) for u. The largest magnitude of AðtÞ occurs for high frequencies, and is equal to 4m. The eigenvalues of the second kind operator for u on a sphere [3] are

3 ð2m þ 1Þð2m þ 3Þ 3 KBu ¼ 1 þ ð2m  1Þð2m þ 1Þ

KAu ¼ 1 

1k for m ¼ 0; 1; 2; 3; . . . 1þk 1K for m ¼ 1; 2; 3; . . . 1þK

ð5:9Þ

The smallest eigenvalue is either 2k/(1 + k) or 2/(1 + k) depending on whether k is less than or greater than 1. In the limit k ? 0 or k ? 1, the double layer operator in Eq. (5.2) has a zero eigenvalue. However, it is common practice to subject this operator to Wielandt deflation [3,4] to remove these zero eigenvalues. These eigenvalues will thus not affect the error amplification. The remaining eigenvalues of the double layer operator are O(1) on a sphere irrespective of the viscosity ratio k, and hence we can write

Au ¼ c;

ð5:10Þ

where c is an O(1) constant. Thus, the amplification factor per iteration, A, can be, at the most

Amax ¼ 4c

ak m; 1þk

ð5:11Þ

For the scheme to be stable, Amax should be less than 1. The maximum allowable frequency is, therefore,

mmax ¼

1þk : 4cak

α = 0.026

ð5:12Þ

The inverse relationship between the slip coefficient and the smallest allowable mesh spacing is not surprising. Note that the velocity Eq. (5.2) can be solved independent of the stress under no-slip conditions when a or k is zero. Therefore, the larger the slip coefficient, the stronger is the coupling between the stress and velocity equations, and the larger is the smallest mesh spacing that can be used. We have verified this behavior in Fig. 13 by plotting the error for different number of points and different combinations of a and k as a function of the iteration number. From these simulations, we were able to estimate the value of c to be approximately 0.5. Thus, if the number of points in the mesh is kept below (1 + k)/2ak, the scheme with the first kind formulation is stable, at least for the

α = 0.024 ( Dots )

Iteration number Fig. 13. The variation of the error in the numerical solution relative to the analytical solution [see Eq. (5.6)] as a function of the iteration number calculated using the first kind equation for the surface stress for different slip coefficients. All the above curves are for a spherical drop described by 128 equally spaced points. The viscosity of the drop relative to the suspending fluid, k, is 0.19. For the above parameters, the simulations are expected from theory [see Eq. (5.12) with c = 0.5] to be numerically stable only up to a = 0.0245, and the above simulations support this prediction. Note that the simulation for a = 0.02 has been carried out only up to 70 iterations to allow the visibility of the a = 0.024 curve.

starting shape of a sphere. In the work of Ramachandran et al. [21], typical non-zero values of the slip coefficient a were between 0.01 and 0.1, while k was varied from 0.19 to 6.8. For the accuracy required in these simulations, it was desirable to have at least 64 points describing the surface of the shape. Convergence was not obtained even for the first time step (steps 2–5 of the scheme described earlier in this section) for (a, k) combinations of (0.01, 6.8), (0.1, 0.19) and (0.1, 6.8), and even though we were able to implement the simulation for the (0.01, 0.19) case for the first few time steps, the addition of grid points at the curved ends of the stretched drop in subsequent time steps in order to maintain accuracy eventually led to numerical instability. Let us now examine the stability of the numerical scheme with the mixed type Eq. (5.5) employed for obtaining the surface stress vector f. Our results from the previous section suggest that the scheme should be well posed except for small values of b(t) = ak, when we expect the tangential amplification factor to be large, particularly for intermediate spatial frequencies. Eqs. (5.2) and (5.5) suggest that the amplification factor per iteration, A, for this scheme is

A¼

k2 a Atangent Au ; 1þk

ð5:13Þ

where Atangent is now the tangential amplification factor for the second kind equation. The dominant behavior of Atangent can be reasonably captured by the approximation 4m/(2b(t)m + 1) = 4m/(2akm + 1) [see Eq. (2.7)], which satisfies the high frequency behavior for both b(t) = 0 and non-zero b(t). Therefore, the largest possible magnitude of the amplification factor per iteration is

Amax ¼

k2 a 4m k 2akm k¼ ; 1 þ k 2akm þ 1 1 þ k 2akm þ 1

ð5:14Þ

where we have used Au ¼ c ¼ 0:5. Eq. (5.14) suggests that Amax is less than 1 for any discretization, and the scheme should be stable irrespective of the values of a,k or m! This is surprising, because AðtÞ grows linearly with the frequency when the tangential coefficient of the operator, b(t) = ak, is small. The key to understanding this is that the stress and velocity equations are coupled weakly when ak is small, as indicated by the prefactor k2a/(1 + k) in Eq. (5.13). We

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

λ = 6.8

α = 0.5

ε α = 0.1 α = 0.05

α = 0.01

Iteration number

(a)

α = 0.05 λ = 0.19

λ =1

ε λ =3

λ = 6.8

λ = 0.01 Iteration number

(b) Fig. 14. The variation of the error in the numerical solution relative to the analytical solution [see Eq. (5.6)] as a function of the iteration number calculated using the second kind equation for the surface stress for different slip coefficients. All the above curves are for a spherical drop described by 128 equally spaced points. In subfigure (a), k is kept constant at 6.8, and a is varied, while in subfigure (b), a is kept constant at 0.05, and k is varied.

confirm this prediction for a variety of combinations of (a, k) in Fig. 14. It is notable that as a is increased at constant k [subfigure (a)] or vice versa [subfigure (b)], the number of iterations required for convergence increases. It must be remembered that the predictions in Eqs. (5.12) and (5.14) are valid when the shape is a sphere. For more elongated shapes, the spatial frequency in the region of highest density of points must be considered in this analysis. However, the trends of ill-posedness related to the tangential/normal coefficients and mesh spacing are quite generally valid for all shapes. 6. Conclusions We have investigated the properties of a general class of second-kind integral equation operators that can arise in the solution of axisymmetric creeping flow problems. This general second-kind operator is characterized by two constants b(n) and b(t), which represent the weights of the tangential and normal components of the double layer potential in the integral equation. We show that the eigenvalues are restricted to a narrow band as long as b(n) and

153

b(t) are nearly equal and O(1) values, and this makes the operator amenable to rapid solution via Richardson iterations (method of successive substitution). For all other combinations of these parameters, the inversion of the operator is not straightforward. Particularly, when the two constants are highly disparate, the condition number is high, and the large number of Richardson iterations required to bring the error down to the desired tolerance can make this approach for inverting the operator prohibitive. A cursory inspection of the operator eigenspectrum suggests that alternative, direct approaches for inverting the operator for these parameter values suffer from problems of ill-posedness. For example, for the slip operator, which has a zero normal coefficient b(n) and a non-zero tangential coefficient b(t), eigenvalues continue to decrease in magnitude as the minimum allowable wavelength, or equivalently, the mesh spacing in the numerical implementation is decreased. This should lead to significant magnification of even the smallest of high-frequency errors in the forcing function. The important point we elucidate in this paper, is that this observation is misleading, and is, in fact, strictly applicable only to the normal component of the inverted field. Examination of the error in the tangential component reveals that in the limit of high frequencies, the amplification factor asymptotes to a constant value. The caveat, however, is that this asymptote varies as 1/b(t), i.e., it increases as b(t) is reduced. Thus, the inversion problem, depending on the mesh and the accuracy of the method employed, can eventually become ill-posed even for the tangential component if the tangential coefficient of the operator, b(t), is small enough. The manifestation of the ill-posedness for small b(t) is more pronounced for problems requiring fine mesh spacing for maintaining accuracy, such as in the thin film regions formed between colliding particles. But as long as b(t) is O(1) or higher, obtaining the tangential component of the solution is a well-posed problem even if the normal coefficient of the operator b(n) is zero. A similar conclusion can be stated for the normal component of the solution when b(t) is set to zero and b(n) is non-zero. This is a significant conclusion for creeping flow problems involving slip, since it allows the use of the boundary integral equation arising naturally in terms of the physical variables, instead of requiring the exploration of alternative formulations with fictitious variables. Acknowledgments We acknowledge funding from the Department of Energy (DOE Award No. DE-FG02-05ER25704, OR# 20080220). We are also grateful to Prof. Hector Ceniceros for helpful discussions. Appendix A. Non self-adjointedness of J± and its nullspace It is notable that the operator J± is not self-adjoint; its adjoint, J T , is

J T qðx0 Þ ¼ bðnÞ nðx0 Þ Z PV 

qðx0 Þ 1

Tðx; x0 Þ  qðxÞdSðxÞ  nðx0 Þ nðx0 Þ  2 8p S Z PV 

qðx0 Þ 1

Tðx; x0 Þ  qðxÞdSðxÞ  nðx0 Þ þ bðtÞ Hðx0 Þ  2 8p S Z 1 Gðx;x0 Þ  qðxÞdSðxÞ ð6:1Þ þ 8p S One can see that the normal vector n(x0) is an eigenfunction of the operator JT irrespective of the parameters b(n) and b(t) of the double layer term; the corresponding eigenvalue is K+ = 0 or K = b(n). Thus, the operator J+ and the operator J with b(n) = 0 have at least one null eigenfunction. In fact, when b(n) = 0, which is consistent

154

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

with operators found in creeping flow problems with slip at interfaces, n(x0) is a null eigenfunction for both operators J± and JT qðx0 Þ. For these special cases, the solution of problem (1.9) is no longer unique. The reason for this null eigenvector is that the single layer integral is invariant to the addition of cn(x0) to the single layer density, which arises from satisfaction of the continuity equation by the free space Green’s function G(x,x0). Physically, this corresponds to the fact that pressure in an incompressible fluid is arbitrary up to a constant, and the addition or subtraction of a constant pressure throughout the fluid does not alter the velocity distribution. The presence of a zero eigenvalue necessitates the following: 1. In order for Eq. (1.10) to have a solution, the RHS of this equation must be orthogonal to the null eigenfunction of the adjoint of the operator J, i.e. the integral

Z

um ðx0 Þ 1 ¼ 2 8p

bðtÞ Hðx0 Þ  qðx0 Þ bðtÞ  2 8p

Z

PV

½HðxÞ  qðxÞ Z 1 ½Gðx; x0 Þ  Tðx; x0 Þ  nðxÞdSðxÞ þ 8p S S

þ cd nðx0 ÞnðxÞ  qðxÞdSðxÞ;

ð6:2Þ

Appendix B. Calculation of the eigenvalues and eigenfunctions of J± In this appendix, we describe the procedure for obtaining the eigenfunctions and eigenvalues of the operator J± defined in Eq. (1.10), on a sphere using a procedure identical to that described in Kim and Karrila. The specific eigenvalue problem that we intend to solve is

bðnÞ ½/ðx0 Þ  nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ  /ðx0 Þ 1  2 8p Z

Tðx; x0 Þ  fbðnÞ ½/ðx0 Þ  nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ Z 1 Gðx; x0 Þ  /ðxÞdSðxÞ  /ðx0 Þg  nðxÞdSðxÞ þ 8p S



S

¼ K/ðx0 Þ

Z

Tðx; x0 Þ

S

ð7:2Þ

^ m ðx0 Þ 1 u ¼ 2 8p

Z S

1 Gðx; x0 Þ  ^f m ðxÞdSðxÞ  8p

Z

^ m ðxÞ Tðx; x0 Þ  u

S

 nðxÞdSðxÞ:

ð7:3Þ

In the above equations, the point x0 is on the surface of the sphere S. Adding Eqs. (7.2) and (7.3), we get



Z ^ m ðx0 Þ um ðx0 Þ þ u 1 ^ m ðxÞ  nðxÞdSðxÞ þ Tðx; x0 Þ  ½um ðxÞ  u 2 8p S Z 1 þ Gðx; x0 Þ  ½^f m ðxÞ  f m ðxÞdSðxÞ ¼ 0 ð7:4Þ 8p S

The eigenfunction Eq. (7.1) may be rearranged slightly to yield

bðnÞ ½/ðx0 Þ  nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ  /ðx0 Þ  K/ðx0 Þ þ 2 Z n o 1 ðnÞ  Tðx; x0 Þ  b ½/ðx0 Þ  nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ  /ðx0 Þ 8p S Z 1  nðxÞdSðxÞ þ Gðx; x0 Þ  /ðxÞdSðxÞ ¼ 0 ð7:5Þ 8p S We can compare the above equation term by term with Eq. (7.4), which yields

where cd is an arbitrary, non-zero scalar constant, usually chosen as an O(1) number. With the above implementation, deflation augments only the zero eigenvalue of the operator JS± to a non-zero value, cdA, where A is the area of the surface S. Otherwise, it retains all the other eigenvalues and all the eigenvectors of JS±. Therefore, the result of solving equation with the operator Jd substituted for JS± is that the resulting normal component of the inverted result q(x0) is affected, but not the tangential components. This is usually acceptable, since it is generally the tangential components of q(x0) that are required in numerical schemes involving slip (e.g. [20–22]).

J  /ðx0 Þ ¼

S

1 8p

while the growing harmonics may be substituted in the boundary integral equation for the internal velocity field to yield

nðx0 Þ  bðx0 ÞdSðx0 Þ

J d qðx0 Þ ¼

Gðx; x0 Þ  f m ðxÞdSðxÞ þ

 um ðxÞ  nðxÞdSðxÞ;

S

must be zero. This is ordinarily true, because the vector b, in most cases, is a velocity vector which obeys the continuity equation over the volume enclosed by the surface S. 2. To obtain a unique solution for Eq. (1.10), we need to deflate the operator JS± in the numerical scheme. The deflated operator Jd is

Z

ð7:1Þ

The general solution to the creeping flow equations, which consists of growing and decaying harmonics (e.g. see [3]) was first presented by Lamb. Using the BI equations [1,3], one can establish a relationship between Lamb’s velocity and stress functions for both growing and decaying harmonics. The decaying harmonics may be substituted in the BI equation for the external velocity field [see Eq. (1.1)] to obtain

^ m ðx0 Þ ¼ ½2KI  bðnÞ nðx0 Þnðx0 Þ  bðtÞ Hðx0 Þ  /ðx0 Þ um ðx0 Þ þ u ^ m ðx0 Þ ¼ ½bðnÞ nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ  /ðx0 Þ um ðx0 Þ  u ^f ðx Þ  f ðx Þ ¼ /ðx Þ m

0

m

0

ð7:6Þ

0

Thus, /(x0) can be an eigenfunction of the operator J only when the following equations are satisfied:

^ m ðx0 Þ ¼ ½2KI  bðnÞ nðx0 Þnðx0 Þ  bðtÞ Hðx0 Þ um ðx0 Þ þ u  ½^f m ðx0 Þ  f m ðx0 Þ and

ð7:7Þ

^ m ðx0 Þ ¼ ½bðnÞ nðx0 Þnðx0 Þ þ bðtÞ Hðx0 Þ um ðx0 Þ  u  ½^f m ðx0 Þ  f m ðx0 Þ Let us consider axisymmetric modes only, as this is the focus of this ^ m , and the stress fields paper. In this case, the velocity fields um and u fm(x0) and ^f m ðx0 Þ may be written as [3]

 um ¼ ir P m ðcos hÞ ðm þ 1ÞBm þ

mþ1 Am 2ð2m  1Þ





dP m ðcos hÞ ðm  2Þ Bm  Am dh 2mð2m  1Þ   m ^ m ¼ ir P m ðcos hÞ mbn þ u an 2ð2m þ 3Þ

dP m ðcos hÞ ðm þ 3Þ bm þ þ ih dh 2ðm þ 1Þð2m þ 3Þ   m2 þ 3v  1 f m ¼ ir Pm ðcos hÞ 2ðm þ 1Þðm þ 2ÞBn  An 2m  1

dP m ðcos hÞ ðm þ 1Þðm  1Þ 2ðm þ 2ÞBm þ Am þ ih dh mð2m  1Þ   2 m  m  3 ^f m ¼ ir Pm ðcos hÞ 2mðm  1Þbm þ am 2m þ 3

dP m ðcos hÞ mðm þ 2Þ 2ðm  1Þbm þ þ ih dh ðm þ 1Þð2m þ 3Þ þ ih

ð7:8Þ

155

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

The problem of determining the eigenspectrum then reduces to the following homogeneous simultaneous equations [radial and tangential components of both equations in (7.7)] in an, bn, An and Bn:

Xy ¼0

ð7:9Þ

where the vector y is

2

am

3

6b 7 6 m7 y¼6 7 4 Am 5

ð7:10Þ

Bm

The last term on the RHS of the above equation can be regarded as a velocity resulting from a single layer density q, and we will henceforth denote this term by the operator Vq(x0). Therefore,

Vqðx0 Þ ¼ ðnÞ

and the matrix X, as an example, for the slip operator JSþ ðb

2 6 6 6 X¼6 6 6 4

Z PV 

qðx0 Þ 1  J Tþ qðx0 Þ ¼ bðnÞ nðx0 Þ  Tðx; x0 Þ  qðxÞdSðxÞ  nðx0 Þ 2 8p S Z PV  qðx 1 0Þ  Tðx; x0 Þ  qðxÞdSðxÞ nðx0 Þ þ bðtÞ Hðx0 Þ  2 8p S Z 1 Gðx; x0 Þ  qðxÞdSðxÞ ð8:2Þ  nðx0 Þg þ 8p S

¼ 0Þ, is

1 8p

m

ðmþ1Þ 2ð2m1Þ

ðmþ3Þ mðmþ2Þ ðtÞ  2ðmþ1Þð2 mþ3Þ  b ðmþ1Þð2mþ3Þ

1  2bðtÞ ðm  1Þ

m2Þ ðtÞ ðmþ1Þðm1Þ  2mðð2 m1Þ þ b mð2m1Þ

2

m

2ð2mþ3Þ ðmþ3Þ 2ðmþ1Þð2mþ3Þ

m3Þ  2K ðmð2mþ3Þ

ðmþ1Þ 2ð2m1Þ

m  4Kmðm  1Þ

dPm ðcos hÞ dh

2

m1Þ  2K ðmð2þ3 m1Þ Þ

Wqðx0 Þ ¼ 

where

ð7:15Þ

dPm ðcos hÞ=dh maxh ½dPm ðcos hÞ=dh

ð7:16Þ

 Tðx; x0 Þ  qðxÞdSðxÞ  nðx0 Þ;

ð8:4Þ

S



qðx0 Þ þ Wqðx0 Þ þ Vqðx0 Þ ¼ kqðx0 Þ: 2

" # qðx0 Þ 1 ðbðnÞ  bðtÞ Þ nðx Þnðx Þ þ Wqðx0 Þ ¼ ðtÞ I  0 0 2 b bðnÞ  ½K qðx0 Þ  Vqðx0 Þ ¼

ð8:1Þ

1 ðtÞ

b

½K qðx0 Þ  Vqðx0 Þ 

ðbðnÞ  bðtÞ Þ bðtÞ bðnÞ

½K qðx0 Þ  nðx0 Þ  nðx0 Þ  Vqðx0 Þnðx0 Þ:

ð8:6Þ

The surface force resulting from the single layer potential Vq is discontinuous across the surface. For a point on the surface S approaching from inside its enclosed volume (see, for e.g., [3,4]), the surface force gI is

qðx0 Þ þ Wqðx0 Þ: 2

ð8:7Þ

The expression for the viscous energy dissipation arising from the flow field Vq(x0) and its corresponding stress field gI(x0) in the volume V enclosed by the integration surface S can be written as

Z

V

In this appendix, we prove that the eigenvalues and the eigenvectors of the operator J+ are real on a sphere. We choose to prove the equivalent result that the eigenvalues of the transpose of operator J T on a sphere are real. The eigenvalue problem of interest is

where K± and q(x0) have, in general, imaginary parts. We recall that the operator J T is

PV

ð8:5Þ

EI ¼ 2l

Appendix C

J Tþ qðx0 Þ ¼ K qðx0 Þ

Z

½bðtÞ I þ ðbðnÞ  bðtÞ Þnðx0 Þnðx0 Þ 

g I ðx0 Þ ¼

^ is The normalized eigenvector /

b m Pm ðcos hÞ þ ih D bm ^ m ¼ ir C /

1 8p

we can write Eq. (8.1) in terms of Eqs. (8.2)–(8.4) as

ð7:14Þ

h i ¼ max C A;B r m P m ðcos hÞ h   dPm ðcos hÞ ¼ max C A;B hm h dh

ð7:11Þ

where V is the single layer operator. If we define W as the adjoint of the double layer operator, i.e.

The normalization factor /0 is chosen such that neither the tangen^m tial nor the radial component of the normalized eigenfunction / can exceed a magnitude of 1 for h 2 [0,p].

/h0m

7 7 7 7 7 ðm þ 1Þ þ 4Kðm þ 1Þðm þ 2Þ 7 5 1  2ð2K  bðtÞ Þðm þ 2Þ 1  2bðtÞ ðm þ 2Þ

The above equation can be inverted to yield

ð7:13Þ

n o ¼ max /r0m ; /h0m

3

ð7:12Þ

The eigenfunction derived above can be normalized as

^ m ¼ /m / /0m

ð8:3Þ

1  m

mðmþ2Þ ðm2Þ ðtÞ ðtÞ ðmþ1Þðm1Þ  ð2K  bðtÞ Þ ðmþ1Þð2 mþ3Þ 1  2ð2K  b Þðm  1Þ  2mð2m1Þ þ ð2K  b Þ mð2m1Þ

/m ¼ ir Pm ðcos hÞC rm þ ih C hm

/r0m

Gðx; x0 Þ  qðxÞdSðxÞ; S

 2ð2mmþ3Þ

The equation system in Eq. (7.9) has a non-trivial solution only when the determinant of X is zero. This yields a quadratic equation in K, and therefore, two roots, KAm and KBm . Analytic solutions for KAm and KBm are possible, but the equations are too long to be presented here. The corresponding null eigenvectors of X give the non-trivial solutions for am,bm,Am and Bm, which can be substituted in Eq. (7.8) to get fm and ^f m , and yield the eigenfunction /m(x0) from Eq. (7.6) as

/0m

Z

e : e dV ¼

Z

½Vqðx0 Þ  g I ðx0 ÞdS:

ð8:8Þ

S

The superscripted asterix denotes the complex conjugate of the preceding variable or expression. Since the energy of dissipation is a real number and, in fact, never a negative quantity, EI P 0. We can write, therefore, using Eqs. (8.6)–(8.8) 2 1 3   Z ðtÞ ½Kþ qðx0 Þ  Vq ðx0 Þ  Vqðx0 Þ  Vq ðx0 Þ 4 b ðnÞ ðtÞ 5dSðx0 Þ ¼ a2 : Þ S ½Kþ qðx0 Þ  nðx0 Þ  nðx0 Þ  Vqðx0 Þnðx0 Þ  Vq ðx0 Þ  ðbbðtÞb bðnÞ ð8:9Þ

156

A. Ramachandran et al. / Computers & Fluids 64 (2012) 141–156

The above equation can be reformulated as ðnÞ

1 bðtÞ

½ Kþ I 1  I 2  

ðb

ðtÞ

b Þ

bðtÞ bðnÞ

½Kþ I3  I4  ¼ a2 ;

ð8:10Þ

where the integrals I1,I2,I3 and I4 are

I1 ¼

Z

qðx0 Þ  Vq ðx0 ÞdS;

ð8:11Þ

Vqðx0 Þ  Vq ðx0 ÞdS;

ð8:12Þ

½qðx0 Þ  nðx0 Þ½nðx0 Þ  Vq ðx0 ÞdS;

ð8:13Þ

½nðx0 Þ  Vqðx0 Þ½nðx0 Þ  Vq ðx0 ÞdS:

ð8:14Þ

S

I2 ¼ I3 ¼

Z

ZS S

and

I4 ¼

Z S

Eq. (8.10) suggests that the eigenvalue K± is a real quantity, provided each of the integrals I1,I2,I3 and I4 are real. Let us examine each of these integrals separately. The integral I1 is a positive number, because the operator Vq is a self-adjoint operator. The integrands in I2 and I4 are products of complex conjugates, and hence, these integrands and the corresponding integrals are real. Let us now consider the integral I3. To determine whether this integral is real or complex, let us first expand the function q in terms of the eigenfunctions w(x0) of the self-adjoint operator V.

qðx0 Þ ¼

1 X

cm wm ðx0 Þ:

ð8:15Þ

m¼1

The integrand of I3 then becomes

½Vq ðx0 Þ  nðx0 Þ½qðx0 Þ  nðx0 Þ 1 X 1 X ¼ cm cp mm ½wm ðx0 Þ  nðx0 Þ½wp ðx0 Þ  nðx0 Þ:

ð8:16Þ

m¼1 p¼1

Here mm is the eigenvalue of the operator V. The integral I3 becomes, therefore,

I3 ¼

Z 1 X 1 X cm cp mm ½wm ðx0 Þ  nðx0 Þ½wp ðx0 Þ  nðx0 ÞdS:

ð8:17Þ

S

m¼1 p

On a sphere, wm(x0)  n(x0) = Pm (cosh) [see Eq. (7.12)], and the integral in Eq. (8.17) becomes

Z

½wm ðx0 Þ  nðx0 Þ½wp ðx0 Þ  nðx0 ÞdS

S

Z

¼ 2p

h¼p

Pm ðcos hÞP p ðcos hÞ sin hdh ¼

h¼0

2p dmp : 2m þ 1

ð8:18Þ

This reduces I3 to

I 3 ¼ 2p

1 X c  c m mm

m¼1

m

2m þ 1

;

ð8:19Þ

which is clearly real. Since the four integrals I1 through I4 are real, the eigenvalues on a sphere, according to Eq. (8.10) must be real. The eigenvalues are also always real for the Hebeker operator [b(n) = b(t)] on any geometry, again from Eq. (8.10). References [1] Leal LG. Advanced transport phenomena: fluid mechanics and convective transport processes. New York: Cambridge University Press; 2007.

[2] Ladyzhenskaya O. The mathematical theory of viscous incompressible flows. New York: Gordon and Breach; 1963. [3] Kim S, Karrila SJ. Microhydrodynamics: principles and selected applications. New York: Dover; 2005. [4] Pozrikidis CJ. Boundary integral and singularity methods for linearized viscous flow. New York: Cambridge University Press; 1992. [5] Youngren GK, Acrivos A. Stokes flow past a particle of arbitrary shape: a numerical method of solution. J Fluid Mech 1975;69(2):377–403. [6] Ingber MI. Numerical simulation of the hydrodynamic interaction between a sedimenting particle and a neutrally buoyant particle. Int J Numer Methods Fluids 1989;9:263–73. [7] Chi BK. The motion of immisicible drops and the stability of annular flow. Pasadena: California Institute of Technology; 1986. [8] Martinez MJ, Udell KS. Axisymmetric creeping motion of drops through circular tubes. J Fluid Mech 1990;210:565–91. [9] Tsai TM, Miksis MJ. Dynamics of a drop in a constricted capillary tube. J Fluid Mech 1994;274:197–217. [10] Maurin AL, Barthes-Biesel D. Motion of a deformable capsule through a hyperbolic constriction. J Fluid Mech 1994;279:135–63. [11] Queguinier C, Barthes-Biesel D. Axisymmetric motion of capsules through cylindrical channels. J Fluid Mech 1997;348:349–76. [12] Ingber MI, Womble DE, Mondy LA. Simulations of hydrodynamics among immersed particles in Stokes flow using a massively parallel computer. Bound Elem 1995;17:515–22. [13] Powers T, Lauga E. The hydrodynamics of swimming microorganisms. Rep Prog Phys 2009;72:096601. [14] Hebeker F. Efficient boundary element methods for three-dimensional exterior viscous flow. Numer Methods Partial Differ Eq 1986;2(4):273–97. [15] Rallison JM, Acrivos A. A numerical study of the deformation and burst of a viscous drop in an extensional flow. J Fluid Mech 1978;89(1):191–200. [16] Power H. On the Rallison and Acrivos solution for the deformation and burst of a viscous drop in an extensional flow. J Fluid Mech 1987;185:547–50. [17] Ratcliffe T, Zinchenko AZ, Davis RH. Buoyancy-induced squeezing of a deformable drop through an axisymmetric ring constriction. Phys Fluids 2010;22(8):082101. [18] Zinchenko AZ, Davis RH. Squeezing of a periodic emulsion through a cubic lattice of spheres. Phys Fluids 2008;20:04080. [19] Zinchenko AZ, Davis RH. Algorithm for direct numerical simulation of emulsion flow through a granular material. J Comput Phys 2008;227: 7841–88. [20] Zinchenko AZ, Davis RH. A boundary integral study of a drop squeezing through interparticle constrictions. J Fluid Mech 2006;564:227–66. [21] Ramachandran A et al. The effect of interfacial slip on the dynamics of a drop in flow: Part I. Stretching, relaxation and breakup. J Rheol 2012;56:45–97. [22] Luo H, Pozrikidis CJ. Interception of two spheres with slip surfaces in linear Stokes flow. J Fluid Mech 2007;581:129–56. [23] Luo H, Pozrikidis CJ. Effect of surface slip on Stokes flow past a spherical particle in infinite fluid and near a plane wall. J Eng Math 2009;62: 1–21. [24] Ajdari A. Slippage at a polymer/polymer interface: entanglements and associated friction. CR Acad Sci Ser II: Mec Phys Chim Sci Terre Univers 1993;317:1159–63. [25] Zhao R, Macosko CW. Slip at polymer-polymer interfaces: rheological measurements on coextruded multilayers. J Rheol 2002;46:145–67. [26] Park CC, Baldessari F, Leal LG. Study of molecular weight effects on coalescence: interface slip layer. J Rheol 2003;47:911–42. [27] Hsu AS, Roy A, Leal LG. Drop-size effects on coalescence of two equal-sized drops in a head-on collision. J Rheol 2008;52:1291–310. [28] Lee PC, Park HE, Morse DC, Macosko CW. Polymer–polymer interfacial slip in multilayered films. J Rheol 2009;53:893–915. [29] Goveas JL, Fredrickson GH. Apparent slip at a polymer–polymer interface. Eur Phys J B 1998;2:79–92. [30] Adhikari NP, Goveas JL. Effects of slip on the viscosity of polymer melts. J Polym Sci B 2004;42(10):1888–904. [31] Schawlabe JT, Vlahovska PM, Miksis MJ. Monolayer slip effects on the dynamics of a lipid bilayer vesicle in a viscous flow. J Fluid Mech 2010;647:403. [32] Richardson LF. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philos Trans Roy Soc Lond Ser A 1911;210: 307–57. [33] Chapra SC, Canale RP. Numerical methods for engineers. New Delhi: McGrawHill; 2000. [34] Yoon Y et al. Coalescence of two equal-sized deformable drops in an axisymmetric flow. Phys Fluids 2007;19:102102. [35] Ramachandran A, Leal LG. The effect of interfacial slip on the rheology of a dilute emulsion of drops for small capillary numbers. J Rheol; submitted for publication.