Resistance Functions for Two Spheres in Axisymmetric Flow—Part I: Stream Function Theory

Abstract

We consider low-Reynolds-number axisymmetric flow about two spheres using a novel, biharmonic stream function. This enables us to calculate analytically not only the forces, but also the dipole moments (stresslets and pressure moments) and the associated resistance functions. In this paper the basics properties of axisymmetric flow and the stream function are discussed. Explicit series expansions, obtained by separation in bispherical coordinates, will be presented in a follow-up paper.

1. Introduction

The grand resistance and mobility tensors describe the hydrodynamic interaction between rigid bodies suspended in a fluid medium and play an all-important role in colloidal science [1–7]. More specifically, they express the linear relationship between the Cartesian force multipole moments exerted by the particles on the fluid and the gradients of the ambient flow velocity taken at the particle centers. Both tensors depend, in general, on the positions and orientations of all the suspended particles. However, in the special case of just two spherical bodies, owing to O(2)-symmetry about the line connecting the particle centers, the full tensors can be reduced to a set of scalar resistance and mobility functions [8–15]. Accurate knowledge of the two-body resistance functions is essential to overcome certain contact singularities—also referred to as lubrication singularities—that dominate the many-body tensors when a pair of particles comes close to touching [16, 17].

The calculation of the 2-body resistance functions is based on the solution of the Stokes boundary value problem for stationary, low-Reynolds number flow about two spheres. In this case there exists a set of curvilinear coordinates that is adapted to the physical boundaries at hand—the bispherical coordinates—and it seems natural to try to solve the problem analytically by separation in these coordinates [18–21]. However, separation of the Stokes equations in bispherical coordinates succeeds only up to a modulation factor [22], which complicates matching of the general solution to the boundary data. In most cases, it turns out that the matching problem can be reduced to a linear recursion problem that involves an infinite tridiagonal matrix. Although this cannot, in general, be solved analytically, the numerical treatment of such a tridiagonal recursion scheme [23–28] is much more efficient than recursions based on the twin-multipole expansion [9, 12–14].

Further progress can be made if the flow is axisymmetric, which happens when the spheres translate along or rotate around the line connecting their centers and when the ambient flow is axially symmetric about that line. A subset of resistance functions, including those with the strongest contact singularities, is determined by this kind of flow problems [13, 15]. It follows from the Stokes equations that for axisymmetric flow, the azimuthal velocity component is decoupled from the remaining flow variables, namely, the pressure and the radial and axial velocity components. Accordingly, there are two types of axisymmetric flow problems, problem I involving only the azimuthal velocity and problem II involving the other flow variables. Problem I can be reduced to a harmonic equation for the azimuthal velocity with Dirichlet boundary conditions. This was solved by separation in bispherical coordinates almost hundred years ago by Jeffery [29], who calculated, in particular, the torques acting on two spheres that rotate with given angular velocities around the line of centers.

Problem II is more difficult, since it involves coupled equations among three variables. However, owing to the incompressibility, the two nonvanishing velocity components, that is, radial and axial velocities, can be expressed as derivatives of a scalar potential, the so-called stream function, whose introduction can be traced back to Stokes [30]. This stream function satisfies a fourth-order differential equation that can again be separated in bispherical coordinates. In this way, Stimson and Jeffrey calculated the forces on two spheres moving with equal velocities along the line of centers [31]. Many years later, Brenner [32] and Maude [33] solved the same problem for spheres moving with opposite velocities along the line of centers. The latter case, where two spheres approach each other, is the one with the strongest contact singularity. It is possible to extract the asymptotic behavior of the forces at close proximity from the series expansion in bispherical coordinates [34, 35], and the results agree with those obtained from a singular perturbation theory in stretched cylindrical coordinates [36, 37]. Singular perturbation theory, valid for two spheres near touching, has also been used to calculate the forces and their dipole moments in the presence of an external, axisymmetric rate-of-strain flow [14, 38, 39]. Moreover, the twin-multipole expansions, valid for large distances, are well-known for these cases [12–14]. Yet, the corresponding series expansions in bispherical coordinates, valid for all distances, have not been published so far.

In this and a follow-up paper [40], we generalize and complete previous studies of axisymmetric two-sphere-flow problems in several ways. First, we suggest a slight conceptual improvement by considering from the beginning an inhomogeneous external flow, as generated, for example, by means of a localized external volume-force density. This allows us to identify easily the contributions from each sphere without being obliged to assume artificial deformations, as in [13, 38, 41]. Second, and more important, we introduce a novel stream function that differs from the classical one [30–33, 37] but has the advantage of being a biharmonic function. Its series expansion in bispherical coordinates matches the corresponding expansions of the flow fields better than the old stream function. In particular, the biharmonic stream function is expressed in associated Legendre functions instead of Gegenbauer polynomials. Third, with the new stream function, we are able to calculate analytically not only the forces, but also their dipole moments (stresslets and pressure moments). And we will show that this can be accomplished without full knowledge of the pressure (which would require solving a tridiagonal recursion scheme). This signifies some progress, as it was claimed earlier in the literature [38] that the stresslets cannot be calculated on the basis of the stream function alone. In this paper, we present the theory of the biharmonic stream function for axisymmetric flow, while in the follow-up paper, we derive the series expansions for the forces and dipole moments in bispherical coordinates.

