Numerical Study of Free Interface Problems Using

Boundary Integral Methods

Thomas Yizhao Hou

Abstract. Numerical study of fluid interfaces is a difficult task due to the presence of high frequency numerical instabilities. Small perturbations even at the round-off error level may experience rapid growth. This makes it very difficult to distinguish the numerical instability from the physical one. Here, we perform a careful numerical stability analysis for both the spatial and time discretization. We found that there is a compatibility condition between the numerical discretizations of the singular integral operators and of the Lagrangian derivative operator. Violation of this compatibility condition will lead to numerical instability. We completely eliminate the numerical instability by enforcing this discrete compatibility condition. The resulting scheme is shown to be stable and convergent in both two and three dimensions. The improved method enables us to perform a careful numerical study of the stabilizing effect of surface tension for fluid interfaces. Several interesting phenomena have been observed. Numerical results will be presented.

1991 Mathematics Subject Classification: 65M12, 76B15. Keywords and Phrases: Free boundary, numerical stability, surface tension, topological singularity.

1 Introduction.

Many physically interesting problems involve propagation of free interfaces. Water waves, boundaries between immiscible fluids, vortex sheets, Hele-Shaw cells, thin-film growth, crystal growth and solidification are some of the better known examples. Numerical simulations for interfacial flows play an increasingly important role in understanding the complex interfacial dynamics, pattern formations, and interfacial instabilities. Many numerical methods have been developed to study these interfacial problems, including phase field models, volume-of-fluid methods, level set methods, front tracking methods, and boundary integral/element methods. Here we will focus on boundary integral methods.

Numerical study of fluid interfaces is a difficult task due to the presence of high frequency numerical instabilities [6, 7]. Small perturbations even at the round-off error level may experience rapid growth. This makes it very difficult to distinguish the numerical instability from the physical one. In our study, we first establish the well-posedness of the linearized motion far from equilibrium. This involves careful analysis of singular integral operators defined on free interfaces [2, 11]. The continuous well-posedness analysis provides a critical guideline for our numerical analysis of the discrete system. We found that there is a corresponding compatibility condition between the discrete singular operators and the discrete derivative operator. Violation of this compatibility condition will lead to numerical instability. We completely eliminate the numerical instability by introducing an effective filtering. The amount of filtering is determined by enforcing this discrete compatibility condition. The resulting scheme is shown to be stable and convergent [3]. The corresponding 3-D problem is considerably more difficult since the singular operators have non-removable branch point singularities and there is no spectrally accurate discretization. A new stabilizing technique is introduced to overcome this difficulty [12, 13]. This technique is very general and effective. It also applies to non-periodic problems with rigid boundaries.

The improved method enables us to perform a careful numerical study of the stabilizing effect of surface tension for fluid interfaces. Water waves with small surface tension are shown to form singular capillary waves dynamically. The mechanism for generating such capillary waves is revealed and the zero surface tension limit is investigated [5]. In another study, surface tension is shown to regularize the early curvature singularity induced by the Rayleigh-Taylor instability in an unstably stratified two-fluid interface. However, a pinching singularity is observed in the late stage of the roll-up. The interface forms a trapped bubble and self-intersects in finite time [4, 9].

2 Stability of Boundary Integral Methods for 2-D Water Waves

In this section, we consider the stability of boundary integral methods for 2-D water waves. The result can be generalized to two-fluid interfaces and Hele-Shaw flows [4, 8]. Consider a 2-D incompressible, inviscid and irrotational fluid below a free interface. We assume the interface is 2π-periodic in the horizontal direction and parametrize the interface by a complex variable, z(α, t) = x(α, t) + iy(α, t), where α is a Lagrangian parameter along the interface. We use the usual convention of choosing the tangential velocity to be that of the fluid. The first boundary integral method for water waves was proposed by Longuet-Higgins and Cokelet [15] who used a single layer representation. Here we will use a double layer representation introduced by Baker-Meiron-Orszag [1]. Following [1], we obtain a system of evolution equations as follows:

