{'fraction_non_alphanumeric': 0.028569251866524457, 'fraction_numerical': 0.013979887246685967, 'mean_word_length': 4.4185758513931885, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Manifold death: the implementation of controlled topological changes in thin sheets by the signature method

A well-known drawback of the Volume-Of-Fluid (VOF) method is that the breakup of thin liquid films or filaments is mainly caused by numerical aspects rather than by physical ones. The rupture of thin films occurs when their thickness reaches the order of the grid size and by refining the grid the breakup events are delayed. When thin filaments rupture, many droplets are generated due to the mass conserving properties of VOF. Thus, the numerical character of the breakup does not allow obtaining the desired convergence of the droplet size distribution under grid refinement. In this work, we present a novel algorithm to detect and perforate thin structures. First, thin films or ligaments are identified by taking quadratic moments of the VOF indicator function. A multiscale approach allows us to choose the critical film thickness independently of the mesh resolution. Then, the breakup is induced by making holes in the films before their thickness reaches the grid size. We show that the method improves the convergence upon grid refinement of the droplets size distribution and of enstrophy.


Introduction
Multiphase fluid mechanics with sharp interfaces involves diverse singularities each with its own difficulties. In this paper we address the change of topology happening when a hole forms in a thin liquid sheet and in particular its numerical modelling in Volume-Of-Fluid (VOF) methods. This topology change has important consequences for fluid fragmentation, and in particular for the droplet size distribution. Several mechanisms for the topology change of thin sheets have been proposed [1,2] and a full description of physical mechanisms is still missing. In some cases, inter-molecular forces may be shown to have a destabilizing effect on thin sheets, while in other cases, as in some soap films, permanent dipoles may form on interfaces having a stabilizing effect [3]. In this paper we are concerned with interfaces that do not display such a stabilizing effect and break before reaching molecular length scales.
Mathematically, the Navier-Stokes equations with the assumption of sharp interfaces between phases are incomplete as they do not describe the piercing of thin sheets. These equations can only describe the indefinite thinning of a sheet when it is compressed or stretched by the flow. Eventually, it may happen that the sheet becomes thinner than molecular scales, at which point the equations cease to be valid. However, before this happens, molecular forces may induce the breakup of films [4,5,6,7,8]. In many cases, however, and in particular in atomization or floating bubble experiments, observations show that sheets are pierced well before their thickness reaches molecular length scales [9,10]. Thus, an ad-hoc prescription for piercing thin sheets at some macroscopic length scale is needed to describe real flows.
In addition to the physical modelling challenge, the numerical simulation of thin sheets and their breakup poses significant challenges. In fact, when the fluid structures become too thin they can not be correctly resolved neither with the Volume-of-Fluid (VOF) or the Level-Set method. The mass conservation property of the VOF method is responsible for the breakup of thin sheets with thickness of the order of the grid size and for the formation of several droplets, while level set schemes continuously evaporate mass [11]. To avoid the occurrence of artificial breakup using VOF, extremely fine grids must be used, such that ∆ < h/6 (where h is the film thickness) and huge computational efforts are required even using codes with adaptive mesh refinement (AMR) capabilities [12,13]. In the other case ∆ > h/6, the artificial breakup of thin sheets does not permit to observe convergence under grid refinement of some meaningful quantities such as droplets size distribution or enstrophy. Automatic reconnection algorithms have been suggested for thin sheets in the Front-Tracking method [14]. Some algorithms have also been suggested to make "numerical" reconnection impossible for the VOF method [15]. We are not aware of any similar mechanisms (either promoting reconnection or blocking it) for the level set method.
In this paper we propose a manifold death algorithm to perforate sheets or ligaments of a given thickness in a controlled way in the VOF method, thus avoiding the purely numerical breakup caused by the grid-size and the associated grid convergence issues. First, the thin regions are detected by computing quadratic forms based on the values of the color function. The sign of the eigenvalues of the quadratic form, known as signature, is then used to identify thin regions. Once their position is known, we create artificial holes that expand and "destroy" the thin regions.
The rest of the paper is organized as follows. The mathematical model is defined in Section 2. The manifold death algorithm is described in Section 3. The numerical results obtained with the manifold death algorithm are discussed in Section 4, followed by our conclusions.

Governing equations
The governing equations for the incompressible two-phase flow with immiscible fluids are written using the one-fluid formulation. When using the VOF method, the volume fraction (or color function) is used to identify the location of the interface such that C = 1 in fluid 1 and C = 0 elsewhere. Then, in the one-fluid formulation, the value of a generic material property p is a function of C, namely p = Cp 1 + (1 − C)p 2 . The volume fraction is obtained by solving the following advection equation The velocity u and pressure p are obtained by solving the Navier-Stokes equations that express the conservation of mass and momentum where ρ and µ are the fluid density and viscosity, respectively. The gravitational force is taken into account with the ρg term. Surface tension is modeled with the term f = σκnδ S , where σ is the surface tension, κ denotes the curvature of the interface, n is the unit vector normal to the interface and δ S is the surface Dirac distribution on the interface (i.e. zero everywhere except on the interface). The open-source code Basilisk is used to solve the two-phase incompressible Navier-Stokes system with quad/octrees on adaptive Cartesian meshes. The efficient adaptive implementation of Basilisk allows to solve with a high resolution only in the relevant parts of the domain, thus reducing the computational cost of the simulation. A momentum-conserving scheme is used for the velocity and a piece-wise linear geometric VOF method is adopted for tracking the interface. More details about the code can be found in [16,17,18] and references therein.

The manifold death algorithm
In this section we describe the manifold death method used to perforate thin structures in a controlled way.

The signature method
The first step of the manifold death algorithm consists in detecting the thin sheets (or ligaments) in the domain. To do that, we use the following signature method : 1. Consider a point x 0 ∈ R 3 and translate the coordinate system so that the new origin is x 0 = 0.
Consider a radius R approximately the same size of the sheet thickness h c one wants to detect, see Figure 1, and the bilinear form The quadratic moments T ij on a sphere S of radius R can be found by integrating where φ = 2C − 1 is a symmetric indicator function. The numerical integration can be done on the volume V obtained by taking a shell of thickness δ around S. 3. After orthonormalization of the quadratic form, one finds a new set of coordinates in which f (X, X) = i X 2 i , where i are the eigenvalues of the operator with matrix T ij . The number of positive, negative and zero values of i is the signature s of the quadratic form. For the given point x 0 , the signature s indicates the average shape of the interface in the vicinity of the point and can be used to determine whether the point is in the bulk of the phase, in a thin sheet (or ligament) or close to the interface, see Figure 1. In Table 1 we reported some of the signatures that can be obtained using this method.

Signature s
Position of x 0 (+, +, +) (+, +) Bulk of the phase The eigenvalues can be found by solving the equation that in two dimensions gives 1,2 = tr(T) ± (tr(T) 2 − 4∆)/2. However, we decided to rely on the GNU Scientific Library for both two and three dimensional cases [19].  For the numerical computation of the quadratic moments in (4), some simplifications and adjustments re necessary to tailor the method to the Basilisk code. For simplicity, we first replace the sphere with a cube centered in x 0 of size L = 2R. Moreover, in order to take advantage of Basilisk's capabilities, we make sure that the cube lies completely within the stencil (i.e. the set of surrounding cells) around x 0 . Since the standard stencil size in Basilisk is five cells (5∆), we can detect thin sheets whose thickness is approximately three cells (3∆). To identify larger sheets (h c > 3∆) we take advantage of the multilevel nature of Basilisk by computing the signature on the appropriate coarser grid. The restriction operator for the indicator function φ is the standard volumetric average, while the prolongation of the signature onto finer grids is done with the injection operator. Therefore, the critical thin sheet thickness can be set independently from the mesh resolution. Figure 2 illustrates the multilevel approach used for detecting thin structures. Using the finer grid (left), the distance between the two interfaces is larger than h c ≈ 3∆ f ine and the shaded cell is marked as "bulk of phase", while on the coarser one a thin structure of thickness h c ≈ 3∆ coarse = 6∆ f ine is detectable.
The thin sheets identified with the signature method are shown in Figure 3. The figure clearly shows that the method is symmetric. In fact, the algorithm detects a thin sheet in the reference phase in the bottom, while, in the two other regions above, the reference phase 1 encloses thin sheets of the phase 0. To summarize, with this method we can identify thin regions in both phases. We would like to remark that the algorithm is compatible with the adaptive-grid approach of Basilisk and then it is suitable for large parallel computation. Our observations suggest that it is much faster than the actual multiphase solver, so it can be used with a negligible computational overhead. Finally, the signature method can be used also for controlled numerical adaptive grid refinement to enforce refinement in thin regions, see [20].

The topology changes
Once we have identified the thin sheets, we want to perforate them before their thickness reaches the cell size and the numerical breakup happens. When the previously detected thin regions is in the reference phase (i.e. C > C * ), we change the topology by placing some cubic "holes" setting the color function C to zero. Otherwise, to force interface reconnections when the thin region is in the other phase (C < C * ), we set the color function C to 1. In our simulations we take C * = 0.3.
An important remark is that the size of the holes (made in both phases) is critical, in fact holes will expand only if their diameter is at least the sheet thickness h c , otherwise they will re-close to minimize the energy of the system [21]. In order to recover the desired grid independency, in addition to the thickness h c , also the size of the holes has to be fixed and independent of the cell size. To this aim, we create cubic holes of the same length h c used for the computation of the signature. By choosing the minimum size for which the expansion criterium is satisfied, we are minimizing the perturbations on the system as well as the amount of mass removed or added. Moreover, the manifold death algorithm is used only at fixed time interval t h to give the holes some time to expand with Taylor-Culick velocity. Finally, we limit the maximum number of holes n h created at every iteration of the algorithm.

Results
In this section we present our results. We apply the manifold death method to a simple single vortex test, to an axisymmetric secondary atomization situation and to a phase inversion problem.

Single vortex
The single vortex test has been designed to check the ability of interface tracking methods when the reference phase undergoes large deformations and thin structures of thickness of the same order of the the grid size are involved, see [22,23]. The divergence-free velocity u = (u x , u y ) is obtained from the following stream function ψ = 1 π sin 2 (πx)sin 2 (πy) , as u x = ∂ψ/∂y and u y = −∂ψ/∂x. The domain is a unit square with the bottom left corner placed in the origin and a circular droplet with radius r = 0.15 is centered at (0.5, 0.75). On the sides of the box, homogeneous Dirichlet boundary conditions are imposed. The rotation ends at T = 3s. This simple test is designed to check, in particular, the first step of the manifold death algorithm, namely the detection of thin structures where the artificial breakup occurs. We first perform the simulation on a coarse mesh, where many artificial breakup events happen. Then, we check if the manifold death method is able to identify and perforate the ligament as soon as it reaches a thickness equal to the coarse grid size. In other words, the critical ligament thickness h c used in the algorithm has to be equal to the grid size ∆ l of the coarse mesh. This requires ∆ l = h c ≈ 3∆ l+1 , where ∆ l+1 indicates the finer grid size. However, since using the standard multilevel ∆ l = 2∆ l+1 holds true and since Basilisk uses quadtree grids, which restricts the resolution to powers of two, the simpler way to satisfy this requirement is to increase the coarse grid size ∆ l by taking a square domain of size 1.5 (note that this would be equivalent to using a 47 2 grid on the unit square). The criterion used for the mesh adaptation is the error on the volume fraction field C with a threshold equal to 1 · 10 −5 . The frequency of the holes creation is t h = 0.05s and the maximum number of holes per manifold death call is n h = 10.
We show in Figure 4 the interface at different grid resolutions with and without manifold death. The four sub-pictures correspond to different times. In the top row of each picture we report the interface with no manifold death for the 64 2 (with larger domain) and 256 2 grids, while in the bottom row the manifold death method is used with h c = 0.023m = 3∆ 128 on the 128 2 and 256 2 grids. Note that, for clarity, the smallest droplets have been removed.
In the standard VOF, the breakup of the thin filament is clearly grid dependent, since with the coarser grid at the end of the simulation the reference phase is made of a set of droplets, while with the finer one the thin filament is still well resolved and only few droplets have been created near the end of the tail.    The validation of the manifold death method, and in particular the detection of thin structures, can be done qualitatively by looking at the position (indicated by the arrows) of the most recent breakup event.
By controlling the breakup, we have replicated on finer grids what happens numerically on a coarser grid, showing that the method accurately detects and perforates thin structures of given thickness.

Secondary atomization
Secondary atomization refers to the physical situation where the droplets originated from the primary atomization of a jet, due to the interaction with the ambient high speed gas flow, deform and fragment. Here, we study the evolution of a single droplet that experiences an impulsive acceleration of the surrounding gas. The droplet is initially at rest, surrounded by fluid with zero velocity. At t = 0 an impulsive velocity U 0 condition is imposed on the left boundary. The drop then stretches and deforms into a film whose shape resembles that of a bag. As the bag inflates, its thickness decreases until holes appear and the bag breaks up bursting the bag into a spray. An extensive analysis of the different fragmentation regimes that the droplets may experience has been carried out in [24,25,26]. The goal of this paper is not to further investigate the secondary atomization problem, instead we want to show the improvements in terms of grid convergence that can be attained using the manifold death method. To this aim, we select a single case from [26] with the intent to show that we are able to control the sheet breakup, while in the standard VOF simulations it happens numerically as soon as the sheet thickness reaches the minimal grid size.
The computational domain is a two-dimensional axisymmetric channel, represented by a square of side 15R 0 , where R 0 is the droplet initial radius used for the adimensionalization. The dynamics of the droplet can be uniquely determined by the Reynolds Re = ρ G U 0 R 0 /µ G and Weber We = ρ G R 0 U 2 0 /σ numbers, and the ratios of densities r d = ρ L /ρ G and viscosities r v = µ L /µ G . For this test we take Re = 1090, We = 7.5, r d = 1110 and r v = 90.9. We perform the simulations on two fine grids 8192 2 and 16384 2 (corresponding to the levels of refinement 13 and 14) with about 500 and 1000 cells per droplet radius, respectively. Adaptation is performed using a threshold of 10 −5 and 5 · 10 −5 on the velocity u and volume fraction C errors, respectively. For this test case we set t h = 0.1 and impose that only one hole is done during the whole simulation. To reduce the effects of the long initial transient where the droplet stretches and forms the forward-facing bag, we start the simulations on the finer grid by using a snapshot obtained at lower resolution. In Figure 5 we reported the evolution of the droplet until its breakup. The droplets t = 0.0s 20 40 60 80 100 115.5 deforms and a forward facing bag is formed. The breakup occurs at t = 115.0 at the lower resolution and at t = 117.5 when using the finer grid, see Table 2. The evolution is similar to the same solution as reported in [26], apart from the bulge near the symmetry axis that was not observed there and is now present. The breakup of the droplet is shown in Figure 6 for the two grid resolutions. The black line indicates the interface, while the gray background is the position of the deformed droplet in the case without controlled perforation and the same mesh resolution. From top to bottom, we report the results obtained reducing the critical thickness h c at which thin sheets are detected and perforated. In the top row we set h c1 = 3∆ 2048 = 0.022R 0 , in the middle one h c2 = 3∆ 4096 and, finally, in the bottom one h c3 = 3∆ 8192 . In Table 2 we report the time at which the first breakup event happens. By reducing the value of h c the bag breaks later no matter the mesh resolution used for the solution of the Navier-Stokes equations. Also, with a given critical thickness and comparing the results for the two resolutions, the breakup occurs in the same spot and almost at the same time, indicating that we successfully managed to control the breakup making it almost grid independent. To summarize, in all the cases with manifold death, the breakup happens earlier than in the reference case without controlled perforation and the anticipation time depends on the critical thickness h c .

Phase inversion problem
The situation of the phase inversion consists of two fluids, initially at rest, with the lighter one placed in the bottom of a box surrounded by the heavier one, see Figure 7. The outer cubic box has size H, while the inner one H/2. On the outer walls we impose a free-slip boundary condition, so that the normal velocity is zero and the tangential components obey a symmetry condition. A π/2 static contact angle is imposed on    the walls and the gravitational acceleration is g = (0, −9.81, 0). Due to gravity the lighter fluid raises and in the end will occupy the top part of the box. In between the simple initial and final configuration, it is possible to observe many breakup and coalescence events, which in standard VOF simulations are due to both numerical and physical aspects. In particular, the role that sheets breakup plays in the lack of convergence of enstrophy has been recently pointed out in [27]. The convergence upon grid refinement of primary moments, such as kinetic and mechanical energy, has already been shown in literature [28,29,30] and therefore here we do not provide a discussion about this aspect. In this work, however, we show that the manifold death method prevents the numerically induced breakup of thin structures and, as a consequence, the grid-convergence of enstrophy and droplets size that plagues many of the literature cases is improved. In this study, we consider moderate values of the Archimedes and Bond numbers, aiming to perform a true DNS multiphase simulation, with an accurate resolution of the fluid flow. The geometrical dimensions and the physical properties of the two fluids are reported in Table 3. Note that the properties of the lighter fluid resembles those of oil and those of the heavier one of water.  We introduce the following definitions for adimensional numbers and quantities The enstrophy in the reference phase (oil) is obtained with where ω = ∇ × u denotes the vorticity. The enstrophy is integrated over the smaller subdomain [0, 0.95H] × [0, 0.95H] × [0, H] to neglect regions where the dynamics is strongly affected by the presence of the walls of the box, see [27].  The simulations are carried out on 512 3 and 1024 3 grids and the critical thickness is h c = 3∆ 512 . The thresholds used for the refinement criteria are 5 · 10 −3 and 5 · 10 −4 for the velocity and volume fraction errors, respectively. Concerning the manifold death parameters, the adimensional time interval between consecutive holes formation is t * h = 0.03 and the maximum number of holes per iteration is n h = 240. In Figure 8 the improvements on the convergence of the oil enstrophy using the manifold death algorithm can be appreciated. In the standard case, the enstrophy peaks later using the 1024 3 grid than with the 512 3 one. Instead, when the manifold death algorithm is used to control the topology changes, the convergence upon grid refinement of enstrophy is obtained. The enstrophy ripples when using the manifold death algorithm are due to the creation of holes and their frequency is clearly 1/t * h . They last for a very short period of time and therefore do not affect the overall convergence. On the right of the same figure, we report the enstrophy difference E 1024 (t * ) − E 512 (t * ) for the uncontrolled case (Ref) and the one with manifold death (MD). At t * = 2.9 a large difference (almost 1/6 of the maximum enstrophy value) can be observed for uncontrolled case, with more enstrophy produced using the coarser grid, while after t * = 3.3 a higher enstrophy production can be observed with the finer grid. This is due to the disruption of the sheet shown in Figure 9 and without controlled perforation its onset depends on the grid size. When using the manifold death method, the difference is much smaller and enstrophy converges.
To confirm that the cause for the non convergence upon grid refinement of enstrophy is the breakup of thin sheets, we report in Figure 9 the interface at the dimensionless time t * = 2.9 for four cases. By comparing the two pictures obtained without controlled perforation (top row) the numerical essence of the breakup is clear. In fact, at lower resolution (left) the sheet is broken and many holes and ligaments can be seen, while at higher resolution (right) the sheet is well resolved since its thickness is larger than the grid size. On the contrary, when the manifold death algorithm is used (bottom row), at the same time the artificial holes have already destroyed the sheet in both cases and the breakup becomes grid independent. A close-up on the expansion of the holes created by the manifold death algorithm and the subsequent ligament formation is shown in Figure 10. Finally, we focus on the droplet size distribution and in particular on the number of larger droplets, since they can be expected to result from both physically dominated breakup and manifold death of thin structures. In Figure 11 we report the droplet diameter frequency N(d) as a function of the droplet diameter in the standard case (left) and with the use of manifold death (right). In the former, we recall that at t * = 3.25 the sheet is well resolved on the 1024 3 grid, while on the coarser one it has been destroyed by the numerical breakup. As a consequence, more droplets have been generated using the 512 3 resolution. At t * = 3.58, the sheet no longer exists also with the 1024 3 grid. The two larger droplets (d > 10 −2 m) have been correctly resolved on both grids, while the frequency of the smaller ones shows similar profiles, with a shift towards the left (smaller droplets) when the finer mesh is used. This confirms the numerical nature of the breakup. With controlled perforation, the number of trustworthy droplets on the right of the dashed line (d = 8∆) is similar for the two grids and at both times. Moreover, when the mesh size decreases, more small droplets (in the untrustworthy region) are again collected.

Conclusions
We have described a new method to detect thin structures and to create holes depending only on the value of the critical thickness h c . We have performed a simple single vortex test to prove the ability to identify thin ligaments of thickness h c independently on the grid resolution. The axisymmetric secondary atomization and the three-dimensional phase inversion tests show the improvements in terms of grid convergence that the method can provide to the study of complex multiphase simulations. The findings of this study suggest that using an algorithm to induce in a controlled way the breakup of thin structures is necessary to obtain grid independence of the droplet size distribution and second order moments such as enstrophy. Future work will concentrate on improving the numerical hole formation, for example by re-distributing the mass near the holes rather than removing it, to better conserve mass and energy.