2. Stokes Equations

We consider two spheres immersed in an unbounded, simple fluid with shear viscosity 𝜂. The spheres—in the following labeled by Greek indices 𝜈,𝜇∈{1,2}—have radii 𝑎_𝜈 and are centered on the 3-axis of a Cartesian frame at positions 𝐗_𝜈=𝑍_𝜈𝐞₃. To exclude overlap, we assume that |𝑍₁−𝑍₂|>𝑎₁+𝑎₂. The dynamics of the fluid is described in terms of the local deviation from the thermal pressure, 𝑝(𝐱), and the flow velocity 𝐮(𝐱) which obey the linear Stokes equations for stationary, low-Reynolds number flow with stick boundary conditions [8, 42] ∇⋅𝐮=0∇𝑝−𝜂Δ𝐮=𝐟^ext𝐱∈𝙱₀,𝑝,𝐮⟶0(‖𝐱‖⟶∞),𝐮=𝐔_𝜈+𝛀_𝜈×𝐱−𝐗_𝜈𝐱∈𝜕𝙱_𝜈,𝜈=1,2. (2.1) Here, 𝙱₀ denotes the fluid region and 𝜕𝙱_𝜈 the surface of particle 𝜈. Moreover, 𝐟^ext is an external force density acting on the fluid, while 𝐔_𝜈 and Ω_𝜈 are the translational and rotational velocities of particle 𝜈, respectively. We assume 𝐟^ext to have a compact support that is not overlapping with the particles. Equations (2.1) then pose a Dirichlet boundary value problem for (𝑝,𝐮) whose solution—for fixed geometry and given sources (𝐟^ext,𝐔,Ω)—is uniquely determined [8, 43].

We are interested here in the forces and dipole moments exerted by the particles on the fluid. These are defined by (we use Latin indices 𝑖,𝑗,…∈{1,2,3} to label Cartesian components and the summation convention over repeated upper and lower Latin indices)ℱ^𝜇_𝑖=_𝜕𝙱𝜇𝑑𝑆𝐹^𝜇_𝑖(𝐱),ℱ^𝜇_𝑖𝑗=𝑎_{𝜇𝜕𝙱𝜇}𝑑𝑆𝐹^𝜇_𝑖(𝐱)𝑁_𝑗(𝐱),𝜇=1,2, (2.2) where the surface-force density [8, 42]𝐅^𝜇=𝐍⋅𝑝𝟏−2𝜂(∇𝐮)_{𝜕𝙱+𝜇}𝐱∈𝜕𝙱_𝜇 (2.3) is determined by the solution of (2.1). Here, 𝐍 is the normal field on 𝜕𝙱_𝜇 (directed outwards), (∇𝐮) is the symmetrized velocity gradient, and the subscript 𝜕𝙱⁺_𝜇 indicates analytic continuation from the fluid regime onto the particle surface. It is useful to split the dipole moments into a trace, a skew-symmetric, and a symmetric-traceless part, according toℱ^𝜇_𝑖𝑗=𝑄^𝜇𝛿_𝑖𝑗−12𝜀^𝑚_𝑖𝑗𝒯^𝜇_𝑚+𝒮^𝜇_𝑖𝑗, (2.4) where (for an arbitrary second-rank tensor with components 𝐵_𝑖𝑗, we use the notations 𝐵_(𝑖𝑗)=(1/2)(𝐵_𝑖𝑗+𝐵_𝑗𝑖),𝐵_[𝑖𝑗]=(1/2)(𝐵_𝑖𝑗−𝐵_𝑗𝑖),𝐵_{⟨𝑖𝑗⟩}=𝐵_(𝑖𝑗)−(1/3)𝛿^𝑘𝑙𝐵_𝑘𝑙𝛿_𝑖𝑗)𝑄^𝜇=13𝛿^𝑚𝑛ℱ^𝜇_𝑚𝑛𝒯^𝜇_𝑖=−𝜀_𝑖^𝑚𝑛ℱ^𝜇_[]𝑚𝑛𝒮^𝜇_𝑖𝑗=ℱ^𝜇_𝑖𝑗 (2.5) denote the pressure moment [14], the torque, and the (deviatoric) stresslet [11, 44], respectively.

