The Saffman-Taylor problem on a sphere

Fernando Parisio, Fernando Moraes, José A. Miranda

Laboratório de Física Teórica e Computacional, Departamento de Física, Universidade Federal de Pernambuco, Recife, PE 50670-901 Brazil

Michael Widom

Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213 USA (July 14, 2000)

Abstract

The Saffman-Taylor problem addresses the morphological instability of an interface separating two immiscible, viscous fluids when they move in a narrow gap between two flat parallel plates (Hele-Shaw cell). In this work, we extend the classic Saffman-Taylor situation, by considering the flow between two curved, closely spaced, concentric spheres (spherical Hele-Shaw cell). We derive the mode-coupling differential equation for the interface perturbation amplitudes and study both linear and nonlinear flow regimes. The effect of the spherical cell (positive) spatial curvature on the shape of the interfacial patterns is investigated. We show that stability properties of the fluid-fluid interface are sensitive to the curvature of the surface. In particular, it is found that positive spatial curvature inhibits finger tip-splitting. Hele-Shaw flow on weakly negative, curved surfaces is briefly discussed.

Typeset using REVTEX

I. INTRODUCTION

Formation and evolution of dynamic structures is an exciting area of nonlinear phenomenology. Of particular practical and theoretical interest is the hydrodynamic pattern formation of the growing interface between two fluids. One of the best studied patternformation systems of this type is the Saffman-Taylor problem [1,2]: it addresses motion of two viscous, immiscible fluids in the narrow space between two parallel, flat plates known as a Hele-Shaw cell. When a fluid of low viscosity displaces a fluid of higher viscosity, the interface between them becomes unstable and deforms. Dynamic competition leads to the formation of fingering structures.

Experiments and theory focus on two principal geometries: (i) rectangular [1,2] and (ii) radial [3–5]. In rectangular cells the unperturbed interface is straight and the unperturbed flow is uniform and parallel to cell walls. In the radial case the unperturbed interface is circular with the less viscous fluid pumped into the more viscous one at a point and the flow radially outward. For both situations, the initial development of the interface instability tracks the predictions of linear stability theory [1–5]. After the initial surface deformation, as the unstable modes of perturbation grow, they become coupled in a weakly nonlinear stage of evolution [6,7]. Finally, the system evolves to a complicated late stage, characterized by formation of fingering structures, in which nonlinear effects dominate [2]. As a result, beautiful finger-bubble undulated structures are formed in rectangular cells, while visually striking, fan-like, branched patterns rise in radial flow. Spreading, splitting and competition are the three basic growth mechanisms of the viscous fingering process [2]. In this work we are particularly interested in tip-splitting events.

Despite the extensive experimental and theoretical work on the Saffman-Taylor problem in both radial and rectangular setups, the majority of the studies focus on flow in flat Hele-Shaw cells. Curiously, the dynamic behavior for flow in flat cells is described by the very same set of equations as those for flow in porous media [2], which is indeed a highly non-flat environment, characterized by voids and curved internal surfaces. Even though viscous fingering is not restricted to occur between flat surfaces, the study of the Saffman-Taylor problem on curved surfaces has been practically neglected. One exception is the experimental work by Zhao and Maher [8], which considers flow in a cylindrical Hele-Shaw cell. They used the cylindrical cell as an experimental realization of a flat rectangular Hele-Shaw cell with periodic boundary conditions. Their experiments showed that instabilities in cells with periodic boundary conditions are qualitatively similar to instabilities in cells with physical sidewalls. The flow is performed parallel to the axis of a cylindrical cell with a large radius of curvature, so the authors did not explore the influence of cell curvature.

On the theoretical side, Entov and Etingof [9] considered the general mathematical problem of viscous flow in non-planar Hele-Shaw cells, in the zero-surface-tension limit. They applied conformal mapping techniques to derive a class of explicit solutions for the shape of the fluid domain under study. As in reference [8], the authors in [9] were not interested in the influence of the cell curvature in the shape of the interfacial patterns.

Spatial curvature has been a relatively overlooked feature in the study of other pattern formation systems as well. Only recently researchers started to investigate how the properties of the patterns could be affected by the curvature of the surface in which such structures evolve [10–15]. In the early 1990’s Levine and collaborators studied coarsening of twodimensional foams [10] and grain growth [11] on curved surfaces. They showed that the stability properties of such froth bubbles and grains depended on the curvature of the surface. Ref. [12] reviews the dependence of many chemical and biochemical surface processes on curvature. Numerical studies of reaction-diffusion systems in curved spaces examine the evolution of spiral waves [13] and the occurence of Turing patterns on a sphere [14]. These numerical simulations [13,14] indicate that curvature imposes geometrical restrictions on the shape of the patterns. Finally, Schoenborn and Desai [15] studied the intra-surface kinetics of phase ordering on curved surfaces.

For the Saffman-Taylor problem, the interplay between Hele-Shaw cell curvature and interfacial pattern formation is largely unexplored in the present literature. However, in a recent mode-coupling analysis of radial flow in flat cells, Miranda and Widom [7] suggested that cell curvature could be used as a control parameter to regulate the tendency towards finger tip-splitting. A thorough investigation of the relationship between cell curvature and the fluid-fluid interface dynamics still needs to be addressed. In this work we begin such investigations, focusing on the Saffman-Taylor problem on a sphere.

