An extension of the multiple marker algorithm for study of phase separation problems at the mesoscale

Multiphase fluid flow is an active field of research in numerous branches of science and technology. An interesting subset of multiphase flow problems involves the dispersion of one phase into another in the form of many small bubbles or droplets, and their subsequent separation back into bulk phases after this has occurred. Phase dispersion may be a desirable effect, for example in the production of emulsions of otherwise immiscible liquids or to increase interfacial surface area for chemical reactions, or an undesirable one, for example in the intermixing of waste and product phases during processing or the generation of foams preventing gas-liquid decoupling. The present paper describes a computational fluid dynamics method based on the multiple marker front-capturing algorithm – itself an extension of the volume-of-fluids method for multiphase flow – which is capable of scaling to mesoscale systems involving thousands of droplets or bubbles. The method includes sub-grid models for solution of the Reynolds equation to account for thin film dynamics and rupture. The method is demonstrated with an implementation in the OpenFOAM R © computational mechanics framework. Comparisons against empirical data are presented, together with a performance benchmarking study and example applications.


Introduction
The problem of separating two immiscible fluids after they have become mixed arises in many branches of process engineering. In particular, phenomena related to contact between two or more immiscible fluids are critical for industries such as petrochemistry, extractive metallurgy, chemical engineering, and many others.
In the oil and gas industry, mixing of immiscible liquid phases is generally undesirable and results in intermixing or contamination of product streams. The contaminant is usually water dispersed in a hydrocarbon phase, and specialised settling equipment is required to separate the two phases to an acceptable degree; this must occur by recombination of droplets of the contaminant phase with each other and the bulk phase by settling and coalescence processes. Improving the performance of such separation equipment may have significant process cost and efficiency benefits, and perhaps unsurprisingly a great deal of both fundamental and experimental research has been conducted in this area (for example [1] and contemporary papers). In the field of pyrometallurgical engineering, phase mixing and separation are similarly associated with undesirable effects -for example slag foaming, in which gas bubbles are generated by or injected into the process and build up at the upper surface of the viscous slag phase, or slag-metal entrainment, in which small droplets of one phase become separated from the phase interface due to mechanical shear or other disturbances [2][3][4]. The latter problem is of particular concern in the operation of smelting furnaces, as entrained product particles discarded in a waste phase may represent a significant economic loss or environmental hazard. Experimental studies have demonstrated that the recombination of molten metal droplets dispersed in a molten slag phase may take a considerable amount of time, even if the droplets are of the order of several millimetres in size [5].
Existing research into the computational fluid dynamics (CFD) modelling of phase separation may be broadly grouped into one of two approaches. Unresolved methods seek to represent droplets of one phase entrained in another using a continuum approximation and model only the bulk motion and large-scale flow structures of the fluids involved. Models of this type are exemplified by the Euler-Euler method for multiphase flow, in which a separate set of Navier-Stokes equations is solved numerically for each phase present and interactions between the phases (such as intermixing, drag, and droplet break-up and coalescensce) must be accounted for using empirical closures such as population balance models. Such methods have great utility for engineering calculations as they are robust and relatively low in computational complexity, but a drawback is that the closures used are generally problemspecific and must be determined for each case of interest. Examples of work using unresolved methods for the study of unit operations in extractive metallurgy can be found in [6] and [7].
By contrast, resolved methods attempt to model the individual entrained droplets and some or all of their interactions with the surrounding phase directly. Within resolved methods, Euler-Lagrange methods treat droplets as discrete particles with trajectories computed from force balances, and pure Eulerian methods model the fluid inside entrained droplets using the same governing equations as the surrounding phase while completely resolving all phase interfaces. The former includes discrete phase models and CFD-coupled discrete element methods and has been widely applied to engineering problems such as the study of bubble plumes in reaction vessels [8], while the latter includes level set and multiple-marker volume of fluids (VOF) methods and is used primarily in the detailed study of individual droplet collisions [9][10][11][12][13][14]. Euler-Lagrange methods require empirical closures between the particle and fluid models similar to those used in unresolved methods, and as a result are subject to many of the same limitations. Pure Eulerian methods are arguably the most accurate, but are prohibitively computationally expensive due to their need to resolve detail across a large range of length scales (typically from nanometres in the thin films between coalescing droplets up to centimetres or metres for the vessel containing the phases).
Although they are not suitable for industrial applications, pure Eulerian methods may be used to perform numerical experiments and determine parameters required for the empirical closures in more scalable methods. Since experimental data on materials for phase separation problems is occasionally unreliable or unavailable (for example in the molten state at high temperatures for pyrometallurgical applications), such bootstrapping of empirically-based simulations using results obtained from more fundamental numerical models may be the only option available to researchers; further development of pure Eulerian methods for phase separation problems is therefore of practical as well as academic value.
In order to address some of the shortcomings of existing pure Eulerian methods, we present an extension to the multiple-marker method which improves computational performance in problems with large numbers of dispersed phase particles present.

