A unified finite volume framework for phase-field simulations of an arbitrary number of fluid phases

While the phase-field methodology is widely adopted for simulating two-phase flows, the simulation of an arbitrary number ( N ≥ 2) of fluid phases at physical fidelity is non-trivial and requires special attention concerning mathematical modelling, numerical discretization, and solution algorithm. We present our most recent work with a focus on validation for multiple immiscible, incompressible, and isothermal phases, enhancing further our library for diffuse interface phase-field interface capturing methods in OpenFOAM (FOAM-extend 4.0/4.1). The phase-field method is an energetic variational formulation based on the work of Cahn and Hilliard where the interface is composed of a physical diffuse layer resembling realistic interfaces. The evolution of the phases is then governed by the minimization of the free energy of the system. The accuracy of the method is demonstrated for a number of test problems, including a floating liquid lens, bubble rise in two stratified layers, and drop impact onto thin


| INTRODUCTION
Numerical simulation of technically relevant applications of fluid dynamics often extends beyond a twophase problem formulation. In numerous processes, more than two phases co-exist with multiple fluid interfaces dynamically evolving. Of particular interest to this study are applications linked to drop-filminteractions in internal combustion engines. In such combustion processes, fuel droplets collide with a thin layer of lubricating oil film on the piston that might weaken/change the oil's lubricating properties and cause combustion of a fuel-oil mixture instead of a pure fuel combustion. A fundamental understanding of such a complex and dynamic process requires detailed numerical simulation of the physics involved during the impact process that includes more than two phases. To this end, we have extended our two-phase diffuse interface modelling library for phase-field interface capturing simulations to cover more than two immiscible, incompressible, and isothermal phases. Here, we refer to a region of matter surrounded by an adjacent medium as a phase, which in the present work can be liquid or gaseous fluid regions, respectively.
The building block of our approach is the phase-field interface capturing method, based on a diffuse-interface model, where the system is described using one field variable in the case of a two-phase flow or multiple field variables in the case of a multiphase flow problem. Unlike interface capturing methods based on sharp-interface modelling, such as volume-of-fluid (VOF) or level-set (LS), the composition profile across the interface in the phase-field method is sigmoid in shape, based on the work of Cahn and Hilliard, [1] that reflects a physically diffuse profile. The ingenuity of this method stems from the fact that the free energy of a non-homogeneous system contains interfacial energy-related effects in the form of composition gradient, and the system evolves by minimization of its free energy. This is different from the aforementioned methods, which underlie a continuum mechanical rationale; see the work of Brackbill et al., [2] for instance. The advantages and challenges of the phasefield methodology are set out and discussed extensively in the work of Feng et al. [3] The fluid dynamics of a system of multiple fluid phases is governed by the coupled Cahn-Hilliard Navier-Stokes equations. For the mathematical foundation for immiscible, incompressible, and isothermal two-phase flows, one can refer to many available treatises in the literature, such as the works of Lowengrub and Truskinovsky, [4] Jacqmin, [5] Jacqmin, [6] Liu and Shen, [7] Yue et al., [8] and Ding et al. [9] Extension from two-phase to ternary-phase systems were proposed for instance in the works of Kim and Lowengrub, [10] Boyer and Lapuerta, [11] Kim, [12] Boyer et al., [13] and He et al. [14] Such models for ternary systems are available in proprietary software, like, for instance, in the microfluidics module of COM-SOL Multiphysics software that is described in the COM-SOL Multiphysics User's guide [15] based on the work of Boyer et al. [13] Open-source code for the simulation of ternary systems has been made available under the numerical simulation platform PELICANS. [16] The model in Boyer et al. [13] was subjected to a study pursued in the work of Řehoř et al. [17] and implemented using the FEniCS Project. [18] Extending the model formulation to systems containing an arbitrary number of fluid phases (i.e., the general case of N > 2 phases) is non-trivial and indeed subject to ongoing research building on continuous efforts during the last decade. [19][20][21][22][23][24][25][26][27][28][29] One central challenge for an Nphase model to be of physical fidelity resides in the requirement to honour thermodynamic and reduction consistency likewise. Based on the pioneering work of Boyer and Minjeaud, [20] Dong [24] presented for the first time a fully reduction-consistent model formulation, that is, an N-phase model, which reduces to the established two-phase model, when only two phases are initialized to co-exist in multiphase setups. Their N-phase model was moreover reported to be thermodynamically consistent, that is, to obey conservation of mass, conservation of momentum, the second law of thermodynamics, and Galilean invariance. [24] Nevertheless, the model of Huang et al. [27] then contains several modifications to fulfil consistency conditions and further improve the work of Dong. [24] Differences are indeed subtle, and these papers have paved the avenue towards the practical deployment of the phase-field method for multiphase flow simulations. On these bases, we have implemented a versatile diffuse-interface model library for multiple phases using the finite volume method (FVM) with support for unstructured meshes of general topology in OpenFOAM (FOAM-extend 4.0/4.1). In future work, we will aim at the deployment for simulating multiphase flows in geometrically complex domains such as porous media or fibre mats. This work shall demonstrate our validation efforts for a broad bandwidth of different literatureknown numerical tests and benchmarks for multiphase flows including a floating liquid lens, bubble rise in two stratified layers, and drop impact onto thin liquid film.