The outline of the work is the following: section II defines Hele-Shaw flow between concentric, thin spherical shells, considering fluid injection (withdrawal) at the north (south) pole. We derive a differential equation describing the early nonlinear evolution of the interface modes. In section III, we interpret results obtained in section II and investigate both the linear and weakly nonlinear evolution of the system. The flat space limit of infinite radius of curfature is examined in section III A. Section III B discusses the linear growth rates of unstable modes and relates these to the degree of lattitude of the unperturbed interface. Nonlinear analysis in section III C concentrates on the effect of cell curvature on finger tip-splitting. We show that positive spatial curvature inhibits splitting. Flow on the northern and southern hemispheres are contrasted, and a symmetry-breaking behavior is detected: for the southern hemisphere tip-splitting is replaced by finger tip-sharpening. Flow on weakly negative curved surfaces is briefly discussed. Section IV presents our final remarks. An appendix derives Darcy’s law for flow between concentric spheres.

II. THE MODE COUPLING DIFFERENTIAL EQUATION

Consider two immiscible, incompressible, viscous fluids, flowing in a narrow gap of thickness b, between two concentric, thin spherical shells (see figure 1). We name this device the spherical Hele-Shaw cell. Assume that b is smaller than any other length scale in the problem, so the system is effectively two-dimensional. The flow takes place in the surface of a two-dimensional sphere, embedded in three-dimensional space, and endowed with the metric [16]

ds^{2} = dρ^{2} + a^{2} sin^{2} (\frac{ρ}{a}) dφ^{2},  (1)

where a is the radius of curvature of the sphere, 0 ≤ φ < 2π denotes the polar angle measured on the sphere and 0 ≤ ρ ≤ πa is the geodesic distance from the radial flow source, located at the sphere’s north pole. The sphere has a constant, positive Gaussian curvature K = 1/a^{2} and a constant mean curvature H = 1/a. We presume the Gaussian curvature is more relevant than the mean curvature for reasons discussed in sections III A and III C. At any stage of our calculation, the “flat-cell” limit a → ∞ (or, equivalently K → 0) gives all the well-known familiar results for flow in flat Hele-Shaw cells.

Denote the viscosities of the upper and lower fluids, respectively as η_{1} and η_{2}. Between the two fluids there exists a surface tension σ (figure 1). The flows are assumed to be irrotational, except at the interface. Fluid 1 is injected into fluid 2 through an inlet located at the sphere’s north pole, at a given flow rate Q, which is the area covered per unit time. Fluid 2 is simultaneously withdrawn, at the same rate, through an outlet placed at the south pole.

During the flow, the fluid-fluid interface has a perturbed shape described as ρ = R ≡ R+ ζ(φ, t). The interface perturbation amplitude is represented by ζ(φ, t), and R = R(t) denotes the time-dependent unperturbed radius. We can calculate R(t) from the time-dependent surface area A(t) = 4πa^{2} sin^{2}(R(t)/2a). For flow rate Q we write A(t) = 4πa^{2} sin^{2}(R_{0}/2a) + Qt then solve for

R(t) = a arccos (C_{0} − \frac{Qt}{2πa^{2}} ) ,  (2)

where C_{0} = cos(R_{0}/a), and R_{0} is the unperturbed radius at t = 0. The unperturbed shape is a polar cap of geodesic radius ρ = R, surface area A and circumference L = 2πa sin (R/a). Note the identity Q = vL where v = dR/dt is the velocity of the unperturbed interface.

We express the net perturbation ζ(φ, t) in the form of a Fourier expansion

ζ(φ, t) = ∑_{n=−∞}^{+∞} ζ_{n}(t) exp (inφ),  (3)

where

ζ_{n}(t) = \frac{1}{2π} ∫_{0}^{2π} ζ(φ, t) exp (−inφ) dφ  (4)

denotes the complex Fourier mode amplitudes and n=0, ±1, ±2, ... is the discrete azimuthal wave number. To keep the area of the perturbed shape independent of the perturbation ζ, we express the Fourier mode corresponding to n = 0 in the expansion (3) as

ζ_{0}(t) = −\frac{1}{2a} cot (\frac{R}{a}) ∑_{n≠0}^{∞} |ζ_{n}(t)|^{2}.  (5)

The constraint (5) is intrinsically a nonlinear concern and is not required in linear stability analysis.

Since we are interested in both linear and early nonlinear behavior of the system, we must derive a differential equation for ζ_{n}, correct to second order. This second-order modecoupling equation considers the presence of a full spectrum of modes. We begin the hydrodynamic study of the system by considering a generalized version of the usual Darcy’s law [1,2], adjusted to describe flow between concentric spheres (see Appendix A)

v_{j} = −\frac{b^{2}}{12η_{j}} ∇p_{j},  (6)

where v_{j} = v_{j}(ρ, φ) and p_{j} = p_{j}(ρ, φ) are, respectively, the velocity and pressure in fluids j = 1 and 2. The gradient in equation (6), associated with the metric (1), is [16]

∇ = \frac{∂}{∂ρ} ρ̂+ \frac{1}{a sin(ρ/a)} \frac{∂}{∂φ} φ̂,  (7)

where the unit vectors ρ̂ and φ̂ point in the direction of increase of ρ and φ, respectively. Equation (7) was obtained from the corresponding three-dimensional expression for the gradient in spherical coordinates (r, θ, φ), by keeping r = a and noting that θ = ρ/a.

At the interface, the pressure difference between the two fluids is governed by the mean curvature (the average of the two principal curvatures) of the fluid-fluid interface [2]. We can identify the directions of the principal curvatures in the limit of b smaller than any other length scale by considering a “tangent Hele-Shaw cell” consisting of two parallel planes tangent, respectively to the inner and outer spheres at some point along the interface. Within the tangent Hele-Shaw cell, one principal curvature is associated with the interface profile in the direction perpendicular to the tangent planes. We call this curvature κ_{⊥} and note that it is of order 1/b with a specific value set by interface contact angles. The remaining direction of principal curvature is parallel to the tangent plane and tangent to the interface. We call this curvature κ_{∥}. Then the pressure jump boundary condition may be written

