Liquid Structure Classification Towards Topology Change Modeling

We present a method that, in combination with a two-plane liquid-gas interface reconstruction in each computational cell, identifies thin features in Eulerian simulations of spray atomization that should undergo imminent breakup. Rapid identification of these features, which include gas sheets between colliding droplets, liquid sheets generated in bag breakup, and ligaments that breakup through Rayleigh–Plateau instabilities, allows for the on-the-fly application of subgridscale (SGS) breakup models for each type of thin feature in order to accurately predict the resulting topology change behavior. Without the feature identification abilities of the proposed method, thin features would experience numerical breakup when their characteristic size falls below that of the computational mesh. Such topology change is mesh-dependent, and therefore, unphysical. The proposed method incorporates the interface normal and centroid into any underlying connected-component labeling (CCL) algorithm to distinguish thin features from other fluid structures. We apply the method on data from simulations of offset binary droplet collision and turbulent air-blast atomization to verify its ability to identify individual thin features and demonstrate its usefulness for calculating relevant properties of the identified features.


Introduction
Atomization involves complex events where the topology of liquid structures abruptly changes, such as the formation and breakup of ligaments and films and the coalescence of droplets. In Eulerian simulations of liquid-gas flows, the computational mesh size controls the occurrence of these topology changes, as the mesh size dictates the minimum size of a liquid or gas structure. As a result, a ligament will experience numerical breakup when its diameter falls below the mesh size, and two droplets will coalesce if the thickness of the gas film between them falls below the mesh size. Due to the multi-scale nature of spray atomization, simply increasing the mesh resolution to delay or prevent topology change becomes computationally intractable. In addition, mesh convergence in topology change behavior should not be expected for simulations using the Navier-Stokes equations alone, as those equations do not predict behavior at the molecular scales relevant to topology change. For numerical simulations to become more relevant to the fundamental research of underlying breakup processes, the collection of accurate spray statistics, and the development of spray control techniques, they should avoid unphysical, mesh-dependent topology change and instead, use a Lagrangian understanding of topology in combination with subgrid-scale (SGS) breakup modeling of films and ligaments to determine the occurrence of topology change events. In geometric volume-of-fluid (VOF) simulations, a Lagrangian representation of thin structures is enabled by the use of two planes to reconstruct the liquid-gas interface in each cell, as in the R2P (reconstruction with 2 planes) method of Chiodi and Desjardins [1,2]. R2P delays or even prevents numerical topology change by enabling the realization of thin films with thickness well below the computational mesh size. Figure 1 shows how R2P captures and retains the gas film formed during a binary droplet collision, while a single-plane interface (PLIC) [3], reconstructed according to ELVIRA [4], fails to resolve the gas film and causes the droplets to coalesce. Sev- The diameter to mesh size ratio is D/∆x = 24.2. Both images are taken at t = 0.1 ms. Single-plane interface elements are colored blue, while two-plane interfaces are colored red.
eral SGS topology change models exist in the literature [5,6], but in order to determine which computational cells within the domain these models should be applied to, it is necessary to identify features that constitute precursors to topology change events, such as liquid sheets, ligaments, and gas films. Such features can be generalized as connected regions of Eulerian field data where each cell in the region meets a criterion, such as a threshold liquid volume fraction. The labeling of such connected regions is related to a class of algorithms known as connected-component labeling (CCL). CCL is an application of graph theory and is used in the computer vision field for analysis of binary digital images [7,8]. In the spray atomization field, CCL has been widely used within experimental research for identifying and tracking liquid structures and their associated topology changes to gather statistics on breakup length and droplet size [9,10]. The first applications of CCL within numerical simulations of atomization identified Eulerian droplets for conversion to Lagrangian point particles in algorithms that couple the VOF method with Lagrangian particle-tracking [11,12]. Recent numerical studies of atomization and bubbly flows have applied CCL for the gathering of droplet and bubble statistics in a similar manner to that of the aforementioned experimental studies [13,14]. This work extends the use of CCL to the identification of liquid and gas films in a VOF simulation using two-plane interface reconstruction. We first verify the algorithm and its performance by performing CCL on prescribed volume fraction fields, then demonstrate its usefulness by applying it on atomization simulation data. Finally, the prospect of using the technique towards SGS modeling is investigated.

