An Edge-based Interface Tracking (EBIT) Method for Multiphase Flows with Phase Change

We present a novel Front-Tracking method, the Edge-Based Interface Tracking (EBIT) method for multiphase flow simulations. In the EBIT method, the markers are located on the grid edges and the interface can be reconstructed without storing the connectivity of the markers. This feature makes the process of marker addition or removal easier than in the traditional Front-Tracking method. The EBIT method also allows almost automatic parallelization due to the lack of explicit connectivity. In a previous journal article we have presented the kinematic part of the EBIT method, that includes the algorithms for piecewise linear reconstruction and advection of the interface. Here, we complete the presentation of the EBIT method and combine the kinematic algorithm with a Navier--Stokes solver. A circle fit is now implemented to improve the accuracy of mass conservation in the reconstruction phase. Furthermore, to identify the reference phase and to distinguish ambiguous topological configurations, we introduce a new feature: the Color Vertex. For the coupling with the Navier--Stokes equations, we first calculate volume fractions from the position of the markers and the Color Vertex, then viscosity and density fields from the computed volume fractions and finally surface tension stresses with the Height-Function method. In addition, an automatic topology change algorithm is implemented into the EBIT method, making it possible the simulation of more complex flows. The two-dimensional version of the EBIT method has been implemented in the free Basilisk platform, and validated with seven standard test cases: stagnation flow, translation with uniform velocity, single vortex, Zalesak's disk, capillary wave, Rayleigh-Taylor instability and rising bubble. The results are compared with those obtained with the Volume-of-Fluid (VOF) method already implemented in Basilisk.


Introduction
Multiphase flows with phase change play a critical role in various industrial applications, including combustion engines [1], electronic cooling systems [2], and nuclear reactors [3].A deeper understanding of the associated mass and heat transfer processes is essential for advancing these technologies [4,5].However, the multiscale nature of phase change flows often makes experimental studies formidable.One example is nucleate boiling in the micro-layer regime [6,7], where a thin film, only a few microns thick, forms between the bubble and the heated wall.In this micro-layer, strong heat transfer occurs and contributes significantly to the overall heat extraction.
Conducting quantitative measurements within this micro-layer proves to be a challenging task [8].Therefore, numerical simulation, with its ability to detail physical processes at microscopic levels [9,10], has become a powerful tool for studying phase change flows and has gained increasing interest in academia.In recent years, various computational methods have been developed to simulate multiphase flows with phase change [11].Generally, based on the description the phase interface, these methods can be classified into Front-Capturing and Front-Tracking methods [12].In Front-Capturing methods, the interface movement is implicitly captured via the integration of a tracer function over time, for example, the Heaviside function in Volume-of-Fluid (VOF) methods [13,14], the signed distance function in Level-Set (LS) methods [15,16], and the phase-field equation in Phase-Field (PF) methods [17,18].These methods can automatically handle topology changes and allow for efficient parallelization, making them popular for simulating phase change flows.For a detailed implementation of phase change models based on these methods, readers are referred to Refs.[19][20][21] for VOF, Refs.[22][23][24] for LS, and Refs.[25,26] for PF.
Compared with Front-Capturing methods, Front-Tracking methods [27,28] explicitly track interface motion by advecting connected marker points.
Their capability to track and control topology makes them superior for complex multiscale problems [29].Simulations of phase change flows using Front-Tracking can be traced back to the work of Juric and Tryggvason [30], which investigates 2D film boiling.In this work, an iterative procedure is employed to impose the correct temperature boundary condition at the interface.Esmaeeli and Tryggvason [31,32] later eliminate this iterative process in their improved algorithm, which is subsequently used to study multi-mode film boiling on a horizontal surface.Furthermore, dendrite solidification [33,34] and film boiling from multiple horizontal cylinders [35] are studied using Front-Tracking and the immersed boundary method [36].The three-phase droplet icing problem, including tri-junction and volume change, is computed in Refs.[37,38], where three types of fronts are tracked.Additionally, Irfan and Muradoglu [39] extend Front-Tracking to droplet evaporation driven by species gradients.All these phase change models are based on the Front-Tracking method of Unverdi and Tryggvason [28], which requires manual handling of topology changes.In contrast, the Level Contour Reconstruction Method (LCRM) [40] and the Local Front Reconstruction Method (LFRM) [41] eliminate the need to store logical connections between neighboring surface elements, thus enabling robust simulations of 3D flows exhibiting topology changes with Front-Tracking.They have been successfully applied to film boiling [40], nucleate boiling [42], and rising bubbles with phase change [43,44].
In this work, using the free software Basilisk [45,46], we aim to extend a novel Front-Tracking method, the Edge-Based Interface Tracking (EBIT) method [47,48], to simulate multiphase flows with phase change.The EBIT method, without storing connectivity information, binds Lagrangian markers to the Eulerian grid to achieve automatic parallelization.Its capability and accuracy for multiphase flow without phase change have been validated through typical benchmark problems [48].To extend the EBIT method to phase change flows, we solve additional energy equations with the temperature boundary condition sharply imposed at the interface.The mass flux is computed from the heat fluxes at the interface and is used to solve the adjusted mass and momentum equations.In the presence of phase change, the velocity is discontinuous across the interface.As we employ a collocated grid, on which the cell-centered velocity is approximately projected, unphysical oscillations will be introduced near the interface.To suppress unphysical oscillations, we adopt the ghost fluid method and solve two separate velocity fields.For each phase, the ghost velocity is determined by the real velocity of the other phase and the jump condition due to phase change, thus removing the discontinuity at the interface.The remainder of this paper is organized as follows: In Section 2, we briefly review the EBIT method for completeness.
Then, we introduce how to solve the energy equation and how to incorporate the mass flux into the solution of the mass and momentum equations via the ghost fluid method.Subsequently, in Section 3, we validate the present method with several benchmark tests, followed by concluding remarks in Section 4.

