Log In Sign Up

Signal processing approach to mesh refinement in simulations of axisymmetric droplet dynamics

by   Kazuki Koga, et al.

We propose a novel mesh refinement scheme based on signal processing for boundary integral simulations of inviscid droplet dynamics with axial symmetry. A key idea is to directly access the Fourier coefficients of a principal curvature as a function of the arclength through a natural change of variable. The trapezoidal rule is applied to those Fourier-type integrals and the resulting formula fits in the framework of the non-uniform fast Fourier transform. This observation enables to efficiently generate smooth guidelines for mesh refinement in two singularity formation scenarios. Applications also include a non-iterative construction of the uniform parametrization for an important class of plane curves, which is used in a convergence study of the time-stepping procedure implemented in the previous work by Nitsche and Steen [J. Comput. Phys. 200 (2004) 299].


page 1

page 2

page 3

page 4


Vico-Greengard-Ferrando quadratures in the tensor solver for integral equations

Convolution with Green's function of a differential operator appears in ...

Digraph Signal Processing with Generalized Boundary Conditions

Signal processing on directed graphs (digraphs) is problematic, since th...

Abstract message passing and distributed graph signal processing

Graph signal processing is a framework to handle graph structured data. ...

NFFT.jl: Generic and Fast Julia Implementation of the Nonequidistant Fast Fourier Transform

The non-equidistant fast Fourier transform (NFFT) is an extension of the...

Understand Slope Limiter – Graphically

In this article, we illustrate how the concept of slope limiter can be i...

Using Frame Theoretic Convolutional Gridding for Robust Synthetic Aperture Sonar Imaging

Recent progress in synthetic aperture sonar (SAS) technology and process...

1 Introduction

A deep understanding of droplet motion under the action of surface tension is essential in various applications such as ink-jet printing, spraying, and atomization. There are a large body of literature on theoretical, experimental, and numerical studies of the physics, and a significant part of those works is devoted to investigating singularity formations observed when a fluid interface self-intersects. Such phenomena have been of great interest for their universality, richness, and elegant characterizations.

To reveal the structures of those singularities, a number of works have performed boundary integral simulations for potential flow models that describe the motion of an interface separating two inviscid fluids. For example, a two-dimensional droplet formation by a roll-up due to the Kelvin-Helmholtz instability is carefully studied by Hou et al. [1] (henceforth HLS97). In an axisymmetric case, Leppinen and Lister [2] investigate the self-similarity of inviscid capillary pinch-off for a range of density jumps across the surface of a droplet. Their results are more firmly confirmed by Nitsche and Steen [3] (henceforth NS04) with additional insights into the nature of the self-similarity. Furthermore, a recent work by Burton and Taborek [4] revisits pinch-off of 2D droplets and finds that the type of the singularity is distinctly different from that of its axisymmetric counterpart.

Besides the dimension of the problem being reduced by one, a major advantage of a boundary integral simulation is its flexibility in controlling the computational mesh. The successful works mentioned above undoubtedly rest on mesh refinement schemes that dynamically redistribute computational points densely in the regions of interest. However, the strategies employed in the previous studies are highly problem-specific. That is, it is required to know “a priori” a few properties of given initial conditions, which precludes applications of those methods to general cases.

As an extension of NS04, this paper presents a novel mesh refinement strategy for resolving singularity formations in the vortex sheet formulation of axisymmetric droplet dynamics driven by surface tension. A vortex sheet is a sharp interface separating two inviscid, incompressible, and irrotational fluids shearing past each other, which has a boundary integral formulation that derives from the Biot-Savart law. In our scheme, guidelines for mesh refinement are generated based on signal processing applied to the Fourier coefficients of the curvature in the symmetry plane as a function of the arclength. These Fourier coefficients are discretized by the trapezoidal rule after a natural change of variable, and the resulting discrete formula is approximated by the non-uniform fast Fourier transform [5] with complexity. It turns out that our method is fully automatic and essentially relies only on the assumption that any singularity will not appear in a very short time. Therefore, it is applicable to a wider class of initial conditions with far less human intervention in comparison to the existing methods. We test the numerical scheme on the pair of initial conditions studied by Nie [6] (henceforth Nie01) that reach two distinct singularity formations at the end, and discuss several aspects of our strategy related to its efficiency and accuracy. Also, applications of our ideas include a non-iterative algorithm for constructing the uniform parametrization of closed plane curves, which is utilized in a convergence study of the time-stepping procedure that has not been sufficiently justified in NS04.

The rest of this paper is organized as follows. Section 2 introduces a boundary integral formulation of axisymmetric vortex sheets with surface tension and two relevant initial conditions. Section 3 describes the theory behind our mesh refinement strategy based on signal processing. Section 4 explains our numerical method and some technicalities in its implementation. Section 5 presents numerical results and validates the implemented scheme. Concluding remarks are given in Section 6.

2 Formulations

In this section, we give a brief introduction to boundary integral equations for axisymmetric vortex sheets with surface tension, which is followed by a short description of finite-time singularities in the formulation reported by other authors. For interested readers, a more detailed derivation can be found in NS04.

2.1 Lagrangian theory

We consider irrotational motion of two inviscid, incompressible, and immiscible fluids in separated by a closed surface . In the following, the inner and outer fluid are denoted by and , respectively, and is assumed to be simply-connected. For simplicity, we also assume that the density of each fluid is uniformly equal to unity and that no external force acts on the system. Then, the fluid velocity and pressure in () satisfy the Euler equation with the boundary conditions


where and are the spatial and time variable, respectively. Also,

is the unit normal vector of

pointing to , is twice the mean curvature that is positive if is a sphere, and is a positive surface tension coefficient. The brackets measure the jump of a function across in the following manner:


