A SOLUTION TO STEFAN PROBLEM USING EULERIAN TWO-FLUID VOF MODEL

A novel approach for the solution of Stefan problem within the framework of the multi-fluid model supplemented with Volume of Fluid (VOF) method, i.e. two-fluid VOF, is presented in this paper. The governing equation set is comprised of mass, momentum and energy conservation equations, written on a per-phase basis and supplemented with closure models via the source terms. In our method, the heat and mass transfer is calculated from the heat transfer coefficient, which has a fictitious function and depends on the local cell-size and the thermal conductivity, and the implementation is straightforward because of the usage of the local value instead of a global parameter. The interface sharpness is ensured by the application of the geometrical reconstruction scheme implemented in VOF. The model is verified for three types of computational meshes including triangular cells, and good agreement was obtained for the interface position and the temperature field. Although the developed method was validated only for Stefan problem, the application of the method to engineering problems is considered to be straightforward since it is implemented to a commercial CFD code only using a local value; especially in the field of naval hydrodynamics wherein the reduction of ship resistance using boiling flow can be computed efficiently since the method handles phase change processes using low-resolution meshes.


Introduction
The sustainable development and design in industrial sectors concerned with high heat outputs, e.g., in power engineering and metallurgy, are strongly dependent on well-optimized processes that involve boiling phenomena. Thus, for example, in power engineering nucleate boiling is a favorable boiling regime because of high heat flow rates that can be achieved during cooling inside pressure vessels like boilers, thereby avoiding dangerous overheating of the heat transfer surface that may occur if the liquid coolant cannot reach that surface. This unfavorable boiling heat transfer regime is known as film boiling and, in contrary, is very common in metallurgy, where achievement of desired mechanical properties of metal material via quenching process is followed by the occurrence of vapor blanked around the hot object due to high heat flow rate released in the quenchant liquid.
In addition to these classical industrial appearances of boiling processes is the application of the process in the field of naval architecture, where it can be addressed to, already well-established, air lubrication systems used in reduction of ship resistance; and whose theoretical foundations are given in Pavlov et al. [1]. To reduce the ship resistance, the bubble injection technique under a vessel has been actively studied in recent years, and it was applied to full-scale ships. For example, Kodama et al. [2] reported approximately 10-15 % of saving of the total energy consumption for an experimental ship, and Mizokami et al. [3] achieved at most 12 % fuel saving. In these systems, blowers are used to inject air below a hull, which requires energy input. From the numerical point of view, two distinct approaches may be found in the literature in this regard: the more fundamental ones, dealing with micro-bubble drag reduction (MBDR); and the macroscopic studies aimed to study the flow features around the bodies. An example of the former is the numerical simulation of an MBDR on a flat plate conducted using the two-fluid model in Mohanarangam et al. [4]. In the latter case, a scaled bulk carrier is a subject of the study by Sindagi et al. [5]. Using the single-fluid VOF method the ship resistance is obtained for two different ship speeds; different air injection rates and different diameters of the air injection holes. The references on previous work carried out in both assessments, micro-and macro-scale can also be found therein. If a heating system is implemented at the fore-bottom of the vessel and boiling takes place there as illustrated in Fig. 1, the reduction of ship resistance due to the bubbles and the low viscosity of water with higher temperature may be achieved. The vapor bubbles condensate in the downstream because of lower temperature, but the non-condensable gas bubble may remain. As the heat source, the heat energy from condenser in a steam power plant on board could be used in case that the plant is equipped. Recent research shows that a few volume percent of vapor bubble can reduce up to 45 % of drag [6]. However, the publications that study the latter effect, that is, the reduction of the ship resistance involving boiling flow using the numerical simulation were not found in the literature.
To optimize such processes efficiently, reliable numerical simulations are required to estimate boiling heat transfer using affordable computational resources and time. Furthermore, the numerical simulations have to be capable to compute various flow situations and to reflect properly the physics underlying the boiling phenomena in the sense that the dominant mass transfer rate occurs at the interface between the phases. The development of computational methods for application in the former case (nucleate boiling) has been subject of extensive studies in the past, simulating thereby the boiling curve from the side of nucleate boiling on. This approach is utilized, for example, in Sato and Ničeno [7]. In the study presented by the Heat source Drag reduction due to vapor and non-condensable gas bubbles Drag reduction due to low viscosity of hot water Resistance reduction due to low viscosity of hot water Resistance reduction due to vapor and non-condensable gas bubbles A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno 143 authors, the boiling curve is being modeled from nucleate to film boiling, with consideration of the Critical Heat Flux (CHF), conjugate heat transfer and bubble micro-layer. However, the reliable computational models are still needed for the latter case (film boiling), where the transition from the initial film boiling to succeeding nucleate boiling regime needs to be modeled. Therefore, the long-term goal of our research is to develop the computational method which can simulate the transition from film boiling to nucleate boiling.
In development of such models, the first step is the verification of Stefan problem: the one-dimensional phase-change problem. This moving boundary problem provides the information about unsteady evolution of liquid-vapor interface and the corresponding temperature field throughout the domain composed of those phases, thus being suitable for verification of the applied interface tracking algorithm and accompanied mass transfer model. Stefan problem has been typically used as a benchmark in the boiling process, which consists of superheated vapor and saturated liquid. The examined numerical methods are the interface tracking methods, such as, for example, single-fluid Volume-of-Fluid (VOF), coupled with Direct Numerical Simulation (DNS) of turbulence, in which the volume fraction of vapor/liquid is tracked by solving the advection equation. Sun et al. [8] proposed the energy jump model, and they obtained solutions for different density ratios, assuming zero thermal conductivity of phase at saturation temperature and equal specific heat capacities of both phases. Stefan problem is used in verification of VOF-based phase change model proposed by Pan et al. [9]. The model is based on fixing the interface temperature to saturation temperature and was verified for three different working fluids: HFE7100, R113 and water. A systematic introduction in the solution of Stefan problem using both, the analytical and numerical approach, was given in Perez-Raya and Kandlikar [10]. The authors presented the necessary steps for obtaining the solution using VOF and presented the theoretical framework that underlies final expressions used for verification purposes. The accuracy of three different mass transfer models, namely: Lee [11], Rattner [12] and Sun model [13], in solving Stefan problem was investigated by Kim et al. [14]. The transient interface evolution was examined therein, using different mesh resolutions, time increment and, in Lee model, model constants. Perez-Raya and Kandlikar [15] used the sharp interface approach for the solution of Stefan problem. In their approach, the saturation temperature is imposed in the interface cells and the estimation of temperature gradient is calculated from the distance between the interface and neighboring cell centers. More recently, the idea of incorporation of a gradient of volume fraction field into a model that relies on empirical constant is proposed in Chen et al. [16]. By considering the cell size, the Chen's approach can be applied to non-uniform quadrilateral mesh. The other studies that utilize interface tracking approach, and some of them apply it in solution of Stefan problem, can be found in the review of Kharangate and Mudawar [17].
Compared to the abovementioned single-fluid models, a two-fluid model is considered to be able to predict more accurate results on coarser mesh owing to the unnecessity of resolving the interface geometry/shape, though two sets of momentum and energy equations must be solved. However, due to the averaging process in the formulation of the governing equations, the shape and area of the liquid-vapor interface are lost, and the dynamics related to the interface, e.g. surface tension and interphase forces, cannot be directly computed. To mitigate this, several approaches to handle interfaces inside the two-fluid model were proposed; a review can be found in Mer et al. [18]. Among the modeling approaches presented in [18], further exploration in the paper is made on the following models, namely: Large Interface Model (LIM), Generalized Large Interface Model (GLIM) and Large Bubble Model (LBM). LIM and LBM were applied in the solution of Stefan problem by Fleau [19]. The description of LIM is given in Coste [20]. LIM identifies the large interface, but does not reconstruct it. In order to determine the distance from interface to the adjacent cell, the method utilizes a three-cell stencil A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno concept: one cell contains the interface, while other two are occupied by single representatives of the involved phases. The general idea behind LBM can be found in Dènefle et al. [21]. LBM is three-or four-field approach in the framework of two-fluid model, that considers two continuous and one or two dispersed phases. The distinction between the continuous phases, that are simulated, and dispersed phases, that are modeled, is made with a spatial cutting length, while the mass exchange between the continuous and dispersed field of the same phase is established via the mass transfer terms. However, to the authors' knowledge, there appears not to exist evaluation of the two-fluid VOF for Stefan problem in the literature. In this paper, the two-fluid VOF model available in CFD package Ansys® Fluent, Release 20.1 as Multi-Fluid VOF model is adopted for the solution of Stefan problem. The advantage of the two-fluid VOF over the single-fluid VOF is in avoidance of averaging of the material properties in the energy equation, which reduces the accuracy in the heat-flux calculation [15].
The novelty of this paper is the application of two-fluid VOF method to Stefan problem and the validation of the model to not only quadrilateral meshes but also triangular meshes. The originality of our model is the usage of local mesh size for the heat transfer coefficient, which is explained in the section of the numerical model. The paper is organized as follows. The theoretical solution of Stefan problem is presented in Chapter 2. The proposed computational method is the content of Chapter 3. Chapter 4 outlines the details related to performed numerical simulations. The obtained results together with the applied error measures are discussed in Chapter 5.