The EBIT method
In this paper, the EBIT method [47,48] is employed for the interface tracking, in which the interfacial markers are restricted to be moved along the grid lines during the interface advection.Note that only one marker is allowed per cell edge in the current work.Let x i denote the position for a given marker, its motion is described by where the interfacial velocity u Γ is obtained from the flow field.By employing the first-order explicit Euler method, the discretization of Eq. ( 1) reads with the superscripts n and n + 1 representing successive discrete time instants.To solve the multi-dimensional interface advection, we adopt the dimension-splitting method [47] and advect the interface along each dimension successively.In each 1D advection, a marker is defined as either an aligned or unaligned marker, depending on whether it is located on the grid line that is aligned with the current advection direction.The example of 1D advection along the x-direction is given in Fig. 1.As shown in Fig. 1(b), the new positions of aligned markers are directly obtained by solving Eq. ( 2).In contrast, an unaligned marker will first be advected by solving Eq. ( 2) (the grey points in Fig. 1(c)), followed by fitting a circle through the surrounding marker points to find the intersection with the unaligned grid line (the red points in Fig. 1(c)).Note that the intermediate unaligned markers that leave the grid lines are finally discarded.
After advection, the connectivity of markers is reconstructed with the help of the Color Vertex method [49], where a color field is adopted to indicate the corresponding fluid phase.As shown in Ref. [48], for the 2D situation, five color vertices (four in the corners and one in the cell center) are sufficient to determine the topological configuration and reconstruct the interface segments without ambiguity.Furthermore, with the Color Vertex method, topology changes are handled automatically.For more details about the EBIT method, readers are referred to Refs.[47,48].

Governing equations
Consider a flow domain occupied by the liquid phase and the vapor phase, which are separated by an evolving interface Γ(t).By assuming the two phases are both incompressible and monocomponent, the governing equations [22] are where u is the fluid velocity, ρ the density, µ the dynamic viscosity, g the gravity acceleration, T the temperature, C p the heat capacity per unit mass A marker is defined as aligned or unaligned according to whether the grid line it is located on is aligned with the current advection direction.
at constant pressure, and λ the thermal conductivity.The source term due to phase change S pc will be given later, and the surface tension term f σ = σκδ s n is computed by using a well-balanced Continuous Surface Force (CSF) method [45,50], where σ is the surface tension coefficient, κ the interface curvature, n the interface normal vector and δ s the Dirac distribution function concentrated at the interface.
Across the interface, fluid properties will change, which can be calculated by where ϕ represents an arbitrary variable, such as density and viscosity.The subscripts liq and vap indicate the physical properties of liquid and vapor, respectively.The Heaviside function H is defined as 1 in the liquid region and 0 in the vapor region.In numerical simulations, Eq. ( 6) is typically computed by adopting the volume fraction [12], which is the integral of H over a given cell.Note that, in this work, volume fractions are directly computed from the reconstructed interface given by the EBIT method.Furthermore, at the interface, the mass flux ṁ in the presence of phase change is defined as With the jump operator [ϕ] = ϕ liq − ϕ vap , Eq. ( 7) can be reformulated and leads to the velocity jump condition and the jump condition for pressure is given by By considering the divergence-free condition in the bulk region of each phase and the velocity jump at the interface, the source term due to phase change S pc [51] becomes where V is the volume of the computational cell and S Γ denotes the area of the interface within it.