The far-field condition (2.1) states that the fluid at infinity is at rest. The kinematic boundary condition (2.2) requires that the normal component of the fluid velocity must be continuous across , while the tangential component can be discontinuous there. The Laplace-Young condition (2.3) represents the effect of surface tension on in terms of the mean curvature and the pressure jump across the interface.

Figure 1: Sketch illustrating: (a) axial symmetry, (b) parametrization in plane.

Throughout the rest of this paper, we assume that the flow is axially symmetric. That is, in some cylindrical coordinates , all functions describing the flow are independent of the azimuthal coordinate (see Fig.1(a)). In this case, for tracking the motion of a free surface , it suffices to solve time evolution of the intersection of and a symmetry plane Thus, we parametrize the curve in the plane by


where and are radial and axial coordinates, respectively, and find its evolution equation


Again, is the unit normal vector pointing to , and is the unit tangent vector in the counterclockwise direction. In particular, we set for all , and direct so that it satisfies (see Fig.1(b)). Here, the subscript denotes differentiation with respect to the variable , and is called the relative spacing that is the derivative of the arclength measured from the starting point . In the present case, the kinematic boundary condition (2.2) implies that , where is the average of the fluid velocities on either side of , whereas it puts no restriction on . In fact, it is well-known that the shape of an evolving plane curve is determined solely by the normal velocity , and the tangential velocity only affects the relative spacing .

To obtain an integral form of , we introduce a dynamical variable called the unnormalized vortex sheet strength that measures the jump of the tangential component of the fluid velocity in a frame-dependent manner:


Then, assuming that the azimuthal component of the fluid velocity is uniformly zero (i.e. without swirl), the Biot-Savart law and the Plemelj’s formula imply that is given by the principal-values (denoted by dashed integrals)


where, for example, and (e.g., see [7] for details). Here, and are the complete elliptic integrals of the first and second kind, respectively. Other functions in (2.8) and (2.9) are defined as


On the other hand, the evolution equation for is a consequence of the unsteady Bernoulli’s theorem and the Laplace-Young condition (2.3):


In axisymmetric geometry, it is convenient to express as the sum of


which are the principal curvatures in the symmetry plane and in a plane normal to it, respectively. Without surface tension (i.e. ), the first term in the right-hand side of (2.11) is absent, and choosing renders time-independent. This special choice for is often referred to as the Lagrangian frame, and it has been studied mainly in the context of the famous Moore’s singularity [8]. With surface tension, however, it is no longer straightforward to keep conserved in time. Moreover, it is known that the Lagrangian frame typically leads to rapid decreases in at points where geometry is not necessarily complex, as demonstrated by HLS97 and NS04. These facts motivate us to choose to control the relative spacing for other practical purposes.

Following Hou et al. [9] (henceforth HLS94), we switch from the Cartesian coordinates to the so-called angle–arclength variables defined via the relation


Having these new dynamical variables, we are able to reconstruct up to a translation by integrating Eq.(2.13). Such a loss of information is not a problem here, for the current interest lies in developing mesh refinement schemes to resolve singularity formations due to self-intersections of the interface . The evolution equations for follow from the orthogonal decomposition in (2.6), the Frenet-Serret formula for plane curves, and the identity :


Eq.(2.15) is central to the present work. As mentioned earlier, because of the tangential velocity being arbitrary, the relative spacing can be controlled by a user-specified while keeping the image of unchanged. In fact, this can be done by deriving a desirable first and calculating accordingly based on Eq.(2.15). Our choice for is described in detail in Section 3.

2.2 Singularity formations

It is well-known that, without any regularization, the dynamics of vortex sheets is subject to the Kelvin-Helmholtz instability and generally leads to the Moore’s singularity [8]. On the other hand, HLS94 points out that surface tension acts as a dispersive regularization of the Kelvin-Helmholtz instability, and also reports that numerical solutions with continue beyond the Moore’s singularity for the initial condition studied by Krasny [10] in the zero surface tension setting. Instead of the Moore’s singularity, however, several works have found that other kinds of singularities appear when a fluid interface self-intersects. Such phenomena are often called topological singularities.

Figure 2: Solution close to singularity for pinch-off problem (2.16): (a) interface position in plane, (b) profile of curvature .

The most famous topological singularity is presumably capillary pinch-off, which is characterized by a fast drop in the radius of a neck due to surface tension acting on the azimuthal curvature . To study this process numerically with the vortex sheet formulation, Nie01 proposes the initial condition


Qualitatively, this problem resembles an experiment by Robinson and Steen [11]. As an intuitive illustration, we show in Fig.2 a snapshot of a numerical solution starting from (2.16) and the corresponding curvature as a function of divided by the total arclength . As seen in Fig.2(a), the solution is symmetric with respect to the –axis and eventually forms two end-droplets at the top and bottom along with an elongated satellite droplet in between. Also, one can easily see that overturning occurs in each pinch-off region, which precludes descriptions of the process by a single-valued function of . NS04 studies the same problem numerically to confirm the self-similarity of inviscid capillary pinch-off


for a range of density jumps across the interface. Here, are the coordinates of a point that attains the minimum neck radius, is the axial coordinate where the pinch-off occurs, and is the time when the radius reaches exactly zero. From a numerical simulation perspective, the initial condition (2.16) is particularly suitable for an accurate numerical study on the scaling law (2.17), for it has an additional line symmetry and the shape of the interface is considerably simple away from the necks (see Fig.2(b)). In this paper, this problem is revisited to assess efficiency and accuracy of our numerical schemes.

Another scenario of topological singularity formations is self-intersections due to inertia resisted by surface tension. In the present formulation, this singularity is found for the initial condition