We decompose the flow (𝑝,𝐮) in the form𝑝(𝐱)=𝑝_ext(𝐱)(𝐱),𝐮(𝐱)+𝑞=𝐮_ext(𝐱)+𝐯(𝐱),𝐱∈𝙱₀, (2.6) where (𝑝_ext,𝐮_ext) is the flow caused by 𝐟^ext in the absence of the particles and called the external flow. The remainder (𝑞,𝐯) will be referred to as the excess flow. It satisfies the Dirichlet problem∇⋅𝐯=0,∇𝑞−𝜂Δ𝐯=0,𝐱∈𝙱₀,𝑞,𝐯⟶0,(‖𝐱‖⟶∞),𝐯=𝐮_𝜈,𝐱∈𝜕𝙱_𝜈,𝜈=1,2, (2.7) with applied velocities𝐮_𝜈(𝐱)=𝐔_𝜈+𝛀_𝜈×𝐱−𝐗_𝜈−𝐮_ext(𝐱),𝐱∈𝜕𝙱_𝜈,𝜈=1,2. (2.8) Inserting (2.6) in (2.3), one obtains a decomposition of the surface-force densities of the form𝐅^𝜇=𝐅^𝜇_ext+𝐆^𝜇,𝐆^𝜇=𝐍⋅𝑞𝟏−2𝜂(∇𝐯)_{𝜕𝙱+𝜇}, (2.9) where 𝐅^𝜇_ext is defined as in (2.3), with (𝑝,𝐮) replaced by (𝑝_ext,𝐮_ext). Notice that the external flow is well defined inside the particles, where it obeys the homogeneous Stokes equations (since the support of 𝐟^ext lies outside the particles). The contribution from 𝐅^𝜇_ext to the force moments (2.2) can therefore be evaluated by converting the surface integrals into volume integrals over the interiors 𝙱_𝜇 of the particles. Expanding then (𝑝_ext,𝐮_ext) about the particle centers, 𝐗_𝜇, and using isotropy, it is easy to calculate these integrals exactly.

Assuming that the external flow varies slowly on the length scale of the particle diameters, we expand the applied velocities (2.8) up to first order about the particle centers. This yields𝑢^𝑘_𝜈(𝐱)=𝒰^𝑘_𝜈+𝑎_𝜈−𝜀_𝑛^𝑘𝑙𝒲^𝑛_𝜈+ℰ_𝜈^𝑘𝑙𝑁_𝑙(𝐱)+⋯𝐱∈𝜕𝙱_𝜈,𝜈=1,2, (2.10) where the expansion coefficients𝒰^𝑘_𝜈=𝑈^𝑘_𝜈−𝑢^𝑘_ext𝐗𝜈𝒲^𝑘_𝜈=Ω^𝑘_𝜈−12𝜀^𝑘_𝑚𝑛𝜕^[𝑚𝑢^𝑛]_ext𝐗𝜈ℰ_𝜈^𝑘𝑙=−𝜕^(𝑙𝑢^𝑘)_ext𝐗𝜈 (2.11) are referred to as the (effective) translational and rotational velocities and the local rate of strain, respectively. We also put𝒫_𝜈=−𝑝_ext𝐗𝜈. (2.12) It follows from incompressibility (cf. the first equation in (2.1)) that the local rates of strain are traceless, 𝛿_𝑚𝑛ℰ_𝜈^𝑚𝑛=0. Notice that by admitting inhomogeneous external flows, we can avoid conceptual problems of earlier approaches which were based on a uniform shear flow and which, therefore, either could not always resolve the individual contributions of the spheres [11], or, in order to be able to do so, had to assume artificial deformations [13, 38, 41].

Upon inserting (2.9) in (2.2) and (2.5), and evaluating the contributions from 𝐅^𝜇_ext in the manner outlined above, one obtains ℱ^𝜇_𝑖=𝒢^𝜇_𝑖,𝑄^𝜇4=−3𝜋𝑎³_𝜇𝒫_𝜇+13𝛿^𝑚𝑛𝒢^𝜇_𝑚𝑛,𝒯^𝜇_𝑖=−𝜀_𝑖^𝑚𝑛𝒢^𝜇_[]𝑚𝑛,𝒮^𝜇_𝑖𝑗=83𝜋𝜂𝑎³_𝜇ℰ_𝜇^𝑖𝑗+𝒢^𝜇_𝑖𝑗, (2.13) where 𝒢^𝜇_𝑖 and 𝒢0𝑥0𝑓106^𝜇_𝑖𝑗 are defined as in (2.2), with 𝐹^𝜇_𝑖(𝐱) replaced by 𝐺^𝜇_𝑖(𝐱). On the right-hand side of 𝒮^𝜇_𝑖𝑗, we have neglected—consistent with the expansion (2.10)—a second-order term proportional to 𝑎²_𝜇Δℰ_𝜇^𝑖𝑗≪ℰ_𝜇^𝑖𝑗.

3. Axisymmetric Flow