Solving the mass and momentum equations
The present work utilizes the free software Basilisk [45,46] and employs a time-staggered approximate projection method to solve the incompressible Navier-Stokes equations.Spatial discretization is achieved using a collocated grid where all variables are stored at the cell center.The advection term is discretized using the Bell-Colella-Glaz (BCG) scheme [52], while the diffusion term is solved implicitly.In cases without phase change, the momentum equation Eq. ( 3) is discretized as follows: where the subscripts c and f denote the cell-centered variable and the facecentered variable, respectively.A second-order accurate scheme [53] is employed to interpolate the variables between the cell center and the cell face, which is denoted by the symbol c → f or f → c.Note that the pressure at the next time step p n+1 in Eq. ( 14) is obtained by solving the Poisson equation with the mass conservation equation is only approximately projected [53], as it is interpolated using Eq. ( 15).
This method works well for flows without phase change [45].However, it may incur numerical instability near the interface in the presence of phase change [53].
When phase change is considered, the mass flux is computed after solving the energy equation, which will be elaborated in the next section.To incorporate the mass flux resulting from phase change into the mass and momentum equations, there are mainly two types of methods: the one fluid method [51] or the ghost fluid method [22].The implementation of the one fluid method is straightforward.During the solution of the Poisson equation, a source term S pc , which is nonzero only in the interface cells (see Eq. ( 10)), is added to the mass conservation equation: The singular source term may lead to numerical oscillations near the interface, especially for the cell-centered velocity u c which is only approximately projected.To avoid such oscillations, Zhao et al. [53] proposed an exact projection method, solving an additional Poisson equation to ensure that u c is exactly projected.
In this work, we employ the more efficient ghost fluid method [22].By setting the ghost velocity, the singularity at the interface is removed.As will be shown in the numerical tests, numerical stability can be effectively improved by employing the ghost fluid method.With the jump condition Eq. ( 9), the ghost velocity is populated by where the color function f C is naturally provided by the Color Vertex technique used in the EBIT method, being 1 in the liquid region and 0 in the vapor region.Eq. ( 17) is only employed in a narrow band region near the interface, which will be discussed in the next section.Note that the ghost velocities will be set for both the cell-centered and the face-centered velocities.
After that, we solve Eqs. ( 11)-( 13) separately for u * f,liq and u * f,vap .Then, for the projection step, Eq. ( 16) is modifed to Note that there is no singular source term on the right-hand side of the mass conservation equation, compared with the one-fluid method, since the velocity jump has been appropriately taken into account by Eq. ( 17) in the ghost fluid method.With p n+1 solved, the face-centered velocity is finally updated by The cell-centered velocity is computed accordingly using Eq. ( 15).To obtain the velocity u Γ used to advect interfacial markers, as illustrated in Fig. 2, we apply the jump condition at the cell center.Then, u Γ is interpolated to the interfacial marker via bilinear interpolation [48].20)), and interpolated to the interfacial marker through bilinear interpolation.

