Some Numerical Methods for the Hele-Shaw Equations ^{1}

Nathaniel Whitaker Department of Mathematics and Statistics University of Massachusetts Amherst, MA 01003

Classification index numbers: 35Q30, 65D07, 65T10, 76C05, 76T05 Keywords: interface, fingering, surface tension, boundary-integral method, vortex method 18 pages, 4 figures, 2 tables Proposed running head: Numerical Methods for the Hele-Shaw Equations Please send proofs to

1This work was partially supported by the National Science Foundation under grant DMA-8913482

Prof. Nathaniel Whitaker

Department of Mathematics and Statistics Lederle Graduate Research Center University of Massachusetts Amherst, MA 01003

Abstract

In ([1, 2]), Tryggvason and Aref use a boundary integral method and the vortex-in-cell method to evolve the interface between two fluids in a Hele-Shaw cell. The method gives excellent results for intermediate values of the nondimensional surface tension parameter. The results are different from the predicted results of McLean and Saffman for small surface tension. For large surface tension, there are some numerical problems. In this paper, we implement the method of Tryggvason and Aref but use the point vortex method instead of the vortex-in-cell method. A parametric spline is used to represent the interface. The finger widths obtained agree well with those predicted by McLean and Saffman. We conclude that the method of Tryggvason and Aref can provide excellent results but that the vortex-in-cell method may not be the method of choice for extreme values of the surface tension parameter.

In a second method, we represent the interface with a Fourier representation. In addition, an alternative way of discretizing the boundary integral is used. Our results are compared to the linearized theory and the results of McLean and Saffman and are shown to be highly accurate.

1 Introduction

In this paper, we examine some numerical methods for evolving the interface between two fluids of different viscosities in a Hele-Shaw cell (two closely placed parallel plates). These fluids are immiscible in the sense that there is a finite surface tension which stabilizes small-scale disturbances at the interface. The equations describing the flow of a fluid in a Hele-Shaw cell are the equation of motion

u⃗ = −\frac{b^{2}}{12μ}∇p  (1)

and the equation of continuity

∇ · u⃗ = 0  (2)

where μ is the fluid viscosity, b is the spacing between the plates in the Hele-Shaw cell. The plates are taken as horizontal, the positive y-axis is in the direction of the flow, the x-axis is parallel to the plates and the z-axis is perpendicular to the plates. The velocity vector u⃗ has two components u and v which are functions of t (time), x and y. These velocities result from taking an average of the three-dimensional velocity field. The pressure p is an averaged pressure and the symbol ∇ represents the vector of partial derivatives. Derivations of the above equations are found in Lamb [3] and Bear [4].

For fluids of different viscosities, we have a set of equations to be solved on each side of the interface. These equations are connected by the following two conditions.

1. The normal component of the velocity is continuous at each point of the interface.

2. There is a prescribed jump in the pressure at each point (x_{0}, y_{0}) on the interface given by

\frac{σ}{R(x_{0}, y_{0})}  (3)

where σ is the surface tension parameter whose value depends on the two fluids involved and R(x_{0}, y_{0}) is the radius of curvature at the point (x_{0}, y_{0}) on the interface. We suppose that a fluid 1 is injected in the Hele-Shaw cell expelling some fluid 2 and that the parameters in fluids 1 and 2 are subscripted 1 and 2 respectively. In [1], it is shown that the solution to the above problem depends on a dimensionless surface tension parameter B given by

B = \frac{σk}{2AW L^{2}μ̄}  (4)

where k is b^{2}/12, (0, W) denotes the velocity of the fluid ahead or behind the interface, μ̄ is the average of the viscosities in the two fluids, L is the dimension of the cell in the x direction and A is defined by

A = \frac{μ_{2} − μ_{1}}{μ_{2} + μ_{1}} .  (5)

In a Hele-Shaw cell, for horizontal flow, an interface becomes unstable when a less viscous fluid displaces a more viscous fluid. Chuoke et al [5] performed experiments in a Hele-Shaw cell and packed bed models. They performed a linearized stability analysis including surface tension. Saffman and Taylor [6] also performed a linearized stability analysis and for zero surface tension, found an analytic solution for the shape of single fingers. McLean and Saffman [7] found an analytic shape for the fingers for nonzero values of the surface tension parameter.