Axisymmetric scalar- and vector-valued fields are of the generic form𝑞(𝐱)(),𝐯(𝐱)=𝑞𝑠,𝑧=𝑣^𝑠(𝑠,𝑧)𝐞_𝑠𝜑+𝑣^𝜑(𝑠,𝑧)𝐞_𝜑𝜑+𝑣^𝑧(𝑠,𝑧)𝐞_𝑧, (3.1) where (𝑠,𝜑,𝑧) denote cylindrical coordinates about the 3-axis and (𝐞_𝑠(𝜑),𝐞_𝜑(𝜑),𝐞_𝑧) the associated unit vectors along the coordinate curves. Explicitly,𝐱=𝑧ℝ𝑠,𝜑,𝑧𝑠cos𝜑𝑠sin𝜑𝑠≥0,𝜑∈,𝐞2𝜋,𝑧∈ℝ_𝑠𝜑=0cos𝜑sin𝜑,𝐞_𝜑𝜑=0−sin𝜑cos𝜑,𝐞_𝑧=001. (3.2) Notice in (3.1) that the scalar and the vector components are independent of the azimuthal angle 𝜑. Since the Stokes operator (𝐮₁,𝐮₂)↦(𝑞,𝐯), defined by the Dirchlet problem (2.7), is invariant under rotations about the 3-axis, it follows [15] that the flow (𝑞,𝐯) is axisymmetric if and only if the applied velocities (2.8) are axisymmetric. Expressing (2.10) in cylindrical coordinates one finds that this is the case if and only if the local velocity gradients (2.11) assume the special form𝒰^𝑗_𝜈=𝒰_𝜈001𝒲^𝑗_𝜈=𝒲_𝜈001ℰ_𝜈^𝑗𝑙=ℰ_𝜈1√6−1000−10002, (3.3) where we have used that tr𝓔_𝜈=0. Inserting (3.3) in (2.10), one obtains𝐮_𝜈(𝐱)=−16𝑠ℰ_𝜈𝐞_𝑠𝜑+𝑠𝒲_𝜈𝐞_𝜑𝜑+𝒰_𝜈+23𝑧−𝑍_𝜈ℰ_𝜈𝐞_𝑧(𝑠,𝑧)∈𝜕𝙱_𝜈, (3.4) where (𝑠,𝑧)∈𝜕𝙱_𝜈 means 𝐱(𝑠,𝜑,𝑧)∈𝜕𝙱_𝜈 for some (and thus all) 𝜑∈[0,2𝜋). Likewise, we will use the notation (𝑠,𝑧)∈𝙱₀.

Upon inserting (3.1) and (3.4) into (2.7), one obtains two Dirichlet problems: one for 𝑣^𝜑 and one for (𝑞,𝑣^𝑠,𝑣^𝑧). For brevity, we refer to these as problem I and problem II, respectively. Introducing the Laplace operatorsΔ_𝑚=1𝑠𝜕_𝑠𝑠𝜕_𝑠−𝑚²𝑠²+𝜕_𝑧𝜕_𝑧(𝑚=0,1), (3.5) problem I readsΔ₁𝑣^𝜑=0,(𝑠,𝑧)∈𝙱₀,𝑣^𝜑⟶0,√𝑠²+𝑧²,𝑣⟶∞_𝜑=𝑠𝒲_𝜈,(𝑠,𝑧)∈𝜕𝙱_𝜈,𝜈=1,2, (3.6) while problem II reads1𝑠𝜕_𝑠𝑠𝑣^𝑠+𝜕_𝑧𝑣^𝑧=0,𝜕_𝑠𝑞−𝜂Δ₁𝑣^𝑠=0,𝜕_𝑧𝑞−𝜂Δ₀𝑣^𝑧=0,(𝑠,𝑧)∈𝙱₀,𝑞,𝑣^𝑠,𝑣^𝑧⟶0,√𝑠²+𝑧²,𝑣⟶∞_𝑠=−16𝑠ℰ_𝜈,𝑣^𝑧=𝒰_𝜈+23𝑧−𝑍_𝜈ℰ_𝜈,(𝑠,𝑧)∈𝜕𝙱_𝜈,𝜈=1,2. (3.7)

These two problems are completely decoupled, since (3.6) is determined by the 𝒲_𝜈, while (3.7) is determined by the 𝒰_𝜈 and ℰ_𝜈. Also, from the general uniqueness theorem [8, 43], it follows that both problems have at most one solution.

