## Abstract

The numerical simulation of the deformation of an inviscid fluid–fluid interface subjected to an axisymmetric impulse in pressure is considered. Using a boundary integral formulation, the interface is evolved for a range of upper-fluid and lower-fluid density ratios under the influence of inertial, interfacial and gravitational forces. The interface is seen to evolve into axisymmetric waves or droplets depending upon the density ratio, level of surface tension and gravity. Moreover, the droplets may be spherical, tear shaped or elongated. These conclusions are expressed in a phase diagram of inverse Weber number *We*^{−1} versus Atwood number *At* at zero gravity, i.e. with the Froude number *Fr*^{−1}→0, and complement the earlier findings of Tjan & Phillips, who present a phase diagram of *We*^{−1} versus *Fr*^{−1} for the case in which the upper fluid has zero density. They too report tear-shaped droplets; however, while, in their paper, they form as a result of gravity, those reported here form as a result of surface tension. It is also found that the pinch-off process which effects drops remains of the power-law type with exponent 2/3 irrespective of the presence of gravity and an upper fluid. However, the constant that relates the necking radius to the time from pinch off, which is universal in the absence of gravity and an upper fluid, is affected by the presence gravity, an upper fluid and the class of drops which form.

## 1. Introduction

In a recent study, Tjan & Phillips (2007; henceforth TP) considered the deformation of a free surface subjected to an axisymmetric impulse in pressure. Depending upon the level of gravity and surface tension, expressed through the inverse Froude, *Fr*^{−1}, and Weber, *We*^{−1}, numbers (later defined), the surface was seen to evolve into axisymmetric jets, waves or drops. The liquid under consideration was inviscid and the density of fluid above it was assumed to be zero, or at least sufficiently small relative to that of the liquid to play no role in the dynamics of the evolving surface. Our purpose here is to relax that restriction by introducing an adjacent layer whose density may differ from zero. We then characterize the density ratio of the adjacent fluid layers with an Atwood number, *At*, and, given *Fr*^{−1} and *We*^{−1}, determine what role *At* has on the evolving interface.

The pressure impulse is assumed Gaussian, this being an approximation for the pressure field induced by focused ultrasound (Kino 1987). Moreover, the impulse is assumed to arise from a burst of ultrasound, the burst being long with respect to that of the period of the ultrasound but small with respect to the time scale of the evolving interface. The impulse is applied normal to the initial interface.

Liquid films excited by ultrasound have long been known to emit droplets, with studies dating from Wood & Loomis (1927) and Lang (1962; in the frequency range of 10–800 kHz) who report the formation of fog on the surface. On the other hand, focused ultrasound can be used to eject isolated droplets from the surface, and, to achieve this, Elrod *et al*. (1989) found bursts of ultrasound (in the range of 5–300 MHz) preferable to continuous forcing. In fact, free surfaces excited by focused ultrasound are observed to exhibit a variety of phenomena, including standing axisymmetric waves, ejected individual drops, multiple drops and a fog of tiny drops. All are evident in the (water to air) photographs of Sakai *et al*. (2003).

Considerable attention has been paid to the formation of drops and most particularly those induced by gravitational or periodic accelerative forcing. These include viscous dripping (Brenner *et al*. 1997; Wikes *et al*. 1999), liquid jet breakup (Lin & Reitz 1998) and parametric excitation (Faraday 1831; James *et al*. 2003), albeit at forcing frequencies well below those of ultrasound. On the other hand, little attention has been given to impulsively generated drops and, since this construct is significantly different, we expect the dynamics which affect droplet formation to likewise differ. We do not expect details of drop pinch off to differ, however, as these are probably ubiquitous in all studies.

Elrod *et al*. (1989) appear to be the first to have studied impulsively generated isolated drops and approach the topic from the viewpoint of ink-jet printing. Their work comprised both experiments and numerical modelling (to predict the time evolution of the surface). But they consider only water into air; there appear to be no studies of ultrasonically excited drops into fluids other than air.

Pinch off, on the other hand, has been studied on its own (e.g. by Eggers 1993; Brenner *et al*. 1996; Monika & Steen 2004), in the context of ink-jet printing (e.g. by Day *et al*. 1998; Basaran 2002; Leppinen & Lister 2003) and in the context of impulsively excited drops (TP).