Method
The proposed method operates in conjunction with an underlying CCL algorithm to identify and label two types of features: detached structures (referred to as structures below) and thin structures such as films and ligaments (referred to as films below). While traditional CCL algorithms assume the existence of a binary or logical array to perform labeling on, our method creates two logical arrays-one for structures and the other for films-and performs distinct labeling operations on each. Since the proposed method only changes the criteria by which computational cells are connected and labeled, it can be applied to any CCL method; however, the ability to perform labeling cheaply for on-the-fly SGS modeling relies on the choice of an optimal CCL method as we discuss in the Section entitled Comparison of Labeling Algorithms. The output of the method is a list of structures and films along with the indices of the cells contained within each structure and film. The film phase is also identified in the output. Along with the application of SGS film models, the output could be used to calculate feature properties such as volume, centroid velocity, and mean film thickness. Specifically, we define a structure in a fluid domain as any connected liquid region bounded by one or more liquid-gas or liquid-solid interfaces. A film is a connected liquid or gas region where the mass contained within the neighborhood of any point is characterized by a diagonalized inertia tensor such that the magnitudes of the principal axes, I 1 , I 2 , and I 3 , are related as or depending on the shape of the film. For instance, a liquid sheet formed during bag breakup or a gas film between two colliding droplets would have inertia tensor components related by (2), while a ligament would have inertia tensor components related by (3). As we only seek to identify ill-resolved regions of fluid in the proposed method, however, we restrict our definition of a film as any liquid or gas connected region with a minimum dimension d < ∆ in any direction, where ∆ is the length spanned by two adjacent computational cells at the location of the fluid region. Note that a fluid region may be labeled as both part of a structure and film. For example, a liquid jet with a diameter and length of several cell lengths may have multiple thin ligaments emanating from its core. The proposed method would label the jet and ligaments as a single liquid structure while also labeling each ligament as a separate liquid film. We now introduce the nomenclature used for the method description and the assumptions necessary for the operation of the algorithm. We assume that a liquid volume fraction α i exists in each computational cell Ω i , and that for mixed-phase cells with 0 < α i < 1, the liquid-gas interface within the cell is represented by N i planes, where 1 ≤ N i ≤ 2. For each plane m in the liquid-gas interface representation of cell Ω i where 1 ≤ m ≤ N i , the normal vector is denoted by n i,m , and the centroid is denoted by c i,m . For simplicity, we assume that two-plane interface reconstruction is only used to resolve features with characteristic length below the grid size such that n i,1 · n i,2 < −0.5, and therefore, we classify all two-plane cells as film cells.
We also assume that no liquid and gas films are adjacent to each other. All interface normal vectors point from the liquid phase towards the gas phase. For a given cell Ω i , we define its cell neighborhood N i as the set of cells Ω j that share a face with Ω i . As an example, for a cell in a structured mesh with indices (i, j, k), its neighborhood N is the set of cells with indices (i ± 1, j, k), (i, j ± 1, k), and (i, j, k ± 1).
Each cell in the computational domain is assigned a structure label s i and a film label f i corresponding to the structure or film that the cell belongs to. If a cell does not belong to a structure or film, then its respective label s i or f i is 0. A structure S with label A is defined as where S is the set of cells that are a part of any structure. Likewise, a film F with label B is defined as where F is the set of cells that are a part of any film.

Structure Classification
The first step in the method is to create the logical array that defines the cells that are connected and labeled as part of a structure. An array value is "true" if and only if the corresponding cell is in the set S. We define S as the set difference where V is the set of all liquid-containing cells and G 2 is the set of two-plane cells containing a gas film as determined by the relative position of the interface centroids and the normal direction: The threshold λ can be set to zero, but is often set to a small non-zero value to avoid the labeling of cells with spurious flotsam and jetsam. We exclude cells in G 2 from structure labeling because two-plane gas film cells may bridge two liquid structures with differing structure labels s and therefore contain liquid volume from both structures. In order to account for the liquid volume inside a two-plane cell in G 2 , it is necessary to identify the liquid structures adjacent to the film along with the corresponding liquid volumes inside the film cell. As it is not strictly necessary for the SGS modeling of thin films, the identification of adjacent liquid structures and recovery of corresponding liquid volumes is considered outside the scope of this paper and will be addressed in future work.

Structure Labeling
After the set of structure cells S has been determined, each cell can be assigned a structure label. As with standard CCL algorithms, any two cells Ω i and Ω j must satisfy the criteria in order for them to share a structure label s i = s j . In addition to the criteria in (9), however, the proposed method requires that at least one of the following be true: where is the set of two-plane cells containing a liquid film.

Film Classification
To identify liquid and gas films, we create another logical array that defines the cells that are connected and labeled as part of a film. Similarly to the logical array for structures, an array value is "true" if and only if the corresponding cell is in the set F . We define F as the union of interfacial cell sets G 1 is union of the set of single-plane cells that neighbor a two-plane gas cell and the set of single-plane cells that neighbor another single-plane cell such that a gas film exists between their respective interfaces: L 1 is the corresponding union of sets for single-plane liquid film cells: Since the structure labeling and film classification steps share common logical tests, the two steps can be performed simultaneously for reduced computational cost.

Film Labeling and Phase Identification
Two cells Ω i and Ω j share the same film label This criterion is unchanged from a standard CCL algorithm. Finally, we obtain the phase of each film by defining the sets of gas film cells G and liquid film cells L: A film F B is a gas film if ∀Ω i ∈ F B , Ω i ∈ G, while it is a liquid film if ∀Ω i ∈ F B , Ω i ∈ L.