Theoretical solution of Stefan problem
The analytical solution of Stefan problem, presented here, refers to expressions for computation of interface position and temperature distribution during planar interface evaporation in the system composed of superheated vapor and saturated liquid, reinterpreting thereby the work of Perez-Raya and Kandlikar [10]. The equations involved in derivation process are the energy balance at the interface and one-dimensional heat conduction equation in a continuum. In what follows, a concise exposition of the theoretical framework underlying the applied analytical expressions is given.

Conditions at the interface
In the analytical solution used in this study, the interface condition that accounts only for occurrence of heat conduction in the system is incorporated. This is because in the case with superheated vapor and saturated liquid the heat transfer takes place only in vapor phase, Hence, the interface condition is given by: where ′( ) is the interface velocity, is the vapor density, ℎ is the latent heat, is the thermal conductivity of the vapor, while the negative signed gradient denotes heat flux in normal direction to the interface. Therefore, the mass generated in phase change process is proportional to the balance of the heat fluxes applied from each phase in normal direction to the interface. In order to derive analytical model for this moving boundary problem, introduction of similarity variable, , is required to relate the spatial coordinate, , and time, , as: Written in terms of the interface position, that is, = ( ) and = ξ, this relation yields the equation for determination of interface displacement: where is the constant that has to be determined using iterative procedure. The derivation of Eq. (4) with respect to time gives the interface velocity: The substitution of Eq. (5) into Eq. (2) and involvement of solution that is obtained from incorporation of similarity variable in one-dimensional unsteady heat conduction equation gives: where is constant that has to be determined using an iterative procedure, is the specific heat capacity of the vapor, is the temperature of the superheated wall and ′ is the saturation temperature. The constant embeds the value of similarity variable at the interface, that is, the constant , and is defined as: where is the thermal diffusivity of the vapor phase. Substitution of Eq. (7) into Eq. (4) yields the instantaneous interface position in the case with superheated vapor and saturated liquid [10]: In the cases where initial vapor layer is needed for computational purposes, the necessary time shift of the simulation time is accomplished using Eq. (8).

