## Appendix A  
Let us consider the planar vapor-liquid interface in a phase change process, as shown in [PERSON] et al. (2013). In addition, let us consider the most general industrial case in which the flow is turbulent, and the molecular thermal conductivity would be appended with the turbulent thermal conductivity in energy balance at the vapor-liquid interface, as demonstrated in [PERSON] and [PERSON] (2015). Assuming that the ratio between the interfacial area and the cell volume is equal to the magnitude of a gradient of a volume fraction, the energy balance at the interface would read\[\dot{m}_{pq}r_{0}=\Big{\{}(\lambda_{l}+\lambda_{t})(\
abla T_{l})\cdot\overset{ \rightarrow}{n}+\Big{[}-(\lambda_{v}+\lambda_{t})(\
abla T_{v})\cdot\overset{ \rightarrow}{n}\Big{]}\Big{\}}|\
abla\alpha_{v}| \tag{1}\]  
By further assumption of laminar flow and the saturated liquid, the above equation, when divided by the specific heat of vaporization, reduces to  
\[\dot{m}_{pq}=\frac{\Big{[}-\lambda_{v}(\
abla T_{v})\cdot\overset{ \rightarrow}{n}\Big{]}|\
abla\alpha_{v}|}{r_{0}} \tag{2}\]  
On the other hand, in the thermal phase change model, the convective heat transfer at the interface is introduced, as shown in [27], and the equation for the mass transfer at the interface in the case of saturated liquid (\(T^{\prime}-T_{\text{Iiq}}\) = 0) would read  
\[\dot{m}_{pq}=-\frac{C_{v}h_{v}|\
abla\alpha_{v}|(T_{sat}-T_{v})}{r_{0}} \tag{3}\]  
where the calibration factor \(C_{v}\) is equal unity because we want to keep the basic principles in the mass transfer modeling, thereby avoiding any fine-tuning of model parameters.  
By the equalization of Equations (2) and (3), and by further assuming that the distance between the interface and neighboring cell's center is half the width of the computational cell, that is, \(dx=\Delta x/2\), and that the temperature gradient projection may be approximated as \(\
abla T_{v}\cdot\overset{\rightarrow}{n}=(T_{sat}-T_{v})/dx=2(T_{sat}-T_{v} )/\Delta x\), we have the interphase heat transfer coefficient that reads  
\[h_{v}=\frac{2\lambda_{v}}{\Delta x} \tag{4}\]  
Bearing in mind that the interfacial heat transfer coefficient is defined as  
\[h_{v}=\frac{\lambda_{v}Nu_{v}}{d_{v}} \tag{5}\]  
We arrive at the expression for the Nusselt number, which has only a fictitious function since it is used to actually approximate the temperature gradient in the interface region, as follows:  
\[Nu_{v}=\frac{2d_{v}}{\Delta x} \tag{6}\]  
The vapor phase diameter, dv, that is present in Equation (6) is the value that is manually inserted into the software. But, due to the heat transfer coefficient definition, e.g., Equation (5), this quantity vanishes, and the resulting heat transfer coefficient is defined with Equation (4).  
The model is implemented within ANSYS Fluent as a UDF in the context of interphase models via the Define-Exchange-Property macro. Equation (6) is thus subjected to the procedure for the calculation of the heat transfer coefficient, and this value has been returned by the Define-Exchange-Property macro.  
Furthermore, it is important to note that, prior to computation, the initial flow fields need to be initialized using the Define-Init macro on a per-phase basis. Hence, due to the interpenetrating continua concept that underlies the two-fluid modeling, one should access all the fields and, in each of them, prescribe the initial values of volume fraction, velocity components in x- and y-directions, fluid temperature, and the turbulent kinetic energy, since the flow is turbulent and the frozen turbulence approach has been used in the study.