which is also suggested by Nie01. This situation can be understood as a uniform flow past a sphere [12] where the solid boundary is instantly dissolved at and evolves under the action of surface tension. Again, Fig.3 plots a numerical solution starting from (2.18) and its curvature .

Figure 3: Solution close to singularity for axisymmetric bag breakup (2.18): (a) interface position in plane, (b) profile of curvature .

As seen in Fig.3(a), the interface collides itself away from the axis of symmetry, and the droplet seems to break into a torus-like component and a simply-connected one. For brevity, we refer to this type of singularities as bag breakup. The nature of bag breakup is expected to be different from that of capillary pinch-off, because in this case surface tension should always tend to restore the droplet to a spherical shape. Regarding this point, motivated by their own experiments [13], Burton and Taborek [4] numerically investigate bag breakup in two-dimensional flows (called 2D pinch-off by the same authors), and they find that its self-similarity is of the second kind. That is, the local behavior of a fluid interface close to a self-intersection is characterized by the memory of an initial condition, and the spatial variables follow scaling laws with different exponents. Furthermore, by using the full model and its “slender” approximation, they find numerical evidences of those self-similarity exponents distinct from that of (2.17). We speculate that a similar analysis may be applicable to the axisymmetric case. In this paper, however, we do not aim to identify the similarity exponents in the axisymmetric bag breakup (2.18), because it forms multiple high-curvature regions away from the self-intersection and is therefore inappropriate for that purpose (see Fig.3(b)). Rather, our interest lies in whether it is possible to develop a mesh refinement scheme that is able to capture singularity formations and such high-curvature regions simultaneously.

3 Mesh refinement

The choice of the tangential velocity , which we have left undefined so far, is crucial in numerical simulations of evolving plane curves. This is particularly true if one needs to resolve small-scale structures of a fluid interface approaching an imminent singularity formation. In the context of potential flow models, several works use the Lagrangian frame in a combination with some regridding techniques. For example, in simulations of capillary pinch-off with axial symmetry, Leppinen and Lister [2] advect computational points by the Lagrangian frame and redistribute them at each time step in order to keep the relative spacing proportional to the distance from the point of the minimum neck radius. A similar strategy is used by Burton and Taborek [4] to study bag breakup in two-dimensional flows. Although this idea provides a simple way of accumulating points to high-curvature regions, it does not naturally extend to general cases because their methods heavily rely on the fact that the initial conditions studied in those works lead to self-intersections right on a targeted coordinate axis.

In the following, we employ a style in mesh refinement that harnesses the arbitrariness of the tangential velocity . Namely, the equation for that realizes a given is derived for a simple open curve identified with an axisymmetric surface whose inner region is simply connected, which implies that


for all . At a glance, this condition might seem to limit the applicability of our method. However, with minor modifications to (3.1), we believe that our ideas easily generalize to the cases of periodic curves and unbounded curves with periodic boundary conditions.

Following HLS94, we define the relative spacing as a product of the total arclength and a strictly positive function whose integral over is equal to 1:


By differentiating both sides of the first equation in (3.2) with respect to and substituting the result into (2.15), we get


Also, by combining (2.15) and (3.1), it is easy to see that


Now, by integrating (3.3) with (3.1), we obtain the following integral equation for :


The uniform parametrization, by which we mean that is independent of the variable , corresponds to [9]. In this case, Eq.(3.5) is explicitly solved for as


The uniform parametrization is a typical choice for various moving boundary problems because, unlike the Lagrangian frame, it prevents computational points from clustering at undesired locations and maintains stable spatial discretization. However, it is almost obvious that this is not the best option for accurate numerical studies of singularity formations in general. From this point of view, for a two-dimensional droplet formation by a Kelvin-Helmholtz roll-up, HLS97 alternatively proposes a variable but time-independent with a few minima, which is empirically constructed from numerical results for the same problem with the uniform parametrization. A major advantage of this method is that it is still possible to solve Eq.(3.5) explicitly for in the same form as (3.6), but such a function is literally problem-specific by definition.

As an extension of HLS97, NS04 suggests a strategy for constructing a variable and time-dependent to solve the pinch-off problem (2.16). Their choice for is initially uniform and develops exactly two symmetric minima at and . More specifically, they first prescribe the function


for some constant , and compose it with a monotonically increasing function that also evolves in time and allows the two minima of to dynamically follow the points of the minimal neck radius. Using this composite function , they define the function as


and apply it to the initial condition (2.16) for a range of density jumps across the interface. Although their computations are successful in many aspects, it has to be said that this definition is still specialized to a class of initial conditions with a line symmetry. Moreover, it requires to know in advance the pinch-off time , the number of singular points, and even a characterization of the point that never causes to be discontinuous in time.

The main difficulty in the design of mesh refinement for singular interfacial problems is that the curvature , a pointwise measure of the complexity of a plane curve, is by itself unlikely to define a promising candidate for . Intuitively, there are two reasons for this consideration. Firstly, as shown in Fig.2(b), a peak in that corresponds to a self-intersection can be significantly sharp and therefore, if is written in terms of , the relative spacing may reflect this sharpness as a dynamical variable. Secondly, it is found in Fig.3(b) that there is an inflection point (i.e. ) in a small interval where the gradient of is very steep, which causes to increase locally and magnifies many times. A similar situation can be seen in Fig.2(b), and it is numerically verified in NS04 that there exists at least one inflection point in the self-similar regime (2.17). In fact, to our knowledge, there is no previous study that succeeds in resolving self-intersections of a fluid interface with mesh refinement based on pointwise evaluations of . Nevertheless, it is no doubt that the graphs of the curvature visually show where singularities are about or at least likely to appear. In some sense, the previous works draw rough sketches of such visual cues as their function with the help of numerical experiments or theoretical considerations.