In forming their model, TP include surface tension and gravity but ignore viscosity, arguing that because viscous terms (in the equation of motion) are small relative to their surface tension counterparts, they play little role until pinch off. They further note that subsequent to the pressure impulse, which gives rise to a vortex sheet on the interface, there is no mechanism to create further vorticity or sufficient time for vorticity in the sheet to diffuse, so that the interior fluid remains irrotational. Here, we follow suit and investigate the evolution of the interface of two unbounded inviscid irrotational fluids of different densities following an axisymmetric impulse in pressure, subject to the interplay of inertial, interfacial and gravitational forces.

Our problem involves solving the Laplace equation (2.1) subject to the Bernoulli equation (2.2) on the interface. For computational purposes, equation (2.1) is best recast into a boundary integral form and we do so in a manner akin to that of Baker *et al*. (1984). This technique has been recently employed by others (e.g. Hou *et al*. 1994, 1997, 2001; Nie & Baker 1998; Nie 2001), albeit on either a finite or a periodic domain. Here, however, the domain is semi-infinite and details of how we handle it are given in TP and outlined in §3. Results are given in §4, followed by a discussion of the finite-time singularity which occurs at pinch-off in §5 and remarks in §6.

## 2. Formulation

Our model is formulated as an initial-value problem in an unbounded axisymmetric domain. The governing equations and boundary conditions are introduced in §2*a* with non-dimensionalization discussed in §2*b*. Conversion of the governing equation to a Fredholm integral equation of the second kind is outlined in §2*c*.

### (a) Governing equations

Consider an axisymmetric domain , where is radial and vertical, composed of two subdomains _{1} and _{2}, which contain inviscid incompressible fluids of densities *ρ*_{1} and *ρ*_{2}, respectively. _{1} and _{2} initially occupy and and are separated for all time by a sharp interface *Γ*, along which the coefficient of surface tension is *σ*.

Assuming irrotationality, velocity potentials *φ*_{i} may be defined in each _{i}, with the velocity vector given by , for *i*=1, 2, while incompressibility demands that *φ*_{i} satisfy the Laplace equation(2.1)in each fluid. Under conservative body forces, the unsteady Bernoulli equation(2.2)then precisely enforces conservation of momentum.

When a pressure field is applied on =0 over a short time scale with respect to the time scale of the evolving surface, the pressure and, from (2.2), the velocity potential *φ*_{1} (since ) register as impulses relative to and thereby provide an initial condition to the problem.

The initial-value problem is completed by boundary conditions which require for all that _{i} vanish at infinity, while the usual kinematic and dynamic boundary conditions apply on the interface *Γ*. The former requires that the velocity normal to the fluid surface be identical to that of the fluid particle at the surface; the latter is the Laplace–Young condition (Young 1805), which requires that the pressure jump across *Γ* be proportional to the surface curvature .

In order to effect the Laplace–Young boundary condition on *Γ*, we introduce the scaled velocity potential *φ* and the dipole strength *μ*, which are defined on *Γ* as(2.3)Accordingly, we note thatwhere is the direction normal to *Γ* pointing from fluid 1 to fluid 2, with the tangential direction along *Γ*. This demands that the velocities of both fluids normal to *Γ* be continuous while permitting a possible jump in tangential velocities, a necessary condition to satisfy the kinematic boundary condition on *Γ*. We then see on *Γ* that , and , where *At* is the Atwood ratio of the densities of the two fluids, *viz*.(2.4)

### (b) Non-dimensionalization

We introduce velocity potential and length scales *Φ* and as per TP, in which *Φ* is a measure of the peak velocity potential and is a measure of the radial width of the impulse. Consequently, time and velocity scales follow as and *Φ*^{−1}, respectively, allowing us to write , , , , and , where *r*, *z*, *κ*, *ϕ*, *μ*, *t* and *u*_{n} are dimensionless quantities.

In terms of the scaled velocity potential and dipole strength, and on substituting the Laplace–Young boundary condition , the unsteady Bernoulli equation (2.2) valid on *Γ* then becomes(2.5)Two additional dimensionless parameters are evident in (2.5): the Weber number, , which measures the relative importance of inertial to interfacial forces and the Froude number, , which measures the relative importance of inertial to gravitational forces. Note that *μ*→2*ϕ* in the limit *At*→0, leaving

### (c) Boundary integral equation

The Laplace equation is recast into a boundary integral equation yielding a Fredholm integral equation of the second kind as (see TP)(2.6)where the directional Green's function to (2.1) isand *E*(*k*) and *K*(*k*) are, respectively, complete elliptic integrals of the first and second kind.

Accordingly, the azimuthal component of the vector potential is given by(2.7)whereFinally, from *B*_{θ}, we define a pseudo-streamfunction, *ψ*, as(2.8)