Algorithm description
The modelling of droplets or bubbles during close interactions is a challenge for standard VOF methods. Once two particles of a dispersed phase approach each other at distances on the order of the local mesh resolution, the particles will immediately merge together. This is physically unrealistic; in reality a thin film of fluid of the surrounding continuum fluid will form between the particles, and the evolution of this thin film delays the coalescence until it is able to rupture. In order to overcome this the multiple marker method [10,12] was developed, in which each droplet or bubble of the dispersed phase is assigned its own unique phase in the VOF algorithm. This method prevents premature coalescence of two colliding particles while maintaining an accurate description of the fluid flow physics throughout the model domain. The multiple marker method is typically combined with semi-analytical subgrid models governing the film behaviour [14], and information from the subgrid models is then used to determine the point of film rupture and particle coalescence.
The main limitation of the multiple marker algorithm is its lack of scalability. It is effective for the study of interactions and collisions between pairs or very small groups of dispersed phase particles, but once the number of particles (and hence marker phases) in the simulation increases, the underlying VOF algorithm on which it is based becomes prohibitively expensive.
In the new approach presented here, which we refer to as the dynamic multiple marker (DMM) method, this is addressed by dynamically reassigning particles of the dispersed phase to a small number of marker phases such that no adjacent particles use the same phase. This is achieved at each time step in a transient VOF calculation by identifying collections of mesh cells belonging to each discrete particle of the dispersed phase, tagging pairs of particles which are in close proximity to one other, and using this information to generate a graph representing the pairwise contacts for the set of all particles in the simulation. The nodes of this graph may then be associated with a set of marker phases (much smaller than the number of particles) using an appropriate graph colouring algorithm. Finally, all mesh cells associated with the dispersed phase are reassigned to their new marker phase before resuming the VOF calculation. This is combined with a refined implementation of a thin film subgrid model describing film drainage which enables prediction of film rupture between any pair of particles in the simulation. When rupture is detected, particle coalescence is effected by immediately assigning the pair of particles involved to the same marker phase and allowing the VOF algorithm to merge them together naturally.
An overview of the main steps in the DMM algorithm and film sub-model are shown in Figure 1. In general these steps are executed after the pressure, velocity, and phase field calculations for each time step in a transient VOF simulation are completed. The major algorithm components are described in more detail in the following sections.

Contact detection and marker phase identification
The contact detection and marker phase assignment part of the DMM algorithm (blue block in Figure 1) begins with a threshold calculation on the marker phase fraction fields. All cells in the mesh are looped over, with each being flagged as either belonging (true) or not belonging (false) to the dispersed phase according to a threshold limit specified by the user; this is typically 0.5.
A second loop is then performed in which all unique, contiguous sets of cells in the dispersed phase are found -in the present implementation this is performed by iterative sweeping of thresholded cells which are directly adjacent to each other, but any appropriate connectivity-detection algorithm could be used. During this procedure each discrete particle of the dispersed phase is assigned a unique integer identifier, stored at mesh cell level.
Once the particles are identified, contact (or near contact) between pairs of particles is determined by expanding each particle's set of cells in a "halo" around the actual particle using a sweep algorithm. The expanded particles' locations are stored at mesh cell level. A loop over all cells in the mesh is then conducted, and if a particular cell has been flagged as belonging to more than one expanded particle (i.e. the expanded particles overlap in space), contact between the two particles is recorded as the integer pair of their identifiers. The size of the contact detection halo is specified by the user, but is typically on the order of five to ten mesh cells.
Using the collection of particle contact pairs, a contact graph is then constructed. In this graph nodes represent particles of the dispersed phase, and edges represent contact between them. The nodes of the graph are then assigned integer tags in such a way that no pair of connected nodes has the same tag, while attempting to minimise the total number of tags in use -this is the classical graph colouring problem, for which a number of effective heuristic algorithms exist. In the present implementation DSatur [15] is used and provides a balance of simplicity, performance, and near-optimality.
Once the contact graph is coloured, the tag of each node in the graph is used to assign the associated particle of the dispersed phase to a parent marker phase. The marker phases are constructed at mesh cell level using the particles' (or expanded particles') cell sets. Care must be taken in this step to preserve mass continuity and phase fraction summation rules, since certain particles and regions may have been reassigned from one marker phase field to another.