z̄_{t} = \frac{1}{4πi} ∫_{−π}^{π} γ(α′)cot(\frac{z(α) − z(α′)}{2})dα′ + \frac{γ(α)}{2z_{α}(α)} ≡ u − iv,  (1)
ϕ_{t} = \frac{1}{2}(u^{2} + v^{2}) − gy ,  (2)
ϕ_{α} = \frac{γ}{2} + Re (\frac{z_{α}}{4πi} ∫_{−π}^{π} γ(α′)cot(\frac{z(α) − z(α′)}{2})dα′) ,  (3)

where ϕ is the potential, γ is the vortex sheet strength, z̄ is the complex conjugate of z. Equations (1)-(3) completely determine the motion of the system. The advantage of using the double layer representation is that the Fredholm integral equation of second kind has a global convergent Neumann series [1]. Thus γ can be solved by fixed point iteration.

The boundary integral formulation of water waves is naturally suited for numerical computation. There are many ways one can discretize the boundary integral equations, depending on how we choose to discretize the singular integrals and the derivatives. These choices affect critically the accuracy and stability of the numerical method. Straightforward numerical discretizations of (1)-(3) may lead to rapid growth in the high wavenumbers. In order to avoid numerical instability, a certain compatibility between the choice of quadrature rule for the singular integral and the discrete derivatives must be satisfied. This compatibility ensures that a delicate balance of terms on the continuous level is preserved on the discrete level. Violation of this compatibility will lead to numerical instability.

Let z_{j}(t) be the numerical approximation of z(α_{j}, t), where α_{j} = jh, h = 2π/N. ϕ_{j}(t), γ_{j}(t) are defined similarly. To approximate the velocity integral, we use the alternating trapezoidal rule:

∫_{−π}^{π} γ(α′) cot (\frac{z(α_{j}) − z(α′)}{2})dα′ ≃ ∑_{\begin{matrix}j=−N/2+1 \\ (j − i) odd \end{matrix}}^{N/2} γ_{k} cot(\frac{z_{j} − z_{k}}{2} )2h .  (4)

The advantage of using this alternating trapezoidal quadrature is that the approximation is spectrally accurate. We denote by D_{h} the discrete derivative operator. In general, we have ˆ(D_{h})_{k} = ikρ(kh) for some nonnegative even function ρ. The specific form of ρ(ξ) depends on the approximation. For example, we have ρ_{c}(kh) = 3 sin(kh)/(kh(2 + cos(kh))) for the cubic spline approximation, and ρ(kh) = 1 for a pseudo-spectral derivative.

Now we can present our numerical algorithm for the water wave equations (1)-(3) as follows:

\frac{dz̄_{j}}{dt} = \frac{1}{4πi} ∑_{(k−j)odd} γ_{k}cot(\frac{z_{j}^{(ρ)} − z_{k}^{(ρ)}}{2} )2h + \frac{γ_{j}}{2D_{h}z_{j}} ≡ u_{j} − iv_{j},  (5)

\frac{dϕ_{j}}{dt} = \frac{1}{2}(u_{j}^{2} + v_{j}^{2}) − gy_{j} ,  (6)

D_{h}ϕ_{j} = \frac{γ_{j}}{2} + Re (\frac{D_{h}z_{j}}{4πi} ∑_{(k−j)odd} γ_{k}cot(\frac{z_{j}^{(ρ)} − z_{k}^{(ρ)}}{2} )2h) ,  (7)

where z^{(ρ)} is a Fourier filtering defined as ˆ(z^{(ρ)})_{k} = ẑ_{k}ρ(kh). The Fourier filtering z^{(ρ)} in (5) and (7) is to balance the high frequency errors introduced by D_{h}. This will become apparent in the discussion of stability below.

Theorem 1. Assume that the water wave problem is well-posed and has a smooth solution in C^{m+2} (m ≥ 3) up to time T. Then if D_{h} corresponds to a r-th order derivative approximation, we have for 0 < h ≤ h_{0}(T)