Then, knowing the value of the triple (*r*, *z*, *ϕ*) at any instant, we can solve the Fredholm integral (2.6) for the dipole strength *μ* and subsequently use (2.7) to evaluate *B*_{θ}(*s*). Then, with the knowledge of the velocities from (2.8), we use the Bernoulli equation (2.5) to evolve *ϕ* forward in time, while enforcing the kinematic boundary condition to evolve the surface forward in time. The process can then be repeated with the updated *ϕ* and new surface profile as input.

Finally, to set the process in motion, we require an initial condition for *ϕ* on *z*=0; this may be any continuous *C*^{2} function that vanishes as *r*→∞, but for comparison with TP we set .

## 3. Numerics

Beyond the introduction of a further parameter, the formulation for handling two fluids is not greatly different from that with one fluid and presented no new numerical challenges. The numerical procedure thus follows exactly that in TP and hence only an outline is presented here. To proceed, we first introduce a mapping from the semi-infinite physical domain *s*[0,∞) to a finite computational domain *η*[−1,+1]. The dependent variables {*r*, *z*, *ϕ*, *ψ*, *μ*} are then expressed in terms of basis functions of order *N* over *η* and coefficients for {*r*, *z*, *ϕ*, *ψ*} are determined using collocation, while those for *μ* are found via a Galerkin method.

In order to follow the evolution of the (initially flat *z*=0) surface, we place on it *N* Lagrangian markers, each located by radial *r*(*s*) and vertical *z*(*s*) coordinates, and monitor them. The markers, together with the velocity potential *ϕ*(*s*), form the triple {*r*, *z*, *ϕ*} that defines the primary dependent variables of the system and these are supplemented by the further dependent variables *ψ* and *μ*. The numerical solution involves marching these variables forward in time.

For example, given the normal and tangential velocities ** u**=(

*u*

_{n},

*u*

_{t}) at the surface, the location of the markers

*r*and

*z*are evolved kinematically as(3.1)We note that while

*u*

_{n}is unambiguously given by (2.8), the choice for

*u*

_{t}remains arbitrary as the Lagrangian markers may be advected with any choice of tangential velocity without changing the surface shape. In other words, while the Lagrangian markers necessarily have the same normal velocities as the fluid on either side of the interface, their tangential velocities will in general be different. The following choice for the tangential velocity is used,

*viz*.as it has the desirable property (see TP) that the distance between markers is preserved. Accordingly, such a choice for

*u*

_{t}is especially useful for a Chebyshev collocation implementation, since now

*s*for each collocation point is fixed in time. Various values for the basis function number, and thus the number of mesh points,

*N*, were explored but

*N*=64 was found to provide an adequate balance between resolution and computational effort. In fact, the results were little changed by higher values of

*N*, suggesting that the calculation had converged. Complete details are given in Tjan (2007).

The evolution of *ϕ* is via the Bernoulli equation (2.5) with the partial derivative written as a material derivative, i.e. , yielding(3.2)which reduces to that of TP when *At*=1. (Note that (3.2) in TP contains a typographical error.)

Values of {*r*, *z*, *ϕ*} are updated using a fourth-order Runge–Kutta scheme, after which they are used to calculate the dipole strength by solving the Fredholm integral (2.6).

## 4. Results

Our problem is in terms of the three-dimensional parameter space (*We*, *Fr*, *At*), which is somewhat unwieldy, hence, to reduce the dimensionality, we hold *Fr*^{−1} constant and vary *At* and *We*. Specifically, we restrict our attention to the interplay between inertial and surface tension forces in the absence of gravity (i.e. *Fr*^{−1}→0), while varying *At* over its full range, namely *At*∈[−1,+1]. TP, on the other hand, set *At*=1 and investigated the role of surface tension versus gravity over a range of Weber and Froude numbers. For reference, their phase diagram for *At*=1 and *Fr*^{−1}>0 is presented in figure 1. It depicts regions where spherical drops, inverted tear-shaped drops and axisymmetric waves form; it also highlights a region, specifically for *We*^{−1}<0.045, where surface tension is too weak to overcome numerical instability. The limit *Fr*^{−1}→0 is not shown, but spherical drops occur in the range 0.045<*We*^{−1}<0.105 and their size scales inversely with *We*^{−1}. Once *We*^{−1}>0.105, however, surface tension is sufficiently dominant to impede the evolving axisymmetric wave that precedes the drop, causing it to collapse upon itself without forming a drop.