Thin film sub-model
When particles of the dispersed phase collide, a thin film of the continuum phase is squeezed between them and begins to drain against viscous resistance. At a critical thickness, intermolecular interactions across the film come into play and it becomes unstable and ruptures [16], leading to rapid coalescence of the particles. This critical thickness is typically on the order of 10 −8 m, and therefore cannot be practically captured within the VOF model.
The Reynolds equation is employed as a sub-grid model governing the film thickness in the present work, and is familiar from lubrication theory [17]. It is a non-linear partial differential equation derived on the assumption that the separation of the particles is much smaller than their radii. A coordinate transform and a volume-based reformulation permit solution on the same mesh as the VOF model, via an iterative coupled approach in which pressure is fed in to the Reynolds equation from the VOF solution. The Reynolds equation is then used to march the film thickness forward in time and additionally to determine a tangential lubrication force which is fed back to the VOF simulation. Interested readers are referred to Musehane et al [14] which contains a full description of this approach.
It should be noted that in many-particle simulations, the volume-based approach for solution of the Reynolds equation has the added benefit of allowing all film equations to be solved simulataneously on the same mesh, and incurs no greater performance penalty than the solution of other CFD variables.

Implementation details
A technology demonstrator for the DMM algorithm was developed as a numerical solver in the OpenFOAM R [18] open source computational mechanics framework.
A standard multiphase VOF solver available in the OpenFOAM R framework, multiphaseInterFoam, was used as a basis for the implementation. multiphaseInterFoam uses a mixed PISO-SIMPLE segregated algorithm for the solution of pressure and velocity fields in the fluids combined with stabilised interface compression schemes for the phase fraction fields. The interfacial tension force calculation was extended by introducing smoothing to assist with control of parasitic currents [19], and modified to correctly model a dispersed phase in a continuum.
Droplet identification, contact detection, and marker phase assignment were implemented as separate functions in a new multiMarkerMixture code module, extended from the original multiphase solver module. Temporary fields for particle identifiers and phase assignment tags were implemented as boolean or integer lists rather than mesh field objects to reduce storage requirements and improve performance. The particle contact graph was generated in a compressed node-based format from the list of detected contact pairs. The graph colouring algorithm was implemented as a separate solver module, graphColouring, to facilitate later modularisation and automatic switching. The new functions in multiMarkerMixture were explicitly parallelised using OpenFOAM R 's Pstream class objects.
To process the sub-grid thin film model, the Reynolds equation is discretised using standard OpenFOAM finite-volume operators and first order backward-difference time integration. Ancillary extrapolation equations for the film thickness and velocity profile are solved for each pair of indicator phases, and the results assembled onto a single domain on which the Reynolds equation is solved in different areas according to an activation function.
Some care is required to deal with numerical noise in the pressure gradient coming from the VOF solution, caused by both the density gradient between the continuum and dispersed phases, and the pressure jump induced by interfacial tension. While both of these contributions to the tangential pressure gradient used in the Reynolds equation are formally zero, numerically their contribution is not perfectly cancelled. To ameliorate this, the surface tension force term in the momentum equation is subtracted from the pressure gradient, and the result is smoothed in a similar manner to the interfacial tension. This also helps to reduce the noisy effects of parasitic velocities at the surface of the particles.

Model testing
In this section, the results of some preliminary testing of the DMM method's performance and accuracy are presented.