In the following, we parameterize the surfaces 𝜕𝙱_𝜈 in terms of the azimuthal angle 𝜑 and a polar angle 𝜗 according to𝜕𝙱_𝜈=𝑠𝐱=^𝜈(𝜗)𝑠cos𝜑_𝜈(𝜗)𝑧sin𝜑_𝜈(𝜗)∣𝜗∈[0,𝜋],𝜑∈[0,2𝜋)(𝜈=1,2), (3.8) where the mapping 𝜗↦(𝑠_𝜈(𝜗),𝑧_𝜈(𝜗))∈𝜕𝙱_𝜈 parameterizes the solutions of the surface constraint𝑠²+𝑧−𝑍_𝜈²=𝑎²_𝜈. (3.9) This mapping is assumed to be smooth and nonsingular, with 𝑧̇_𝜈(𝜗)≠0 for all 𝜗∈(0,𝜋). Then, sgn𝑧̇_𝜈=±1 denotes the orientation of the parameterization, and 𝑑𝑆=𝑎|𝑧̇_𝜈(𝜗)|𝑑𝜗𝑑𝜑 is the surface element. The parameterization will be made explicit in the following paper [40], where we introduce bispherical coordinates; for the present purpose, it is sufficient stay general.

Since the mapping (𝑞,𝐯)↦(𝐆¹,𝐆²), defined by the constitutive equation (2.9), is invariant under rotations about the 3-axis, it follows [15] that the surface-force densities 𝐆^𝜇(𝐱) caused by axisymmetric flow (𝑞,𝐯) are again axisymmetric𝐆^𝜇(𝐱)=𝐺^𝜇_𝑠(𝜗)𝐞_𝑠𝜑+𝐺^𝜇_𝜑(𝜗)𝐞_𝜑𝜑+𝐺^𝜇_𝑧(𝜗)𝐞_𝑧. (3.10)

Calculating the associated moments, defined as in (2.2), one can easily perform the azimuthal integrations to obtain the generic form𝒢^𝜇_𝑖=00𝒢^𝜇₃𝒢^𝜇_𝑖𝑗=𝒢^𝜇₁₁−𝒢^𝜇₂₁0𝒢^𝜇₂₁𝒢^𝜇₁₁000𝒢^𝜇₃₃, (3.11) where the nonvanishing elements involve polar integrations over the components of (3.10). Inserting these in (2.13), one findsℱ^𝜇_𝑖=ℱ^𝜇001,𝒬^𝜇4=−3𝜋𝑎³_𝜇𝒫_𝜇+Δ𝒬^𝜇,𝒮^𝜇_𝑖𝑗=𝒮^𝜇1√6,𝒯−1000−10002_𝜇𝑖=𝒯^𝜇001,𝒮^𝜇=83𝜋𝜂𝑎³_𝜇ℰ_𝜇+Δ𝒮^𝜇, (3.12) withℱ^𝜇=𝒢^𝜇₃=2𝜋𝑎_𝜇sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝐺^𝜇_𝑧,𝒯^𝜇=2𝒢^𝜇₂₁=2𝜋𝑎_𝜇sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝑠_𝜇𝐺^𝜇_𝜑,Δ𝒬^𝜇=132𝒢^𝜇₁₁+𝒢^𝜇₃₃=23𝜋𝑎_𝜇sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝑧_𝜇−𝑍_𝜇𝐺^𝜇_𝑧+𝑠_𝜇𝐺^𝜇_𝑠,Δ𝒮^𝜇=23𝒢^𝜇₃₃−𝒢^𝜇₁₁=83𝜋𝑎_𝜇sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝑧_𝜇−𝑍_𝜇𝐺^𝜇_𝑧−12𝑠_𝜇𝐺^𝜇_𝑠. (3.13) To determine the components of (3.10), we insert (3.1) in (2.9) and use the incompressibility condition. This yields𝐺^𝜇_𝑠(𝜗)=1𝑎_𝜇()−𝜂𝑠𝑞𝑠,𝑧𝑎_𝜇𝜕_N−1𝑣^𝑠(𝑠,𝑧)+𝜕_T𝑣^𝑧(,𝐺𝑠,𝑧)_𝜇𝜑(𝜗)𝜂=−𝑎_𝜇𝜕_N−1𝑣^𝜑(),𝐺𝑠,𝑧_𝜇𝑧(𝜗)=1𝑎_𝜇𝑧−𝑍_𝜇𝑞()+𝜂𝑠,𝑧𝑎_𝜇𝜕_T+𝑧−𝑍_𝜇𝑠𝑣^𝑠(𝑠,𝑧)−𝜕_N𝑣^𝑧(𝑠,𝑧), (3.14) where𝜕_N=𝑠𝜕_𝑠+𝑧−𝑍_𝜇𝜕_𝑧,𝜕_T=𝑧−𝑍_𝜇𝜕_𝑠−𝑠𝜕_𝑧 (3.15) denote the derivatives in the normal and tangential directions, respectively, and (𝑠,𝑧) is to be put equal to (𝑠_𝜇(𝜗),𝑧_𝜇(𝜗)) after the differentiations have been carried out. The expressions (3.14) must now be inserted in (3.13). To simplify the integrals we use the boundary conditions on 𝜕𝙱_𝜇, as specified in (3.6) and (3.7), and the relations𝑑=̇𝑑𝜗𝑠𝜕_𝑠+𝑧𝜕̇_𝑧𝑠𝑠̇+𝑧−𝑍_𝜇𝑧̇=0𝑧𝜕̇_T𝑑=−𝑠𝑑𝜗, (3.16) which follow from (3.9) and (3.15). Employing the notation𝑞_𝜇(𝜗)=𝑞(𝑠,𝑧)_{(𝑠,𝑧)=(𝑠𝜇(𝜗),𝑧𝜇(𝜗))}𝜕_N𝑣^𝑧_𝜇(𝜗)=𝜕_N𝑣^𝑧(𝑠,𝑧)_{(𝑠,𝑧)=(𝑠𝜇(𝜗),𝑧𝜇(𝜗))}, (3.17) and so forth, and using that𝑠_𝜇(0),𝑧_𝜇(0)=0,𝑍_𝜇−𝑎sgn𝑧̇_𝜇𝑠_𝜇(𝜋),𝑧_𝜇(𝜋)=0,𝑍_𝜇+𝑎sgn𝑧̇_𝜇, (3.18) we finally obtainℱ^𝜇=2𝜋𝜂sgn𝑧̇_𝜇^𝜋₀1𝑑𝜗𝑠2𝜂^2𝜇𝑞̇_𝜇−𝑧̇_𝜇𝜕_N𝑣^𝑧_𝜇,𝒯^𝜇=83𝜋𝜂𝑎³_𝜇𝒲_𝜇−2𝜋𝜂sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝑠_𝜇𝜕_N𝑣^𝜑_𝜇,Δ𝒬^𝜇=23𝜋𝑎³_𝜇𝑞_𝜇(0)+𝑞_𝜇(𝜋−2)3𝜋𝜂sgn𝑧̇_𝜇^𝜋₀𝑎𝑑𝜗^2𝜇𝜂𝑧_𝜇−𝑍_𝜇𝑞̇_𝜇+𝑧̇_𝜇𝑧_𝜇−𝑍_𝜇𝜕_N𝑣^𝑧_𝜇+𝑠_𝜇𝜕_N𝑣^𝑠_𝜇,Δ𝒮^𝜇4=−3𝜋𝜂𝑎³_𝜇ℰ_𝜇+83𝜋𝜂sgn𝑧̇_𝜇×^𝜋₀1𝑑𝜗𝑧2𝜂^𝜇−𝑍_𝜇𝑠²_𝜇𝑞̇_𝜇−𝑧̇_𝜇𝑧_𝜇−𝑍_𝜇𝜕_N𝑣^𝑧_𝜇−12𝑠_𝜇𝜕_N𝑣^𝑠_𝜇. (3.19)