Solving the energy equation
In this section we elaborate the solution of the energy equation.We use the finite difference method to compute the temperature field separately in the liquid domain and in the vapor domain: where k = liq or k = vap denotes the current region of interest.The advection term is discretized by using a third-order accurate WENO scheme [54], and the diffusion term is solved implicitly.When solving the energy equation, we need to correctly impose the Dirichlet boundary condition [51] T Γ = T sat (23) at the interface, where T sat is the saturation temperature.This is implemented in a similar manner as in Ref. [55].The main principle is to fill the uncompleted discretization stencil with an extrapolated value by considering T Γ = T sat .For example, by using the standard second-order accurate central difference scheme for a cell (i, j) which is far from the interface, T xx in ∇ 2 T can be written as When the cell is near the interface, sometimes the stencil cell may across the interface, as shown in Fig. 3(a).The ghost temperature T ghost for the cell (i + 1, j) in Fig. 3(a) can be extrapolated by where θ = (x Γ − x i )/∆x, ranging from 0 to 1. Substituting Eq. ( 25) into Eq.(24) gives leading to a symmetric linear system [56].
As shown in Eq. ( 21) and Eq. ( 22), the solution of the energy equation is split into two steps, utilizing an intermediate temperature field T * for the contribution of the advection terms.This strategy is employed because, during one time step, the cell centers may be swept by the interface, leading to unknown temperature values inside the newly created vapor/liquid cells [22,55].By solving the advection equation Eq. ( 21) with T * and extrapolating it to the ghost region, the unknown values are determined.The extrapolation [22] is achieved by iteratively solving where τ is a pseudo time, and f the variable to be extended.Instead of solving Eq. ( 27) in the whole domain, we solve it within a narrow band near the interface.A cell is defined as being in the narrow band if it or any of its neighbours is cut by the interface.
After solving the energy equation, the mass flux ṁ is calculated as where h lg is the latent heat, and q liq/vap = λ liq/vap ∇T liq/vap • n represents the heat flux at each side of the interface.As the temperature gradient is not continuous across the interface [53], to evaluate the heat flux accurately and robustly, we compute the normal temperature gradients in the two phases separately by adopting an Embeded Boundary Method [57].The computation in the liquid side is sketched in Fig. 3(b).Two points are found along the normal direction at first, where the temperature values are obtained from neighboring cells using bi-quadratic interpolation.Then the normal temperature gradient at the interface is interpolated by where d 1 and d 2 denote the distances from the interface centroid to the two points, respectively.It is noteworthy that the boundary condition at the interface, as specified in Eq. ( 23), is directly imposed in Eq. ( 29).Once the heat fluxes at the interface are determined, the mass flux ṁ in the cells cut by the interface is obtained.To compute the ghost fluid velocities, the mass flux is then extended to the entire narrow band [22] by solving Eq. ( 27).

Overall numerical process
For clarity, the numerical procedures for each time iteration are outlined as below: 1. Solve the advection part of the energy equation Eq. ( 21).
2. Advect the interface using the EBIT method, with the interfacial velocity obtained from Eq. ( 20).
3. Solve the diffusion part of the energy equation Eq. ( 22).
4. Compute the mass flux and accordingly set the ghost velocities using Eq. ( 17).

Numerical Results
In this section, the proposed method is validated through a series of numerical benchmark tests.The codes developed for this work, as well as the configurations for all simulations, are freely available in the Basilisk sandbox [58].A quad/octree based adaptive mesh refinement (AMR) technique [45] is employed to improve computational efficiency.For a given grid level L, the corresponding number of cells in one direction for a cubic domain is In cases involving a non-cubic domain, it is important to note that 2 L always corresponds to the number of cells employed in the direction with the maximum length.