Performance
In order to evaluate the performance of the DMM method against traditional multiphase approaches for resolved dispersed-phase problems, simulations were performed using a planar two-dimensional metallurgical system in which droplets of dense molten metal settle through a lighter molten slag continuum phase. The parameters of the test case are given in Table 1. A range of initial droplet counts between 1 and 30 were considered -in each case a simulation was executed with both the standard multiple-marker method, in which a unique marker phase must be assigned to each droplet a priori, and the new dynamic method. The standard multiple-marker cases were run using OpenFOAM R 's reference multiphaseInterFoam solver. In the DMM solver, the thin film sub-model was active but coalescence was disabled to maintain a constant droplet count. Tests were performed on an Ubuntu Linux 18.04 LTS platform equipped with an AMD Ryzen 7 1700X eight-core processor and 16GB RAM. All cases were parallelised on four CPU cores. The performance of the two methods in terms of execution time required per unit of simulated time (and the resulting performance improvement factor) is shown in Figure 2. It can be seen that the DMM method confers significant advantages when more than three or four particles of the dispersed phase are present. Moreover the scaling of the DMM method at larger droplet counts is essentially flat, implying that much larger numbers of droplets may be simulated while retaining practical run times.

Droplet collisions
To test the thin film sub-model and dynamic coalescence between particles of the marker phase, collisions between droplets of tetradecane in a nitrogen atmosphere were simulated using the DMM solver. This problem has been used as a validation test case in a number of previous publications. Here, we use the setup and initial conditions of Musehane et al [14] in which the droplets are modelled on an axisymmetric two-dimensional mesh. Parameters for the model and cases are given in Table 2. Domain mesh resolution in all cases was 360 x 80 cells, resulting in a mesh element size of approximately 3 × 10 -6 m.
In experimental tests [9] coalescence after impact between droplets was observed in cases 1 and 4, but not in cases 2 and 3. Following the approach of Musehane et al, a critical film thickness of 6.8 × 10 -7 m was estimated in order to reproduce this behaviour in the DMM  solver. The evolution of the film thickness between the droplets in each simulation as captured by the film sub-model is shown in Figure 3. Selected visualisations of the droplet impacts from case 3 and 4 are shown in Figures  4 and 5 respectively. Film thickness and droplet shape behaviour are in reasonable qualitative agreement with the data and simulation results reported by [9] and [14], however more rigorous validation is advised.

Hindered settling
In the final test case the veracity of the DMM method for predicting particle motion in groups of droplets is examined. A planar two-dimensional region containing silicone oil is simulated, and initialised with a number of static water droplets of a fixed diameter. As the simulation proceeds, the droplets settle under gravity toward the bottom of the region. Interactions between nearby droplets and counterflow of the continuum phase progressively reduces the terminal settling velocity as the number of droplets in the region increases -this is referred to as hindered settling. The Richardson-Zaki equation (1) is commonly used to characterise this problem, where u t is the terminal velocity of a single free droplet in the continuum fluid, u h is the hindered settling velocity, φ is the volume fraction of droplets present, and n is an empirical exponent which is related to the Reynolds number of the flow. n typically takes on values between 2.4 and 4.65 [20]. The parameters for the hindered settling test case are given in Table 3. Droplet volume fractions between 0.075 and 0.35 were simulated using the DMM solver. Periodic boundary conditions were applied on the left and right walls of the region. The thin film sub-model was deactivated since droplet motion, not coalescence, is of interest here. Runs were performed on the South African Centre for High Performance Computing's Lengau HPC cluster. All cases were parallelised on three nodes (72 CPU cores) and took between four and six hours to complete.
An example of the results obtained from the hindered settling simulations is shown in Figure 6, in which the evolution of the velocity field and droplet positions in time can be seen. The structure of the settling band becomes complex over time, with evidence of droplet clumping and bed channelling.
Upon completion each case was post-processed to obtain an average settling velocity of the droplets at each time step. The results were filtered by time (>0.5 s, to allow initial conditions to decay) and vertical position of the droplets (>0.1 m, to prevent settled droplets from skewing the average). The overall average settling velocity, as well as the minimum and    Figure 7 and show reasonable agreement with (1).