We attempt to further extend the idea of NS04 by developing a new tool that automatically generates guideline functions, denoted by , in place of the composite function . In particular, we automate the “sketching” process by applying filtering methods in signal processing to the curvature as a periodic function of the arclength . For this purpose, the first step is to extend the parametric curve to the interval with the symmetry


as illustrated in Fig.1(b), and consider the Fourier coefficients


where denotes the total arclength of the periodic curve . This quantity can be accessed via the natural change of variable :


and is expected to be invariant under a change of parametrization that preserves the starting point and the direction of increasing . This property is essential as an element in mesh refinement, because we cannot construct an appropriate from any frame-dependent information.

Next, we introduce the so-called analytic signal [14]


where is the Hilbert transform, which is the Fourier multiplier of the form


In harmonic analysis, this operator plays important roles in the study of the convergence of the Fourier series for -integrable functions with [15]. On the other hand, in signal processing, the analytic signal is often called a pre-envelope, because, by rewriting as


one can extract the analytic envelope [16]


Here, the angle is called the instantaneous phase, which is also defined in terms of and its Hilbert transform. To illustrate this concept schematically, the analytic envelope is demonstrated in Fig.4 for the following functions:


where , , and .

Figure 4: Signals and corresponding analytic envelopes: (a) Gaussian pulses , (b) product of two sine functions .

The signal is a superposition of two Gaussian pulses [17] that have compact numerical supports because of the fast decaying prefactors. We show the raw and its analytic envelope in Fig.4(a). Here, note that the left pulse has one zero point between two peaks with different signs, while the right pulse has one sharp peak surrounded by multiple zero points. In both cases, the operator ignores all these zero points and transforms the original pulses to simple and strictly positive peaks. Therefore, it is expected that the analytic envelope can remove the difficulty caused by inflection points near self-intersections. On the other hand, the signal is a product of sine functions with a high and low frequency, as shown in Fig.4(b). In this case, the operator fails to detect any local complexity and simply returns a function similar to , the true envelope of . These results imply that the analytic envelope is effective only if high-curvature regions are localized and sufficiently isolated from each other. Moreover, the second example shows that the analytic envelope can be non-differentiable at zero points of the original signal. To avoid this problem, we use a regularized envelope


Adding a constant inside the square root renders the final form of slightly less sharp, but this regularization is clearly inevitable for the sake of smoothness. Also, it should be noted that the analytic envelope of a nearly-singular function is inherently very sharp. In this paper, to obtain smoother than the original , we apply the Gaussian filter to the analytic envelope . That is, the function is convolved with the heat kernel


where the parameter controls the smoothing effect. To summarize, we suggest the following guideline function as a “rough” measure of the complexity of an evolving plane curve:

Figure 5: Curvature and generated guideline function at times close to singularity: (a) pinch-off problem at , (b) axisymmetric bag breakup at .

In Fig.5, the functions with are shown for the curvature in Fig.2(b) and Fig.3(b). In both cases, our technique is successful in generating simple and strictly positive functions that outline relevant features of the original signals .

Finally, we define the relative spacing in terms of our guideline function . In our definition, the function is split into a few components using a parameter , which determines the maximal amount of refinement, and another Gaussian as a function of the time :


Here, and are defined as


and these definitions trivially guarantee that


Moreover, it is easy to show that the definition (3.21) satisfies the relations


Hence, at , the parametrization is uniform and the tangential velocity is given by the formula (3.6). At , however, the equation (3.5) is no longer solvable explicitly for . To circumvent this problem, NS04 suggests the backward difference approximation


for some . This formula allows to find in the same way as for the time-independent , although errors of are introduced to the relative spacing . We specify the values of in Section 4, and a validation of Eq.(3.25) is given in Section 5.

4 Numerical methods

We numerically solve the initial value problem of Eqs.(2.11)(2.14)(2.15) with the normal velocity and the tangential velocity specified in Section 3. Our numerical method consists mainly of spatial discretization, evaluations of the guideline function , and temporal discretization compatible with the approximation (3.25).

4.1 Spatial discretization

In our computations, the dynamical variables are sampled at points equally spaced in the extended parameter space :


and all the derivatives and anti-derivatives are evaluated using the standard fast Fourier transform (FFT). For the principal-values (2.8) and (2.9), we employ the third-order modified dBM approximation [18] with 8193 quadrature nodes. The values of between the grid points (4.1

) are obtained by the Fourier interpolation after the reconstruction from

. In our code, evaluations of the integrands in (2.8) and (2.9), including the iterative algorithms for the complete elliptic integrals and [19], are parallelized on a single NVIDIA Tesla P100, a Graphical Processing Unit (GPU) for scientific computing in the double-precision arithmetics. Also, the direct summation after the integrand evaluations is performed with the warp-shuffle reduction algorithm [20] on the same GPU. Here, we are aware that higher-order quadrature rules for the axisymmetric Biot-Savart integrals have been suggested [21, 22] and successfully implemented by Nie01 and NS04 for the case with surface tension. However, those methods are originally designed for studies on an axisymmetric analogue of the Moore’s singularity, which requires to approximate the principal values beyond the double-precision accuracy. In the present work, it is found that the method described above yields satisfactory results in validating our mesh refinement scheme and is therefore useful in practice.

4.2 Evaluating

Next, we apply the trapezoidal rule to the Fourier coefficients (3.10) after the change of variable (3.11):


Unfortunately, this discrete sum cannot be directly computed using the standard FFT, because the nodes are redistributed non-uniformly in by mesh refinement. Instead, we can use the non-uniform fast Fourier transform (NUFFT), which is an algorithm approximating the following types of sums to a prescribed relative accuracy :