Temperature field
The analytical solution of the temperature field inside the system composed of two phases, vapor and liquid, distinguished by a sharp interface, is derived from the partial differential equation (PDE) for one-dimensional unsteady heat conduction in a single continuum. The time and space coordinates present therein are transformed using the similarity variable introduced in the previous section. Two boundary conditions are required to obtain the A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno solution of the PDE: at the interface and at the wall. The temperature distribution at certain time instance, , is computed as:

Numerical method
Within the framework of two-fluid model adopted here for solution of laminar two-phase flow with phase change, two separate sets of conservation equations are solved in each computational cell; one for each phase present in the domain. The theoretical foundations underlying the two-fluid model are given in Drew and Passmann [22], while the model equations available in Fluent and applied in this study are taken from [23]. The communication between the two phases is established via source terms in the equations, referred to as phase interaction terms. In addition, a single pressure field is shared between the phases. The temperature in the liquid phase is assumed to be constant at the saturation temperature.

Conservation equations
The evolution of planar liquid-vapor interface due to evaporation, studied here, is described with three laws of classical physics, namely: the conservations of mass, momentum and energy: ( ℎ ) + ∇ ⋅ ( ⃗⃗ ℎ ) = + ̿ : ∇⃗⃗ + +̇ℎ (12) where subscripts and refer to particular phase and the phase pair involved in certain interfacial transfer process, respectively. Appropriate closure models are required for the interphase transfer terms that appear at the right-hand side in these equations, those are, interfacial mass transfer, ̇, the momentum transfer, ⃗⃗ , and the heat transfer, , terms. These terms stem from averaging process used in derivation of the set of conservation equations and more details regarding the closures applied therein is given in the forthcoming section.