∥z(t) − z(·, t)∥_{l^{2}} ≤ C(T)h^{r} .  (8)

Similar convergent results hold for ϕ_{j} and γ_{j}. Here ∥z∥_{l^{2}}^{2} = ∑_{j=1}^{N} |z_{j}|^{2}h .

2.1 Discussion of stability analysis

Here we discuss some of the main ingredients in the stability analysis of the scheme given by (5)-(7). We will mainly focus on the linear stability. Once linear stability is established, nonlinear stability can be obtained relatively easily by using the smallness of the error and an induction argument. The reader is referred to [3] for details.

To analyze linear stability, we first derive evolution equations for the errors ż_{j}(t) ≡ z_{j}(t)−z(α_{j}, t), etc., and try to estimate their growth in time. If we take the difference between the sum in (5) for the discrete velocity and the corresponding sum for the exact solution, the linear terms in ż_{j}, γ̇_{j} for the difference are

\frac{h}{πi} ∑_{(k−j)odd} \frac{γ̇_{k}}{z(α_{j})^{(ρ)} − z(α_{k})^{(ρ)}} − \frac{h}{πi} ∑_{(k−j)odd} \frac{γ(α_{k})(ż_{j}^{(ρ)} − ż_{k}^{(ρ)})}{(z(α_{j})^{(ρ)} − z(α_{k})^{(ρ)})^{2}} ,  (9)

where we have expanded the periodic sum, with k now unbounded. To identify the most singular terms, we use the Taylor expansion to obtain the most singular symbols

\frac{1}{z(α_{j}) − z(α_{k})} = \frac{1}{z_{α}(α_{j})(α_{j} − α_{k})} + f(α_{j}, α_{k}) ,

where f is a smooth function. Thus, the most important contribution to the first term in (9) is (2iz_{α})^{−1}H_{h}γ̇_{j}, where H_{h} is the discrete Hilbert transform

H_{h}(γ̇_{j}) ≡ \frac{1}{π} ∑_{(k − j)odd} \frac{γ̇_{k}}{α_{j} − α_{k}} 2h .  (10)

Similarly, the most important contribution to the second term in (9) is −γ(2iz_{α}^{2})^{−1}Λ_{h}(z_{j}^{(ρ)}), where Λ_{h} is defined as follows:

Λ_{h}(ḟ_{j}) ≡ \frac{1}{π} ∑_{(k − j)odd} \frac{ḟ_{j} − ḟ_{k}}{(α_{j} − α_{k})^{2}} 2h .  (11)

Let H and Λ be the corresponding continuous operators for H_{h} and Λ_{h} respectively, i.e. replacing the discrete sums by the continuous integrals. In the continuous level, it is easy to show by integration by parts that

Λ(f) = H(D_{α}f) ,  (12)

where D_{α} is the continuous derivative operator. It turns out that in order to maintain numerical stability of the boundary integral method, the quadrature rule for the singular integral and the discrete derivative operator D_{h} must satisfy a compatibility condition similar to (12). That is, given a quadrature rule, which defines a corresponding discrete operators H_{h} and Λ_{h}, and a discrete derivative D_{h}, they must satisfy the following compatibility condition:

Λ_{h}(ż_{i}) = H_{h}D_{h}(ż_{i}) ,  (13)

for ż satisfying ˆ ż_{0} = ˆ ż_{N/2} = 0. If (13) is violated, it would cause a mismatch of a singular operator of the form (Λ_{h} − H_{h}D_{h})(ż) in the error equations. This will generate numerical instability.

By performing appropriate Fourier filtering in the approximations of the velocity integral, we can ensure a variant of the compatibility condition (13) is satisfied,

Λ_{h}(ż_{j}^{(ρ)}) = H_{h}D_{h}(ż_{j}) .  (14)

This can be verified from the spectrum properties of H_{h} and Λ_{h} and the definition of the ρ filtering. This modified compatibility condition is sufficient to ensure stability of our modified boundary integral method. This explains why we need to filter z in (5) and (7) when we approximate the velocity integral. The modified algorithm also allows use of non-spectral derivative operators.