It is easy to see that the formula (4.2) naturally fits in the framework of the Type-1 NUFFT (4.3). On the other hand, once the Fourier coefficients of are obtained, it can be interpolated back to the non-uniform grid using the Type-2 NUFFT (4.4). Here, we note that there are several implementations of the NUFFT algorithms for general purposes, including CMCL NUFFT [23], NFFT3 [24], and FINUFFT [5]. In our code, we employ FINUFFT with a few cores of Intel Xeon Gold 6136, for, to our knowledge, it is currently the fastest library parallelized with OpenMP. Also, it should be noted that the Type-1 NUFFT itself is not related to any quadrature rules and simply approximates sums of a particular form. In our scheme, the factors serve as quadrature weights compatible with the non-uniform nodes that evolve in time, which automatically solves this issue.

Our practical procedure for evaluating is as follows. Firstly, we decide the largest wavenumber of the Fourier coefficients computed for the mesh refinement purpose. This is equivalent to choosing a target spatial resolution in terms of the uniform grid that we aim to achieve with a non-uniform grid realized by our mesh refinement. Therefore, the number should be related to the parameter in (3.21) so that the target resolution will always be no finer than the best possible local resolution determined by . Secondly, if we assume that our scheme is working ideally, the locally refined grid resolves successfully, whereas complex exponentials, which oscillate uniformly in , can be underresolved for sufficiently large . To resolve both and those complex exponentials, the functions , , and are interpolated by the truncated Fourier serires to an upsampled grid with points in , and then the coefficients are computed using the Type-1 NUFFT. Thirdly, once the coefficients are obtained for , the curvature and its Hilbert transform are computed by inverting the corresponding Fourier coefficients with the standard FFT, and the analytic envelope is evaluated in . Lastly, we transform back to the Fourier space, apply the Gaussian filter, and finally interpolate the resulting to the non-uniform grid using the Type-2 NUFFT. Since each step involves and operations only, additional costs for our mesh refinement strategy are .

4.3 Temporal discretization and filtering

We evolve the dynamical variables using the classical fourth-order Runge-Kutta method. With surface tension, it is well-known that high-order terms in introduce stiffness into the governing equations and the stepsize is required to satisfy a severe stability constraint. In the case of axisymmetric vortex sheets, this constraint is known to be , where is the minimum distance between computational points in the symmetry plane. In this paper, we follow NS04 and try to satisfy the condition with as long as possible. For a general time-independent , a semi-implicit scheme that remove the stiffness has been developed by HLS97, and an extension of this method to more general cases is a direction of our future works.

For the approximation (3.25), we choose at the -th Runge-Kutta stage with , , and . Combining these coefficients with the property (3.24), we can perform time-stepping by the classical Runge-Kutta method without any interpolation technique in the temporal direction. It is not clear that the coefficients used in NS04 are the same as our choice, and there may be some alternative values. A drawback of this procedure is that it amounts to introducing a highly oscillatory delay inevitably associated with and may affect the convergence of our temporal discretization. This issue is carefully discussed in Section 5.

Also, the spectral filter introduced by Krasny [10] is applied to numerical solutions at every time step. This technique removes all Fourier coefficients whose magnitudes are below a cut-off level , and it is originally used to control unphysical growths of round-off errors due to the Kelvin-Helmholtz instability in the zero surface tension case. In the present case, another source of noises is the principal curvature , which is a quotient of two functions that vanish on the axis of symmetry. Since Eq.(2.11) involves the derivative of this function, large noises can be introduced to numerical solutions and grow under, for example, the Kelvin-Helmholtz instability. The Krasny’s filter is also expected to control the latter problem, as demonstrated in NS04.

5 Numerical results

Now we show some numerical results to validate the numerical method described in Section 4. In the following computations, the parameters for our mesh refinement scheme are given by


with the relative accuracy for the NUFFT algorithms. The other parameters , and are specified in each subsection.

5.1 Construction of uniform parametrization

As a limitation of the definition (3.21), it is required to numerically construct the uniform parametrization from a given initial condition before any time-stepping is performed. Fortunately, this initialization can be done by an immediate application of the ideas described in Section 4. Remember that the formula (4.2) is the trapezoidal rule applied to a smooth periodic function of the variable , which is known to be spectrally accurate. On the other hand, the arclength parametrization is a special case of the uniform ones whose domain is . Hence, by inverting the approximate coefficients (4.2) via the standard FFT, we expect to obtain a function represented in the uniform parametrization with spectral accuracy. To our knowledge, this procedure is distinctly different from the existing methods designed for the same purpose. For instance, Baker and Shelley [25] suggest an iterative scheme based on the Newton’s method, which searches for a set of values corresponding to equally spaced in . In another direction, Mikula and Ševčovič [26] propose a “time-evolving” method based on a concept called asymptotically uniform parametrization. The main idea is to fix and evolve a plane curve solely by the tangential velocity with which the relative spacing converges to a constant. Our method needs neither the Newton’s method nor temporal discretization, and it is completely non-iterative.

To justify the argument above, our initialization algorithm is tested on the following example in polar coordinates :


Here, is the polar angle, is the radius as a function of , and is the Legendre polynomial of order two. This special choice, which is presented in [27], corresponds to the zero-velocity state of a periodic solution to the linearized problem in [28]. The representation in the polar coordinates allows to compute the “inverse” mapping and check the accuracy of numerical results against the analytical data (5.2). Following HLS94, we first compute the curvature and the relative spacing from (5.2), and find the representation of in the uniform parametrization by our initialization algorithm. Then, with , the derivative is integrated to obtain , which is followed by the reconstruction of . In Fig.6, the plane curve (5.2) and the convergence of the scheme are shown for .