### (a) *Fr*^{−1}→0, *At*∈[−1,+1]

Here, we find that spherical drops also form when *At*<1 for some *We*^{−1}, as shown in a sequence of snapshots (for *At*=0.75) in figure 2. Accordingly, the size of the drop continues to scale inversely with the Weber number, although the maximum height attained by the drop is affected by *At*. But spherical drops do not form for all *At* and certainly not for all *We*^{−1} investigated. Indeed, as *At* is reduced from unity, which means that the density of the upper fluid is increasing relative to that of the lower fluid, topologies beyond those reported by TP are observed. A phase diagram is given in figure 3.

First, we observe for all *At* that an upper bound in *We*^{−1} exists beyond which no pinch off to drops is observed. Rather, in this range of *We*^{−1} and *At*, the interface oscillates as an axisymmetric wave about the mean elevation and is eventually damped. Snapshots of the free surface of a typical case are shown in figure 4. Such behaviour is reminiscent of a situation identified by TP where gravity acted to damp the axisymmetric wave, only here gravity is absent and oscillations are damped over time by surface tension. Thus, while TP observed axisymmetric gravity waves, we observed axisymmetric capillary waves. Furthermore, as the density differential is progressively reduced to the point that *At* becomes negative, i.e. the heavier fluid on top, the upper bound in *We*^{−1} is reduced. This suggests that although the level of surface tension required to prevent drop formation has reduced, there is, as we might expect physically, an increased inertial restraint from the upper fluid.

On lowering *We*^{−1} to just below the no-drop upper bound, we find for *At*∈(−1,+1) that tear-shaped, rather than spherical, drops form. Tear-shaped drops were also observed by TP, but the dynamics which gave rise to them is there different. Specifically, the evolving interface there collapsed onto itself under gravitational forces, whereas here, in the absence of gravity, the evolving interface is slowed jointly by inertial effects from the upper fluid coupled with surface tension, which eventually causes necking at the base of the tear-shaped drop. Snapshots of the free surface of such a case are depicted in figure 5. Here, there is no collapse, so we might say that while TP observed gravity tears, we observed surface tension or capillary tears.

Further reductions in *We*^{−1} expose, mainly for *At*>1, the formation of spherical drops and ultimately another class of drops not previously reported. For reasons evident in figure 6, we refer to this class as ‘pancake’ or ‘elongated’ drops. Here, surface tension would appear to have less influence than inertia to the extent that the upper-layer fluid grossly distorts the shape of the drop.

Eventually, we encounter a lower bound in *We*^{−1} for all *At*∈[−1,1], below which surface tension is insufficient to overcome numerical instabilities and the calculation breaks down.

Finally, in viewing figure 3, we should point out that although the boundary between the formation of axisymmetric waves and tear-shaped drops is clearly defined, the interface between other classes is less clear. For example, tear-shaped drops become more spherical (or elongated at some *At*) in the vicinity of the boundary and there is a range of *We*^{−1} near the boundary in which they are not distinctly one class or the other. To that end, the boundaries (other than the axisymmetric wave one) indicated on figure 3 should be viewed as indicators of where transition occurs rather than rigid demarcations.

## 5. Scaling

Pinch-off is an example of a singularity which is formed in a finite time, the bifurcation being a topological singularity where the surface self-intersects at time *t*=*t*_{0}. Of particular interest is how the system approaches this singularity and efforts to understand it began with an experimental study of a spherical pendant drop evolving from the end of a nozzle by Peregrine *et al*. (1990).

Mathematically, the self-intersection of the surface may be precisely phrased as *r*(*s*≠0)=0, where the radial coordinate vanishes for some location *s*≠0. To proceed, therefore, we define a radius *r*_{min} and track it as the surface evolves. Thus, let *r*_{min} be the radial coordinate of the point on the surface satisfying the two conditions d*r*/d*s*=0 and d^{2}*r*/d*s*^{2}>0, which decree that there will be a solution for *r*_{min} only after such time that the surface becomes vertical at some location.

TP found that *r*_{min} exhibits a power-law type singularity of the form (*t*_{0}−*t*)^{γ}, where *γ*≈2/3. They further note, as do Keller & Miksis (1983), that such behaviour can be argued on dimensional grounds, on the assumption that the behaviour should, near the singularity time, be independent of initial conditions and thus that the relevant (dimensional) parameters are *σ*, *ρ*_{1}, and _{0}−. From this set of parameters, the Buckingham pi-theorem dictates that only one dimensionless group can form and must itself be a universal constant, say . When further parameters play a role, a second density *ρ*_{2} and gravity *g*, for example, further dimensionless groups arise rendering a functional as (say) . Then, withwe have for all *g* that(5.1)Equation (5.1) thus recovers TP's expression in the double limit , at which they find . So of interest here is how non-zero values of affect . To that end, we study, in §5*a*,*b*, respectively, the limits , and , .

