Temperature changes around interface cells in a one-dimensional Stefan condensation problem using four well-known phase-change models

International Journal of Thermal Sciences xxx (xxxx) xxx

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts

Temperature changes around interface cells in a one-dimensional Stefan condensation problem using four well-known phase-change models Jong Hyeon Son, Il Seouk Park * School of Mechanical Engineering, Kyungpook National University, 80 Daehakro, Bukgu, Daegu, 41566, Republic of Korea

A R T I C L E I N F O

A B S T R A C T

Keywords: Phase-change heat transfer Phase-change model Stefan problem Volume of fluid (VOF)

Discontinuous changes in thermodynamic properties across a very thin-phase interface and irregular de­ formations in the interfacial shape are accompanied by phase changes. Therefore, to simulate the phase-change phenomenon with a finite-sized mesh system is very challenging. Various numerical phase-change models have been developed, and several of them have been embedded in commercial computational fluid dynamics codes. However, the temperature reproducibility has not been dealt with carefully. In this study, we focused on the fact that most of the numerical phase-change models treat phase changes from a volumetric perspective, even though the phase change is definitely an interfacial phenomenon. We solved the one-dimensional (1D) Stefan conden­ sation problem using four well-known numerical phase-change models (two temperature-difference phasechange models and two heat-flux phase-change models). Since the 1D Stefan problem has no convective effects, it is appropriate for investigating the inherent features of models exclusively. The temperature changes and interface movement with time were compared according to the applied phase-change model. Non-physical stepwise temperature changes over time were observed in models based on a volumetric perspective. The error analysis for different grid and time-step sizes was presented for the four different phase-change models.

1. Introduction Phase-change heat transfer is a physical phenomenon that occurs in a variety of industrial applications, including power generation, refriger­ ation, and seawater desalination and in electronic systems. During a phase-change process, enormous amounts of thermal energy (latent heat) are exchanged, while the pressure and temperature remain con­ stant. At the same time, the desired heat transfer function can be per­ formed with a relatively smaller device or a very small mass flow rate. It is important to note that the phase-change process does not require extreme temperatures or pressures to enhance the heat transfer rate, and therefore, it is a highly reliable process with a low possibility of leakage. In general, phase-change heat transfers from gas to liquid or vice versa are characterized by complex convective flows caused by the interaction between buoyancy, surface tension, and viscous forces. Although the phase-change phenomenon seems irregular and chaotic, many experimental researches have shown that the phenomena are sensitive to operating conditions but fairly predictable. Therefore, the focus of most previous studies [1–8] was to provide the resulting heat-flux data for various operating conditions. In recent years,

experimental devices and measurement techniques have significantly advanced, and devices such as high-speed cameras and technologies such as infrared thermometry and transparent electric heating, are now used to investigate the phase-change phenomena in detail [8–13]. However, despite these advances in measurements, there is still need for improving the accuracy of numerical analysis methodologies to better understand the underlying mechanisms of the phase-change phenomena. In numerical analyses for the phase-change heat transfer based on finite volume or finite element methods, mass exchanges across the thinphase interface should be resolved using finite-sized mesh systems. Furthermore, the abrupt changes in thermodynamic properties across phase interfaces requires more careful consideration; for example, the specific volume of the gas phase in a phase-change process is several thousand times that of the liquid phase. Early numerical studies on phase-change heat transfers did not simultaneously solve for both the gas and liquid phases, but instead solved for a single phase of a gas or liquid (usually the unsaturated phase) by assuming a fixed phase interface and a constant interfacial temperature at saturation. Later, moving-grid techniques were employed in numerical approaches to reflect the time progression of the

* Corresponding author. E-mail address: [email protected] (I.S. Park). https://doi.org/10.1016/j.ijthermalsci.2020.106718 Received 11 June 2020; Received in revised form 6 October 2020; Accepted 28 October 2020 1290-0729/© 2020 Elsevier Masson SAS. All rights reserved.

Please cite this article as: Jong Hyeon https://doi.org/10.1016/j.ijthermalsci.2020.106718

Son,

Il

Seouk

Park,

International

Journal

of

Thermal

Sciences,

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Nomenclature C cp hLG k m˙ L n Q ST t Δt T V

ΔV Δx

Cell index specific heat at constant pressure latent heat thermal conductivity phase-change rate time-step index sum of the heat flows entering or exiting through the cell faces heat source time time-step size temperature velocity vector

cell volume size of grid cell

Greek symbols αL liquid volume fraction γ numerical coefficient δ interface position for Stefan problem ρ density σ surface tension coefficient χ solution of Eqn. (11) Subscripts G gas phase L liquid phase sat saturation

phase interface [14–16]. In recent years, numerical techniques for interface tracking, such as the volume of fluid (VOF) [17,18], marker and cell (MAC) [19–21], and level set methods [22,23], have enabled numerical simulations for very complex multi-phase flows in a compu­ tational domain involving both gas and liquid phases [24–27]. Meanwhile, the numerical phase-change model to determine the rate of phase-change is very important for the successful numerical simula­ tion of phase-change phenomena. Various models were developed for accurately predicting the phase-change rate based on the heat transfer information obtained from the solution procedures. Lee proposed a temperature-difference phase-change model in which the degree of subcooling or superheating of the interface cell was considered [28]. Many researchers employed this model for analyzing various phase-change problems because of its simplicity for calculating the phase-change rates. However, the numerical results were reported to be significantly affected by the numerical coefficient (γ) employed in this model; moreover, the value of γ ranges over a very wide band from 10− 3 to 1012 depending on the sizes of the grid and time step. The recently proposed Rattner–Garimella temperature-difference phase-change model [29] is similar to the Lee model in that it also predicts the phase-change rate based on the degree of subcooling or superheating of the interface cell; however, it significantly improves the numerical stability by adjusting the allowable time-step size and adopting a selective method for phase-change rate, considering the fluid amount remaining in an interface cell. Thus, the ambiguity about the numerical coefficient was cleared by removing the numerical coefficient from the model. Phase-change heat transfer is an interfacial phenomenon that occurs at a phase interface. However, the aforementioned temperaturedifference phase-change models predict the phase-change amount using the difference between the average temperature of the gas-liquid mixture in an interface cell and the saturation temperature. That is, these models tend to regard the phase-change heat transfer as a volu­ metric phenomenon. In addition, these models commonly apply the heat sources that represent the latent heat to the energy conservation