Interfacial area
The interphase transfer terms in the presented set of conservation equations for a continuum composed of interpenetrating phases have to be supplemented with suitable correlations. A common issue, however, in all the terms is related to modeling of the interfacial area, a flow regime dependent ratio of dispersed phase's projected area and cell volume. When free surface flow modeling is considered, the interfacial area per unit volume is computed as: This idea is also inherited for free surface modeling within Algebraic Interfacial Area Density (AIAD) modeling approach for handling different interface scales inside two-fluid A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno 147 model, presented by Höhne and Vallée [24]. Since the application of Eq. (13) yields the application of closure models only within the interface zone, the closures used in this study can be referred as "locally imposed closures".

Interphase momentum transfer
In general, the interphase momentum transfer consists of the drag force, lift force, wall lubrication force, turbulent dispersion force and virtual mass force. Since in the solution of Stefan problem there is a distinction between the phases with a free surface, the incorporation of drag force is necessary according to the study of stratified two-fluid modeling by Štrubelj et al. [25]. The role of drag force, as noted by the authors, is to ensure equal velocities of the phases within the interface region, and the physical basis of the correlation applied for computation of the drag force coefficient is not mandatory since the drag force term has no real physical meaning in a stratified free-surface flow as it has in a dispersed flowits presence in the momentum equation is due to averaging process used in derivation of two-fluid model. Thus, the drag force is only taken into account in this paper. The general formulation of interphase drag is given by: where is drag force coefficient and ⃗⃗ − ⃗⃗ is the relative velocity between the phases. In this study, the drag force coefficient is modeled using the default anisotropic drag model available in Fluent. This model allows for higher drag in the normal direction to the interface and lower drag in the direction tangential to the interface, with the default values of normal and tangential interfacial drag friction factors equal to 1e+6 and 1e+3, respectively. In addition, application of this model contributes to the stability of the computation [26].