By using properties of the discrete Hilbert transform: (i) H_{h}^{2} = −I , (ii) Λ_{h}(z^{(ρ)}) = H_{h}D_{h}(z), (iii) the commutator, [H_{h}, f], is a smoothing operator, i.e. [H_{h}, f](ż^{(ρ)}) = A_{−1}(ż) for smooth f, we can derive an error equation for ż_{j} which is similar to the continuum counterpart in the linear well-posedness study [2, 3]

\frac{dż_{j}}{dt} = z_{α}^{−1}(I − iH_{h})D_{h}Ḟ + A_{0}(ż) + A_{−1}(ϕ̇) + O(h^{r}),

where Ḟ = ϕ̇ − uẋ − vẏ, A_{0} is a bounded operator from l^{p} to l^{p}, and A_{−1} is a smoothing operator of order one, i.e. D_{h}A_{−1} = A_{0} and A_{−1}D_{h} = A_{0}. The leading order error equation suggests that we project the error equation into the local tangential and normal coordinate system. In this local coordinate, the stability property of the error equations becomes apparent. Let ż^{N}, ż^{T} be the normal and tangential components of ż, with respect to the underlying curve z(α), N being the outward unit normal, and δ̇ = ż^{T} +H_{h}ż^{N}. We obtain after some simplification

δ̇_{t} = A_{−1}(Ḟ) + A_{0}(ż),  (15)
ż_{t}^{N} = \frac{1}{|z_{α}|}H_{h}D_{h}Ḟ + A_{−1}(Ḟ) + A_{0}(ż),  (16)
Ḟ_{t} = −c(α, t)ż^{N} + A_{−1}(ż),  c(α, t) = (u_{t}, v_{t} + g) · N,  (17)

where equation (17) is obtained by performing error analysis on Bernoulli’s equation and using the Euler equations. In this form it is clear that only the normal component of ż is important. This is consistent with the physical property of interfacial dynamics. Now it is a trivial matter to establish an energy estimate for the error equations. Note that H_{h}D_{h} is a positive operator with a Fourier symbol ρ(kh)|k|. The discretization is stable if the water wave problem is well-posed, i.e. the sign condition, c(α, t) > 0, is satisfied. We refer to [3] for details.

3 Generalization to 3-D Water Waves.

Numerical stability of 3-D boundary integral methods is much more difficult. Let z(α_{1}, α_{2}) be a parametrization of a 3-D surface. Recall that the 3-D free space Green function for the Laplace equation is given by G(z) = −1/(4π|z|). The corresponding velocity integral is given by

\frac{dz}{dt} = ∫ Ω(α′) × ∇_{z′}G(z(α) − z(α′))dα′ + w_{loc}(α) ,

where Ω = μ_{α_{1}}z_{α_{2}} − μ_{α_{2}}z_{α_{1}}, w_{loc} = Ω(α) × (z_{α_{1}} × z_{α_{2}}) /|z_{α_{1}} × z_{α_{2}}|^{2} is the local velocity, and μ is the dipole strength.

One of the main difficulties for 3-D boundary integral methods is that the velocity integral has a branch point singularity which is not removable by desingularization. Also, unlike the 2-D case, we cannot express the leading order contribution of the singular operator as an integral operator defined on a flat surface. Now, the leading order singular operators depend on the free surface and have variable coefficients (assuming the tangent vectors are orthogonal for simplicity):

H_{l}(f) = \frac{1}{2π} ∫ \frac{(α_{l} − α′_{l})f(α′)dα′}{(|z_{α_{1}}(α)|^{2}(α_{1} − α′_{1})^{2} + |z_{α_{2}}(α)|^{2}(α_{2} − α′_{2})^{2})^{3/2}} , l = 1, 2,