| Two-phase Cahn-Hilliard Navier-Stokes equations
The flow of two immiscible, incompressible, and isothermal Newtonian fluid phases is governed by the coupled Cahn-Hilliard Navier-Stokes equations. Here, only a single phase-field variable c [À1, 1] (volume fraction contrast) has to be used to discriminate the presence of the phases. The concentration in the bulk phases is represented by c = ±1 and the position of the fluid interface can be represented as the zero iso-surface c = 0. Following the work of Cahn and Hilliard, [1] the total free energy (Helmholtz free energy functional) of this system consists of the sum of the bulk and the gradient energy contributions, which can be written as [8] : where Ω is the fluid domain, λ is the interfacial mixing energy density parameter, ε is the capillary width, and Ψ(c) is the so-called double-well potential which can be phenomenologically modelled according to Ginzburg and Landau, [30] that is, Ψ c ð Þ ¼ 1 4 c 2 À 1 ð Þ 2 . Assuming a planar interface and equilibrium conditions, the interfacial mixing energy density parameter λ can be shown to be λ ¼ 3 2 ffiffi 2 p σε where σ is the interfacial tension coefficient. The evolution/kinetics of the system is described by the minimization of free energy through a variational procedure, [1] which leads to the chemical potential Φ of the system: where Ψ 0 (c) is the first derivative of the double-well potential with respect to the order parameter c. Using the Ginzburg and Landau double-well potential, Ψ 0 (c) = c 3 À c. From this, the equilibrium interface profile can be obtained when requiring Φ = 0, which yields: where x is the local coordinate normal to the interface. The fluid dynamics of two immiscible incompressible phases are then governed by the Cahn-Hilliard equation being coupled with the Navier-Stokes equation following the work of Jacqmin, [5] viz.: where t is the time, u is the divergence-free velocity field, and ρ c and μ c are the volumetric average density and dynamic viscosity of the two fluids, given as The term e p denotes a modified pressure, since parts of the known Korteweg tensor term accounting for interfacial capillarity have been absorbed into the pressure gradient term. For Newtonian fluids, the viscous stress tensor is τ = μ c (ru + (ru) T ). The gravitational acceleration is g and J is the phase-field flux that generalizes the Fick's law as J = ÀM rΦ, with M being the mobility parameter. The term f s = Φ rc models the capillarity of the diffusive interfaces on the basis of the Korteweg stress tensor. [31] Note that (Equation (4a)) is a fourth-order non-linear parabolic partial differential equation with respect to c, which renders it challenging to solve numerically. [32] 2.2 | N-phase Cahn-Hilliard Navier-Stokes equations The evolution of the order parameters c p , 1 ⩽ p ⩽ N, for N immiscible fluids is governed by the Cahn-Hilliard transport equations, as [20,33] : Due to isochoric conditions (phase volume conservation), it applies for multiple phases that P N p¼1 The phase-field flux in Equation (5) is: where is the chemical potential of phase p, and is the scalar mobility between phases p and q, where M 0 is a non-negative mobility constant. In recognition of the intuitive statement that the mobility should be tensorial, [34] we devise Q p,q in Equation (6) as a modified projection operator I À n p n p , viz. [34] : Note that in the sharp-interface-limit ε ! 0 ð Þ the operator r Á M p,q r in Equation (5) reduces to the Laplace-Beltrami operator Δ s in the asymptotic limit. Details about the exponent in Equation (6), which we choose to m = 4, are given in the work of Gugenberger et al. [34] Using the potential according to Ginzburg and Landau, [30] the phenomenological potential functions are The interfacial mixing energy density parameters λ p,q in Equation (7) are λ p,q ¼ 3 2 ffiffi 2 p σ p,q ε, as a function of pairwise interfacial tension coefficients and capillary widths.
Here, the matrix λ p,q È É N p,q¼1 is symmetric and zero-diagonal. It should be noted that zero gradient (homogeneous Neumann) boundary conditions are used for both the order parameters and the chemical potentials on the computational boundaries in our work.
The coupling with the Navier-Stokes equations, is accomplished via constitutive models, which hold the order parameter and/or the chemical potential. The average density and dynamic viscosity are 2 , respectively. The diffusive mass flux is where which is the total free energy density for N immiscible, incompressible, and isothermal phases. The chemical potentials Φ p are determined through the variational derivative of the total free energy with respect to the order parameter. This results in the well-known Korteweg stress term, but in the formulation for a multiphase system. This term is transformed and partially absorbed into the pressure gradient term: 1