Interphase heat transfer
The relation for computation of heat transfer across the interface between the phases reads: where ℎ is the heat transfer coefficient, is the interfacial area and − is the temperature difference between the interacting phases. The heat transfer coefficient at the interface is given by: where is the thermal conductivity of the liquid, while and are Nusselt number in the vapor phase and the user-prescribed characteristic dimension, respectively. In the present study, the heat transfer coefficient, ℎ , is modeled using two-resistance model; wherein the different correlation or expression of the heat transfer coefficient is applied at each side of the vapor-liquid interface.
At the vapor side, the heat transfer coefficient is computed using cell-size dependent Nusselt number correlation, which we proposed in this paper: where is the user-prescribed vapor phase characteristic dimension and Δ is the approximated cell size. The derivation of Eq. (17) is given in Appendix. Note that the vapor phase diameter, , defined in Eq. (17) is cancelled, when it is substituted into Eq. (16). Thus, an arbitral value can be given as . The cell size is calculated as: where is the cell volume.
At the liquid side, the saturation temperature is maintained at the interface using zeroresistance model. This model refers to application of a large heat transfer coefficient, hl → +∞, to ensure the saturation temperature at the interface throughout the computation process. This idea corresponds to application of large coefficient method to fix the saturation temperature at the interface, as discussed in the work from Perez-Raya and Kandlikar [15].

Interphase mass transfer
The mass transfer at liquid-vapor interface is modeled using the thermal phase change model implemented in Fluent that reads: where is the interfacial area computed with Eq. (13), ′ is saturation temperature, ℎ is the latent heat of vaporization, and are, respectively, vapor and liquid phase's scaling factors, ℎ and ℎ are, respectively, the heat transfer coefficients at vapor and liquid side of the interface, while and are temperatures of vapor and liquid phase, respectively. In Eq. (19), Cv and Cl equal to unity theoretically, although they are treated as the tuning parameters in Fluent. Since we intend to obtain a solution of the governing equations without introducing any tuning parameters except computational mesh, they are set to be unity in this paper.

Discretization of equations
The interpolation near the interface is accomplished using Geo-Reconstruct scheme. Furthermore, the momentum and energy equations are discretized using second-order upwind interpolation scheme. The pressure equation is solved with PRESTO! method [23], while the gradient discretization is carried out using least square cell based method. The coupling between pressure and velocity fields is realized using phase coupled SIMPLE algorithm. The governing equations are discretized in time using the implicit first-order Euler method.
A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno

Material properties
Since our final objective of the simulations is an application to engineering problems, we select the water as the working fluid instead of an artificial material. The properties of the water at 1.013 bar are listed in Table 1. The saturation temperature is 373.15 K at this system pressure.

Computational domain and boundary condition
The computational domain is designed as a beam divided into finite number of cells in the x-direction of the two-dimensional (2-D) Cartesian coordinate system and with one cell layer in the y-direction, i.e. a special case of 2-D domain that is also known as 1.5-D mesh. A schematic view of the computational domain together with the boundary condition and initial phase distribution is shown in Fig. 2. The left boundary of the computational domain is the wall, where the constant superheat of 10 K, ΔT, is prescribed for all the simulation cases. The symmetry condition is imposed on the top and bottom boundaries to ensure the one-dimensionality. The right boundary is defined as the pressure outlet, with prescribed value of volume fraction and the bulk liquid's temperature.

Initial condition
The vapor volume fraction is set to unity in the first cell adjacent to the superheated wall, while in the remainder of the domain is set to zero; a vice versa approach is applied for liquid volume fraction as depicted in Fig. 3a. Here Δ is the cell length of the wall adjacent cell. The velocity fields of both phases are initially set to zero in all the cells throughout the domain. The initial vapor temperature inside the vapor phase is defined by linear profile changing from T´+ΔT to T´ as illustrated in Fig. 3b. The initial liquid temperature is, on the other hand, set to T´ in the whole domain.

