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.  (henceforth HLS97). In an axisymmetric case, Leppinen and Lister  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  (henceforth NS04) with additional insights into the nature of the self-similarity. Furthermore, a recent work by Burton and Taborek  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  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  (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.
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 ofpointing 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.
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)
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 . 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.  (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 . 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  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.
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 . 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  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 .
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 , Burton and Taborek  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  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  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:
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 
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 . 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 
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 .
The signal is a superposition of two Gaussian pulses  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:
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  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 , 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  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.
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 , NFFT3 , and FINUFFT . 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  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  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č  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 , corresponds to the zero-velocity state of a periodic solution to the linearized problem in . 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 .
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 .
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  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.
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.
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 .
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.
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 .
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 . 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.
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 . 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  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).
-  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.
-  D. Leppinen and J.R. Lister. Capillary pinch-off in inviscid fluids. Phys. Fluids, 15:568–578, 2003.
-  M. Nitsche and P.H. Steen. Numerical simulations of inviscid capillary pinchoff. J. Comput. Phys., 200:299–324, 2004.
-  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.
-  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.
-  Q. Nie. The nonlinear evolution of vortex sheets with surface tension in axisymmetric flows. J. Comput. Phys., 174:438–459, 2001.
-  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.
-  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.
-  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.
-  R. Krasny. A study of singularity formation in a vortex sheet by the point-vortex approximation. J. Fluid Mech., 167:65–93, 1986.
-  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.
-  G.K. Batchelor. An Introduction to Fluid Mechanics. Cambridge University Press, Cambridge, 1967.
-  J.C. Burton and P. Taborek. Role of dimensionality and axisymmetry in fluid pinch-off and coalescence. Phys. Rev. Lett., 98:224502, 2007.
-  D. Gabor. Theory of communications. J. Inst. Elect. Engrs. (Pt. III), 93:429–457, 1946.
-  C. Muscalu and W. Schlag. Classical and Multilinear Harmonic Analysis, Volume 1. Cambridge University Press, Cambridge, 2013.
-  S.O. Rice. Envelopes of narrow-band signals. Proc. IEEE, 70(7):692–699, 1982.
-  C. Bauer, R. Freeman, T. Frenkiel, J. Keeler, and A.J. Shaka. Gaussian pulses. J. Magn. Reson., 58:442–457, 1984.
-  M. Nitsche. Axisymmetric vortex sheet motion: accurate evaluation of the principal value integral. SIAM J. Sci. Comput., 21(3):1066–1084, 1999.
-  M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions. Dover Publications, New York, 1965.
-  J. Cheng, M. Grossman, and T. McKercher. Professional CUDA C Programming. John Wiley & Sons, Inc., 2014.
-  Q. Nie and G. Baker. Application of adaptive quadrature to axi-symmetric vortex sheet motion. J. Comput. Phys., 143:49–69, 1998.
-  M. Nitsche. Singularity formation in a cylindrical and a spherical vortex sheet. J. Comput. Phys., 173:208–230, 2001.
-  L. Greengard and J. Lee. Accelerating the nonuniform fast Fourier transform. SIAM Review, 46(3):443–454, 2004.
-  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.
-  G.R. Baker and M.J. Shelley. On the connection between thin vortex layers and vortex sheets. J. Fluid Mech., 215:161–194, 1990.
-  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.
-  H. Lamb. Hydrodynamics. Cambridge University Press, Cambridge, 1932.
-  J.W.S. Rayleigh. On the capillary phenomena of jets. Proc. R. Soc. Lond., 29:71–97, 1879.
-  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.
-  D.M. Ambrose. Singularity formation in a model for the vortex sheet with surface tension. Math. Comput. Simul., 80:102–111, 2009.