Figure 6: (a) Interface position of test problem (5.2) in plane, (b) exponential convergence of initialization algorithm based on inverting approximate coefficients (4.2).

As shown in Fig.6(a), it has the variable radius and therefore its parametrization is non-uniform. Next, in Fig.6(b), numerical errors in and are plotted versus in terms of the –norm divided by the maximal value of each function. The exponential convergence of the trapezoidal rule is clearly visible, and each graph reaches a level slightly above for . In fact, we do not perform a construction of the uniform parametrization for the initialization purpose, because the initial conditions (2.16) and (2.18) are uniform by definition. Instead, our initialization algorithm is effectively used in the convergence study of the time-stepping with the approximation (3.25), which requires to remove effects of mesh refinement for comparisons between two solutions with different .

5.2 Capillary pinch-off

Our numerical method is first tested on the pinch-off problem (2.16) with , , and .

Figure 7: Snapshots of interface position at indicated times for pinch-off problem (2.16).

Fig.7 plots the numerical solutions in the plane at indicated times, including the symmetric image in . To clearly show arrangements of computational points, in each figure the solid line is drawn at the the full resolution and every 16th point is plotted again as a dot. At , the droplet starting from a sphere is axially deformed and acquires a capsule-like shape. Since it is still relatively close to a sphere, the effect of mesh refinement is only slightly visible as a sparseness near . At , the capsule-like droplet is further elongated along the

–axis, and now it has developed two round ends connected by a thick rod. At this moment, the rod in the middle is almost straight near

, while geometry is relatively simple around the end-points and . For these simplicities, computational points are beginning to accumulate near the junctions between the rod and the round ends, as clearly seen in the corresponding figure. Almost the same argument applies to the solution at . Here, the radius of the rod is no longer uniform, and two isolated points of the minimum neck radius are formed near the junctions, implying where pinch-off is about to occur. At , which is just before noises become visible in the numerical solution, the minimum neck radius has rapidly decreased and the droplet is forming two end-droplets at the top and bottom along with a sharply elongated satellite droplet in between. Also, as discussed in Section 2, overturning is found in each pinch-off region, and it develops the so-called cone-crater structure carefully studied by NS04.

For reference, this GPU-accelerated simulation took approximately 1.25 hours with a single CPU core and one hour with two CPU cores. We could not find significant speed gain from further increasing the numbers of cores, which is presumably because the orders of and were relatively small compared to the main target of the parallelized NUFFT algorithms (see [5] for details). The same simulation was performed with the uniform parametrization and took roughly 10 minutes, although the numerical solution at around was contaminated by noises at the highest wavenumber.

Figure 8: Interface profile for pinch-off problem (2.16) at : (a) closeup of pinch-off region with every 4th point plotted as dot, (b) plots of and .

In Fig.8(a), we show a closeup of the upper pinch-off region at . As easily seen in the figure, a small region containing the shrinking neck is densely resolved by our mesh refinement scheme, and the minimum value of is six times smaller than that of the uniform parametrization. This refinement is roughly 1.3 times weaker than that of HLS97, and two times weaker than that of NS04. Nevertheless, considering that our method is more general and fully automatic, we believe that this result is successful. In Fig.8(b), the curvature and the corresponding are shown for the same . As expected, the procedure outlined in Section 4 successfully generates a simple and strictly positive from the curvature with two significantly sharp peaks. Here, note that the smoothness of can be controlled by the parameter of the heat kernel (3.19), and the allowable range of actually depends on the number . Relations among the parameters in our scheme are a subject of further investigation.

Figure 9: Confirmation of scaling law (2.17) for pinch-off problem (2.16): (a) vs , (b) vs .

Next, we confirm the scaling law (2.17) using our numerical results. The procedure is essentially the same as in NS04. Firstly, in Fig.9(a), we plot versus and fit a straight line to the data points. As one can easily see, the fitted line clearly explains the linear behavior of , which strongly suggests that the solution is in the self-similar regime. Secondly, by extrapolating the line to

, we obtain an estimate

, which is slightly smaller than obtained in NS04. Thirdly, the coordinate is plotted in Fig.9(b) as a function of using the value of estimated at the previous step, and then another straight line is fitted to the data points. Although it is not as clear as for , it is still reasonable to conclude that the results are well explained by the fitted line. Again, by extending the line to , we obtain an estimate , while the estimate by NS04 is . The deviations of our estimates from those in NS04 are arguably due to the gaps between the last data points and the end-points of the fitted lines. However, we believe that it is safe to claim that our results are slightly better than the computation by NS04 with the uniform parametrization and .

Figure 10: Convergence of time-stepping with approximation (3.25): (a) plots of errors measured as relative –norm, (b) correspondence between stepsize and filter level .

5.3 Convergence study

NS04 justifies the controversial approximation (3.25) by showing that changing the resolution or the parameter in (3.7) does not affect the results in their figures similar to Fig.(9). In this paper, we aim to verify in a more direct way that the time-stepping involving (3.25) actually converges faster than . Before doing any special trick, we first perform a standard convergence study with the pinch-off problem (2.16) and the time (see Fig.7), in which solutions with various are compared in against an “exact” solution obtained with a sufficiently small stepsize. Here, to weaken the stiffness , the number is reduced to 256, and the “exact” solution is computed with and . For large that still violate the constraint, the instability typically appears as a rapid growth in the Fourier coefficients of numerical solutions at high wavenumbers. This fact is numerically demonstrated by HLS94 for two-dimensional vortex sheets with surface tension, and the same phenomenon should be seen in the axisymmetric case. To continue simulations for such large , we control the growth of the Fourier coefficients by adjusting the filter level . In Fig.10(a), we plot errors in the mapping measured as the –norm divided by the maximum value of , and the correspondence between and is shown in Fig.10(b). As seen in Fig.10(a), the error curve for (denoted by “Variable” in the legend) is parallel to a linear function (blue line) and the convergence rate is apparently . From this result, one may be tempted to immediately conclude that the order of the classical Runge-Kutta method is completely lost. However, noting that a change in the stepsize can alternate the relative spacing linearly through the approximation (3.25) with , we still hesitate to eliminate the possibility that the true convergence rate may be hidden behind the behavior of (3.25).