equation instead of forcing the interfacial temperature to remain at the level of the saturation temperature. Meanwhile, the sharp-interface model computes the rate of phase change using the temperature gradient on the phase interface [30,31]. In this model, the phase-change rate is determined by the heat fluxes at the phase interface toward the saturated and unsaturated phases. For the model, the temperature gradient should be precisely computed at the phase interface and for the interfacial position to be accurately traced first. In contrast to the previously described temperature-difference phase-change models, the sharp-interface model can be classified as a heat-flux phase-change model, and it strictly follows the interfacial-phenomenon characteristics. Sun et al. [32] proposed a new heat-flux phase-change model that computes the phase-change rate using the finite-volume-based heat balance in an interface cell rather than using the heat fluxes at the exact phase interface. Therefore, the Sun model can be regarded as the finite-volume-based sharp-interface model. Note that though both the sharp-interface and Sun models predict the phase-change rate using the heat balance around the phase interface; however, the Sun model de­ scribes the phase interface using a control volume with a finite thick­ ness; thus, strictly speaking, it is not an interfacial-perspective model. In this study, we compare the time variations in the temperature around the phase interface obtained using the four aforementioned phase-change models. To evaluate possible distortions in the tempera­ ture and interface movement by the applied numerical phase-change models in a finite volume method solution procedure, a onedimensional (1D) Stefan condensation problem without any convec­ tive effect is solved in this study. We have confirmed that some models provide a non-physical stepwise temperature change with time because of the differences in perspectives with regard to the phase-change phe­ nomena. Incorrect temperatures may affect the whole phase-change process and amplify errors over time; although this effect is expected to be more critical in the phase-change heat transfers accompanying a fluid flow, it was not recognized well since the effect was simply hidden by strong convective effects. The present study is the very first study, to the best of our knowledge, that delves into the time progression characteristics of temperature when numerical phase-change models are applied. Although numerical phase-change models are widely used, stepwise temperature variations have not been considered. The contribution of this study will be helpful in the appropriate selection of phase-change models and evaluation of their numerical results. The rest of the paper is structured as follows. The Stefan problem is presented in Section 2. The advantages and disad­ vantages of the aforementioned four phase-change models, as well as the temperature changes and interface movements for different phasechange models, are presented in Section 3. Section 4 presents the con­ clusions of this study.

Fig. 1. Schematic of the Stefan condensation problem. 2

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

2. Stefan problem

Table 1 Thermal properties of isobutane.

In this study, a 1D Stefan condensation problem is simulated using four numerical phase-change models. The Stefan problem does not include convective effects; the movements in the phase interface are only caused by conductive heat transfer, and large discontinuities in the thermodynamic properties appear across the phase interface. A sche­ matic of the Stefan condensation problem is shown in Fig. 1. Initially, the 0.2-mm-long computational domain is filled with saturated gas, and the right and left boundaries are maintained at the saturation temper­ ature and a constant subcooled temperature, respectively. Here, the saturation temperature is 300 K, and the subcooled temperature is 295 K. The condensation of gas starts from the subcooled left boundary, and the phase interface grows over time. While the temperature of the gas on the right side of the phase interface is maintained at the saturation temperature, that of the liquid on the left side of the phase interface decreases almost linearly from the saturation temperature at the phase interface to the subcooled temperature at the left boundary. The phase change at the interface continues to progress because of the heat sink at the left boundary; thus, the interfacial position moves steadily toward the right. The 1D Stefan condensation problem can be described by the following conservation equations for volume fraction and energy:

∂αL m˙ L = , ∂ t ρL

Properly 3

Density [kg/m ] Thermal conductivity [W/m⋅K] Specific heat [J/kg⋅K] Latent heat [J/kg]

αn+1 − αnL L Δt