Several numerical methods have been presented for the evolution of an interface in a Hele-Shaw cell. These include the methods of Tryggvason and Aref in [1, 2], Degregoria and Schwartz in [8, 9], Meiburg and Homsy in [10] and Whitaker [11]. Some of the most impressive results are due to the method of Meiberg and Homsy. Their predicted finger widths are extremely close to the predictions of McLean and Saffman. They also compare their results with an interface which is predicted to form a singularity in a finite time(see Aitchison and Howison [12]). The numerical results of Tryggvason and Aref seem to agree well with the predictions of the linearized theory. Their results also agree, in general, with the McLean and Saffman fingers. However, their method predicts a much narrower finger for small values of the surface tension parameter B. The numerical method also gives fingers slightly wider than the McLean and Saffman fingers for very large values of the surface tension parameter but it is believed that the finger may have not evolved fully. The method uses a boundary integral formulation along with Christiansen’s vortex-in-cell method [13].

In this paper, we present several modifications to the method of Tryggvason and Aref. In particular, a method is implemented using the boundary integral method presented by Tryggvason and Aref along with the point vortex method. In this method, the interface is represented by a parametric cubic spline. The purpose of this method is to compare the vortex-in-cell method and the point vortex method for extreme values of the nondimensional surface tension parameter. The point vortex method is known to be unstable when vortices lie too close to each other. It is shown here that it can be used to high accuracy for the Hele-Shaw equations with a modest number of vortices. The results of McLean and Saffman are compared with those produced by the vortex-in-cell method and the point vortex method. The point vortex method does much better for extreme values of the surface tension parameter. In addition, instabilities using the point vortex method are noticed only when time steps are taken too large. Chorin and Bernard [14] point out that the point vortex approximation will be reasonable when the velocity field of the point vortices is smoothed out and made bounded. In particular the vortex-in-cell method and the vortex blob [15] method are examples of methods which achieve this. In this problem, it appears that surface tension already stabilizes the vortex sheet(see Birkhoff([16]).

A second method is also presented which is based on interpolating the point vortices by a complex polynomial instead of a cubic spline and using an improved quadrature formula. This method produces excellent results when compared with the linearized theory and the fingers of McLean and Saffman. The results compare well with the results from the method of Meiberg and Homsy. The remainder of this paper is organized as follows. In section 2, the method of Tryggvason and Aref is described along with the vortex-in-cell method. In section 3, it is described how the point vortex method is implemented, and in section 4 the results from using the point vortex method are compared with the results from the method of Tryggvason and Aref and the McLean and Saffman fingers. In section 5, we describe the second method where the interface is represented by a complex polynomial and the boundary integral is discretized in a more accurate way. In section 6, the results from the presented methods are compared with the linearized theory. In section 7 the results from the second method are compared with the McLean and Saffman fingers.

2 Method of Tryggvason and Aref

In this method, the interface is represented by a vortex sheet. A vortex sheet in two dimensions is a curve along which vorticity, ω, is concentrated as a delta function. The vorticity is related to the velocity field by

ω = ∇ × u⃗.  (6)

Using equation (1), we see that ω = 0 on both sides of the interface. The sheet is characterized by a strength, γ, at each point of the interface which has the dimensions vorticity per unit area. Suppose that the fluids in the Hele-Shaw cell meet at an interface with unit tangent vector τ⃗. It can be shown then that the vortex sheet strength is exactly the jump in the tangential component of the velocity, i.e.,

γ = u⃗_{1} · τ⃗ − u⃗_{2} · τ⃗.  (7)

Equation (2) implies that there exists a stream function,ψ, such that

\frac{∂ψ}{∂x} = −v,  (8)

and

\frac{∂ψ}{∂y} = u.  (9)

Equations (6), (8) and (9) are then combined to give

△ψ = −ω.  (10)

The Green’s function associated with equation (10) can be used to formally write down the velocity field associated with ψ, for an arbitrary ω. For a singular vorticity concentrated along a curve, the velocity field for the sheet can be written in terms of γ, i.e.,

U⃗(s, t) = \frac{1}{2π} ∫ K(x⃗(s, t) − x⃗(ŝ, t))γ(ŝ, t)dŝ,  (11)

see Birkhoff [17]. The integral is taken over the sheet. x⃗ = (x, y) denotes a point on the interface being parametrized by time t and arclength s, and

K(x) = \frac{1}{2π} \frac{1}{|x⃗|^{2}} (−y, x).  (12)

Equation (11) could be used as an evolution equation for the interface assuming that γ(s, t) is known. This is the idea behind the point vortex method.

We now give the an integral equation for the appropriate γ for the Hele-Shaw equations. An equation for γ is derived by taking the dot product of (1) with τ⃗ in fluid 1 and fluid 2, subtracting one from the other. After solving for γ one obtains the equation

γ = \frac{Δμ}{μ̄} [\frac{1}{2}(u⃗_{1} + u⃗_{2})] · τ⃗ + \frac{b^{2}}{12μ̄}∇(Δp) · τ⃗,  (13)

where Δp = p_{2} − p_{1}, Δμ = μ_{2} − μ_{1} and μ̄ = \frac{1}{2}(μ⃗_{1} + μ⃗_{2}). It can be shown that the velocity at the sheet, given by (11), is the average of the limiting values of the velocities approaching from each side, i.e.,

U⃗ = \frac{1}{2}(u⃗_{1} + u⃗_{2}).  (14)

8 This fact is used and surface tension (3) is introduced in order to write (13) in the form

γ(s, t) = \frac{Δμ}{μ̄}U⃗ · τ⃗ + \frac{Δμ}{μ̄}W⃗ · τ⃗ + \frac{σb^{2}}{12μ̄} \frac{∂}{∂s} \frac{1}{R(s, t)},  (15)

where W⃗ is the velocity far ahead and far behind the interface. Equation (15) is an integral equation for the strength of the vortex sheet since U⃗ is a function of γ. The interface is discretized into a set of point vortices located at arclengths s_{i} with coordinates (x(s_{i}, t), y(s_{i}, t)) denoted (x_{i}, y_{i}) where i = 1, …, N. Equation (11) is replaced by a numerical integration formula and substituted into a discretized version of (15). γ is then solved for at each point vortex through iteration. The point vortices are then advected using a modified version of the vortex-in-cell method. The vortex-in-cell method is a method proposed by Christiansen [13] to find the velocity field resulting from some arbitrary distribution of vorticity. A grid is superposed over the (singular) vorticity distribution. Using an interpolation method presented by Meng and Thompson [18], the strength is used to approximate the vorticity at the grid points. The stream function and hence the velocity field is then obtained by solving (10), using a fast Poisson solver. The point vortices are then advected according to this velocity field. The interface is evolved by applying this algorithm iteratively. The interface is assumed periodic in the x-direction.

As the vortex sheet evolves it stretches unevenly which could lead to a large segment of the interface represented by only a few point vortices. This is addressed through a redistribution of the vortices after each time step. In this redistribution, the interface is assumed to be a curve consisting of piecewise linear segments connecting the point vortices. New points are then redistributed evenly onto this curve.

3 The Point Vortex Method

Equation (11) can be written in the complex form

\frac{dz(s, t)}{dt} = \frac{i}{2π} ∫_{−∞}^{∞} \frac{γ(ŝ, t)}{x(s, t) − x(ŝ, t) − i(y(s, t) − y(ŝ, t))}dŝ  (16)

where z = x + iy where i = \sqrt{(} − 1). As in the formulation of Tryggvason and Aref, let us assume that the interface is periodic in the x-direction. Without loss of generality, let us assume that this period is 1. It can then be written as

\frac{dz(s, t)}{dt} = −\frac{1}{2π} ∑_{j=−∞}^{∞} ∫_{0}^{S} \frac{γ(ŝ, t)}{x(s, t) − x(ŝ, t) − i(y(s, t) − y(ŝ, t)) − j} dŝ  (17)

where S is the total arclength taken over one period of x. We interchange the integral and summation, evaluate the complex sum using the well known formula

πcoth(πz) = ∑_{j=−∞}^{∞} \frac{1}{z − ij}  (18)

to obtain

\frac{dz(s, t)}{dt} = −\frac{1}{2} ∫_{0}^{S} γ(ŝ, t)coth(π(Δy + iΔx))dŝ  (19)

where Δx = x(s, t) − x(ŝ, t) and Δy = y(s, t) − y(ŝ, t). After separating real and imaginary parts we have

\frac{dx(s, t)}{dt} = −\frac{1}{2} ∫_{0}^{S} \frac{sinh(2πΔy)}{cosh(2πΔx) − cos(2πΔx)}γ(ŝ, t)dŝ  (20)

and

\frac{dy(s, t)}{dt} = \frac{1}{2} ∫_{0}^{S} \frac{sin(2πΔy)}{cosh(2πΔx) − cos(2πΔx)}γ(ŝ, t)dŝ  (21)

For 0 ≤ s ≤ S, our interface (x(s, t), y(s, t)) is discretized into a finite number of points, (x(s_{i}, t), y(s_{i}, t), denoted (x_{i}, y_{i}) for i = 1, …, N. We assume a uniform s-mesh with spacing Δs. Equations (20) and (21) then become

\frac{dx_{i}}{dt} = −\frac{S}{2N} ∑_{j=1,j≠i}^{N} \frac{sinh(2π(y_{i} − y_{j}))}{cosh(2π(x_{i} − x_{j})) − cos(2π(x_{i} − x_{j}))}γ_{j}Δs  (22)

and

\frac{dy_{i}}{dt} = \frac{S}{2N} ∑_{j=1,j≠i}^{N} \frac{sin(2π(y_{i} − y_{j}))}{cosh(2π(x_{i} − x_{j})) − cos(2π(x_{i} − x_{j}))}γ_{j}Δs.  (23)

The above discretization is the trapezoidal rule applied to the 2N ordinary differential equations (20) and (21). In the above, we use a uniform mesh only for convenience. Equations (22) and (23) are solved using the fourth order Runge-Kutta method. γ_{i} represents the strength of the vortex sheet at point (x_{i}, y_{i}) and is approximated using Tryggvason and Aref’s boundary integral method at each time step as described in the previous section. Because of (22) and (23), it would appear that the method for approximating the arclength between vortices would be important. In addition, one needs approximations for the tangent vector and the curvature at each vortex in the discretization of (15). We choose here to represent the interface by a parametric cubic spline to obtain approximations for arclength, the tangent vector and the curvature. The points along the interface are redistributed equally with respect to arclength at each time step. The interface is assumed to be symmetric about the centerline for this version of our code.

4 A Comparison of Results with the McLean and Saffman Fin-

gers

In this section, we compare results from using the point vortex method with those obtained by Tryggvason and Aref using the vortex-in-cell method and the McLean and Saffman solutions. The results of McLean and Saffman assume that A = 1, therefore all our computations shown here are for this value of A. This implies that the viscosity of the injected fluid is negligible.

Our interface is initially a cosine wave of small amplitude. This wave is discretized into particles and the numerical method is applied to evolve this wave which evolves into a long finger. By the linearized theory, there is a wavelength λ_{m} which grows faster than all other wavelengths.

Tryggvason and Aref plotted their results in terms of λ_{m} divided by λ_{i} where λ_{i} is the wavelength that an initially flat interface is perturbed by, i.e., the wavelength of the cosine wave. For fixed A, B is related to this ratio by

\frac{λ_{m}}{λ_{i}} = 2π\sqrt{(}3B).  (24)

In the remainder of this paper, \frac{λ_{m}}{λ_{i}} will be denoted by p.

In figure 1, three fingers are shown for the values p = .41,1.3 and 1.68. The fingers equilibrate faster for small values of p. This is why the fingers are shorter from left to right. The first finger to the left represents an extremely small value of p(p = .41). The middle finger represents a moderate value of p (p = 1.3). The last finger from the left represents an extremely large value of p. This value of p = 1.68 is extremely close to values of p where any initial perturbation would be linearly stable and would therefore not grow. The fingers develop a more pronounced neck as p increases. This is observed by other methods also.

In figure 2, we show the predicted finger widths of McLean and Saffman which correspond to the solid curve. The ”+” in figure 2 represents finger widths obtained by Tryggvason and Aref using their finest resolution. For intermediate values of p, Tryggvason and Aref’s finger widths are in excellent agreement with the results of McLean and Saffman, but for small values of p, the method produces a different finger width. Also for large values of p the fingers are slightly wider.

In figure 3, we compare the results of the point vortex method with the McLean and Saffman fingers. The symbols ”x” and ”o” in figure 3 correspond to 8 and 16 vortices per most unstable wavelength(λ_{m}) respectively. This number of vortices, per most unstable wavelength, is maintained as the interface evolves. The results clearly converge to the McLean and Saffman fingers. We do not observe a narrower finger for small values of p as Tryggvason and Aref conjecture. The method also seems to do a much better job for the large values of the surface tension parameter. The method becomes unstable if the time step is taken too large and the instability varies depending on p. The method becomes unstable in the sense that the wave begins to exhibit a chaotic behavior and the integral equation for the strength requires more and more iterations to converge.

The number of vortices used to compute the equilibrated fingers for different values of p ranged from 178 to 322, using 16 vortices per most unstable wavelength. The 322 corresponds to the smallest value of p in figure 3 and 178 corresponds to the largest value of the p given in figure 3. The fingers grown from smaller surface tensions initially need more vortices since λ_{m} is smaller but these interfaces settle into equilibrated fingers earlier. The fingers grown from larger surface tension require less vortices initially and require longer to settle down.

5 Some Modifications to the Method

The previous method will be denoted method I and the method described below will be referred to later as method II. Suppose that the interface is represented by the N point vortices (x_{0}, y_{0}), …, (x_{N}, y_{N}). The interface is assumed to be periodic in the x-direction with period κ where x_{N} = x_{0} + κ and y_{0} = y_{N}. A parametric curve (x(r), y(r)) is sought, parametrized by 0 ≤ r ≤ 1, which interpolates the points (x_{0}, y_{0}), …, (x_{N}, y_{N}). The curve is sought in complex form, i.e., x(r) + iy(r) where i = \sqrt{−1}. In principle, we want a phase polynomial of the form

∑_{j=0}^{N−1} β_{j}e^{ijr}  (25)

which represents the periodic part of x(r) + iy(r), i.e., x(r) − rκ + iy(r). We seek β_{0}, …, β_{N−1} such that (25) interpolates the data (0, x_{0} + iy_{0}), (1/N, x_{1} − (κ/N) + iy_{1}), …, (1, x_{N} − κ + iy_{N}).

The interface, x(r) + iy(r), can then be written as

rκ + ∑_{j=0}^{N−1} β_{j}e^{ijr}.  (26)

Krasny [19] uses this representation to implement filtering. The phase polynomial is used to compute the tangent vector given by

x′(r) + iy′(r) = κ + ∑_{j=0}^{N−1} ijβ_{j}e^{ijr}.  (27)

Higher derivatives and curvature are also computed from the phase polynomial by differentiation. The quantities are easily computed explicitly at the data points using the Fast Fourier Transform.

Our curve is parametrized according to the parameter r given through the Fourier representation above. Equation (17) can be written according to this reparametrization, in the form,

\frac{dz(r, t)}{dt} = −\frac{1}{2π} ∑_{j=−∞}^{∞} ∫_{0}^{1} \frac{γ(r̂, t)}{x(r, t) − x(r̂, t) − i(y(r, t) − y(r̂, t)) − jκ} \frac{ds}{dr̂}dr̂  (28)
= −\frac{1}{2π} ∑_{j=−∞}^{∞} ∫_{0}^{1} \frac{γ(r̂, t)}{x(r, t) − x(r̂, t) − i(y(r, t) − y(r̂, t)) − jκ} \frac{ds}{dr̂}dr̂  (29)

where \frac{ds}{dr̂} is the magnitude of the tangent vector, |x′(r̂) + iy′(r̂)|. This expression, as before, simplifies into

\frac{dx(r, t)}{dt} = −\frac{1}{2} ∫_{0}^{1} \frac{sinh(2πΔy)}{cosh(2πΔx) − cos(2πΔx)}γ(r̂, t)\frac{ds}{dr̂}dr̂  (30)

and

\frac{dy(r, t)}{dt} = \frac{1}{2} ∫_{0}^{1} \frac{sin(2πΔy)}{cosh(2πΔx) − cos(2πΔx)}γ(r̂, t)\frac{ds}{dr̂}dr̂  (31)

where Δx = x(r, t) − x(r̂, t) and Δy = y(r, t) − y(r̂, t).

In the discretization of equations (20) and (21), the trapezoidal rule is used with a point excluded where the integral is singular. Since the grid size is Δs, the approximation will have error O(Δs). For method II, equations (30) and (31) are approximated using an idea presented by Sidi and Israeli [20] and used by Shelley [21]. The basic idea is to use Richardson’s extrapolation to eliminate the O(Δs) term from the error in the trapezoidal rule. This results in the formulas

\frac{dx_{k}}{dt} = −\frac{1}{2N} N−1 ; ∑ ; j = 0_{(j+k)odd} \frac{sinh(2π(y_{k} − y_{j}))}{cosh(2π(x_{k} − x_{j})) − cos(2π(x_{k} − x_{j}))}γ_{j}|x′_{j} + iy′_{j}|Δr  (32)

and

\frac{dy_{k}}{dt} = \frac{1}{2N} N−1 ; ∑ ; j = 0_{(j+k)odd} \frac{sin(2π(y_{k} − y_{j}))}{cosh(2π(x_{k} − x_{j})) − cos(2π(x_{k} − x_{j}))}γ_{j}|x′_{j} + iy′_{j}|Δr.  (33)

The error in using these formulas is O((Δr)^{2}). The vortices are redistributed evenly after every time step using a parametric cubic spline.

As it was mentioned earlier, the integral equation given by (15) is solved by iteration in method I. This is equivalent to using Jacobi’s method for a matrix system. For method II, the quadrature formula for the integral is changed. The Jacobi type iteration does not converge, in general, for method II, due to the new discretization. It is probably possible to find a similar inexpensive iterative algorithm to compute the strength of the sheet. This, however, is not our focus so we simply use gaussian elimination.

6 Linearized Results

In [5] and in [6], it is shown that when the Hele-Shaw equations are satisfied on both sides of an interface, given by the equation

y = −ϵcos(nx),  (34)

then the interface, after a short time Δt, evolves into an interface given by the equation

y = ϵe^{ρΔt}cos(nx).  (35)

The analysis assumes that ϵe^{ρΔt} is small enough. The rate of growth ρ is given by

ρ = \frac{(μ_{2} − μ_{1})W n − τkn^{3}}{μ_{2} + μ_{1}} .  (36)

n is the wave number, ϵ is a small parameter, τ is the surface tension coefficient, (0, W) is the velocity far ahead and far behind the interface and k = \frac{b^{2}}{12}.

Our interface is initially a wave of the form (34). It is evolved numerically over a time step, Δt, and the numerical ρ is computed using a method given in [11]. This is done for methods I and II and the results are compared with the exact ρ given by equation (36). For A = 1, these results are presented in tables I and II. 5 vortices are used to represent the entire interface for all results in method II. 102 vortices are used for the results of method I. Rates of growth ρ are shown for values of B = 0.0, 1.0 and 1.68. The results in table I are for Δt = .000001 and ϵ = .001. The results in table II are for Δt = .000001, and ϵ = .001 and .000001. The results of method II are excellent with only 5 vortices. In method II for ϵ = .001, we get at least 3 digits of accuracy. Using ϵ = .000001, we get at least 10 digits of accuracy.

7 Comparison with the McLean-Saffman fingers

In this section, the results of the method II are compared with the predicted finger widths of McLean and Saffman [7]. In figure 4, we compare the results of method II with the McLean and Saffman fingers. The symbols ”x” and ”o” in figure 4 correspond to 4 and 8 vortices per most unstable wavelength(λ_{m}) respectively. As in method I, this number of vortices, per most unstable wavelength, is maintained as the interface evolves. The results clearly converge to the McLean and Saffman fingers. The method also seems to do an excellent job for the extreme values of the surface tension parameter used. The results are comparable with those of Meiberg and Homsy.

The number of vortices used to compute the equilibrated fingers for different values of p ranged from 57 to 121, using 8 vortices per most unstable wavelength. The 121 corresponds to the smallest value of p in figure 4 and 57 corresponds to the largest value of the p given in figure 4. The fingers grown from smaller surface tensions initially need more vortices since λ_{m} is smaller but these interfaces settle into equilibrated fingers earlier. The fingers grown from larger surface tension require less vortices initially and require longer to settle down.

The finger widths were approximated by 1 divided by the speed of the finger as in [9] and in [10]. This value oscillates slightly as the finger evolves(see [10]), so the width shown in figure 4 represents an average value. This variance in the oscillation is greater as the nondimensional surface tension parameter increases.

Acknowledgements. I would like to thank Professors Tom Buttke, Gretar Tryggvason and Robert Krasny for several helpful discussions. These computations were done at the GANG lab at the University of Massachusetts at Amherst.

References

[1] G. Tryggvason and H. Aref, J. Fluid Mech. 136, 1 (1983).

[2] G. Tryggvason and H. Aref, J. Fluid Mech. 154, 287 (1985).

[3] H. Lamb, Hydrodynamics, Dover, New York, 1945.

[4] J. Bear, Dynamics of Fluids in Porous Media, Elsevier, New York, 1972.

[5] R. Chuoke, P. van Meurs, and C. van der Poel, Trans. Amer. Inst. Min. Eng. AIME-216, 188 (1959).

[6] P. Saffman and G. Taylor, Proc. Roy. Soc. London 245, 312 (1958).

[7] J. McLean and P. Saffman, J. Fluid Mech. 102, 455 (1981).

[8] A. Degregoria and L. Schwartz, Phys. Fluids 28, 2313 (1985).

[9] A. Degregoria and L. Schwartz, J. Fluid Mech. 164, 383 (1986).

[10] E. Meiburg and G. Homsy, Phys. Fluids 31, 429 (1988).

[11] N. Whitaker, J. Comput. Phys. 90, 176 (1990).

[12] J. Aitchison and S. Howison, J. Comput. Phys. 60, 376 (1985).

[13] J. P. Christiansen, J. Comput. Phys. 13, 363 (1973).

[14] A. J. Chorin and P. Bernard, J. Comput. Phys. 13, 423 (1973).

[15] A. J. Chorin, J. Fluid Mech. 57, 785 (1973).

[16] G. Birkhoff, Helmholtz and Taylor instability, in Proceedings of Symposia in Applied Mathematics, Vol. XIII, edited by G. Birkhoff, R. Bellman, and C. Lin, p. 55, American Mathematical Society, Providence, 1962.

[17] G. Birkhoff, Taylor instablity and laminar mixing, Technical Report 1862, Los Alamos National Laboratory, Los Alamos, N.M., 1954.

[18] J. Meng and J. Thomson, J. Fluid Mech. 84, 433 (1978).

[19] R. Krasny, J. Comput. Phys. 65, 292 (1986).

[20] A. Sidi and M. Israeli, J. Sci. Comp. 3, 201 (1988).

[21] M. J. Shelley, J. Fluid Mech. 244, 493 (1992).

Figure 1. Three evolved fingers for the nondimensional surface tension parameter p = .41, 1.3 and 1.68. Figure 2. Finger widths from the method of Tryggvason and Aref(+) versus predicted widths from McLean and Saffman(solid line). Figure 3. Finger widths from the point vortex method(Method I) versus the predicted widths from McLean and Saffman(solid line). The symbols ”x” and ”o” correspond to 8 and 16 vortices per unstable wavelength respectively. Figure 4. Finger widths from the modified point vortex method(Method II) versus the predicted widths from McLean and Saffman(solid line). The symbols ”x” and ”o” correspond to 4 and 8 vortices per unstable wavelength respectively.