Stefan problem
Here, the 1D Stefan problem is considered, which is widely used to verify phase change models [19,51].As shown in Fig. 4(a), a thin vapor layer is placed between a liquid at saturation temperature T sat and a heated wall with temperature T w .The left boundary (at x = 0) is a solid wall, while the right boundary (at x = l) is set as an outlet.In this case, the vapor layer T !vapor liquid  grows over time due to boiling caused by the heat flux from the vapor side, pushing the liquid out of the domain.An analytical solution is available [51], and the interface position X Γ (t) is given by where α v = λ vap /ρ vap C p,vap represents the thermal diffusivity of the vapor, and χ is the solution of a transcendental equation: During the simulation, the liquid remains at the saturation temperature, while the theoretical temperature distribution within the vapor layer is given by This problem can be characterized by the Jacob number Ja, given by Ja = where ∆T = T w − T sat is the superheat.
Following Ref. [19], the domain length l for this problem is set to 10 mm.
With the Jakob number being Ja = 29.84, the physical properties are set as Initially, the thickness of the vapor layer is 322.5 µm, corresponding to the theoretical solution Eq. ( 30) at t = 0.282 s.The simulation is performed up to t = 10.282 s, and three grid levels, ranging from 4 to 6, are used.In Figs.
4(b) and 4(c), the interface positions and the final temperature distributions, obtained with different grid levels, are compared with the theoretical solutions.Excellent agreement is observed for all grid resolutions.The relative error of the final interface position can be calculated by where the subscripts theo and num represent the theoretical solution and the numerical result, respectively.Fig. 4(d) shows the relative errors for different grid levels, indicating a second-order convergence rate.
Additionally, we also examine the velocity distribution to further validate our method.The velocity distributions obtained by the one fluid method and the ghost fluid method are plotted in Fig. 5.It is shown that a sharp velocity jump appears at the interface due to phase change.With the one fluid method, unphysical oscillations are observed for u c as it is approximately projected [53].In contrast, by using the ghost fluid method, the sharp velocity jump is accurately captured without introducing oscillations near the interface.As discussed earlier, the discontinuity at the interface is removed by setting the ghost velocity according to the jump condition Eq. ( 17), thus improving numerical stability.

Sucking problem
Proposed by Welch and Wilson [59], the 1D sucking problem is another T !"# vapor liquid physical properties for this problem are set as With the superheat ∆T = T ∞ −T sat = 10 K, the Jakob number is Ja = 29.95.
In this case, the domain length l is 10 mm.Initially the interface is located at x = 0.476 mm, which corresponds to the theoretical solution at t = 0.1 s.
We initialize the temperature field according to the theoretical solution at t = 0.1 s and conduct the simulation up to t = 1.1 s.The interface positions, obtained with different grid levels ranging from 6 to 8, are plotted against the simulation time in Fig. 6(b).It is shown that the results converge to the analytical solution with increasing grid refinement.The relative errors of the final interface position at t = 1.1 s are shown in Fig. 6(d), from which a second-order convergence rate can be observed.Furthermore, in Fig. 6(c), the temperature distributions at the end of the simulation are compared against the theoretical solution.The thermal boundary layer is well resolved for the simulation at grid level 8.

Bubble evaporation/condensation with a constant mass flux
A 2D bubble growing with a constant mass flux is considered in Ref. [22] to verify the solution of the mass and momentum equations in the presence of phase change.Here, we also simulate this problem, additionally including a condensation setup.The computational domain is a square of (a) (b) Note that there is no need to solve the energy conservation equation, since the mass flux is manually imposed.We use ṁ = 0.1 kg/m 2 • s for the evaporation case and ṁ = −0.1 kg/m 2 • s for the condensation case.Theoretically, the bubble radius will evolve linearly with time [22]: where R 0 is the initial radius of the bubble, set to 0.1 m and 0.

Scriven problem
The growth of a static bubble at saturation temperature T sat in a superheated liquid with temperature T ∞ is considered to validate the present method.This setup is known as the Scriven problem [22,53], named after  the theoretical solution obtained by Scriven [60].As illustrated in Fig. 9, it is computed in an axisymmetric configuration.The computational domain is a square with length 160 µm.The left boundary is set as the axis of symmetry, and an outflow boundary condition is applied to the top and right boundaries.
To improve computational efficiency, a symmetry boundary condition is employed on the bottom boundary, allowing the computation of only a quarter of the bubble.With the Jakob number Ja = 5.99, the physical properties are set as For this problem, the theoretical evolution of the bubble radius [60] is given by where α l = λ liq /ρ liq C p,liq is the thermal diffusivity of the liquid, and the growth constant β is obtained from solving the equation During the bubble growth, the temperature inside the bubble remains at the saturation temperature T sat while the temperature distribution within the liquid can be theoretically calculated by In the simulations, the bubble radius is initially set to 50 µm, corresponding to the theoretical value at t = 94.7 µs.The temperature is accordingly initialized based on the theoretical solution Eq. ( 42).This problem is computed with four different grid levels, ranging from 6 to 9, and the final time is set as t = 378.8µs.The time histories of the bubble radius, obtained with different grid levels, are presented in Fig. 9(b).It is observed that the numerical results converge to the theoretical solution.Quantitatively, the relative errors of the bubble radius at the final time are plotted against the grid level in Fig. 9(d), indicating a second-order convergence rate.Furthermore, Fig. 9(c) presents the temperature distributions along the z-axis at the end of the simulation for different grid levels, as well as the theoretical solution.It is evident that a grid level of 8 is adequate for obtaining a converged result, which agrees well with the theoretical solution.The velocity vectors and the superheat distributions obtained by the one fluid method and the ghost fluid method at the final time are compared in Fig. 10.For the simulation with the one fluid method, the spurious currents are so strong that the spherical shape of the thermal boundary layer cannot be maintained during the bubble growth.In contrast, weaker spurious currents are observed for the result obtained using the ghost fluid method, and the spherical shape of the thermal boundary layer is also well preserved.