(p_{1} − p_{2})|_{R} = σ (κ_{∥} + κ_{⊥}) |_{R}.  (8)

As was the case for flow in flat Hele-Shaw cells, κ_{⊥} is much larger than κ_{∥} but is is nearly constant [17,18]. This curvature does not significantly affect the motion in our problem, because its gradient is nearly zero.

Since the closed boundary describing the fluid-fluid interface is itself on the top of a curved surface (sphere), the calculation of the “intra-surface” interface curvature κ_{∥} is not as simple as it was in the flat-cell case [4,16]. Taking into consideration the fact that the interface evolves in the surface of a sphere of radius a, we derive a slightly involved expression for the fluid-fluid interface curvature

κ_{∥} = \frac{[cos(ρ/a) a^{2} sin^{2}(ρ/a) + 2 cos(ρ/a) (\frac{∂ρ}{∂φ})^{2} − a sin(ρ/a) \frac{∂^{2}ρ}{∂φ^{2}} ]}{[a^{2} sin^{2}(ρ/a) + (\frac{∂ρ}{∂φ})^{2}]^{3/2}}.  (9)

The sign convention for the curvature κ_{∥} is such that a circular interface above the equator has positive curvature, whereas it has negative curvature below the equator. Keeping terms up to second order in the perturbation amplitude ζ, we rewrite the interfacial curvature as

(κ_{∥})|_{R} = \{\frac{C}{aS} − \frac{1}{a^{2}S^{2}} (ζ + \frac{∂^{2}ζ}{∂φ^{2}} ) + \frac{C}{a^{3}S^{3}} [ζ^{2} + \frac{1}{2} (\frac{∂ζ}{∂φ})^{2} + 2ζ \frac{∂^{2}ζ}{∂φ^{2}} ]\} ,  (10)

where we have introduced the shorthand notation S = sin(R/a) and C = cos(R/a). Exactly at the equator (C = 0), only the term in (10) that is linear in ζ survives.

Taking advantage of the irrotational and incompressible flow conditions, we define the velocity potential ϕ_{j} in each of the fluids, where v_{j} = −∇ϕ_{j}. The velocity potential satisfies Laplace’s equation ∇^{2}ϕ_{j} = 0, where ∇^{2} is the two-dimensional Laplacian defined on the surface of the sphere. Combining the velocity potential with equations (8) and (9) for the pressure difference and the generalized Darcy’s law (6), we write the equation of motion

A (\frac{ϕ_{1}|_{R} + ϕ_{2}|_{R}}{2} ) − (\frac{ϕ_{1}|_{R} − ϕ_{2}|_{R}}{2} ) = −α (κ_{∥})|_{R},  (11)

where

A = \frac{η_{2} − η_{1}}{η_{2} + η_{1}}  (12)

is the viscosity contrast and

α = \frac{b^{2}σ}{12(η_{1} + η_{2})},  (13)

contains the surface tension.

Now define Fourier expansions for the velocity potentials ϕ_{j}. Far from the interface the velocity field should approach the unperturbed steady flow with a circular interface of radius R. Thus for ρ → 0 and ρ → πa the velocity potentials ϕ_{j} approach ϕ_{j}^{0}, the velocity potentials for purely radial (ρ̂ direction) flow, satisfying Laplace’s equation

ϕ_{j}^{0} = −\frac{Q}{2π} log [\frac{tan(ρ/2a)}{tan(R/2a)}] + D_{j},  (14)

where D_{j} are independent of ρ and φ. The general velocity potentials obeying all these requirements are

ϕ_{j} = ϕ_{j}^{0} + ∑_{n≠0} ϕ_{jn}(t) [\frac{tan(R/2a)}{tan(ρ/2a)}]^{(−1)^{j} |n|} exp(inφ).  (15)

The trigonometric dependence on ρ transforms, in the flat-cell limit a → ∞, into the ratio R/ρ. In order to calculate the mode coupling differential equation for the system, we substitute expansions (14) and (15) into the equation of motion (11), keep second order terms in the perturbation amplitudes, and Fourier transform them.

To conclude our derivation we need additional relations expressing the velocity potentials in terms of the perturbation amplitudes. To find these, consider the kinematic boundary condition which states that the normal components of each fluid’s velocity at the interface equals the velocity of the interface itself [19]. Using the gradient (7) we write the kinematic boundary condition for flow in a sphere as

\frac{∂R}{∂t} = [\frac{1}{a^{2} sin^{2}(ρ/a)} \frac{∂ρ}{∂φ} \frac{∂ϕ_{j}}{∂φ} ]_{ρ=R} − (\frac{∂ϕ_{j}}{∂ρ} )_{ρ=R} .  (16)

Inserting expression R = R + ζ(φ, t) and equation (15) for ϕ_{j} into the kinematic boundary condition (16), we solve for ϕ_{jn}(t) consistently to second order in ζ to find

ϕ_{1n}(t) = −\frac{aS}{|n|}ζ̇_{n} − \frac{QC}{2πaS|n|}ζ_{n}
+ ∑_{n′≠0} (sgn(nn′) − \frac{C}{|n|}) ζ̇_{n′}ζ_{n−n′} + \frac{QC}{2πa^{2}S^{2}} ∑_{n′≠0} (sgn(nn′) + \frac{S^{2}}{2C|n|}) ζ_{n′}ζ_{n−n′},

(17)

and

ϕ_{2n}(t) = \frac{aS}{|n|}ζ̇_{n} + \frac{QC}{2πaS|n|}ζ_{n}
+ ∑_{n′≠0} (sgn(nn′) + \frac{C}{|n|}) ζ̇_{n′}ζ_{n−n′} + \frac{QC}{2πa^{2}S^{2}} ∑_{n′≠0} (sgn(nn′) − \frac{S^{2}}{2C|n|}) ζ_{n′}ζ_{n−n′}.