Computational mesh
The presented numerical model for solution of Stefan problem using two-fluid VOF model has been assessed to three different types of computational meshes: (a) uniform quadrilateral, (b) stretched quadrilateral and (c) the hybrid of stretched quadrilateral and triangular cells (hereinafter also named hybrid mesh), as shown in Fig. 4, respectively. The mesh is generated by using the GMSH code [28]. The first mesh type is referred as the standard mesh, while the last two types serve to prove the capability of the proposed model under more realistic conditions regarding computational domain, that is, when the gradients of the dependent variables are computed using graded mesh in wall adjacent region. With respect to mesh resolution, the model verification is carried out using meshes with 50, 100 and 200 cells, hereinafter also named coarse, medium and fine, respectively. The graded mesh, where included, is composed of ten quadrilateral cells with expansion ratio of 1.2, while the number of cells with the dominant cellshape has been varied.
The time steps applied in the simulations are mesh resolution dependent and written in Table 2. The total simulation time is divided in two parts. Firstly, the simulation is run for 0.1 ms to obtain more accurate initial distribution of interface displacements; otherwise, an initial discrepancy from the analytical solution is observed that shifts the normal interface evolution to later time. Thus, two-orders of magnitude lower time step size than in the A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno remainder of simulation is used together with the number of time steps that is 10 % of the number of time steps in the second part of the simulation. Then, the simulation is run for 0.1 s with increased time step size. Furthermore, due to presence of initial vapor cell in the simulation, the simulation time has to be shifted for a constant value obtained from Eq. (8). Therefore, the results at specific time instances are obtained using linear interpolation from neighboring values. The simulations carried out using medium and fine hybrid meshes required, however, finer time step sizes, listed in Table 3.

Implementation into Fluent
Since the two-fluid modeling approach utilized the concept of interpenetrating continua, in each cell of the computational mesh the relevant dependent variable's has to be assigned for both phases. Thus, three field variables are initialized on a per phase basis using Define-Init User-Defined Function (UDF) macro: volume fraction, velocity and temperature. Furthermore, the proposed correlation for computation of the heat transfer coefficient, Eq. (16) is incorporated into the Fluent code using Define-Exchange-Property UDF macro. Finally, in the cases where the hybrid mesh is used, the interface position is calculated using Define-Report-Definition-FN UDF macro. In all the computations, the default relaxation factors for the iterative calculations implemented in Fluent were used.

Definition of the interface position
In the computations performed with the meshes composed of quadrilateral cells, the interface is considered as an iso-surface of vapor volume fraction equal to 0.5, that is, = 0.5. Thus, the information on interface position is obtained by tracking the x-coordinate of the isosurface. In the case of hybrid mesh, inside the mesh region composed of triangular cells, this definition of interface, however, leads to interface distributed over two or more cells in a A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno 152 piecewise-linear segments. Therefore, a cell-shape independent approach is used, with the interface position given by: where , is the vapor volume fraction in i-th cell, , is i-th cell's volume and is the length of the domain. The sum in denominator of Eq. (20) is the volume of the domain.

Definition of computational errors
We define the computational errors of the interface position and the temperature field here. The relative error of the interface position at time t is given by: where , ( ) and , ( ) refer to the x-coordinate of the interface position obtained from the analytical solution and the numerical simulation, respectively. In this study, the relative error of the interface position is calculated for two time instances: in the middle of the simulation, = /2, and at the end of the simulation, = . We evaluate the relative error in the temperature field at the end of the simulation time which is defined as: where ( , ) is the error of the temperature at the cell center of ci-th cell calculated by: where ( , ) and ( , ) are the temperatures at the location of the cell center, x, for the analytical solution and the numerical simulation, respectively, and Δ is wall superheat. The temperature of the vapor is used for Tn, although in the two-fluid modeling approach two energy equations are solved and two temperature fields are present in the domain. This is because the temperature field calculated from the energy equation of the liquid phase is constant at the saturation temperature.