Λ(f) = \frac{1}{2π} ∫ \frac{(f(α) − f(α′))dα′}{(|z_{α_{1}}(α)|^{2}(α_{1} − α′_{1})^{2} + |z_{α_{2}}(α)|^{2}(α_{2} − α′_{2})^{2})^{3/2}} .

As in 2-D, there are certain compatibility conditions among singular operators and the derivative operator. For example, we have Λ = H_{1}D_{1} +H_{2}D_{2}. Stability of the boundary integral method requires a similar compatibility condition to hold:

Λ_{h}(z) = (H_{1}^{h}D_{1}^{h} + H_{2}^{h}D_{2}^{h})z,

which, unfortunately, is generically violated by almost all discretizations. Although this compatibility condition can be imposed by applying a Fourier filtering as in the 2-D case, such filtering can no longer be evaluated efficiently by Fast Fourier Transform (FFT) since the singular operator, H_{l}^{h} or Λ_{h}, is not a convolution operator. The kernel depends on a variable coefficient.

In additional to the above compatibility condition, there are several other compatibility conditions that need to be satisfied for 3-D surfaces. Since there are no spectrally accurate approximations to the singular integrals in 3-D, it is almost impossible to enforce all the other compatibility conditions by using Fourier filtering alone.

To overcome this difficulty, we introduce a new stabilizing method without using the Fourier filtering. This new technique can be illustrated more clearly for the 2-D point vortex method [12]. Let us illustrate how we enforce the compatibility condition Λ_{h} = H_{h}D_{h} indirectly by adding a stabilizing term. The modified point vortex method approximation for 2-D water waves is given by

\frac{dz̄_{j}}{dt} = \frac{1}{2πi} ∑_{k≠j} \frac{γ_{k}h}{z_{j} − z_{k}} + \frac{γ_{j}}{2D_{h}z_{j}} + C_{j}^{I},

where

C_{j}^{I} = \frac{γ_{j}}{2i(D_{h}z_{j})^{2}} (Λ_{h} − H_{h}D_{h})z_{j}.

This method is clearly consistent since the point vortex method gives a first order approximation to the singular integral: (H_{h}D_{h} − Λ_{h})z(α_{j}) = O(h). Let ż_{j} = z_{j} − z(α_{j}), γ̇_{j} = γ_{j} − γ(α_{j}) be the errors in z_{j} and γ_{j}. Let E_{i} = \frac{1}{2πi} ∑_{j≠i} \frac{γ_{j}}{z_{i}−z_{j}} h,

and define Ė_{i} = E_{i} − E(α_{i}), Ċ_{i}^{I} = C_{i}^{I} − C^{I}(α_{i}). Using the same argument as before, we can show that the linear variation in E_{i} is given by

Ė_{i} = \frac{1}{2iz_{α}(α_{i})}H_{h}(γ̇_{i}) − \frac{γ(α_{i})}{2iz_{α}(α_{i})^{2}} Λ_{h}ż_{i} + A_{0}(ż_{i}) + A_{−1}(γ̇_{i}) .

Similarly, we have

Ċ_{i}^{I} = \frac{γ(α_{i})}{2iz_{α}(α_{i})^{2}} (Λ_{h} − H_{h}D_{h})ż_{i} + A_{0}(ż_{i}) + A_{−1}(γ̇_{i}),

where we have used the fact that (Λ_{h} − H_{h}D_{h})z(α_{i}) = O(h). Now combining Ė_{i} with Ċ_{i}^{I}, we obtain

Ė_{i} + Ċ_{i}^{I} = \frac{1}{2iz_{α}(α_{i})}H_{h}(γ̇_{i}) − \frac{γ(α_{i})}{2iz_{α}(α_{i})^{2}} Λ_{h}ż_{i}

+\frac{γ(α_{i})}{2iz_{α}(α_{i})^{2}} (Λ_{h} − H_{h}D_{h})ż_{i} + A_{0}(ż_{i}) + A_{−1}(γ̇_{i})