(18)

The overdot denotes total time derivative. The sign function sgn(nn′) = 1 if (nn′) > 0 and sgn(nn′) = −1 if (nn′) < 0.

We can use relations (17) and (18) to replace the velocity potentials ϕ_{j} in the equation of motion (11) with the perturbation ζ and its time derivative ζ̇. Keeping only quadratic terms in the perturbation amplitude, and equating Fourier modes n on each side of equation (11), leads to the differential equation for perturbation amplitudes ζ_{n}. For n ≠ 0,

ζ̇_{n} = λ(n) ζ_{n} + ∑_{n′≠0} [F (n, n′) ζ_{n′}ζ_{n−n′} + G(n, n′) ζ̇_{n′}ζ_{n−n′}] ,  (19)

where

λ(n) = [\frac{Q}{2πa^{2}S^{2}} (A|n| − C) − \frac{α}{a^{3}S^{3}} |n|(n^{2} − 1)]  (20)

is the linear growth rate, and

F (n, n′) = \frac{|n|}{aS} \{\frac{QAC}{2πa^{2}S^{2}} [\frac{1}{2} − sgn(nn′)] − \frac{αC}{a^{3}S^{3}} [1 − \frac{n′}{2} (3n′ + n)]\} + \frac{Q}{4πa^{3}S} ,  (21)
G(n, n′) = \frac{1}{aS} \{A|n|[1 − sgn(nn′)] − C\}  (22)

are the second-order mode coupling terms. Equation (19) is the mode coupling equation of the Saffman-Taylor problem in a spherical Hele-Shaw cell. It gives us the time evolution of the perturbation amplitudes ζ_{n} accurate to second order. In the following sections we study equation (19) in more detail, and investigate the role played by cell geometry in the interface dynamics.

III. DISCUSSION

We use the mode coupling equation (19) to investigate the linear instability of individual modes and the coupling of a small number of modes. The most noteworthy effect of curvature we identify concerns its influence on finger tip-splitting. Tip-splitting is related to the influence of a fundamental mode on the growth of its harmonic [7]. We take n as the fundamental and 2n as the harmonic. To observe interfacial instability of the fundamental mode n, we must have λ(n) > 0. This occurs if the destabilizing contribution QA in eq. (20) is positive and sufficiently large compared with the stabilizing surface tension term proportional to α. To observe growth of the harmonic mode 2n, we presume that QA is sufficiently large that λ(2n) is non-negative.

To simplify our discussion it is convenient to rewrite the net perturbation (3) in terms of cosine and sine modes

ζ(θ, t) = ζ_{0} + ∑_{n=1}^{∞} [a_{n}(t) cos(nθ) + b_{n}(t) sin(nθ)] ,  (23)

where a_{n} = ζ_{n} + ζ_{−n} and b_{n} = i (ζ_{n} − ζ_{−n}) are real-valued. Without loss of generality we may choose the phase of the fundamental mode so that a_{n} > 0 and b_{n} = 0. We replace the time derivative terms ȧ_{n} and ḃ_{n} by λ(n) a_{n} and λ(n) b_{n}, respectively, for consistent second order expressions. Under these circumstances equations of motion become

ȧ_{2n} = λ(2n) a_{2n} + \frac{1}{2} T (2n, n) a_{n}^{2}  (24)
ḃ_{2n} = λ(2n) b_{2n},  (25)

where the tip-splitting function is defined as

T (2n, n) = [F (2n, n) + λ(n) G(2n, n)] .  (26)

Note that the sign of T (2n, n) dictates if finger tip-splitting is favored or not by the dynamics. If T (2n, n) < 0, at second order the result is a driving term of order a_{n}^{2} forcing growth of a_{2n} < 0. With this particular phase of the harmonic forced by the dynamics, the n outwards-pointing fingers of the fundamental mode n tend to split. In this case the driving term in equation of motion (24) spontaneously generates the harmonic mode. In contrast, if T (2n, n) > 0 growth of a_{2n} > 0 would be favored, leading to outwards-pointing finger tip-sharpening. Note that mode b_{2n}, whose growth is uninfluenced by a_{n}, skews the fingers of mode n. In the presence of a_{2n} < 0, the role of b_{2n} is to favor one of the two split fingers over the other.

A. The flat-cell limit

We begin our discussion by analyzing the flat-cell limit of the mode-coupling expression (19). We hold fixed the unperturbed interface velocity v = Q/2πL. Three distinct flat space limits can be taken: (i) (north pole) let a → ∞ holding R finite so that C → 1 and aS → R; (ii) (south pole) let a → ∞ holding R′ ≡ πa − R finite; so that C → −1 and aS → R′; (iii) (tropical) let R/a be constant and hold fixed k ≡ 2πn/L. These three limits correspond to three physically distinct flat space flow problems. Our goal in this section is to verify that the linear and nonlinear terms in the equation of motion (19) reduce to their expected forms in the flat space limit. The actual evolution of interfaces according to the equations of motion is then discussed in greater detail in sections III B and III C.

At the north pole (case (i), C → 1), we recover the mode-coupling equations [7] of flat, radial divergent flow [20], related to outward radial motion in which fluid 1 pushes fluid 2:

λ(n) = [\frac{Q}{2πR^{2}} (A|n| − 1) − \frac{α}{R^{3}} |n|(n^{2} − 1)]  (27)

F (n, n′) = \frac{|n|}{R} \{\frac{QA}{2πR^{2}} [\frac{1}{2} − sgn(nn′)] − \frac{α}{R^{3}} [1 − \frac{n′}{2} (3n′ + n)]\}

G(n, n′) = \frac{1}{R} \{A|n|[1 − sgn(nn′)] − 1\} .