Notice from (3.7) and (3.16) that1𝜂𝑞̇_𝜇=𝑠_𝜇̇Δ₁𝑣^𝑠_𝜇+𝑧̇_𝜇Δ₀𝑣^𝑧_𝜇. (3.20)

From (3.19), it follows immediately that for axisymmetric flow, the 𝒯^𝜇 are determined via 𝑣^𝜑 by the 𝒲_𝜈 (problem I, (3.6)), while the (ℱ^𝜇,𝒮^𝜇,Δ𝒬^𝜇) are determined via (𝑞,𝑣^𝑠,𝑣^𝑧) by the (𝒰_𝜈,ℰ_𝜈) (problem II, (3.7)). It is convenient to quantify the resulting linear relationships between the force moments and the applied velocity gradients in the scaled form [15]𝒯^𝜇=_𝜈=1,2𝜎^𝜇𝜈_𝑇𝑊𝒲_𝜈,ℱ^𝜇𝒮^𝜇Δ𝑄^𝜇=_𝜈=1,2𝜎^𝜇𝜈_𝐹𝑈𝜎^𝜇𝜈_𝑆𝑈𝜎^𝜇𝜈_𝑄𝑈𝜎^𝜇𝜈_𝐹𝐸𝜎^𝜇𝜈_𝑆𝐸𝜎^𝜇𝜈_𝑄𝐸𝒰_𝜈ℰ_𝜈𝜇=1,2, (3.21) where the 𝜎^𝜇𝜈_⋅⋅ are dimensionless resistance functions andℱ^𝜇=6𝜋𝜂𝑎_𝜇^−1/2ℱ^𝜇𝒰_𝜈=6𝜋𝜂𝑎_𝜈^1/2𝒰_𝜈,𝒯^𝜇=8𝜋𝜂𝑎³_𝜇^−1/2𝒯^𝜇𝒲_𝜈=8𝜋𝜂𝑎³_𝜈^1/2𝒲_𝜈,𝒮^𝜇=203𝜋𝜂𝑎³_𝜇^−1/2𝒮^𝜇ℰ_𝜈=203𝜋𝜂𝑎³_𝜈^1/2ℰ_𝜈,Δ𝑄^𝜇=203𝜋𝜂𝑎³_𝜇^−1/2Δ𝑄^𝜇. (3.22) The scaling (3.22) is based on the one-sphere results [13]ℱ_𝑖=6𝜋𝜂𝑎𝒰_𝑖,𝒯_𝑖=8𝜋𝜂𝑎³𝒲_𝑖,ℱ_𝑖𝑗=203𝜋𝜂𝑎³ℰ_𝑖𝑗, (3.23) and differs from earlier scalings [11–13]. An advantage of the present choice is that our resistance functions satisfy 𝜎^𝜇𝜈_𝐹𝑈,𝜎^𝜇𝜈_𝑇𝑊,𝜎^𝜇𝜈_𝑆𝐸→𝛿^𝜇𝜈 as the distance between the particles goes to infinity. From Lorentz-reciprocity [8, 9], it follows that 𝜎^𝜇𝜈_𝑆𝑈=𝜎^𝜈𝜇_𝐹𝐸.