= \frac{1}{2iz_{α}(α_{i})}H_{h}(γ̇_{i}) − \frac{γ(α_{i})}{2iz_{α}(α_{i})^{2}} H_{h}D_{h}ż_{i} + A_{0}(ż_{i}) + A_{−1}(γ̇_{i}).

Note that the two Λ_{h}ż_{i} terms cancel each other in the above equation, and only the H_{h}D_{h}ż_{i} term survives in place of Λ_{h}ż_{i}. This in effect enforces the compatibility condition Λ_{h} = H_{h}D_{h}. This stabilizing technique is very general, and it applies to 3-D water waves. For 3-D water waves, we have four more compatibility conditions that need to be satisfied. We need to handle each one of them by adding a corresponding stabilizing term just as we outlined above. This will give a stable discretization for 3-D water waves. Moreover, by using a generalized arclength frame which enforces |z_{α_{1}}|^{2} = λ_{1}(t)|z_{α_{2}}|^{2} and (z_{α_{1}}, z_{α_{2}}) = λ_{2}(t)|z_{α_{2}}|^{2}, these correction terms can be evaluated efficiently by FFT, see [13].

4 Stabilizing Effect of Surface Tension

Surface tension plays an important role in understanding fluid phenomena such as pattern formation in Hele-Shaw cells, the motion of capillary waves on free surfaces, and the formation of fluid droplets. On the other hand, surface tension also introduces high order spatial derivatives into the interface motion through local curvature which couples to the interface equation in a nonlinear and nonlocal manner. These terms induce strong stability constraints on the time step if an explicit time integration method is used. These stability constraints are generally time dependent, and become more severe by the differential clustering of points along the interface.

Hou, Lowengrub, and Shelley [8] proposed to remove the stiffness of surface tension for 2-D fluid interfaces by using the Small Scale Decomposition technique and reformulating the problem in the tangent angle θ and arclength metric s_{α}. Curvature has a very simple expression in these variables, κ = θ_{α}/s_{α}. One important observation is that the stiffness only enters at small scales. The leading order contribution of these singular operators at small scales can be expressed in terms of the Hilbert transform, which is diagonalizable using Fourier Transform. By treating the leading order terms implicitly, but treating the lower order terms explicitly, we obtain a semi-implicit discretization which can be inverted efficiently using FFT. This reformulation greatly improves the time step stability constraint. Many interfacial problems that were previously not amenable are now solvable using this method. This idea has been subsequently generalized to 3-D filaments by Hou, Klapper, and Si who use curvature and arclength metric as the new dynamic variables [10]. Applications to Hele-Shaw flows, 3-D vortex filaments, and the Kirchhoff rod model for protein folding all give very impressive results.

In the following, we would like to present two numerical calculations using our numerical methods. In Fig. 1, we show that water waves with small surface tension generate singular capillary waves dynamically. Our study shows that the dynamic generation of capillary waves is a result of the competition between convection and dispersion. The capillary waves originate near the crest in a neighborhood where both the curvature and its derivative are maximum. For fixed but small surface tension, the maximum of curvature increases in time and the interface develops oscillatory capillary waves in the forward front of the crest. The minimum distance between adjacent capillary crests appears to approach zero, suggesting the formation of trapped bubbles as observed in Koga’s experiments of breaking waves [14]. On the other hand, for a fixed time, as the surface tension coefficient τ is reduced, both the capillary wavelength and its amplitude decreases nonlinearly. The interface converges strongly to the zero surface tension profile [5].

We study the stabilizing effect of surface tension for an unstably stratified two-fluid interface in Fig. 2. This problem was first investigated by Pullin in [16]. Due to the numerical instability, Pullin’s calculations were not conclusive. Using our improved method, we do not observe any numerical instability and we are able to perform well-resolved calculations to study the stabilizing effect of surface tension. Our study shows that surface tension indeed regularizes the early curvature singularity induced by the Rayleigh-Taylor instability. The interface rolls up into two spirals as time evolves. Note that the tips of the fingers broaden as they continue to roll, and that the interface bends towards the tip of the fingers. At around t = 1.785, the interface forms a trapped bubble and self-intersects. The minimum distance between the neck of the bubble is approximately 5 × 10^{−4} [4]. This process of bubble formation through self-intersection of a fluid interface has been observed in [8, 9] for a vortex sheet. In both cases, we found a convincing evidence that the minimum distance between the neck of the bubble scales like (t_{c} − t)^{2/3}, providing a partial agreement with the self-similar scaling.