Provided the viscosity contrast A > 0, so that the less viscous fluid pushes the more viscous fluid, the interface is linearly unstable and exhibits finger growth. As the interface perturbation grows, the nonlinear mode coupling broadens and splits the outward-pointing fingers and sharpens the inward-pointing fingers. This occurs because T (2n, n) < 0 when λ(2n) ≥ 0.

In contrast, at the south pole (case (ii), C → −1), we obtain flat, radial convergent flow [20] equivalent to the inward radial motion corresponding to withdrawal of fluid 2 surrounded by fluid 1. Provided that QA > 0, so that the less viscous fluid displaces the more viscous fluid, the interface remains linearly unstable. However, now the outwardspointing fingers sharpen because for negative C we find T (2n, n) > 0 when λ(2n) ≥ 0. The asymmetry between north and south pole behaviors occurs primarily in the nonlinear term and comes from the terms proportional to C. In contrast, the linear growth rate λ is nearly symmetric between the north and south pole limits. This can be understood because interchanging the north and south poles is equivalent to reversing the sign of Q (the direction of flow) and the sign of A (interchanging the fluid viscosities) while holding the surface tension α unchanged. However, the term in λ proportional to C breaks this symmetry slightly.

A special example of the tropical case (iii) is the equatorial limit, for which C = 0 and S = 1. This case should be compared with the cylindrical flow geometry, taking the cylinder tangent to the sphere at the equator. The flat space limit of this problem indeed reduces to the problem of rectangular flow in flat space with periodic boundary conditions. In particular,

λ(n) → |k| [Av − αk^{2}]  (28)
F (n, n′) → 0
G(n, n′) → A|k| [1 − sgn(kk′)] .

As expected, there is no tendency for tip-splitting on the cylinder. This suggests our suggestion that tip-splitting is controlled by Gaussian curvature K rather than mean curvature H because variation of the radius of curvature of the cylinder alters the mean curvature while the Gaussian curvature remains zero. Tip-splitting is absent for any value of the mean curvature.

Before the limit is reached, there are small differences between the spherical and cylindrical cases. For example, the surface tension contribution to λ vanishes for k = 0 on a cylinder corresponding to translation invariance of the interface length. On a sphere the corresponding displacement moves the interface from the equator to the tropics, shortening the interface and lowering the surface energy. However, on a sphere the modes n = ±1 correspond to a global off-center shift of the circular interface preserving circular shape and perimeter. Thus the surface tension term vanishes for n = ±1 on the the sphere while this mode increases the perimeter and raises the energy on a cylinder.

B. Linear growth

Consider the purely linear contribution, which appears as the first term on the right hand side of equation (19). Since R varies with time, the linear growth rate λ(n) is time dependent as well. This implies that the actual relaxation or growth of mode n is not proportional to the factor exp[λ(n)t], but rather

ζ_{n}(t) = ζ_{n}(0) exp [∫_{0}^{t} λ(n)dt′] .  (29)

If ∫_{0}^{t} λ(n)dt′ > 0 the disturbance grows, indicating instability. Two relevant facts can be extracted from the linear growth rate: (i) the existence of a series of critical radii R_{c}(n)(defined by setting λ(n) = 0) at which the interface becomes unstable for a given mode n; (ii) the presence of a fastest growing mode n^{∗}, given by the closest integer to the maximum of equation (20) with respect to n (defined by setting dλ(n)/dn = 0). In view of equation (29) n^{∗} is not simply related to the number of fingers present, even in the early stages of pattern formation. Furthermore, in the nonlinear regime the subsequent tip-splitting process and mode competition result in a final number of fingers which can differ from the number present in the linear regime.

Cardoso and Woods [23] analyze flat space radial flow assuming the presence of a constant low level of noise during the whole evolution of the interface (their “model B”, also see Ref. [7]). The sources of noise may come from inhomogeneities on the surface of the Hele-Shaw cell, irregularities in the gap thickness b, or from thermal or pressure fluctuations [24]. The predictions of this model are in qualitative agreement with experimental observations within the linear regime [23] and the nonlinear regime [7].

Suppose that we begin with an initially circular interface that is steadily expanding. During the interface expansion each mode n is perturbed with a constant (in time) random complex amplitude ζ_{n}(0). This noise amplitude contains an n dependent random phase but its magnitude |ζ_{n}(0)| is independent of n by assumption. As the interface continues to expand, it progressively reaches critical radii R_{c}(n) for n = 2 , 3 , ... . Once a particular R_{c}(n) is reached, the perturbation amplitude ζ_{n} starts to vary with time. Within this model, the first order (linear) solution of equation (19) can be written as