Bubble rising in superheated liquid
Subsequently, the bubble growth in a superheated liquid under gravity is considered to validate the present method for problems involving coupled effects of phase change and buoyancy.During the rising, the bubble expands due to phase change at the interface, and deforms due to the interaction between surface tension, buoyancy, and advection effect.As a result, the bubble growth rate and the temperature distribution within the flow region cannot be described by the Scriven solution used for a spherical symmetric system [22,51].For the simulations, we follow the experiment setup of Florschuetz et al. [61], which is also simulated by Bureš and Sato [51,62] to verify their phase change model on axisymmetric grids.We consider an ethanol system at atmospheric pressure, with the following physical properties: solution is a very good approximation as the bubble is so small that the effect of buoyancy is negligible and that its spherical shape is well preserved.
Hence, initially, a bubble with a diameter of 420 µm, which corresponds to the theoretical solution at t = 0.0056 s, is placed 1 mm above the bottom boundary.The temperature field is initialized based on the theoretical solution described by Eq. ( 42).The gravitational acceleration, acting in the negative z-direction, is set to 9.81 m/s 2 .
The simulations are conducted with different grid levels, ranging from 10 to 12, up to t = 0.0856 s.Fig. 13 shows the interface shape and superheat distribution at different time instants.As the bubble rises, its initially spherical interface deforms into an ellipsoidal shape, and a cooler region forms at the wake of the bubble.During the bubble growth, the dimensionless Péclet number can be used to compare the heat transfer resulting from convection and diffusion: where R is the effective radius and u r is the rising velocity, which are computed by and respectively.Fig. 12(a) shows the evolution of the Péclet number as a function of time for different grid levels, confirming the convergence of our results.
It can be observed that the Péclet number increases from the initial zero value to about 3000 at the end of the simulation.This trend suggests that the early stage of the process is dominated by heat diffusion, which is reasonable given that the bubble is released from rest with a small initial radius.As time progresses, the bubble grows larger and its velocity increases, leading to an increasing influence of advection on the heat transfer [62].The time history of the normalized effective radius is compared with the experimental result of Florschuetz et al. [61] and the numerical result of Bureš and Sato [62].
The normalization factor is 2β √ α l , with β being the Scriven growth constant defined by Eq. ( 41).It is shown that our results agree with the numerical solution of Bureš and Sato [62], which was confirmed to be converged in their study.The finest grid resolution used by Bureš and Sato [62]   The last case is the film boiling problem, in which the Rayleigh-Taylor instability at the interface will be triggered by buoyancy [63,64].As shown