4. Stream Function

Since problem I is already a scalar one and needs no further treatment at this point, we now focus on problem II and the calculation of (ℱ^𝜇,𝒮^𝜇,Δ𝒬^𝜇). To solve (3.7), we make the ansatz𝑣^𝑠=𝜕_𝑧𝑣Ψ,_𝑧1=−𝑠𝜕_𝑠𝑠Ψ, (4.1) where Ψ(𝑠,𝑧) is a scalar stream function. Obviously, (4.1) satisfies the condition of incompressibility, (1/𝑠)𝜕_𝑠𝑠𝑣^𝑠+𝜕_𝑧𝑣^𝑧=0. Our stream function differs by a factor 𝑠 from the classical stream function introduced by Stokes [30] and used hitherto [31–33, 36, 37] but has the advantage of being a biharmonic function. In fact, using the identity𝜕_𝑠1𝑠𝜕_𝑠1𝑠=𝑠𝜕_𝑠𝑠𝜕_𝑠−1𝑠², (4.2) one obtains from (3.5), (3.7), and(4.1)𝜕_𝑠𝑞=𝜂𝜕_𝑧Δ₁𝜕Ψ,_𝑧𝑞=−𝜂Δ₀1𝑠𝜕_𝑠1𝑠Ψ=−𝜂𝑠𝜕_𝑠𝑠Δ₁Ψ, (4.3) which leads to0=𝜕_𝑧𝜕_𝑠𝑞−𝜕_𝑠𝜕_𝑧𝑞=𝜂𝜕_𝑧𝜕_𝑧+𝜕_𝑠1𝑠𝜕_𝑠𝑠Δ₁Ψ=𝜂Δ₁Δ₁Ψ. (4.4) To satisfy the boundary condition (𝑞,𝑣^𝑠,𝑣^𝑧)→0 as √𝑠²+𝑧²→∞, it is sufficient to require that Ψ stay bounded at infinity.

Next, we turn to the boundary conditions on the particle surfaces. Since the applied velocities 𝐮_𝜈 satisfy the homogenous Stokes equations inside particles, they can be expressed in terms of stream functions, too. We denote these applied stream functions by Φ_𝜈, (𝜈=1,2), and determine them from the conditions𝑢^𝑠_𝜈=𝜕_𝑧Φ_𝜈=−16𝑠ℰ_𝜈,𝑢^𝑧_𝜈1=−𝑠𝜕_𝑠𝑠Φ_𝜈=𝒰_𝜈+23𝑧−𝑍_𝜈ℰ_𝜈(𝑠,𝑧)∈𝙱_𝜈,Φ,𝜈=1,2_𝜈regular(𝑠⟶0), (4.5) which follow from (3.4) and (4.1). The solutions areΦ_𝜈()1𝑠,𝑧=−2𝑠𝒰_𝜈−16𝑠𝑧−𝑍_𝜈ℰ_𝜈(𝑠,𝑧)∈𝙱_𝜈,𝜈=1,2. (4.6)

It is easy to show that the boundary condition 𝐯=𝐮_𝜈 on 𝜕𝙱_𝜈 is satisfied, ifΨ=Φ_𝜈𝜕_NΨ=𝜕_NΦ_𝜈𝐱∈𝜕𝙱_𝜈. (4.7)

The proof is based on the identities𝑑𝑑𝜗Ψ−Φ_𝜈=−𝑠̇𝑣^𝑧−𝑢^𝑧_𝜈+𝑧̇𝑣^𝑠−𝑢^𝑠_𝜈−𝑠̇𝑠Ψ−Φ_𝜈𝜕_NΨ−Φ_𝜈=𝑧−𝑍_𝜈𝑣^𝑠−𝑢^𝑠_𝜈−𝑠𝑣^𝑧−𝑢^𝑧_𝜈−Ψ−Φ_𝜈(𝜗∈(0,𝜋)), (4.8) which follow from (3.15), (3.16), (4.1), and (4.5). Now, if (4.7) holds, these identities reduce to the homogeneous 2×2-system𝑧̇−𝑠̇𝑧−𝑍_𝜈𝑣−𝑠^𝑠−𝑢^𝑠_𝜈𝑣^𝑧−𝑢^𝑧_𝜈=00, (4.9) which has a nonzero determinant. Hence, (4.9) admits only the trivial solution, 𝐯−𝐮_𝜈=0.