ζ_{n}^{lin}(t) = \{ \begin{matrix}ζ_{n}(0) & if R < R_{c}(n) \\ ζ_{n}(0) \{(\frac{1+C_{c}}{1+C} ) [\frac{1}{T} ]^{A|n|−1} exp [C(A|n| − C) (T − 1)]\} & if R ≥ R_{c}(n). \end{matrix}  (30)

where C_{c} = cos [R_{c}(n)/a] and T = tan [R_{c}(n)/2a]/ tan (R/2a).

To see the overall effect of eq. (30), we plot the time evolution of the interface using the experimental parameters given in Paterson’s classical experiment [5]. Paterson observed the rapid growth of fingers, as air (η_{1} ≈ 0) was blown at a relatively high injection rate, Q = 9.3 cm^{2}/s, into glycerine (η_{2} ≈ 5.21 g/(cm s)) in a radial source flow Hele-Shaw cell. The thickness of the cell b = 0.15 cm and the surface tension σ = 63 dyne/cm. We take into account modes n ranging from n = 2 up to 20. We evolve from initial radius R_{0} = 0.05 cm. The noise amplitude |ζ_{n}(0)| = R_{0}/1000. Figure 2 depicts the evolution of the interface, for a random choice of phases, up to time t = 20 s. We set the radius of the sphere a = 10. We encourage the reader to compare the resulting interface with equivalent figures in flat space. Ref. [7] contains a figure in which growth conditions are similar. In particular, the phases of the initial perturbations are identical, although in the present case we choose smaller initial amplitudes. It is also of interest to compare with flat space experimental patterns found in the literature [5,25,20,3].

C. Nonlinear behavior and tip-splitting

To visualize the consequences of the second order term in the equations of motion, we solve equation (19) to second order accuracy. If we substitute the linear solution given in equation (29) into the second-order terms on the right hand side of equation (19), we obtain the differential equation

ζ̇_{n} = λ(n)ζ_{n} + W (n, t),  (31)

where

W (n, t) = ∑_{n′≠0} [F (n, n′)ζ_{n′}^{lin}ζ_{n−n′}^{lin} + G(n, n′)ζ̇_{n′}^{lin}ζ_{n−n′}^{lin}]  (32)

acts as a driving force in the linearized equation of motion (31). Equation (31) is a standard first order linear differential equation [26] with the solution

ζ_{n}(t) = \{ \begin{matrix}ζ_{n}(0) & if R < R_{c}(n) \\ ζ_{n}^{lin}(t) \{1 + ∫_{t_{c(n)}}^{t} [\frac{W (n,t′)}{ζ_{n}^{lin}(t′)}] dt′\} & if R ≥ R_{c}(n). \end{matrix}  (33)

Here t_{c}(n) is the time required for the unperturbed growth to reach radius R_{c}(n) and can be easily calculated from equation (2). This solution describes the weakly nonlinear evolution, where the dominant modes just become coupled by nonlinear effects.

We use the second order solution (33) to investigate the nonlinear coupling among various modes n. In figure 3, we plot the interface for a certain time (t = 20 s), considering the same random choice of initial phases as was employed in figure 2, and coupling all modes with 2 ≤ n ≤ 20. The nonlinear evolution leads to wider fingers and their tips become more blunt. These fingers spread and some of them start to bifurcate, by splitting at the tip. As was shown in eq. (24), tip splitting is caused by fingers of mode n driving the growth of their own harmonic 2n.

Now we turn our attention to the investigation of how curvature influences finger tipsplitting for Hele-Shaw flow on a sphere. It has been shown in reference [7] that, for flat, divergent radial flow under the condition that λ(2n) ≥ 0, we have T (2n, n) < 0. Thus, when the harmonic is able to grow, finger tip-splitting is favored. In section III A we showed that in the limit of small cell curvature K we recover the flat space limit where we know T (2n, n) is negative. To carry out our current analysis, we consider the case in which λ(2n) = 0 and see how the curvature influences T (2n, n).

Since spherical Hele-Shaw flow involves many independent parameters we change only one relevant quantity at a time to see what each one does. We hold the unperturbed interface velocity v and the unperturbed interface contour length L fixed, while varying the curvature K. This isolates the influence of spatial curvature K from the effect of variations in v and L. Fixed L means that fixed mode number n corresponds to a fixed wavelength.

We adopt an instantaneous approach: we look at the linear growth rate and mode coupling at an instant in time, ignoring the past history of how a given interface arose from some initial condition followed by growth. Such an instantaneous approach, at which L and v have a particular value, enables us to compare the behavior of interfaces evolving in distinct background curvatures, but under dynamically equivalent circumstances. Moreover, at the instant when the interface circumference L in curved space equals the circumference in flat space, if the two velocities also match, the identity Q = vL shows that the Q value in flat space equals the Q value in curved space. Therefore, it is advantageous to look at the instantaneous tendency towards tip-splitting for interfaces of identical unperturbed v and L.

Consider a particular v and L combination at the onset of growth of mode 2n (using the condition λ(2n) = 0) in the limit of flat space, where it is known that T (2n, n) < 0. To illustrate how tip-splitting varies with curvature for flow on a sphere, we plot in figure 4 the tip-splitting function T (2n, n) as a function of K. The solid(dashed) curves refer to the behavior in the northern (southern) hemisphere. By inspecting figure 4 we notice that, in the northern hemisphere there is a suppression of tip-splitting for increasingly larger curvature because the magnitude of T (2n, n) is maximum for the flat, divergent radial case (K = 0) and decreases as curvature is increased. For each n there is some value of K for which the circle of circumference L hits the equator. This is precisely the point at which northern and southern hemisphere branches of the curves meet in figure 4.

For the southern hemisphere T (2n, n) is always positive. Consequently, we should not expect finger tip-splitting of outward-pointing fingers in this hemisphere. Actually, finger tipsplitting is replaced by tip-narrowing, along with a splitting of the inward-pointing fingers. Finger tip-narrowing is regulated by the curvature K.

Another noteworthy point about figure 4 is the evident symmetry-breaking in T (2n, n) between northern and southern hemispheres. The justification for this asymmetry is similar to the one that explains why the interface is always in-out asymmetric in flat radial flow: if we are located at either pole we can always distinguish regions that are inside and outside the interface. That is why convergent and divergent flat, radial flows are not equivalent [20]. From this point of view, the asymmetry observed in figure 4 should be somehow expected. We conclude this section presenting the lowest order curvature expansion for T (2n, n) and studying how curvature K influences tip-splitting for small curvature values. Using the instantaneous approach described above, series expansion of equation (26) for the northern hemisphere yields

T (2n, n) = \frac{Q^{2} (2n^{2} + 1)}{32π^{2}α n(2n + 1)^{2}} [−\frac{Q^{2}}{8π^{2}α n^{2}(2n + 1)^{2}} + K] + O (K^{2}).  (34)

The perturbative expansion for small K (34) explicitly shows a linear correction to the flat space limit (K = 0) that reduces the magnitude of T (2n, n), inhibiting finger tip-splitting. A similar kind of expansion can be done for the southern hemisphere, resulting in enhanced tip-sharpening. Those results are in agreement with the small curvature behavior depicted in figure 4.

Equation (34) allows us to make predictions about tip-splitting behavior for Hele-Shaw flow in (weakly) negatively curved backgrounds. This type of flow could happen, for instance, between two saddle-like surfaces. Actually, flow in porous media seems to be somehow linked to flow on negatively curved surfaces. Porous materials define multiply-connected surfaces, presenting negative average Gausssian curvature [16]

K̄ = \frac{∫ KdS}{∫ dS} = \frac{4π (1 − g)}{∫ dS},  (35)

where dS denotes an infinitesimal area element and g is the so-called genus [16], which denotes the number of holes present in a given surface. Expression (35) relates the integral of the Gaussian curvature K of a given surface to its topological properties, by virtue of the Gauss-Bonnet theorem [16]. From (35) we verify that K̄ becomes progressively negative when the number of holes is increased: that is why K̄ > 0 for a sphere (g = 0), K̄ = 0 for a torus (g = 1) and K̄ < 0 for a g-torus (g ≥ 2). In this sense, a medium which is rich in pores (holes) would present negative curvature features. However, flow on surfaces of constant negative Gaussian curvature is complex and to treat the problem rigorously would require an interesting generalization of our formalism beyond the scope of the present paper.

Here we simply point out what would be the behavior for flow on surfaces that are slightly (negatively) curved. It is easy to see, by performing the substitution K → −K in (34), that negative curvature should enhance finger splitting in comparison to the flat case. This is in striking contrast to the problem of flow in positively curved surfaces, such as the sphere, where curvature leads to suppression of tip-splitting. This sets an important distinction between Hele-Shaw flows on surfaces with negative curvature and those on surfaces with positive K. Our results confirm explicitly the predictions made in reference [7] about the role of cell curvature on finger splitting.

The fact that the correction is linear in K supports our suggestion that Gaussian curvature is more relevant than mean curvature. The same expansion (eq. (34) could be written in terms of the square of mean curvature H. However, then we would predict supression of tip-splitting for both positive and negatively curved surfaces. We believe tip-splitting is enhanced for negative Gaussian curvature because the metric creates ample space for the split fingers to penetrate without mutual competition.

IV. CONCLUDING REMARKS

In this paper we generalized the traditional Saffman-Taylor problem by studying viscous flow on curved surfaces. Our main purpose was to investigate the influence of spatial curvature on viscous fingering pattern formation, when fluid flow takes place on a sphere. By deriving the equation of motion for the interface pertubation amplitudes, using a modecoupling approach, a study of both linear and weakly nonlinear stages of evolution could be carried out.

We have shown that cell curvature can be used as a control parameter to discipline splitting of the viscous fingers. The fluid-fluid interface can be more stable or unstable, with respect to tip-splitting, depending on the curvature of the surface to which the flow is confined. We also detected an asymmetry on tip-splitting behavior depending where the interface evolves: while tip-splitting may be still present on the northtern hemisphere, it is completely replaced by finger tip-sharpening on the southern hemisphere. We also found evidence that Hele-Shaw flows on negatively curved surfaces would present enhanced tendency to tip-splitting, so highly branched patterns may be expected. In summary, we have explicitly verified that interfacial behavior is coupled to the geometry of the Hele-Shaw cell, so that curvature has important consequences for flow dynamics.

ACKNOWLEDGMENTS

F.P., F.M. and J.A.M. thank CNPq and FINEP (through its PRONEX Program) for financial support. Work of M.W. was supported in part by the National Science Foundation grant No. DMR-9732567.

APPENDIX A: DARCY’S LAW ON A SPHERE

This appendix derives Darcy’s law for viscous flow between concentric spheres. The derivation is based upon results presented in Bird, et al. [21]. We begin with a coordinatefree representation of the continuity equation for an incompressible fluid

∇ · u = 0  (A1)

and the Navier-Stokes equation

ρ [\frac{∂u}{∂t} + (u·∇)u] = −∇p + η∇^{2}u  (A2)

where u denotes the three-dimensional fluid velocity and we neglect the acceleration due to gravity. Neglecting the inertial terms on the left-hand side of eq. (A2), transforming to spherical coordinates (r, θ, φ), and specializing to the case of polar flow (u_{r} = u_{φ} = 0) [21,22], we rewrite the Navier-Stokes equation

\frac{1}{r} \frac{∂p}{∂θ} = η (∇^{2}u_{θ} − \frac{u_{θ}}{r^{2} sin^{2} θ} ) .  (A3)

The noteworthy term in eq. (A3) is the second term multiplying viscosity which enters as a result of the curvilinear coordinate system.

To determine u_{θ}(r, θ) note that the continuity equation (A1) adapted for polar flow

\frac{1}{r sin θ} \frac{∂(u_{θ} sin θ)}{∂θ} = 0  (A4)

demands a solution of the form

u_{θ}(r, θ) = \frac{u(r)}{sin θ}.  (A5)

Insert this form into the Navier-Stokes equation (A3) and multiply by r sin θ to separate variables:

sin θ\frac{∂p}{∂θ} = \frac{η}{r} \frac{∂}{∂r} (r^{2} \frac{∂u}{∂r}) .  (A6)

Because the left-hand side involves only angle θ and the right-hand side is radial, involving only r, each side must be constant sharing a common value B. The solution of the radial equation subject to no-slip boundary conditions u(a) = u(a + b) = 0 is

u(r) = \frac{au_{0}}{b^{2}} \frac{(r − a)(a + b − r)}{r}  (A7)

where B = −2ηau_{0}/b^{2}.

Darcy’s law for Hele-Shaw cells is obtained by averaging the three-dimensional velocity u over the gap width b. We choose to average in the radial direction r and obtain v̄_{θ} = ū/ sin θ with

ū = \frac{u_{0}}{6} F (\frac{b}{a}) ≈ \frac{u_{0}}{6} (1 − \frac{b}{2a} + ⋯) ,  (A8)

where

F (\frac{b}{a}) = 6 (\frac{a}{b})^{3} [\frac{b}{a} (1 + \frac{b}{2a}) − (1 + \frac{b}{a}) log (1 + \frac{b}{a})] .  (A9)

Equation (A9) could also be written, in a more compact form, as a hypergeometric function F(b/a) = F (1, 2; 4; −b/a).

Equation (A8) generalizes the usual flat-cell velocity average. Darcy’s law becomes

v̄_{θ} = −\frac{b^{2}F (\frac{b}{a})}{12η} (\frac{1}{a} \frac{∂p}{∂θ}) .  (A10)

The effect of curvature can thus be incorporated entirely into a reduced gap width b or an enhanced viscosity η. Equation (A10) recovers the usual Darcy’s law for flow in flat Hele-Shaw cells with corrections of higher order in b/a. In this work we will be interested in the case b ≪ ζ ≪ ρ, with ρ ∼ a.

REFERENCES

[1] P. G. Saffman and G. I. Taylor, Proc. R. Soc. London Ser. A 245, 312 (1958).

[2] D. Bensimon, L. P. Kadanoff, S. Liang, B. I. Shraiman, and C. Tang, Rev. Mod. Phys. 58, 977 (1986); G. Homsy, Ann. Rev. Fluid Mech. 19, 271 (1987); K. V. McCloud and J. V. Maher, Phys. Rep. 260, 139 (1995).

[3] J. Bataille, Revue Inst. Pétrole 23, 1349 (1968).

[4] S. D. R. Wilson, J. Colloid Interface Sci. 51, 532 (1975).

[5] L. Paterson, J. Fluid Mech. 113, 513 (1981).

[6] J. A. Miranda and M. Widom, Int. J. Mod. Phys. B 12, 931 (1998).

[7] J. A. Miranda and M. Widom, Physica D 120, 315 (1998).

[8] H. Zhao and J. V. Maher, Phys. Rev. A 42, 5894 (1990).

[9] V. M. Entov and P. I. Etingof, Eur. J. Appl. Math. 8, 23 (1997).

[10] J. E. Avron and D. Levine, Phys. Rev. Lett. 69, 208 (1992).

[11] P. Peczak, G. S. Grest, and D. Levine, Phys. Rev. E 48, 4470 (1993).

[12] S. Hyde, S. Anderson, K. Larsson, Z. Blum, T. Landh, S. Lidin, and B. W. Ninham, The Language of Shape. The Role of Curvature in Condensed Matter: Physics, Chemistry and Biology (Elsevier, Amsterdam, 1997).

[13] V. S. Zykov, A. S. Mikhailov, and S. C. Müller, Phys. Rev. Lett. 78, 3398 (1997).

[14] C. Varea, J. L. Aragón, and R. A. Barrio, Phys. Rev. E 60, 4588 (1999).

[15] O. Schoenborn and R. Desai, J. Stat. Phys. 95, 949 (1999).

[16] B. A. Dubrovin, A. T. Fomenko, and S. P. Novikov, Modern Geometry - Methods and

Applications, Parts 1 and 2 (Springer-Verlag, New York, 1984).

[17] J. W. McLean and P. G. Saffman, J. Fluid Mech. 102, 455 (1981); P. G. Saffman, in Macroscopic Properties of Disordered Media, 1982, Lecture Notes in Physics, vol. 154, p. 208, edited by R. Burridge, S. Childress and G. Papanicolaou (Springer-Verlag, New York, 1982).

[18] C. -W. Park and G. M. Homsy, J. Fluid Mech. 139, 291 (1984); D. A. Reinelt, J. Fluid Mech. 183, 219 (1987).

[19] R. E. Rosensweig, Ferrohydrodynamics (Cambridge University Press, Cambridge, 1985).

[20] H. Thomé, M. Rabaud, V. Hakim and Y. Couder, Phys. Fluids A1, 224 (1989).

[21] R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena (Wiley, New York, 1960).

[22] G. H. Vatistas and W. Ghaly, Computers & Fluids 26, 683 (1997).

[23] S. S. S. Cardoso and A. W. Woods, J. Fluid Mech. 289, 351 (1995).

[24] M. J. P. Gingras and Z. Rácz, Phys. Rev. A 40, 5960 (1989).

[25] J.-D. Chen, J. Fluid Mech. 201, 223 (1989); J. -D. Chen, Exp. Fluids 5, 363 (1987).

[26] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products (Academic Press, New York, 1994).

FIGURE CAPTIONS

FIG. 1: Schematic configuration of the flow in a spherical Hele-Shaw cell.

FIG. 2: Time evolution of the fluids according to eq. (30), including modes 2 ≤ n ≤ 20. The initial perturbation amplitudes |ζ_{n}(0) = R_{0}/1000 and R_{0} = 0.05cm. Other parameters are given in the text. We show the interface at t = 10, 20 and 30.

FIG. 3: Nonlinear evolution. All physical parameters and initial conditions are the same as those used in Fig. 2. interfaces.

FIG. 4: Variation of T (2n, n) as a function of the spherical Hele-Shaw cell curvature K, for modes n = 8, 10, 12. The solid (dashed) curves describe behavior in the Northtern (Southern) hemisphere.