Uniform mesh
The advancement of vapor-liquid interface in the performed numerical simulation reproduces the analytical solution with a commensurate level of accuracy. This is shown in Fig. 5a. Additionally, the obtained result is supplemented with ± 2 % error band in order to make a stringent estimation of solution accuracy, as indicated in Figs. 5b and 5c. The comparison between Fig. 5b and Fig. 5c, together with exact values given in Table 4, points out that increased mesh resolution does not necessarily provide the result that is more accurate. This is in accordance with the study in Gauss et al. [29], where is shown that, when two-fluid model is used in typical interface tracking case, the solution is not fully dependent on the mesh size. Here, close to the start of the computation, within the period from 10 % to 15 % of the total simulation time (Fig. 5b), the solution obtained on fine mesh exhibits the highest accuracy. However, when the end of the simulation is being approached, that is, 85 % to 90 % of the prescribed total time (Fig. 5c), the result obtained using medium mesh is found as the most accurate. The temperature field inside the domain at t = 0.1 s exhibit distributions shown in Fig. 6. Near the interface a discrepancy from exact, linear, temperature profile is observed. In this case the increased mesh resolution leads to more accurate result, as quantified in Table 5.  Since the solution of Stefan problem with the approach presented in this paper differs from standard solutions obtained with single-fluid approach, it is noteworthy to outline some additional features of the model.
The consideration of liquid phase at saturation temperature imposes as the primary goal to study the vapor phase's temperature field. However, due to application of two-fluid model, there is also a temperature field associated with liquid phase that is distributed throughout the domain, as shown in Fig. 7. Fixing the liquid phase's temperature at the interface to saturation temperature by application of zero-resistance heat transfer model ensures stagnant temperature profile of the liquid, present also in the cell adjacent to interface zone with = 0 that is preceded by a jump in temperature decrease.
A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno The tendency to equalize the velocities of the phases by application of the appropriate drag model, mentioned before, is shown in Fig. 8. The close agreement of the maximum velocity magnitude of the vapor (square marked) and liquid phase's velocity magnitude is outlined in the interface zone.
The liquid phase's temperature distribution and the velocity magnitude field of both phases are general features of the model and, therefore, are not considered further in this study.

Stretched quadrilateral mesh
The incorporation of stretched mesh in near wall region of a quadrilateral mesh also results in accurate estimation of interface displacement, as shown in Fig. 9a. Although shortly after the simulation start only the result obtained on the fine mesh enters the prescribed error band, as displayed in Fig. 9b, finally, as shown in Fig. 9c, it is accompanied with the coarse mesh solution.
According to exact measures, at two selected time instances, as indicated in Table 6, the interface displacement solution obtained with fine stretched quadrilateral mesh is more accurate than the solution obtained with quadrilateral uniform mesh of same number of cells. The medium mesh that yields most accurate interface displacement in the selected time instances in the case of uniform cell distribution, however, in this case has highest discrepancy, but still in reasonable limit, that is, lower than 5 %.   The temperature field obtained using stretched quadrilateral meshes, shown in Fig. 10, has the maximum relative error in computation of temperature field lower than 10 %, when all the applied mesh resolutions are considered (Table 7). In the case of medium mesh, however, the better results are obtained in the case without stretched cells in wall adjacent zone.

Hybrid mesh
The interface displacement in the case of hybrid mesh is calculated using Eq. (20) and exhibits the distribution shown in Fig. 11. Although detailed view in Fig. 11b suggests the solution obtained using the fine mesh as the most promising one, the result in Fig. 11c  with the exact values of the interface displacement error, given in Table 8, indicate the reliability of medium size mesh.  Sampling the temperature field data from cell's centers yields a "zig-zag" temperature field distribution inside the zone composed of triangular cells in the hybrid mesh, as shown in Fig. 12. The error in computation of temperature field for coarse and medium resolution mesh is in the range of the errors already reported for other types of meshes with these resolutions, that is, less than 10 %. However, the maximum relative temperature field error obtained using fine mesh at the simulation end time is, compared to fine mesh calculations with other mesh types, the largest ( Table 9).
The distribution of the volume fraction of the vapor phase and the temperature for the hybrid mesh at t = 0.1 s are shown in Fig. 13. The volume fraction sharply changed from zero to one, and the temperature field features linear change from the wall temperature to the saturation temperature in the vapor phase.