To summarize the flow velocity (𝑣^𝑠,𝑣^𝑧) can be expressed by (4.1) in terms of a scalar stream function Ψ(𝑠,𝑧) that satisfies the biharmonic two-body Dirichlet problem Δ₁Δ₁Ψ=0(𝑠,𝑧)∈𝙱₀,Ψbounded√𝑠²+𝑧²⟶∞,Ψ=Φ_𝜈𝜕_NΨ=𝜕_NΦ_𝜈(𝑠,𝑧)∈𝜕𝙱_𝜈,𝜈=1,2, (4.10)where the applied stream functions Φ_𝜈 are given by (4.6).

The pressure at a given point (𝑠,𝑧)∈𝙱₀ can, in principle, be calculated from the stream function by integration of (4.3) along any curve in 𝙱₀ that starts at infinity and ends in the point (𝑠,𝑧). However, for all but a few special points, these integrals cannot be performed analytically [40].

To establish the existence of a stream function as specified by (4.10), we consider two harmonic functions, 𝜓₁ and 𝜓₂, that vanish at infinity,Δ₁𝜓_𝑖𝜓=0_𝑖1→𝒪√𝑠²+𝑧²√𝑠²+𝑧²()⟶∞𝑖=1,2. (4.11) It is known from potential theory that such functions exist in bispherical geometry, and this is, of course, explicitly confirmed by the solution of problem I. Following [31], we putΨ(𝑠,𝑧)=𝜓₁(𝑠,𝑧)+𝑧𝜓₂(𝑠,𝑧), (4.12) and obtainΔ₁Ψ=2𝜕_𝑧𝜓₂Δ₁Δ₁𝑧Ψ=0,Ψ⟶𝒪√𝑠²+𝑧²√𝑠²+𝑧²⟶∞. (4.13) Thus Ψ, as given by (4.12), is biharmonic and stays bounded at infinity. Moreover, the ansatz (4.12), with two independent harmonic functions, is general enough to allow specification of Ψ  and 𝜕_NΨ on the boundaries, as will be confirmed explicitly in [40].

Finally, we express the force moments (ℱ^𝜇,𝒮^𝜇,Δ𝒬^𝜇) in the stream function by inserting (3.20) and (4.1) into (3.19). We quote the auxiliary formulae, valid on 𝜕𝙱_𝜇,1𝜂𝑑𝑞𝑧̇1𝑑𝜗=−𝑠𝜕_N+1Δ₁Ψ𝑧𝜕̇_N𝑣^𝑧=−𝑧𝑠̇Δ₁𝑑Ψ+𝑑𝜗𝑠𝑣^𝑠,𝑧̇𝑧−𝑍_𝜇𝜕_N𝑣^𝑧=−𝑧̇𝑧−𝑍_𝜇𝑠Δ₁Ψ−𝑧𝑠̇𝑣^𝑠+𝑑(⋯),̇𝑑𝜗𝑧𝑠𝜕_N𝑣^𝑠=𝑧̇𝑧−𝑍_𝜇𝑠Δ₁Ψ+𝑧𝑠̇𝑣^𝑠+𝑑(⋯),𝑑𝜗 (4.14) where (𝑠,𝑧)=(𝑠_𝜇(𝜗),𝑧_𝜇(𝜗)) and the unquoted expressions (⋯) tend to zero for 𝑠→0 and, therefore, do not contribute to the integrals in (3.19). The final results areℱ^𝜇=−𝜋𝜂sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑧̇_𝜇𝑠_𝜇(𝜕_N−1)Δ₁Ψ_𝜇,𝒮^𝜇=23𝜋𝜂sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑠_𝜇̇𝑠²_𝜇𝜕_N−2Δ₁Ψ_𝜇,Δ𝒬^𝜇=23𝜋𝑎³_𝜇𝑞_𝜇(0)+𝑞_𝜇(𝜋)−2𝜋𝜂3𝑎²_𝜇sgn𝑧̇_𝜇^𝜋₀𝑑𝜗𝑠_𝜇̇𝜕_N+1Δ₁Ψ_𝜇. (4.15) With (4.15), we have succeeded to express the forces and stresslets in the stream function alone. The pressure moments require, in addition, the quantities 𝑞_𝜇(0) and 𝑞_𝜇(𝜋), that is, the pressure at the poles of the spheres. As will be shown [40], the poles belong to the few points, where the pressure can be determined analytically in terms of the stream function.

The further evaluations in bispherical coordinates and the associated series expansions for the resistance functions will be presented in [40].

Acknowledgment

T. El Naqeeb expresses her gratitude to RWTH Aachen University for its hospitality.