ρcp Δx ( Δt

(

(2)

where αL is the liquid volume fraction; t, time; m˙ L , the phase-change rate (the mass generation rate of the created phase (i.e., the liquid phase) per unit volume); ρ, the density; cp, the specific heat at constant pressure; k, the thermal conductivity; T, the temperature; ST, the heat source per unit volume; and the subscript L, the liquid phase. In the present study, the VOF method was used to numerically pro­ duce movement of the phase interface. The liquid volume fraction αL in Eqn. (1) is defined as the ratio of the volume of liquid in a cell to the cell volume. The right-hand side of Eqn. (1) represents the source of liquid volume fraction by phase changes. The source terms m˙ L / ρL and ST in Eqns. (1) and (2) are obtained differently depending on the applied phase-change model (the details will be described in the following section). In this study, the interface cell is defined as a cell of volume fraction between 0 and 1, and the phase interface is the position where the volume fraction is equal to 0.5. The position of the phase interface can be computed via linear interpolation among the volume fractions of neighboring cells. The density and thermal conductivity of the interface cell where gas and liquid coexist can be computed using the volumefraction-weighted average: (3)

k = kG + (kL − kG )αL ,

(4)

αG ρG cp,G + αL ρL cp,L , αG ρG + αL ρL

550.6 0.0892 2446

9.12 0.0169 1816

329,400

(6)

,

) T n+1 − T n = Dn+1 − Dn+1 + ST Δx, e w

T(x, t) = Twall +

(7)

( ) ∂T Dw = k , ∂x w

(8)

( ) ( √̅̅̅̅̅̅̅̅̅̅̅̅̅ ) Tsat − Twall x kL erf , 2 ρL cp,L t erf(χ )

(10)

where erf is the error function, and χ is a solution to the following transcendental equation: (

)

χ exp χ 2 erf(χ ) =

cp,L (Tsat − Twall ) √̅̅̅ , hLG π

(11)

where hLG is the latent heat. In the present study, isobutane was used as a working fluid, and the thermal properties of the liquid and gas phases of this fluid are listed in Table 1. Note that the thermal properties are assumed constant regardless of changes in temperature. 3. Numerical solutions obtained using different phase-change models 3.1. Lee model The Lee model employs the difference between the temperature of the interface cell and the saturation temperature to compute the phasechange mass generation rate and the resulting latent heat source per unit volume according to the following equations:

and the specific heat of the interface cell can be computed using the mass-weighted average: cp =

ρL

Gas

where Δt is the time-step size; Δx, the mesh size; and superscript n, the time-step index. The diffusive heat-flux terms in Eqns. (7) and (8) were discretized by the 2nd order central differencing scheme. The two con­ servation equations were coupled via the source terms. The properties that depend on a volume fraction were applied. Several iterations were required in every time step until temperature and volume fraction converge. The discretized algebraic equations were solved using the biconjugate gradient stabilized (BICGSTAB) method [33]. Since nonlinear convective terms are not included in the Stefan problem, the exact solutions for the interface position and temperature profile can be determined as follows [34,35]: √̅̅̅̅̅̅̅̅̅̅̅̅ kL t δ(t) = 2χ , (9) ρL cp,L

)

ρ = ρG + (ρL − ρG )αL ,

m˙ L Δx

( ) ∂T De = k , ∂x e