| SOFTWARE DESCRIPTION, ALGORITHM, AND NUMERICAL SCHEMES
The diffuse-interface phase-field method for an arbitrary number of fluid phases, as described in Section 2.2 is implemented in OpenFOAM (FOAM-extend 4.1) as model library used by a top-level solver referred to as phaseFieldFoam. OpenFOAM is a comprehensive open-source C++ library for computational continuum physics including computational fluid dynamics (CFD). [35][36][37] The code structure follows a rigorous objectoriented and physics-guided development paradigm, which results in a versatile code design that is extendable towards coupled (multi-)physics applications. The deprecated two-phase version of phaseField-Foam has been extensively validated for different physics in the works of Cai et al., [38] Cai et al., [39] Fink et al., [40] Wörner et al., [41] Samkhaniani et al., [42] and Bagheri et al. [43] and its solution algorithm is described in the work of Jamshidi et al. [44] The algorithms for pressurevelocity coupling lean on other interface capturing approaches implemented in OpenFOAM and described in great detail in the work of Deshpande et al. [45] and have been massively extended in our work to cope with high-fidelity solution of the coupled Cahn-Hilliard-Navier-Stokes equations for N immiscible, incompressible, and isothermal phases. The viscous stress tensor τ implementation is generic in the sense that it allows the use of any fluid rheology model available in OpenFOAM. The interfacial Korteweg tensor term f s is implemented as a volumetric density, whereupon in the context of phase-field method it substantially reduces the parasitic/ spurious currents by orders of magnitude when compared to methods based on sharp-interface models. [5,44,46] The basic structure of our diffuse-interface model library is depicted in Figure 1. It consists of two main base classes, namely, multiphaseSystem that includes constitutive models and phaseFieldEquations that includes the transport equations for the order parameters.
The pressure-velocity coupling is carried out using the PIMPLE algorithm, which is a combination of the Pressure Implicit with Splitting of Operator (PISO) [47] and the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [48] or the SIMPLE-Consistent (SIMPLEC) [49] algorithms, respectively, which can be selected.
The N-phase system of equations described in Section 2.2 are discretized using the FVM with support for unstructured meshes of general topology. Moreover, the code is equipped with dynamic load balancing (DLB) along with adaptive mesh refinement (AMR) techniques. Deploying this solution-adaptive approach allows to increase the mesh resolution dynamically where needed (at the interface) while keeping the computational load approximately balanced on parallel computations on distributed-memory computer architectures where domain decomposition is applied.
The system of equations described in Section 2.2 are applied to a series of validation cases in Section 4 to illustrate the capabilities of the software in the simulation of N > 2 immiscible phases. The discretization of the equations is performed using second order accurate schemes for both temporal and spatial terms. The advection terms (divSchemes), that is, the divergence terms for the scalar transport of the order parameters, are discretized using the high-resolution scheme SuperBee, [50] while the convection term within the momentum equation is discretized using the second-order accurate Gam-maV scheme. The temporal discretization scheme (ddtSchemes) is set to a second-order accurate backward scheme. The gradient terms (gradSchemes) are discretized using the second order accurate Gauss linear scheme and the diffusion terms (lapla-cianSchemes) are discretized using the second order accurate Gauss linear orthogonal scheme.