From this point of view, we repeat the same convergence study after eliminating the effects of our mesh refinement scheme. That is, the computed solution for each is transformed to its arclength parametrization by the initialization technique based on inverting (4.2), and then errors are measured again as the relative –norm in the same manner as above. This convergence study is geometric in the sense that the difference between two plane curves is considered in a canonical form, which is the arclength parametrization in the present case. Again in Fig.10(a), we show some error curves obtained after removing potential artifacts due to the approximation (3.25). Firstly, the error curve computed from the regridded solutions (denoted by “Regrid”) is plotted with a theoretical curve (purple curve). As opposed to the result above, the convergence rate here is clearly superlinear and estimated to be at least , although some order reduction may be occurring presumably due to the fast oscillation of . Also, to check that the numerical solutions with our mesh refinement converge to the same curve as with the uniform parametrization, we plot the error curve (denoted by “TrueUni”) obtained from comparisons between the regridded solutions and an “exact” solution with the tangential velocity (3.6). The two dotted curves are indistinguishable except for the left-most data point, which strongly suggests that the numerical method described in Section 4 converges to a correct solution faster than the order of the approximation (3.25).

5.4 Bag breakup

Next, the performance of our numerical method is further examined in its application to the axisymmetric bag breakup (2.18) with the same , , and as for the pinch-off problem. Again, Fig.11 plots the numerical solutions in the plane (including ) at indicated times, and every 16th point is plotted as a dot after drawing the solid line at the the full resolution.

Figure 11: Snapshots of interface position at indicated times for axisymmetric bag breakup (2.18).

Also, since we do not evolve the end-point in this simulation, the interface position in each figure is reconstructed with . At , the droplet, which is initially a sphere, has developed a small indentation at the top and the inner region is no longer convex. The effect of mesh refinement is not evident here because no complex geometry is found in the interface. At , as the indentation further develops, the droplet begins to acquire a bowl-like shape whose edge is bent toward the axis of symmetry. Now, computational points weakly cluster near the two round corners of the edge, and consequently some sparseness is found away from those small regions. At , a large volume of the inner fluid has flowed into the edge region and the bottom of the bowl has become thinner. At , which is just before the numerical solution is contaminated by noises, the fluid interface is about to collide itself at the point connecting the edge and the bottom. As one can see, most of computational points are now devoted to resolving the two corners of the edge and the small region around the self-intersection.

The topological change shown in Fig.11 should be considered qualitatively different from that in Fig.7. In the case of bag breakup, the droplet is about to break into one torus-like component and one simply-connected component, whereas the symmetric pinch-off seems to form three simply-connected ones. In Fig.12(a), we show a closeup of the region near the point of the self-intersection at .

Figure 12: Interface profile for axisymmetric bag breakup (2.18) at : (a) closeup of self-intersection region with every 4th point plotted as dot, (b) plots of and .

As opposed to capillary pinch-off, overturning or a cone-crater structure is not observed in this case, and the fluid interface has two corners facing each other. This structure is more similar to those found in, for example, HLS97 and Burton and Taborek [4]. For the same , the curvature and the corresponding are shown in Fig.12(b). In this case, the solution forms two high-curvature regions away from the point of the self-intersection, which precludes applications of mesh refinement based on the prescribed guideline (3.7). As expected, our method succeeds in generating a simple and strictly positive from the curvature with four peaks with different magnitudes. Here, the relative spacing is roughly two times smaller for the smallest peak than that of the uniform parametrization, and 3.5 times smaller for the others three. On the other hand, a drawback of mesh refinement based on the analytic envelope is found in the interval between the two right-most peaks in Fig.12(b). In the function , the two adjacent peaks are “bridged” and our mesh refinement scheme overestimates the complexity of a region where the curvature is small and smooth (see Fig.11). This phenomenon is not due to the effect of the Gaussian filter. In fact, this bridged function is found “before” applying any smoothing filter, and therefore inherent to the use of the analytic envelope.

6 Conclusion

In this paper, we have proposed a new mesh refinement scheme for the vortex sheet formulation of inviscid droplet dynamics with axial symmetry. The method is based on signal processing using the Fourier coefficients of the curvature in the plane as a function of the arclength. A combination of the analytic envelope and the subsequent Gaussian filter enables to automatically detect two kinds of singularity formations, and additional computational costs for our mesh refinement are with the help of the non-uniform fast Fourier transform. We have also proposed a non-iterative construction of the uniform parametrization for closed plane curves. It is used to confirm, at least in a geometric sense, a convergence of the time-stepping procedure with the approximation (3.25) introduced by NS04.

There are two important directions for future developments. Firstly, as discussed in Section 5, a hanging bridge between two peaks of the guideline function , which is found in the case of bag breakup, results in some inefficiency in our mesh refinement scheme. This drawback motivates us to seek an alternative technique in signal processing that replaces the analytic envelope employed in this paper. Also, for the smoothing step, it may be interesting to switch from the heat kernel to other window functions such as the Kaiser‐Bessel family [29]. Secondly, an adaptive stepsize control is desirable in simulations of vortex sheets with surface tension. In addition to the small time-scales near self-intersections, HLS94 reports that the “shadow” of the Moore’s singularity and its saturation will lead to fast changes in solutions and can cause sudden drops in overall accuracy. Moreover, Ambrose [30] points out using their model problem that a nearly-singular behavior unrelated to any self-intersection can appear even in the non-zero surface tension case. Unfortunately, it is not trivial to change the stepsize adaptively in our current method because of the formula (3.25) with . To overcome this difficulty, it may be a key to gain more insights into the approximation as a delay term in the formulation. These projects are currently in progress.