### (a) *Fr*^{−1}→0 with variable *At* and *We*

As is evident from (5.1), the pinch-off scaling law is not affected by the density ratio , although can affect the magnitude of . Our purpose here is to determine as a function of density ratio and, to expose that relationship, we exclude by working in the *Fr*^{−1}→0 limit. In addition, because , we prefer, with no loss of generality, to use *At* as the independent variable (because *At*∈[−1,1] (2.4)) and plot , which we do in figure 7.

For each class of drops, we see that decreases monotonically with increasing *At* and recovers TP's universal constant only for spherical drops when *At*=1. Observe also that is lowest for spherical drops and successively higher, at a given *At*, for tear-shaped and elongated drops, although only tear-shaped drops occur over the full range of *At* (figure 3). Thus, in applying (5.1), we must be aware that is affected not only by the density ratio but also by the topology of the emerging drop.

### (b) *At*=1 with variable *We*^{−1} and *Fr*^{−1}

We turn now to the case and note that although TP explored the phase space for with variable *We*^{−1} and *Fr*^{−1}, they explored the near singularity scaling relationship only in the limit. Here, we are interested in the role of gravity, so we set and explore . To ascertain this dependence, we plot against in figure 8 and for clarity depict spherical drops and tear-shaped drops separately, in figure 8*a*,*b*, respectively. As in TP's case, the r.h.s. of figure 8*a*,*b* indicates that is essentially independent of , indicating that in spite of gravity being present it plays no dynamical role in pinch off for either spherical or tear-shaped drops. However, as we found with the density ratio, gravity does influence the magnitude of , causing a significant departure from its zero-gravity single-layer universal value; indeed, we here find . The value of is also affected by the topology of the drops which form, being at the lower end for spherical drops and higher end for tear-shaped ones.

For , on the other hand, there would appear from figure 8 to be a noticeable gravitational influence, but our conclusion is that there is not. Rather, numerical uncertainties in deducing *r*_{min} and *t*_{0}−*t* as each approaches zero are compounded in evaluating *r*_{min}(*t*_{0}−*t*)^{−2/3}, leading to specious results in the double limit. To pursue this further, we avoid evaluating the double limit and instead plot *Π*_{1} against *Π*_{2}, as shown in figure 9. Here, there is less sensitivity as *r*_{min} and *t*_{0}−*t* approach zero and the data collapse about distinct straight lines of slope representative of the class of drops which form.

## 6. Remarks

We continue, in this work, a study by TP who consider the evolution of a free surface subject to an axisymmetric impulse in pressure. Their research was part of a study of haemorrhage in the lung caused by ultrasonic imaging and, in particular, to explore a non-thermal, non-cavitational damage mechanism. The mechanism they consider is built around the notion that ultrasound focused near a tissue–liquid interface acts to expel tiny droplets of blood or other fluids which then puncture the soft bubble-wrap-like air-filled sacs (alveolar) of the lung pleural surface. Moreover, they further show that droplets can indeed be ejected over the range of intensity levels, measured by a dimensional quantity denoted the mechanical index, employed in clinical ultrasonography. Of course, the interface is not always liquid–gas, rather the plural surface may at times be coated by a layer of mucus, and for that reason we here extend the work to include fluid–fluid interfaces with a broad range of density ratios. With that as background, however, we preferred not to restrict ourselves solely to the conditions associated with lung haemorrhage, but rather to look in general at impulsively generated waves and drops. Indeed, the phenomenon has the tantalizing prospect that it could be related to other physical scenarios such as ink-jet printing. Our findings also shed light on what might well be a ubiquitous finite-time singularity at drop pinch off, as we find that parameters other than those evident in zero gravity, liquid–air pinch-off, affect only the constant relating the necking radius to the time from pinch off, not the scaling law itself.

## Acknowledgments

This work was supported by the National Institutes of Health grant R01-EB02641 (formerly HL58218). Our thanks to Jonathan Freund, Gustavo Gioia and Jonathan Higdon for their interest and helpful comments.

## Footnotes

- Received November 6, 2007.
- Accepted January 7, 2008.

- © 2008 The Royal Society