(1)

) ∂ ∂( ∂T ρc T = k + ST , ∂t p ∂x ∂x

=

Liquid

(5)

m˙ L = γ αG ρG

where subscript G in Eqns. (3)–(5) denotes the gas phase. Uniform mesh systems of sizes ranging from 1.25 μm to 10 μm were used in the simulations. The unsteady terms in governing equations were discretized using the first-order backward differencing scheme, and the governing equations were spatially integrated via the finite volume method:

Tsat − T Tsat

ST = m˙ L hLG ,

if ​ T < Tsat ,

(12) (13)

where γ is a numerical coefficient. As mentioned earlier the value of γ can be changed over a very wide range depending on the grid size and time-step size. In the Lee model, the latent heat corresponding to the mass 3

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 2. Temperature and volume fraction distributions around the interface cell for the Lee model. Fig. 4. Results for the solution dependence of the numerical coefficient in the Lee model for the time-step size of 10− 6 s and grid size of 2.5 μm.

stepwise change in temperature occurs over time. This is regarded as an inevitable physical defect of the Lee model because the subcooledness ( = Tsat − T) should be inversely proportional to the gas volume fraction (αG = 1 − αL ) in order to secure a continuous and nearly constant phasechange rate in a cell according to Eqn. (12). Further, as shown in the graph in the inset of Fig. 3, the liquid volume fraction in the 9th cell starts increasing from zero before the 8th cell is fully filled with the liquid. This is because the temperature of the neighboring saturated-phase cell (the 9th cell) falls below the saturation temperature level by the cooling effect even before the interface cell (the 8th cell) is completely filled with the liquid; in other words, two consecutive cells can be flagged as an interface cell. In the Lee model, these defects cannot be controlled. Fig. 4 shows the changes in the interface position over time when applying different values of the numerical coefficient in the Lee model. In the present 1D Stefan problem, the thickness of the liquid film can be obtained as follows: ∑ δ(t) = αL,i Δxi , (14)

Fig. 3. Variations of the cell temperature and volume fraction over time in the 8th to 10th cells for the Lee model.

i

generation of the created phase is used as a source term in the energy equation; thus, the cell temperature of the interface cell is maintained around the saturation temperature: the energy balance in an interface cell can be completely explained with the heat flows passing through the cell faces and the latent heat, and the phase-change rate is determined by only the degree of subcooling (or superheating) of the interface cells in which both phases coexist. The Lee model presumes the phase-change process, which was originally a surface phenomenon, to be a volumetric phenomenon in a finite-sized control volume containing the phase interface. Fig. 2 shows a schematic of the temperature and volume fraction distributions around the interface cell when the Lee model is applied to the Stefan condensation problem. In the Lee model, condensation occurs only when the gas volume fraction is non-zero and the cell temperature is lower than the saturation temperature. This is captured in Eqn. (12). Therefore, if the temperature spreads via numerical diffusion, phase change can simultaneously occur in multiple cells that do not include the actual phase interface, as shown in the figure. Fig. 3 depicts the variations in the cell temperature and volume fraction over time in the 8th to 10th cells from the left wall. As shown in the figure, the cells sequentially turn to an interface cell; the interface cell shows a nonlinear decreasing pattern of temperature while the liquid volume fraction shows a nearly linear increase. However, sur­ prisingly, this nonlinear temperature decrease in the cell continues even after the phase interface completely moves to the next cell; thus, a

where Δx is the size of the grid cell. For a constant time-step size of 10− 6 s and grid size of 2.5 μm, when the numerical coefficient is 103, the interface growth rate differs significantly from the exact solution. As the numerical coefficient gradually increases to 108, the film thickness converges to the exact solution. However, if a numerical coefficient of 1012 is chosen, a totally false result is obtained. This characteristic of the Lee model has been previously reported [27,36–38]. The errors depicted in the graph legend in Fig. 4 are the relative errors to the exact solution at time 0.2 s. For a constant time-step size of 10− 6 s, the results of the grid de­ pendency test for numerical coefficients of 104 and 106 are depicted in Fig. 5 (a) and (b). When the numerical coefficient is 104, for all tested grid sizes, the interface position shows somewhat large differences from the exact solution; in particular, the error slightly increases as the grid size decreases since the numerical coefficient of 104 is inappropriate for the present time-step size of 10− 6 s. On the other hand, for a numerical coefficient of 106, the predicted interface position is in good agreement with the exact solution for all tested grid systems, and the error reduces as the grid size decreases. Note that the film thickness in a single grid cell increases linearly in all cases, as shown in the magnified view of Fig. 5 (b). This implies that the Lee model needs much smaller grids and correspondingly small time-step sizes to precisely reproduce nonlinear interfacial growth in phase-change problems.

4

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 6. Temperature and volume fraction around the interface cell for the Rattner–Garimella model.

Fig. 5. Results of grid dependency test for the numerical coefficient of 104 and 106 with the time-step size of 10− 6 s in the Lee model.

Fig. 7. Variations of the cell values in the 8th to 10th cells for the Ratt­ ner–Garimella model: (a) volume fraction and cell temperature, and (b) heat capacity.

5

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

3.2. Rattner–Garimella model In the temperature-difference phase-change model proposed by Rattner and Garimella [29], the mass generation rate of the created phase can be computed as follows: ( ) m˙ L = − max m˙ 1 , m˙ 2 , m˙ 3 , (15)

m˙ 1 =

ρcp (T − Tsat )

m˙ 2 = −

m˙ 3 = −

hLG Δt

αG ρG Δt

,

,

)− 1 ( 1 1 1 − , Δt ρG ρL

ST = m˙ L hLG .

(16) (17) (18) (19) Fig. 8. Changes in the interfacial position of Rattner–Garimella model for different time-step sizes with a grid size of 2.5 μm.

Eqn. (16) can be used to compute the phase-change rate based on the difference between the cell temperature and saturation temperature; in this regard, this model is similar to the Lee model. However, unlikely the Lee model, the time-step size influences the phase-change rate. The density and specific heat of the fluid in Eqn. (16) are those of the gasliquid mixture in the interface cell. Eqn. (17) gives the phase-change rate when all of the remaining saturated phase in an interface cell changes, thus ensuring that the generated mass does not exceed the mass of the saturated phase remaining in the cell. Eqn. (18) limits the volume change of the created phase for a single time step to be less than the cell volume. Similar to the Lee model, the Rattner–Garimella model con­ siders the phase-change process as a volumetric phenomenon, and the latent heat corresponding to the mass generation rate of the created phase is considered as a heat source in the energy conservation equation. Fig. 6 depicts a schematic for the temperature and volume fraction around the interface cell when the Rattner–Garimella model is applied. Similar to the Lee model, the temperature of the saturated-phase cell adjacent to the interface cell may be less than the saturation temperature because of the normal heat transfer with the interface cell and the nu­ merical diffusion. Hence, the cell temperature in multiple consecutive cells (the i, i+1, and i+2th cells in the figure) around the cell (the ith cell) containing the actual phase interface may be lower than the saturation temperature, as shown in Fig. 6. However, this model forces only two consecutive cells to be flagged as an interface cell, as shown in the figure. Therefore, the subcooled temperature in the i+2th cell did not activate phase change in this model. The changes in the volume fraction and cell temperature in the 8th to 10th cells are depicted in Fig. 7 (a). Similar to the Lee model, the cells sequentially change to an interface cell in the Rattner–Garimella model. The increasing rate of the liquid volume fraction is nearly constant in an interface cell, but decreases slightly as the interface cell moves to the next cell. Whenever the phase interface moves to the next cell, the cell temperatures decreases abruptly and then rapidly recovers. This abnormal temperature fluctuation was observed throughout the whole computing domain regardless of the fluid phase: Even when the 8th to 10th cells were in the saturated phase (i.e., earlier than 0.06 s or when the phase interface was to the left of the 8th cell), the cell temperature of these cells showed such a fluctuation. Further, even when the cells completely changed to an unsaturated phase, the cell temperatures showed a fluctuating pattern whenever the phase interface moved to the next cell, but the fluctuation was weaker than that observed when the cell was in a saturated phase. This is clearly a non-physical result; this result is obtained because the Rattner–Garimella model uses the heat capacity (ρcp ) for the mixture in an interface cell to predict the phase-change rate, as shown in Eqn. (16). Fig. 7 (b) shows the time variation in the heat capacity of the 8th to 10th cells; the heat capacity of each cell abruptly increases from that for

Fig. 9. Results of the grid dependency test for the Rattner–Garimella model with the time-step size of 10− 6 s.

the gas phase as soon as the cell turns to an interface cell. Despite such an increase in heat capacity, to realize a continuous change in the phasechange rate, a non-physical temperature fluctuation such as that shown in Fig. 7 (a) is inevitable when using Eqn. (16). We think that this fluctuation arises mainly because this model retains the volumetricphenomenon perspective; the representative temperature for the mixture in a control volume is not sufficient to guarantee a continuous change in phase-change rate and physically reasonable temperature changes simultaneously. Furthermore, the non-physical temperature fluctuation in the saturated-phase cells (the 8th to 10th cells before 0.06 s, as shown in Fig. 7 (a)), may trigger phase changes at the cells according to Eqn. (16); therefore, as mentioned earlier, this model forces the phase change to occur only in two consecutive cells to prevent non-physical phase changes caused by the numerical fluctuation in temperature; however, this is a clearly artificial effect. Nevertheless, the Rattner–Garimella model is an innovative advancement from the Lee model: the Rattner–Garimella model did not employ the numerical coefficient (γ), which has been controversial among researchers for a long time, and instead used the time-step size to predict the phase-change rate. Therefore, to additionally ensure nu­ merical stability, the selection method for phase-change rate as per Eqns. (15)–(18) was used, and the following limitation for the time-step 6

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 11. Variations of the volume fraction and cell temperature in the 8th to 10th cells for the Sun model.

Fig. 10. Temperature and volume fraction distributions for the Sun model.

size was proposed: [ /( ) ] 1 k Δt ≤ min Δx2 . 4 ρcp eff

(20)

i

Fig. 8 shows the changes in the interfacial position for three different time-step sizes and a constant grid size of 2.5 μm. When a time-step size of 10− 4 s is applied, the result shows a large difference from the exact solutions. However, the difference decreases as the time-step size is reduced. The maximum of the recommended time-step size as per Eqn. (20) is about 1.531 × 10− 6 s. When the time-step sizes are smaller than the recommended maximum, the results are in good agreement with the exact solution, with an error of less than 1%. Fig. 9 shows the results of the grid-size dependency test for the Rattner–Garimella model with a constant time-step size of 10− 6 s. Even with a somewhat coarse grid system of 5 or 10 μm, the changes in the interfacial position matched well with the exact solution, and the error was reduced as the grid size decreased. However, a linear growth was observed in the interfacial position in a computing cell similar to that observed in the Lee model. The recommended maximum for the timestep size for a grid size of 1.25 μm was less than 10− 6 s; therefore, the error was slightly higher than those for the coarser grids.

Fig. 12. Results of the grid dependency test for the Sun model with time-step size of 10− 6 s.

(21), and thus, the heat flow between the both cells is simply offset. In the interface cell, the latent heat corresponding to the phase change is applied as a source term in the energy conservation equation. Meanwhile, the cell temperature in the neighboring saturated-phase cell is forcibly maintained at the saturation temperature using the following source term linearization technique: { / at ​ Cint QCint ΔVCint ST = . (22) 30 10 Tsat − 1030 TCsat at ​ Csat

3.3. Sun model In the Sun model [32], the phase-change amount is calculated using the heat flows through the cell faces of the interface cell, which is defined as a cell with a volume fraction between 0 and 1, as shown in Fig. 10. The sum of the heat flows for the interface cell and its neigh­ boring saturated-phase cell are used to calculate the phase-change rate as follows: QCint QCsat m˙ L = + , hLG ΔVCint hLG ΔVCsat

Thus, the Sun model resolves the non-physical problems that the temperature of a saturated-phase cell does not remain at the saturation temperature and that multiple interface cells are flagged. In addition, since this model estimates the phase-change based on the heat balance in an interface cell, it can be regarded as a heat-flux-based phase-change model. A schematic of the temperature and volume frac­ tion distributions near the interface cell in the Sun model is shown in Fig. 10. If the energy source of Eqn. (22) is applied in an interface cell, the right-hand-side terms in Eqn. (7) for energy conservation are perfectly self-balanced, i.e., the temperature of the interface cell in the left-hand side of Eqn. (7) does not change and remains constant at the level of the saturation temperature if there are no numerical diffusions during the phase change. Recall that the interfacial temperature remains

(21)

where Q is the sum of the heat flows entering or exiting through the cell faces; ΔV is the cell volume; and the subscripts Cint and Csat denote the interface cell and its neighboring saturated-phase cell, respectively. The Sun model presumes the phase-change rate to be determined by thermal communication with only the unsaturated phase; as the phase change progresses in a cell, the cell temperature decreases by numerical diffu­ sion, and the heat exchange with the neighboring saturated-phase cell occurs naturally; in order to exclude the effect of heat flow from or to the saturated phase, the Sun model uses the heat balances for both the interface cell and its neighboring saturated-phase cell as shown in Eqn. 7

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 14. Variations of the volume fraction and cell temperature in the 8th to 10th cells for the sharp-interface model.

Fig. 13. Temperature and volume fraction distributions near the interface cell for the sharp-interface model.

at the saturation temperature; that is, the situation is similar to the occurrence of a phase interface of finite thickness of the cell size. Therefore, the stepwise temperature changes are expected to be inevi­ table in this model too. Owing to these characteristics, the Sun model seems like a volumetric model even if this model employs the heat-flux information to compute the phase-change rate. Furthermore, the discontinuous changes in temperature in all volumetric models are ex­ pected to require much higher grid density for problems with highly complicated changes in the interfacial shape. Fig. 11 shows the changes in the volume fraction and cell tempera­ ture in the 8th to 10th cells. As mentioned earlier, in the Sun model, the cell temperature of the interface cell continues to remain at the satu­ ration temperature throughout the phase change. Therefore, the tem­ perature in the unsaturated-phase area also remains almost unchanged unless the phase interface moves to the next cell. Then, the rate of phase change is constant when the phase interface stays in one cell because the phase-change rate is determined by the heat balance in the interface cell as shown in Eqn. (21). However, in a very short moment when the interface cell changes to the next cell, there are abrupt changes in temperature, producing a stepwise variation in temperature, as shown in Fig. 11. Fig. 12 shows the results of the grid dependency test for the Sun model with a constant time-step size of 10− 6 s. The Sun model shows a linear growth rate of the interfacial position in a computing cell, similar to other volumetric models, as shown in the figure. However, even with coarse grid systems with a grid size of 5 or 10 μm, the interfacial position matched well with the exact solution, and the error reduces gradually as the grid size decreases.

Fig. 15. Changes in the interfacial position of sharp-interface model for various grid sizes with the time-step size of 10− 6 s.

temperature at the phase interface is set to remain at the saturation temperature perfectly. Thus, the phase-change rate can be calculated using only the heat-flux at the interface toward the unsaturated phase by using Eqns. (23) and (24): ΔT , h−

(23)

q˙′′int Af , hLG ΔV

(24)

q˙′′int = k

3.4. Sharp-interface model

m˙ i =

The sharp-interface model computes the phase-change rate using the heat-flux data of the exact phase interface. Therefore, this model can be naturally classified as a heat-flux phase-change model. It treats the phase change as an interfacial phenomenon perfectly. In this model, to obtain the heat fluxes at the phase interface, first, the interfacial position should be traced accurately. This tracing includes determining the location and orientation of the phase interface. Orientation is not an issue in the present 1D problem, and the phase interface could be accurately located using the cell volume fraction. Next, the temperatures in the cells adjacent to the interface cell should be calculated accurately; for this purpose, in this model, the cell center of the interface cell is relocated to the center of the phase interface as shown in Fig. 13, and the

where h− is the distance between the center of an adjacent unsaturatedphase cell and the phase interface (the relocated center of the interface cell), and ΔV is the cell volume of the interface cell as shown in Fig. 13. The sharp-interface model does not employ the latent heat source in the energy conservation equation (i.e., ST = 0 in Eqn. (7)); instead, it sets the temperature at the exact phase interface as the saturation temper­ ature. Accordingly, the temperatures in the saturated-phase cells are strictly maintained at the saturation temperature, and there is no nu­ merical diffusion in the temperature and volume fraction. However, even with these advantages, the interfacial position must still be accu­ rately known in advance, especially, for multidimensional problems. 8

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

variation. In particular, the Rattner–Garimella model showed nonphysical temperature fluctuation. Although this unrealistic tempera­ ture variation is treated so that it does not influence the local phase change in the model, the false temperature distribution inevitably ac­ cumulates over time and creates an error that is not negligible. We tested four well-known phase-change models in this study. Non-physical tem­ perature change is commonly observed in the models based on a volu­ metric perspective. Table 2 lists the relative errors in temperature compared with the exact solution over time for the four phase-change models. The time-step size was 10− 6 s, and the grid size was 2.5 μm. The data in the table show that sharp-interface model has a significantly lower relative error comparing with other models. Table 3 summarizes the relative errors in the film thickness for various grids and time-step sizes for the four phase-change models. For the Lee model, the errors obtained for two different numerical co­ efficients (104 and 106) are presented together; the numerical coefficient in this model should be determined carefully considering the grid size and time-step size; using the numerical coefficient of 106 yielded totally different results from the exact solution for the time-step size of 10− 4 s for the grid sizes of 10, 5, 2.5, and 1.25 μm, and the error rapidly increased as the grid size decreased. To analyze the errors in the different phase-change models, we plotted the errors on the graphs with horizontal axes in the grid size as shown in Fig. 17. First, the solutions of the Lee model are significantly influenced by the numerical coefficient as shown in Fig. 17 (a) and (b); however, for the lower numerical coefficient of 104, despite somewhat large errors, the solution showed no dependency on the grid size and time-step sizes as shown in Fig. 17 (a). For a higher numerical coefficient of 106, the errors reduce with decreasing grid size as shown in Fig. 17 (b), and very large errors are observed for the time-step size of 10− 4 s; these errors are caused by the incorrect choice of the numerical coefficient for the given

Fig. 16. Temperature changes at the point of 18.75 μm over time for different phase-change models.

The changes in the volume fraction and cell temperature in the 8th to 10 cells are shown in Fig. 14. In this model, the interface cell center continues to move with the phase interface progression; thus, the cell temperature remains at the saturation temperature until the cell be­ comes a completely unsaturated-phase cell, as shown in the figure. Further, the cell temperature of the saturated-phase cells well remains constant at the saturation temperature level. Once the phase interface moves to the next cell, the cell temperatures begin to decrease contin­ uously. The stepwise temperature fluctuations no longer appears in this model. The rate of increase in the volume fraction in an interface cell is almost constant and the change in the volume fraction in the saturatedphase cells does not appear at all. All of these results of no-stepwise temperature change and no multiple-interface-cells flagging can be attributed to that this model well preserves the interfacial-phenomenon perspective. Fig. 15 shows the growth in the interfacial position for different grid sizes when the sharp-interface model is applied. Even for the largest grid size of 10 μm, the error is very small (less than 0.15%). This is a big advance in comparison to the errors in the other phase-change models. Except when the phase interface resides in the 1st cell from the left boundary, the predicted interfacial position traces the exaction solution well, as shown in the figure. The linear change in the interfacial position in an interface cell, which is commonly observed in all volumetric models, completely disappears in the sharp-interface model. Thus, the sharp-interface model is more advantageous when tracking a complex change at the phase interface with coarser grids. th

Table 3 Results for the relative errors in the film thickness at t = 0.1 s by four different phase-change models. Grid Δx [mm] 10

5

2.5

3.5. Comparison of phase-change models 1.25

Fig. 16 shows the temperature changes over time at a point 18.75 μm from the left boundary; the results obtained using the aforementioned four phase-change models are presented in a graph. The Lee, Ratt­ ner–Garimella, and Sun models produced a stepwise temperature

Time step size [s]

Error (t = 0.1 s) Lee (γ = 104)

Lee (γ = 106)

RattnerGarimella

Sun

Sharpinterface

2.329% 2.465% 2.465% 2.465% 2.536% 2.593% 2.593% 2.593% 3.582% 3.584% 3.584% 3.584% 4.681% 4.691% 4.692% 4.692%

43.65% 1.626% 1.612% 1.612% 44.78% 0.503% 0.647% 0.647% 124.7% 0.104% 0.101% 0.101% 173.1% 0.114% 0.104% 0.109%

11.99% 2.819% 1.728% 1.617% 23.13% 3.038% 0.875% 0.654% 38.40% 4.727% 0.533% 0.104% 59.34% 90.89% 0.942% 0.123%

1.695% 1.612% 1.605% 1.605% 0.803% 0.637% 0.629% 0.629% 0.383% 0.072% 0.055% 0.055% 0.877% 0.050% 0.029% 0.030%

0.139% 0.207% 0.219% 0.220% 0.035% 0.093% 0.107% 0.107% 0.176% 0.022% 0.049% 0.050% 0.425% 0.025% 0.021% 0.025%

10–4 10–5 10–6 10–7 10–4 10–5 10–6 10–7 10–4 10–5 10–6 10–7 10–4 10–5 10–6 10–7

Table 2 Results for the temperature change at a point 18.75 μm by four different phase-change models for a time-step size of 10− Time [s] 0.05 0.10 0.15 0.20 0.25 0.30

6

s and grid size of 2.5 μm.

T - Twall [K] 6

Error 6

Exact

Lee (γ = 10 )

Rattner- Garimella

Sun

Sharp- interface

Lee (γ = 10 )

Rattner- Garimella

Sun

Sharp- interface

5.000 4.260 3.484 3.019 2.702 2.467

5.000 4.400 3.565 2.998 2.769 2.421

4.956 4.405 3.566 2.993 2.775 2.412

5.000 4.412 3.571 3.000 2.778 2.421

5.000 4.266 3.486 3.021 2.703 2.468 Avg.

0.000% 3.286% 2.325% 0.696% 2.480% 1.865% 1.775%

0.880% 3.404% 2.354% 0.861% 2.701% 2.230% 2.072%

0.000% 3.568% 2.497% 0.630% 2.813% 1.865% 1.896%

0.000% 0.110% 0.057% 0.066% 0.037% 0.040% 0.057%

9

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 17. Relative errors in the film thickness for the different phase-change models with the grid sizes: (a) Lee model, γ = 104, (b) Lee model, γ = 106, (c) Ratt­ ner–Garimella model, (d) Sun model, and (e) Sharp-interface model.

time-step and grid sizes. The Rattner–Garimella model in Fig. 17 (c) showed a strong de­ pendency on the time-step size as this model directly uses the time-step size in calculating the rate of phase-change as shown in Eqns. (16)–(18). For small time-step sizes of 10− 6 and 10− 7 s, the error decreases as the grid size is reduced; however, for larger time-step sizes, the error in­ creases as the grid size is reduced. The error level achieved in the best choice of the grid and time-step sizes is similar to that of the Lee model with a numerical coefficient of 106. The Sun model shown in Fig. 17 (d) exhibits weak dependency on the time-step size if the adequately small time-step size of less than 10− 4 s is adopted. Except for the time-step size of 10− 4 s, the errors decrease as the grid size is reduced. The Sun model shows a lower error level than the Lee and Rattner–Garimella models. The sharp-interface model shown in Fig. 17 (e) exhibits almost the same error tendency as the Sun model; however, the achieved error level

is smaller than that of the Sun model. In particular, for large grid sizes of 10 and 5 μm, a much smaller error level was achieved compared with that of the other models. Generally, the sharp-interface model showed the least error among the tested phase-change models. For the benefit of the readers, we replotted the errors in the graphs with the horizontal axes showing the time-step size in Fig. 18. 4. Conclusion Although the phase-change phenomena have been widely examined using numerical phase-change models, the merits and demerits of phasechange models have not been dealt with systematically. We examined the characteristics of four well-known phase-change models by solving the 1D Stefan condensation problem in this study. Many phase-change models treat phase changes with a volumetric perspective, though this is originally an interfacial phenomenon. We found that the phase10

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

Fig. 18. Relative errors in the film thickness for the different phase-change models with the time-step sizes: (a) Lee model, γ = 104, (b) Lee model, γ = 106, (c) Rattner–Garimella model, (d) Sun model, and (e) Sharp-interface model.

change models based on a volumetric perspective produce stepwise changes in temperature with time; this non-physical temperature vari­ ation has not been importantly dealt with because the time variation of the local temperature is not a crucial matter and has not been explicitly observed because of convective effects. Although each numerical phasechange model has its own countermeasure to treat such non-physical effects, it is inevitable that the false temperature spreads over the whole domain and that errors are accumulated. Meanwhile, the sharpinterface model that strictly takes an interfacial perspective produced a continuous change in temperature and showed lower error than the other three models in volumetric perspective; however, accurately tracking the interfacial position is more important in the sharp-interface model, and this requirement is expected to be much more crucial in multidimensional problems. An error analysis was conducted for each phase-change model for various grid and time-step sizes; we expect the results will be helpful for researchers to choose the sizes of grid and time step as well the appropriate phase-change model for their study.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgements This research was supported by the National Research Foundation of Korea (NRF) funded by the Korea Government (Ministry of Science and ICT - Grant No. 2017M2A8A4017283, and No. 2019R1A2C3003890).

11

J.H. Son and I.S. Park

International Journal of Thermal Sciences xxx (xxxx) xxx

References

[19] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182–2189. [20] B.J. Daly, A technique for including surface tension effects in hydrodynamics calculations, J. Comput. Phys. 4 (1969) 97–117. [21] R.K.C. Chan, R.L. Street, A computer study of finite-amplitude water waves, J. Comput. Phys. 6 (1970) 68–94. [22] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [23] Y.C. Chang, T.Y. Hou, B. Merriman, S. Osher, A level set formulation of eulerian interface capturing methods for incompressible fluid flows, J. Comput. Phys. 124 (1996) 449–464. [24] C.H. Sohn, J.H. Son, I.S. Park, Numerical analysis of vortex core phenomenon during draining from cylinder tank for various Initial swirling speeds and various tank and drain port sizes, J. Hydrodyn. B 25 (2013) 183–195. [25] J.H. Son, C.H. Sohn, I.S. Park, Numerical study of 3-D air core phenomenon during liquid draining, J. Mech. Sci. Technol. 29 (2015) 4247–4257. [26] M. Magnini, J.R. Thome, Computational study of saturated flow boiling within a microchannel in the slug flow regime, J. Heat Tran. 138 (2016), 021502. [27] D.G. Kim, C.H. Jeon, I.S. Park, Comparison of numerical phase-change models of Stefan vaporizing problem, Int. Commun. Heat Mass Tran. 87 (2017) 228–236. [28] W.H. Lee, A pressure iteration scheme for two-phase flow modeling, in: T. N. Veziroglu (Ed.), Multiphase Transport Fundamentals, Reactor Safety, Applications, vol. 1, Hemisphere Publishing, Washington DC, 1980. [29] A.S. Rattner, S. Garimella, Simple mechanistically consistent formulation for volume-of-fluid based computations of condensing flows, J. Heat Tran. 136 (2014), 071501. [30] Y. Sato, B. Niˇceno, A sharp-interface phase change model for a mass-conservative interface tracking method, J. Comput. Phys. 249 (2013) 127–161. [31] I. Perez-Raya, S.G. Kandlikar, Numerical modeling of interfacial heat and mass transport phenomena during a phase change using ANSYS-Fluent, Numer. Heat Tran. B 70 (2016) 332–339. [32] D. Sun, J. Xu, Q. Chen, Modeling of the evaporation and condensation phasechange problems with FLUENT, Numer. Heat Tran. B 66 (2014) 326–342. [33] H.A. Van der Vorst, Bi-CGSTAB, A fast and smoothly converging variant to Bi-CG for the solution of Nonsymmetric systems, SIAM J. Sci. Stat. Comput. 13 (1992) 631–644. [34] V. Alexiades, A.D. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere Pub. Corp., Washington D.C., 1993. [35] S.W.J. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, J. Comput. Phys. 160 (2000) 662–682. [36] Z. Yang, X.F. Peng, P. Ye, Numerical and experimental investigation of two phase flow during boiling in a coiled tube, Int. J. Heat Mass Tran. 51 (2008) 1003–1016. [37] J. Wei, L. Pan, D. Chen, H. Zhang, J. Xu, Y. Huang, Numerical simulation of bubble behaviors in subcooled flow boiling under swing motion, Nucl. Eng. Des. 241 (2011) 2898–2908. [38] M. Bahreini, A. Ramiar, A.A. Ranjbar, Numerical simulation of bubble behavior in subcooled flow boiling under velocity and temperature gradient, Nucl. Eng. Des. 293 (2015) 238–248.

[1] M.K. Dobson, J.C. Chato, Condensation in smooth horizontal tubes, J. Heat Tran. 120 (1998) 193–213. [2] M. Zhang, R.L. Webb, Correlation of two-phase friction for refrigerants in smalldiameter tubes, Exp. Therm. Fluid Sci. 25 (2001) 131–139. [3] P.S. Lee, S.V. Garimella, Saturated flow boiling heat transfer and pressure drop in silicon microchannel arrays, Int. J. Heat Mass Tran. 51 (2008) 789–806. [4] W. Li, Z. Wu, A general correlation for evaporative heat transfer in micro/minichannels, Int. J. Heat Mass Tran. 53 (2010) 1778–1787. [5] K. Balasubramanian, P.S. Lee, L.W. Jin, S.K. Chou, C.J. Teo, S. Gao, Experimental investigations of flow boiling heat transfer and pressure drop in straight and expanding microchannels – a comparative study, Int. J. Therm. Sci. 50 (2011) 2413–2421. [6] R. Ali, B. Palm, C. Martin-Callizo, M.H. Maqbool, Study of flow boiling characteristics of a microchannel using high speed visualization, J. Heat Tran. 135 (2013), 081501. [7] H. Lee, I. Mudawar, M.M. Hasan, Flow condensation in horizontal tubes, Int. J. Heat Mass Tran. 66 (2013) 31–45. [8] H. Huang, J.R. Thome, An experimental study on flow boiling pressure drop in multi-microchannel evaporators with different refrigerants, Exp. Therm. Fluid Sci. 80 (2017) 391–407. [9] A. Kalani, S.G. Kandlikar, Flow patterns and heat transfer mechanisms during flow boiling over open microchannels in tapered manifold (OMM), Int. J. Heat Mass Tran. 89 (2015) 494–504. [10] G.M. Carlmagno, G. Cardone, Infrared thermography for convective heat transfer measurements, Exp. Fluid 49 (2010) 1187–1218. [11] H. Kim, J. Buongiorno, Detection of liquid-vapor-solid triple contact line in twophase heat transfer phenomena using high-speed infrared thermometry, Int. J. Multiphas. Flow 37 (2011) 166–172. [12] T. Oka, Y. Abe, Y.H. Mori, A. Nagashima, Pool boiling of n-pentane, CFC-113, and water under reduced gravity: parabolic flight experiments with a transparent heater, J. Heat Tran. 117 (1995) 408–417. [13] J.P. McHale, S.V. Garimella, Nucleate boiling from smooth and rough surface-Part 1: fabrication and characterization of an optically transparent heater-sensor substrate with controlled surface roughness, Exp. Therm. Fluid Sci. 44 (2013) 456–467. [14] I.S. Park, D.H. Choi, Heat- and mass-transfer analysis for the condensing film flow along a vertical grooved tube, Int. J. Heat Mass Tran. 44 (2001) 4277–4285. [15] I.S. Park, M.Y. Kim, Numerical investigation of the heat and mass transfer in a vertical tube evaporator with the three-zone analysis, Int. J. Heat Mass Tran. 52 (2009) 2599–2606. [16] I.S. Park, Numerical analysis for flow, heat and mass transfer in film flow along a vertical fluted tube, Int. J. Heat Mass Tran. 53 (2010) 309–319. [17] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [18] I.S. Partom, Application of the VOF method to the sloshing of a fluid in a partially filled cylindrical container, Int. J. Heat Mass Tran. 7 (1987) 535–550.

12