We would like to thank Jon Wilkening for providing us with helpful comments and computational resources at Berkeley. We are also grateful to Toshio Aoyagi for providing us with computational environments at Kyoto. This research was partly conducted in the semester program “Singularities and Waves in Incompressible Fluids” at The Institute for Computational and Experimental Research in Mathematics (ICERM). We acknowledge the support from JSPS Overseas Challenge Program for Young Researchers (No. 201880131).


  • [1] T.Y. Hou, J.S. Lowengrub, and M.J. Shelley. The long-time motion of vortex sheets with surface tension. Phys. Fluids, 9:1933–1954, 1997.
  • [2] D. Leppinen and J.R. Lister. Capillary pinch-off in inviscid fluids. Phys. Fluids, 15:568–578, 2003.
  • [3] M. Nitsche and P.H. Steen. Numerical simulations of inviscid capillary pinchoff. J. Comput. Phys., 200:299–324, 2004.
  • [4] J.C. Burton and P. Taborek. Two-dimensional inviscid pinch-off: An example of self-similarity of the second kind. Phys. Fluids, 19:102109, 2007.
  • [5] A.H. Barnett, J. Magland, and L. af Klinteberg. A parallel non-uniform fast Fourier transform library based on an “exponential of semicircle” kernel. SIAM J. Sci. Comput., 2019. in press.
  • [6] Q. Nie. The nonlinear evolution of vortex sheets with surface tension in axisymmetric flows. J. Comput. Phys., 174:438–459, 2001.
  • [7] R.E. Caflisch and X. Li. Lagrangian theory for 3D vortex sheets with axial or helical symmetry. Transport Theory Statist. Phys., 21:559–578, 1992.
  • [8] D.W. Moore. The spontaneous appearance of a singularity in the shape of an evolving vortex sheet. Proc. R. Soc. Lond. A., 365:105–119, 1979.
  • [9] T.Y. Hou, J.S. Lowengrub, and M.J. Shelley. Removing the stiffness from interfacial flows with surface tension. J. Comput. Phys., 114:312–338, 1994.
  • [10] R. Krasny. A study of singularity formation in a vortex sheet by the point-vortex approximation. J. Fluid Mech., 167:65–93, 1986.
  • [11] N.D. Robinson and P.H. Steen. Observations of singularity formation during the capillary collapse and bubble pinch-off of a soap film bridge. J. Colloid and Interface Science, 241:448–458, 2001.
  • [12] G.K. Batchelor. An Introduction to Fluid Mechanics. Cambridge University Press, Cambridge, 1967.
  • [13] J.C. Burton and P. Taborek. Role of dimensionality and axisymmetry in fluid pinch-off and coalescence. Phys. Rev. Lett., 98:224502, 2007.
  • [14] D. Gabor. Theory of communications. J. Inst. Elect. Engrs. (Pt. III), 93:429–457, 1946.
  • [15] C. Muscalu and W. Schlag. Classical and Multilinear Harmonic Analysis, Volume 1. Cambridge University Press, Cambridge, 2013.
  • [16] S.O. Rice. Envelopes of narrow-band signals. Proc. IEEE, 70(7):692–699, 1982.
  • [17] C. Bauer, R. Freeman, T. Frenkiel, J. Keeler, and A.J. Shaka. Gaussian pulses. J. Magn. Reson., 58:442–457, 1984.
  • [18] M. Nitsche. Axisymmetric vortex sheet motion: accurate evaluation of the principal value integral. SIAM J. Sci. Comput., 21(3):1066–1084, 1999.
  • [19] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions. Dover Publications, New York, 1965.
  • [20] J. Cheng, M. Grossman, and T. McKercher. Professional CUDA C Programming. John Wiley & Sons, Inc., 2014.
  • [21] Q. Nie and G. Baker. Application of adaptive quadrature to axi-symmetric vortex sheet motion. J. Comput. Phys., 143:49–69, 1998.
  • [22] M. Nitsche. Singularity formation in a cylindrical and a spherical vortex sheet. J. Comput. Phys., 173:208–230, 2001.
  • [23] L. Greengard and J. Lee. Accelerating the nonuniform fast Fourier transform. SIAM Review, 46(3):443–454, 2004.
  • [24] J. Keiner, S. Kunis, and D. Potts. Using NFFT3—a software library for various nonequispaced fast Fourier transforms. ACM Trans. Math. Software, 36(4):19, 2009.
  • [25] G.R. Baker and M.J. Shelley. On the connection between thin vortex layers and vortex sheets. J. Fluid Mech., 215:161–194, 1990.
  • [26] K. Mikula and D. Ševčovič. Computational and qualitative aspects of evolution of curves driven by curvature and external force. Comput. Vis. Sci,, 6:211–225, 2004.
  • [27] H. Lamb. Hydrodynamics. Cambridge University Press, Cambridge, 1932.
  • [28] J.W.S. Rayleigh. On the capillary phenomena of jets. Proc. R. Soc. Lond., 29:71–97, 1879.
  • [29] J.F. Kaiser. Digital filters. In F.F. Kuo and J.F. Kaiser, editors, System Analysis by Digital Computer, pages 218–285. New York: Wiley, 1966.
  • [30] D.M. Ambrose. Singularity formation in a model for the vortex sheet with surface tension. Math. Comput. Simul., 80:102–111, 2009.