| Floating liquid lens
The floating liquid lens problem, which is considered a three-phase benchmark problem, is simulated here. This two-dimensional problem has been also studied in the works of Dong, Huang et al., and Yuan et al. [21,22,26,27,51] In this problem, an initially circular liquid oil drop is floating on the air-water interface (see Figure 2) and will spread depending on the magnitude of the gravitational acceleration and interfacial tension coefficients. The simulation is performed until the oil drop reaches an equilibrium configuration. The final thickness of the drop will then be compared to an exact equilibrium solution given in the work of Huang et al. [27] for the case when jgj = 0 m/s 2 and an asymptotic analysis provided by Langmuir and De Gennes et al., [52,53] viz.: for all other values of the gravitational acceleration. The transport properties of the fluids are shown in Table 1.
The computational domain is [ÀL, L] Â [0, 0.8 L] with the length scale L = 0.04 m. The wall (no-slip) boundary condition is applied to the upper and lower boundaries, the cyclic (periodic) boundary condition is applied at the left and right boundaries, and the empty boundary condition is applied to the front and back boundaries, given that the computational domain is 2D.
The capillary width is ε = 0.01 L where the Cahn number is selected to be equal to 0.01 (Cn Similar to the work of Huang et al. [27] we have also considered three different magnitudes for the gravitational acceleration, which are (0, 5, 9.8) m/s 2 in the negative Z-direction. Figure 3 depicts the final configuration of the threephase system after reaching an equilibrium state corresponding to different magnitudes of the gravitational acceleration. The final shape of the air-oil-water system depends strongly on the interaction of the magnitude of the gravitational acceleration, pair-wise interfacial tension coefficient, and also the influence of fluid densities. [22] As the value of the gravitational acceleration increases further, the final thickness of the drop at the centre point in the longitudinal direction reduces. In the event where the interfacial tension dominates, that is, jgj = 0 m/s 2 , the oil drop resembles two circular 'caps' in the upper and lower regions, which is in agreement with the work of De Gennes et al. [53] and in the event where the gravity dominates, the oil drop resembles a 'puddle'. Figure 4 shows the normalized drop/puddle thickness versus the normalized gravity. The drop/puddle thickness is normalized with respect to the length scale L and the gravity is normalized with respect to jgj = 1 m/s 2 . The final maximum thickness is measured between the upper and lower interfaces of the equilibrium configuration in the vertical direction. It is observed in Figure 4 that our numerical results are in excellent quantitative agreement with the exact and asymptotic solutions.

| Rising bubble in two stratified layers
Here we consider the three-phase rising bubble in two stratified layers problem studied in the works of Boyer and Lapuerta, Boyer et al., and Fontes [11,13,54] where a single air bubble of a specified volume rises initially in a heavy fluid and attempts to cross the interface between a heavy and a light fluid. Depending on the bubble volume,  [52,53] Huang et al. [27] F I G U R E 4 Normalized thickness of the floating liquid lens versus normalized gravity it can either cross the fluid-fluid interface or remain captured at the interface. A criterion is proposed in the works of Greene et al. [55,56] that predicts the state of the bubble at the interface depending on a critical volume, viz.
For V bubble > V* the air bubble will penetrate across the interface, while for V bubble < V* it will be trapped at the interface. The criterion is based on 'macroscopic balance between buoyancy and interfacial tension forces' [13] and has been validated experimentally. It should be noted that the criterion only considers Archimede's and interfacial tension forces and ignores hydrodynamic effects. [11] A schematic diagram illustrating the geometrical details of the rising bubble case setup is shown in Figure 5. The computational domain considered in this study is axisymmetric and has a dimension of [0, D] Â [0, 10D], where D is the bubble initial diameter and the gravity is applied in the negative Z-direction. The initial barycenter position of the air bubble is at [0, 1.25D]. The wall (no-slip) boundary condition is applied at the lower, upper, and the right boundaries, and the wedge boundary condition is applied at the left boundary.
The capillary width is ε = 0.025D and the computational domain is discretized by [80 Â 800] cells, resulting in N c ≈ 8 interfacial cells. The mobility parameter is set to be M = χ Â ε 2 with the constant χ = 0.01 m Á s/kg. An adaptive time stepping criteria is utilized with the maximum Courant number Co max = 0.25 and initial time step of Δt = 0.001 s. The computational domain is initialized with zero initial velocity at t = 0.
The transport properties of the fluids are shown in Table 2.
Based on the transport properties of the fluids given in Table 2, Relation (14) gives the critical bubble radius r* ≈ 2.7664 mm. In other words, an air bubble of radius r bubble > r* can penetrate across the interface and an air bubble of radius r bubble < r* remains trapped in the interface.
Based on the calculated r*, we have performed numerical simulation for four different bubble radii, which are r bubble = 2.5, 3.0, 5.0, and 7.0 mm. The numerical results are depicted in Figures 6-9, respectively. It should be pointed out that since the computational domain is axisymmetric, a rotational extrusion is performed for illustration of the bubble during the postprocessing phase.
Both bubble penetration and entrapment are predicted successfully in our numerical simulations. It can be seen from Figure 6 that since r bubble < r*, the bubble is entrapped at the liquid-liquid interface. Figures 7-9 show bubble penetration through the liquid-liquid interface since r bubble > r*. It can also be seen that the penetrating bubble carries part of the heavy fluid into the light fluid domain in its wake (trailing edge), which is also seen in the works of Boyer and Lapuerta, Boyer et al., and Fontes [11,13,54] and as the radius of the air bubble increases further, the length of the 'heavy fluid tail' becomes longer.

| A four phase problem
To illustrate the capabilities of the N-phase solver, the four phase fluid mixture problem proposed in the work of Dong [22] has been simulated and the results are reported in this section. Here, the dynamics of four immiscible incompressible phases with transport properties as given in Table 3 is considered.
A schematic of the computational domain is illustrated in Figure 10. The computational domain is [ÀL/2,  with an initial time step of Δt = 0.25 μs. The system is initialized with no initial velocity at t = 0 s. At t = 0 s, the upper half of the domain is filled with air, which contains the F2 drop, and the lower half of the domain is filled with water, which contains an air bubble and the F1 drop. When the system is released, the F2 drop free falls into the water through air and the air bubble and F1 drop move upward through the water. Figure 11 depicts sequential time instances of the dynamics of the four phase mixture problem. In the upper half of the computational domain, gravity causes the F2 drop to descend rapidly through air, while maintaining its original circular shape, and impact the water surface. In the lower half of the computational domain, buoyancy causes the air bubble and the F1 drop to slowly ascend through water while deformation of the original circular shape of both phases is observed. At the moment of the impact of the F2 drop onto the water surface, an air pocket (trapped air) is formed beneath the F2 fluid. The air pocket remains trapped during the entire span of the simulation. The impact also generates outward-moving waves or ripples on the water surface. As a result of the small density contrast, fluid F2 remains floating but mostly submerged in water. The air bubble rises through water trying to join the air side of the computational domain. As the air bubble approaches the surface, it touches the F2 fluid, disconnects a portion of it from the water, escapes the water interface, and joins the air side of the domain. At the same time, the F1 fluid rises to the water surface and stays floating. By the end of the computation, the surface of the water is mainly covered by floating fluids of F1 and F2, which includes an air pocket. It is important to note that the characteristic features described in this section have been observed in the work of Dong. [22] But a close inspection also reveals minor differences, such as the dynamics of the air bubble rupture at the free-surface at t = 125 ms ( Figure 11E), which might be due to different numerical schemes used and subtle differences in the underlying models related to reduction-consistency.

| Drop impact onto thin liquid film
The drop impact process of identical liquids is widely studied both experimentally [43,[57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] and numerically. [43,[72][73][74][75][76][77][78][79][80][81][82][83][84][85][86][87] On the contrary, experiments conducted to describe the drop impact process onto thin liquid film of non-identical immiscible liquids are limited. [88][89][90][91][92] By non-identical immiscible liquids, we explicitly mean that the film liquid is different from that of the drop and a homogeneous mixture is not formed when the liquids are mixed. The majority of the experimental studies focus on the impact morphology of miscible systems [92][93][94][95][96][97][98] and typically the impact parameters are selected such that they lead to rim instabilities and, consequently, the generation of secondary droplets. To the best of our knowledge, even fewer numerical investigations of this complex process have been performed so far, such as the works of Yeganehdoust et al. and Wang et al. [99,100] This can be due to complex impact dynamics of interfaces and the presence of triple-junction points. [99,101] The work of Yeganehdoust et al. [99] is concerned with a numerical study of the entrapped air layer when a water droplet of D = 2 mm impinges (Dupont Krytox103) film layer at three different Weber numbers and film thickness parameters. Wang et al. [100] consider numerical computations of impact dynamics of a microsized water droplet falling onto an oil layer. Both studies utilize the VOF method for capturing the dynamics of interfaces. The above discussion clearly shows the lack of numerical analysis covering drop impact process onto thin liquid films of immiscible liquids. Therefore, in this section we perform numerical computations of a silicone oil drop impacting onto a thin liquid film of water dyed with fuchsine (maroon coloured). Fuchsine does not bias the transport properties of water. The numerical results are compared to the in-house experimental results qualitatively. This study is a direct continuation of our previous work on drop-film interaction. [43] The experimental setup is shown in Figure 12. The setup consists of a drop generator system, an impact substrate wetted by a liquid film, as well as an optical system for recording the drop impact. The drop generated by the drop generator system is accelerated by gravity and impinges onto the liquid film. The substrate is a sapphire plate. A foil made of polyvinyl chloride with a recess of 50 mm in diameter and 0.6 mm in height is applied to the plate to contain the liquid film. The film thickness is precisely controlled by a confocal-chromatic film thickness sensor (Micro Epsilon confocalDT 2421 IFS 2405-1).
F I G U R E 1 0 Schematic of the geometrical details of the four phase mixture case The optical system consists of a Photron SA-X2 chromatic high-speed camera and two high-performance LEDs (Constellation 120E). The LEDs in combination with a white screen in the background provide uniform illumination. The drop impact is captured with a frame rate of 20 000 fps and a resolution of 16.5 μm per pixel. The surface tension between maroon coloured water and air σ 1,2 = 0.0728 ± 0.003 kg/s 2 as well as the interfacial F I G U R E 1 1 Contours of order parameter-four phase mixture problem tension between maroon dyed water and silicone oil σ 2,3 = 0.027 ± 0.003 kg/s 2 are measured with the DCAT 25 tensiometer from Dataphysics.
The transport properties of the liquids are shown in Table 4. In our study, the film height is h = 0.  changes of the process are captured. The impacting drop is initialized very close to the liquid film. Justifications regarding the selected domain size and initial position of the impacting drop can be found in the work of Bagheri et al. [43] The capillary width is ε = 0.01D and the mobility parameter is set to be M = χ Â ε 2 with the constant χ = 0.001 m Á s/kg. The computational grid is generated with refinement regions that resolve the interfaces with approximately 16 interfacial cells during the entire time span of the simulation. This is vital in phase-field simulation of the drop impact onto thin liquid films since the interface profile affects calculation of the gradients of the order parameters and, thereupon, computation of the interfacial energies. The boundary conditions applied and the time stepping criteria employed are similar to those of our previous work. [43] Figure 14 depicts sequential time instances of the impact process where the experiment (left) is compared side by side to the simulation (right). Rotational extrusion is used for the post-processing of the numerical results and better visualization of the impact process. Moreover, the liquid film (water-fuchsine) is coloured maroon to represent the experiment realistically and to be distinguished from the grey impacting oil drop. We also reduced the opacity of the liquid film to 45% so that the topological changes that the oil drop undergoes can be visualized clearly in Figure 14.
The topological changes of the drop impact process for identical liquids include the coalescence of the impacting drop and the liquid film in the early stages of the impact phenomenon. This is completely different for the case of non-identical liquids. Since the liquids are immiscible, the interfaces experience compression at the bottom half of the impacting drop and expansion at the upper half of the impacting drop. This clearly causes alteration of the initially specified capillary width ε and directly changes the interfacial energy associated with the interfaces. To address the accurate calculation of the interfacial mixing energy density parameter, we utilized our relaxation methodology developed in previous work [43] to estimate an out-of-equilibrium mixing energy density parameter. Through numerical investigation, we found that the equilibrium formulation of λ p,q results in the incorrect prediction of the dynamics of the impact process.
As depicted in Figure 14, the crown rises as a result of a kinematic discontinuity caused by the jump in both the liquid film thickness and the local velocity field. [102,103] The crown wall expansion is dominated by both the outer lamella sheet, which consists of water, and the inner lamella sheet, which consists of silicone oil. It can be seen from the simulation results that the crater widening is mainly due to the existence of the inner oil lamella sheet. The free rim, which restricts the crown wall expansion from above, consists of two layers of both liquids. In other words, there exist two free rims in this process. The crater descends as the interfacial tension forces become dominant compared to the inertial force. At the end of the impact process a dome-like structure consisting of the impacting drop liquid is formed. This structure can change to a central (Worthington) jet at higher Weber numbers. The final configuration of this impact process is an oil disk floating on the water surface. It can be seen from Figure 14 that the entire process is captured very well by the numerical method described in this paper.

| SUMMARY AND OUTLOOK
In this contribution, we have performed exhaustive validation of our new phase-field solver for an arbitrary number (N ≥ 2) of fluid phases against a broad bandwidth of different literature-known numerical tests and benchmarks for multiphase flows including a floating liquid lens and bubble rise in two stratified layers. The accuracy of the implemented diffuse-interface model has been shown for N = 3 phases against the aforementioned problems. The numerical method also showed a very good qualitative agreement against in-house experimental study of drop impact onto thin liquid film for nonidentical immiscible liquids. We have demonstrated the potential of the method for simulation of N = 4 phases. The main objective has been to substantiate the physical fidelity of our new diffuse-interface model library and solver, demonstrating suitability for flows of multiple fluid phases, using the FVM with support for unstructured meshes of general topology in OpenFOAM (FOAM-extend 4.0/4.1).
In future work, we will build on this development aiming at simulations of multiphase flows in geometrically complex domains such as porous media or fibre mats. We will also extend our methodology to simulate miscible systems and enhance the computational efficiency by proper utilization of load-balanced AMR techniques.