Figure 1: Comparison of the zero surface tension interface profile with the corresponding ones for decreasing surface tension τ at t = 0.45, N = 2048. (a)τ = 2.5 × 10^{−4}. (b)τ = 1.25 × 10^{−4}. (c)τ = 6.25 × 10^{−5}. (d)τ = 3.125 × 10^{−5}.

Figure 2: Rayleigh-Taylor instability: Atwood ratio A = −0.1, surface tension τ = 0.005, N = 2048 and Δt = 1.25 × 10^{−4}.

References

[1] G. Baker, D. Meiron, and S. Orszag, Generalized vortex methods for freesurface flow problems, J. Fluid Mech. 123, 477-501 (1982).

[2] J. T. Beale, T. Y. Hou and J. Lowengrub, Growth Rates for Linearized Motion of Fluid Interfaces Away from Equilibrium, Comm. Pure Appl. Math, 46, 1269-1301 (1993).

[3] J. T. Beale, T. Y. Hou and J. Lowengrub, Convergence of a Boundary Integral Method for Water Waves, SIAM J. Numer. Anal., 33, 1797-1843 (1996).

[4] H. Ceniceros and T. Y. Hou, Convergence of a Non-Stiff Boundary Integral Method for Interfacial Flows with Surface Tension, Math. Comp., 67, 137-182 (1998).

[5] H. Ceniceros and T. Y. Hou, Dynamic Generation of Capillary Waves, submitted to Phys. Fluid. A, 1998.

[6] J. W. Dold, An Efficient Surface-Integral Algorithm Applied to Unsteady Gravity Waves, J. Comput. Phys., 103, 90-115 (1992).

[7] T. Y. Hou, Numerical Solutions to Free Boundary Problems, Acta Numerica, 1995, pp. 335-415.

[8] T. Y. Hou, J. Lowengrub and M. Shelley, Removing the Stiffness from Interfacial Flows with Surface Tension, J. Comput. Phys, 114, 312-338 (1994).

[9] T. Y. Hou, J. Lowengrub and M. Shelley, The Long-Time Motion of Vortex Sheets with Surface Tension, Phys. of Fluid, A, 9, 1933-1954 (1997).

[10] T. Y. Hou, I. Klapper, and H. Si, Removing the Stiffness of Three Dimensional Interfacial Flows with Surface Tension, to appear in J. Comput. Phys., 1998.

[11] T. Y. Hou, Z.-H. Teng, and P. Zhang, Well-Posedness of Linearized Motion for 3-D Water Waves Far From Equilibrium, Comm. in PDE’s, 21, 1551-1586 (1996).

[12] T. Y. Hou and P. Zhang, Stability of the Point Vortex Method for 2-D Water Waves, submitted to SIAM J. Numer. Anal., 1998.

[13] T. Y. Hou and P. Zhang, Stability of the Point Vortex Method for 3-D Water Waves, submitted to SIAM J. Numer. Anal., 1998.

[14] M. Koga, Bubble entrainment in breaking wind waves, Tellus, 34, 481-489 (1982).

[15] M. S. Longuet-Higgins and E. D. Cokelet, The deformation of steep surface waves on water, I. A numerical method of computation, Proc. Roy. Soc. London A, 350, 1-26 (1976).

[16] D.I. Pullin, Numerical studies of surface tension effects in nonlinear Kelvin-Helmholtz and Rayleigh-Taylor instability, J. Fluid Mech., 119, 507-532 (1982).

Applied Mathematics, 217-50 California Institute of Technology Pasadena, CA 91125 USA Email: hou@ama.caltech.edu