Results and Discussion
In this section, we demonstrate the usefulness of the proposed method for the identification of film features and the extraction of relevant quantities from these features. We first use the method to track the collision of binary offset droplets and the gas film formed between them. We then test the method on a turbulent air-blast atomization simulation. We implement the labeling algorithm and perform all simulations inside the NGA flow solver [15], and use the discretely conservative, unsplit geometric VOF method of Owkes and Desjardins [16] for advection of the liquid-gas interface. All meshes used are Cartesian with uniform grid spacing.

Collision of Offset Droplets
We simulate the collision of binary offset droplets following the experiments of Qian and Law, specifically case (r) of Figure 4 in their paper [17]. For binary offset droplet collisions, the additional relevant nondimensional parameter aside from W e, Re, and R is the impact factor B which is defined as where r ⊥ is the initial distance between the droplet centroids in the direction perpendicular to their initial velocities. Figure 2 shows the interface reconstruction and feature identification at times t = 0.00 and 0.25 ms along with a top-down view of the cells corresponding to the gas film, colored by film thickness. The use of R2P allows the gas film between the colliding droplets to be resolved for the majority of the collision time, after which R2P predicts coalescence due to the lack of an SGS model for the gas film, whereas the Qian and Law experiments show that the droplets bounce. While the gas film is resolved, the proposed method identifies each droplet with a unique structure label s and identifies the gas film cells with a single film label f . For each labeled gas film cell, the local film thickness h is calculated as where N i is the set of all cells that share a face, edge, or vertex with cell Ω i , V i is the cell volume, and A i is the surface area of the liquid-gas interface in the cell. In combination with a SGS model for the draining rarefied gas flow between two droplets, the local film thickness can be used to determine if the gas film should break, thereby coalescing the droplets.

Air-Blast Atomization
We next demonstrate the proposed method on a turbulent air-blast atomization simulation using co-flowing streams of low-speed water and high-speed air emanating from a canonical nozzle, shown in Figure 3a. The gas Reynolds number, Re g ≡ ρ g U g (d o − d i )/µ g , is 16709, while the liquid Reynolds number, Re l ≡ ρ l U l d l /µ l , is 1000, and the momentum flux ratio, M ≡ (ρ g U 2 g )/(ρ l U 2 l ), is 6.4. Here, ρ is the fluid density, U is the fluid velocity, µ is the fluid dynamic viscosity, d o is the gas outer diameter, d i is the gas inner diameter, d o − d i is the hydraulic diameter, d l is the liquid inner diameter, and the subscripts of g and l refer to the gas and liquid phases, respectively. The simulation is initially run using a single-plane PLIC-ELVIRA reconstruction, then switched to R2P after the flow has developed. A snapshot of the flow immediately before a bag breakup event is shown in Figure 3b, where each of the detected films, one liquid sheet and two ligaments, is displayed with a unique color corresponding to its film label f . Although the purple and orange films are a part of the same liquid structure that includes the jet from the nozzle, they are identified as distinct features from the jet because of their limited thicknesses. As with the gas film between colliding droplets, the method allows for calculations to be performed on individual film features. For example, the local thicknesses of the liquid sheet and ligaments can be used, in combination with an SGS breakup model, to determine the moment of breakup for each thin feature, as well as the number, size, position, and velocity of the resulting droplets.

Comparison of Labeling Algorithms
The computational cost of the underlying CCL algorithm determines its feasibility for the SGS modeling of films and ligaments. For SGS modeling to be feasible at every timestep, the computational cost of the labeling step should be much less than that of the pressure Poisson step  in the flow solver. Our implementation of the proposed method uses the union-find tree data structure [18] found in many state-of-the-art CCL algorithms and can be classified as a labelequivalence method [8]. In contrast, the method of Herrmann [12], used in [19,20,21], can be classified as a label-propagation method. Label-equivalence methods and label-propagation methods have similar cost when the domain only contains a small number of structures, but label-equivalence methods are faster when numerous long structures, such as the ligaments and films produced in primary atomization, are present in the computational domain. As there are only two structures and three films in the snapshot in Figure 3, both label-equivalence and label-propagation methods are fast enough to be used at every timestep. In other cases where there is a higher concentration of fluid features, label-propagation methods remain quick enough, relative to the cost of the pressure solver, to be used at every timestep. In practice, the computational cost of the proposed labeling method is less than 1% of the cost of the underlying flow solver.

Conclusions
Through testing of the proposed method on simulations of binary colliding droplets and turbulent air-blast atomization, we demonstrate that thin fluid features can be quickly and accurately identified through the addition of interface normal and centroid information to a connectedcomponent labeling (CCL) algorithm. In order to employ this method to model topology change in spray atomization simulations, our current focus is on the development of methods for the classification of thin features by shape, such as ligament or sheet, and the development and integration of subgrid-scale breakup models in a multiphase flow solver in a manner that allows for their localized, per-feature application.