Film boiling
t = 0.42 s t = 0.48 s in Fig. 14(b), a superheated vapor layer is placed between a heated wall with temperature T w and a liquid at saturation temperature T sat .An axisymmetric configuration is simulated within a domain with a size of [0, where λ 0 is the unstable Taylor wavelength and is calculated by The bottom boundary is a no-slip wall, while the outflow boundary condition and the symmetry boundary condition are applied to the top and the right boundries, respectively.The initial interface shape is given by and the physical properties are set as Initially, the velocities of the two fluids are zero, and a linear temperature distribution increasing from T sat on the liquid-vapor interface to T w on the wall is imposed inside the vapor layer.The gravitational acceleration, acting in the negative z-direction, is 9.81 m/s 2 .Note that the current setup results in a Jakob number Ja of 8 and a Taylor wavelength λ 0 of 0.0787 m.
The simulations are performed up to t = 1.6 s, with increasing grid levels from 8 to 10. Fig. 15 shows the interface shapes and superheat distributions at different time instants for grid level 10.It is observed that, the vapor bubble periodically forms and grows due to the boiling and Rayleigh-Taylor instability, and eventually detaches from the vapor layer.This observation is consistent with the conclusion from previous studies that film boiling is a quasi-steady phase change phenomenon [63,64].The interfaces and time histories of the vapor volume, normalized to be dimensionless, are plotted in Figs.16(a), (b), and (c) for different grid levels, confirming the convergence of our results.The dimensionless vapor volume V is computed as follows: where V 0 is the initial vapor volume and V bc represents the vapor volume flowing out from the top boundary.Additionally, we emphasize that, unlike the pure 2D case simulated in Ref. [20], the bubble will pinch off in the axisymmetric configuration due to surface tension, as seen in Fig. 16(b).Interested readers can refer to Ref. [55] for a detailed analysis.The quantitative study is performed using the space-average Nusselt number Nu where λ ′ = σ (ρ liq −ρvap)g is the characteristic length.In Fig. 16(d), the time history of the Nusselt number obtained at grid level 10 is compared with the value of 1.91 predicted by the Klimenko correlation [65].The time-averaged Nusselt number from the numerical simulation is 1.695, indicating an 11.2% deviation from the theoretical prediction.

Conclusion
We have extended the Edge-Based Interface Tracking (EBIT) method [47,48] to the simulations of multiphase flows with phase change.By using the EBIT method, the interfacial markers are restricted to move along the grid lines, so that automatic parallelization can be achieved.Based on the EBIT method, we additionally solve the energy equations for each phase to include phase change effects.The Dirichlet boundary condition for the temperature at the interface [51] is sharply imposed using the geometric information provided by the EBIT method.Moreover, the mass flux is determined according to the heat fluxes at the interface, which is then incorporated in the solution of the mass and momentum equations.In the presence of phase change, the velocity is not continuous across the interface, leading to numerical instability.This issue is worse on the collocated grids as the cell-centered velocity is approximately projected [45].Instead of solving a second Poisson equation [53] to suppress such oscillations, we have shown that this issue can be addressed by using the ghost fluid method [22].In the ghost fluid method, two velocities are employed for each phase, and the coupling between two phases is achieved by populating the ghost velocities.
As the ghost velocity is set according to the jump condition, the discontinuity at the interface is removed, thus the numerical stability of the present method is improved.Several benchmark problems have been simulated to validate the present method.It is shown that the present method agrees very well with the theoretical solutions and the experimental results.The present method has been implemented in the free software Basilisk [46].The developed codes, along with the configurations for all tests, are freely available in the Basilisk sandbox [58].Currently, the 3D EBIT method is under development.The 3D simulations of phase change flows with the EBIT method is one of our future works.In addition, as Front-Tracking methods track the interface without diffusion and accurately provide the exact position of the interface, they are very suitable for coupling with multiscale models on the interface.For instance, in Ref. [66], a mass boundary layer model is employed on the interface, improving the computation of mass transfer from buoyant bubbles.We are interested in coupling a subgrid model with the EBIT method to more effectively resolve the thin thermal boundary layer in the flows with high Jakob numbers.Applying the present method to the study of nucleate boiling is also one of our main goals.

Figure 1 :
Figure 1: Illustration of the 1D advection of the EBIT method along the x-direction: (a) Initial interface.(b) Aligned markers after advection (blue points).(c) Unaligned markers after advection (grey points) and the newly created markers intersecting with grid lines (red points).Note that the grey points will be discarded.(d) Interface after 1D advection.

Figure 2 :
Figure 2: Schematic of the computation of interfacial velocity: The green and red points represent the vapor and liquid cell centers, respectively, while the yellow point denotes the interfacial marker.The fluid velocities (indicated by the green and red arrows) are obtained by solving the momentum equation.These velocities are then converted to the interfacial velocity (denoted by the yellow arrow) using the jump condition (Eq.(20)),

Figure 3 :
Figure 3: (a) The discretization stencil of ∇ 2 T for a cell near the interface.(b) The computation of the normal temperature gradient at the liquid side for an interfacial cell.

Figure 4 :
Figure 4: Stefan problem: (a) Schematic of the 1D Stefan problem.(b) Time history of the interface position.(c) Temperature distribution at t = 10.282s.(d) Relative error of the interface position on different grid resolutions.Grid levels 4 to 6 correspond to grid resolutions ranging from 16 × 1 to 64 × 1.

Figure 5 :
Figure 5: Stefan problem: The cell-centered velocity uc and face-centered velocity uf distributions across the interface at t = 10.282 s, obtained by the one fluid method (a) and the ghost fluid method (b).
popular benchmark test used to validate phase change models under superheated conditions.As shown in Fig.6(a), a superheated liquid with temperature T ∞ is adjacent to a vapor layer at saturation temperature T sat .The left boundary (at x = 0) is a solid wall with the temperature fixed at T sat , and the right boundary (at x = l) is an outlet.Due to the superheat, the liquid starts to boil, causing the volume of vapor to expand and push the interface rightwards.Compared with the Stefan problem, the sucking problem is more challenging because of the formation of a thin thermal boundary layer close to the interface on the liquid side.For this problem, a theoretical solution can be obtained from the similarity solution[59].Following Ref.[53], the

Figure 6 :
Figure 6: Sucking problem: (a) Schematic of the 1D sucking problem.(b) Time history of the interface position.(c) Temperature distribution at t = 1.1 s.(d) Relative error of the interface position on different grid resolutions.Grid levels 6 to 8 correspond to grid resolutions ranging from 64 × 1 to 256 × 1.

Figure 7 :
Figure 7: Bubble evaporation: (a) Comparison of the theoretical interface shape with the numerical results obtained with different grid levels.(b) Relative error of the bubble radius on different grid resolutions.

Figure 8 :
Figure 8: Bubble condensation: (a) Comparison of the theoretical interface shape with the numerical results obtained with different grid levels.(b) Relative error of the bubble radius on different grid resolutions.
2 m for the evaporation and condensation cases, respectively.The simulations are performed with increasing grid levels from 5 to 7, and the final time is t = 0.01 s.Figs.7(a) and 8(a) present the interfaces at the end of the simulation for different grid levels, which converge well to the theoretical solution.The circular shape of the bubble is well preserved for all grid levels.Additionally, the relative errors of the final bubble radius at different grid resolutions are plotted in Figs.7(b) and 8(b), confirming the second-order convergence rate of the present method.

Figure 9 :
Figure 9: Scriven problem: (a) Schematic of the Scriven problem.(b) Time history of the bubble radius.(c) Temperature distribution along the z-axis at t = 10.282s.(d) Relative error of the bubble radius on different grid resolutions.

Figure 10 :
Figure 10: Scriven problem: The velocity vectors and the superheat distributions (T −T sat ) at the end of the simulation, obtained by the one fluid method (Figs.(a) and (b)) and the ghost fluid method (Figs.(c) and (d)) at the grid level L = 6.Only a quarter of the bubble is simulated in an axisymmetric configuration, and the results are mirrored about the axes r = 0 and z = 0 for better visualization.

Figure 12 :
Figure 12: Bubble rising in superheated liquid: (a) The Péclet number for the bubble as a function of time obtained at different grid levels.(b) Time history of the normalized bubble radius, compared with the experimental data of Florschuetz et al. [61] and the numerical result of Bureš and Sato [62].

Figure 13 :
Figure 13: Bubble rising in superheated liquid: The interface shapes and superheat distributions at different time instants, obtained with the grid level 12.The results are mirrored about the axis r = 0 for better visualization.

Figure 14 :
Figure 14: Schematic of the film boiling.

Figure 15 :
Figure 15: Film Boiling: The interface shapes and superheat distributions at different time instants, obtained with the grid level 10.The results are mirrored about the axis r = 0 for better visualization.

Figure 16 :
Figure 16: Film Boiling: (a) and (b) are the interface shapes obtained at different time instants with different grid levels.Note that the interfaces are mirrored about the axis r = 0 for better visualization.(c) The time variation of the dimensionless vapor volume, which is normalized by the initial vapor volume.(d) The comparison of the numerically computed Nusselt number at level 10 with the Klimenko correlation.