Conclusions
The reduction of ship resistance via air injection beneath the ship has been extensively studied in the past, but the corresponding idea of ship resistance reduction using boiling flow is relatively new in the literature, and the numerical assessments in this regard are still missing. To overcome this, a novel method for computation of mass transfer is proposed herein. In this paper, a novel method for computation of Stefan problem was proposed and the following conclusions can be drawn: 1. Stefan problem was accurately solved using the two-fluid VOF model. To this end, the local-parameter based, heat and mass closure was derived, under the assumption of the constant liquid temperature at the saturation. 2. The sharp distribution of the volume fraction of liquid/vapor was achieved by using the geometrical reconstruction scheme implemented in the VOF model. 3. The velocity jump was captured at the interface owing to the sharp distribution of the interfacial area density calculated as the volume fraction gradient. 4. The method is capable to solve Stefan problem on the hybrid mesh predominantly composed of triangular cells with a commensurate level of accuracy. The accurate solutions are also obtained for the stretched quadrilateral mesh. Although the proposed method was validated only for Stefan problem, the application of the method to engineering problems is considered to be straightforward because the method only uses a local value: the cell volume. Hence, the problem of ship resistance reduction using boiling flow can now be studied in the field of naval hydrodynamics CFD modeling.

ACKNOWLEDGEMENT
The help received from Mr. Nithin Moran Narayan (TU Bremen) and Dr. Amine Ben Hadj Ali (ANSYS, via Ansys Student Community forum) is gratefully acknowledged. The valuable suggestions from Dr. Mijail Febres Soria and Mr. Lubomir Bureš (both Paul Scherrer Institut) are also sincerely appreciated

APPENDIX: The derivation of correlation for calculation of local Nusselt number
The expression for local, cell-size dependent, Nusselt number stems from the balance between thermal phase change and energy jump mass transfer models. It is thereby assumed that liquid is at saturation temperature, = ′. Other, model related, assumptions are also incorporated.
Thus, in thermal phase change model it is also assumed that vapor and liquid phase's scaling factors are equal to unity, that is, = = 1.0 and that interfacial area is computed as magnitude of volume fraction field, Eq. (13). Hence, Eq. (19), reduces to: where hfg denotes specific heat of vaporization.
As stated in the introduction, one of the goals of the proposed model is to ensure the occurrence of dominant mass transfer rate in the interface zone. Therefore, in order to define energy jump model, consider the phase change process in a three-cell stencil, as shown in Fig. 14. Let the interface is thereby located between positions i and i + 1 at x-coordinate of right sided Cartesian coordinate system. Furthermore, it is assumed that the superheated vapor phase occupies the zone on the left side of the interface, while the remainder of the domain is filled with liquid at saturation temperature.
The temperature gradient at the right side of the central cell's center is calculated as [30]: where ℎ denotes the distance between the interface and the center of adjacent vapor cell, see Fig. 14.

Fig. 14
The mass transfer process in a three-cell stencil.
The mass transfer rate calculated using energy jump model is defined as [30]: where indexes and denote, respectively, liquid and vapor phase, ∇ is the temperature gradient, ⃗⃗ is the interface normal directed from vapor to liquid phase, is the interface area and is the cell volume.
Since the liquid is at saturation temperature, ∇ = 0, and Eq. (A 3) is reduced to: The consideration of the problem as one-dimensional allows for the following simplification:

|∇ | ≈
Hence, the energy jump model reads: The heat transfer is present only in vapor phase and, therefore, the heat transfer coefficient is defined as: Further approximations refer to i-th cell temperature and distance h, and are given by: Thus, the heat transfer coefficient is a fictitious quantity since it is used to express the approximate distance from the cell center of the control volume that contains the interface to A solution to Stefan problem Alen Cukrov, Yohei Sato, using Eulerian two-fluid model Ivanka Boras, Bojan Ničeno the interface itself and to involve the thermal conductivity of the vapor phase in the mass transfer model. Furthermore, its inclusion in the mass transfer model yields an approximation of the temperature gradient at the interface, that appears from energy balance at the interface, i.e., in the heat flux in the normal direction to the interface.