Application example
As a demonstration of the full capabilities of the DMM method, an example application in metallurgical phase separation was selected in which a number of small droplets of molten metal must settle through a viscous molten slag layer under gravity, and coalesce into a bulk metal phase. This is a prototypical problem in pyrometallurgical engineering, where undesirable emulsions between liquid phases can easily form during handling and processing of furnace products in the molten state. Of particular interest is the effect of thermophysical properties of the materials on the phase separation behaviour; since these are strongly dependent on chemical composition of the phases as well as their temperature, they are able to vary across wide ranges in pyrometallurgical processes. A further complication is that the properties of the phases affect both the settling rate of the droplets [3] as well as film rupture and droplet coalescence [21]. Meso-scale methods which are able to account for both effects in a single framework are therefore particularly useful in this case. The parameters chosen for the application study are shown in Table 4, and represent a generic slag-metal separation problem. The film sub-model was active, with a critical film thickness of 1 × 10 -5 m.
Periodic boundary conditions were imposed on the left and right sides of the domain, with no-slip walls applied at the upper and lower surfaces. The model region was initialised with 1000 randomly-sized and positioned droplets, and a small layer of bulk metal phase was included at the bottom of the domain. As before calculation runs were performed on the South African Centre for High Performance Computing's Lengau HPC cluster, using three nodes (72 CPU cores) for each case and taking between 18 and 24 hours to complete. Selected visualisations of the velocity field and phase boundaries from the base case simulation are shown in Figure 8. As time progresses it can be seen that droplets initially separate into a dispersed settling band and a densely-packed layer which accumulates above the metal pool. Within the dense layer adjacent droplets begin to coalesce as they move downward, increasing the average droplet size. The largest droplets at the bottom of the dense layer come into contact with the pool and coalesce with it to join the bulk metal phase. The bulk phase  The effect of changing the physical properties of the system was then examined by adjusting their values relative to those given in Table 4 and re-running the model. The total number of droplets remaining in the simulation was recorded at every time step as an indication of how the phase separation is impacted. The results are shown in Figure 9 for slag viscosity and interfacial tension between the phases. It can be seen that in this case the slag viscosity has a moderate effect and influence both settling and film drainage, with higher values causing slower separation. The interfacial tension (which affects only the film drainage) has a marked effect, with lower values significantly increasing the phase separation time; this is to be expected in metallurgical systems which typically have very high interfacial tensions.

Conclusions
The dynamic multiple marker method for resolved modelling of multiphase fluid flow problems with dispersed phases was successfully developed and demonstrated. The method extends existing multiple marker algorithms to dynamic situations involving dispersed-phase particle coalescence and breakup, and improves performance to the point where it is practical to consider systems up to thousands of particles. Sub-grid modelling of the behaviour of thin films between dispersed phase particles was implemented via numerical solution of the Reynolds equation, enabling physically realistic control of particle coalescence. Test cases in hindered settling and droplet collisions were seen to be in good agreement with previous work. It is anticipated that such methods will find some utility in exploring the behaviour of liquid-liquid emulsion and gas-liquid mixing problems which arise in numerous branches of process engineering. They are especially applicable to problems for which experimental data is difficult or impossible to obtain, and may provide a supplemental or alternative route to developing empirical closures for unresolved models. The exploration of a pyrometallurgical slag-metal system demonstrated the new method's capabilities in studying scenarios involving complex phase separation and decoupling.
Naturally much work still remains to be done. The current implementation of the DMM method is limited to two-phase systems with one dispersed phase and one continuum phase, and although the method can in principle be extended to an arbitrary number of each, this will add logical complexity and impact performance. A careful analysis and optimisation of the contact detection and graph colouring algorithms used for marker phase selection is advisable; for heavily parallelised systems with many thousands of dispersed-phase particles, these components start to become significant performance barriers. The accuracy of the thin film sub-model as well as VOF methods in general when simulating materials with extreme physical properties (very high viscosities and surface tensions) also remains an open question, and it is likely that further algorithm development and refinement will be required here together with reliable validation data for problems and regimes of interest. Finally, extension of the thin film sub-models to include additional physics such as local electrostatic effects, non-Newtonian behaviour, and other phenomena is possible in principle and is currently being explored.