The VOF model proposed by Hirt and Nichols (1981) is one of the most widely used two-phase CFD models for capturing the interface. In the VOF model, two-phase flows are recognized as homogeneous mixtures consisting of immiscible fluids. The continuous and dispersed phases can be lumped into a single continuum in which they share the same pressure and velocity fields. Therefore, the VOF model is not capable of describing dispersed phase flow in which the interface length is smaller than the mesh size.

The mixture model is applicable to interpenetrated phases which are treated as a mixture in a cell. The mixture model is similar to the VOF model in that it solves a single set of transport equations. However, the model cannot be expected to resolve a sharp interface even on a fine grid, as the VOF does. In addition, the phase slip velocity is considered in the mixture model; hence, the behavior of the dispersed phase, which cannot be captured with a coarse grid, can be properly simulated.

A commercial CFD code, STAR-CCM + v13.02, introduced the VOF–slip model as a hybrid model combining the VOF and the mixture models. This hybrid model can treat a large-scale interface and small-scale interface simultaneously. In the VOF–slip model, the volume fraction transport equation and momentum equation can be expressed as follows: ∂ α i ∂ t + ∇ · ( α i v ) = - 1 ρ i ∇ · ( α i ρ i v d , i ) (1) ∂ ( ρ v ) ∂ t + ∇ · ( ρ v v ) = − ∇ p + ρ g + ∇ · T + f s + ∇ · Σ α i ρ i v d , i v d , i (2) where α i is the volume fraction of phase i , ρ i is the density of phase i , v is the mixture velocity, v d , i is the diffusion velocity, p is the pressure, T is the stress tensor, and f s is the surface tension force. The above two equations are the same as those that result from adding a diffusion velocity term to the existing transport equations in VOF. The diffusion velocity is the difference between the phase velocity and the mixture velocity, and it is defined as follows: v d , i = v i − v (3) The diffusion velocity and slip velocity, v p s are related as follows: v d , p = ( α p ρ p ρ − 1 ) v p s (4) v p s = v s − v p (5) where p and s denote the primary and secondary phase, respectively. Here, the slip velocity that determines the diffusion velocity can be modeled with a drag coefficient C D and specific body force b , as shown below. v p s = C D b (6) b = g − ( v · ∇ ) v − ∂ v ∂ t (7) In STAR-CCM+, the Schiller–Naumann drag model (Schiller and Naumann, 1935) can be used to obtain C D , which assumes that the primary phase (referred to as the gas phase in this paper) is continuous and that the secondary phase (referred to as the liquid phase in this paper) is dispersed. The drag coefficient is modeled as follows: C D = C D p s = τ s f d r a g p s   ( ρ s − ρ ) ρ s (8) τ s = ρ s d s 2 18 μ p (9) f d r a g p s = { 1 + 0.15 R e p s 0.687   i f   R e p s ≤ 1000 0.0183 R e p s                   i f   R e p s > 1000 , R e p s = ρ p | v p s | d s μ p (10) where d s is the droplet diameter and μ p is the viscosity of the primary phase. Equations 8–10 show that the drag coefficient increases as the droplet diameter increases, which means that the slip velocity increases.

The wettability of the wall is reflected as a contact angle that is required in VOF-slip model. Because the liquid film is formed at the acrylic test section and the film sensor which surface is made of the polyimide film and gold layer, it was difficult to determine the value of the contact angle. Besides, there are anisotropic surface irregularities due to the electrode arrangement in the sensor, which makes the distortion of the wall friction inevitable. For this reason, in the case of no air injection, the contact angle that can yield a reasonable film spreading width was found to be 0°. However, when the air velocity became large enough to cause an entrainment, it was confirmed that the effect of the contact angle became minor and the gravitational force and the interfacial friction force had a dominant influence on the simulation results.

The diffusion velocity terms in the transport equations of the VOF–slip model (Eqs 1, 2) should not be activated in the gas–liquid interface region, and they should only work effectively in the unresolved dispersed phase. Therefore, a different value of the droplet diameter was used, based on the γ term defined as follows: γ = | ∇ α |   d x (11) d s = { d s ∗                                                                                                                   i f   γ ≤ 0.33   − d s ∗ 0.17 ( γ − 0.5 )                                                   i f   0.33 < γ < 0.5 1 × 10 − 7                                                                                               i f   γ ≥ 0.5 (12) where d x is the cell size and d s ∗ is the input droplet diameter.

The variable γ plays a role in determining whether the regime in the cell is a free-surface or unresolved dispersed phase. If its value is greater than 0.5, the corresponding region is interpreted as an interface. For this case, a very small droplet diameter of 0.1 μm is utilized in the slip velocity estimation so that the VOF approach can be applied, i.e., without phase slip. If γ is less than 0.5, it is interpreted as a dispersed-phase region. In this case, the droplet diameter set by the user is utilized in the slip velocity estimation for the interpolation region.

Another required condition is that the droplet diameter should be used only when the liquid volume fraction is reasonably large. This condition can be expressed as follows: d s = { d s ∗                                                                                   i f       6 α l V g r i d π 3 ≥ d s ∗ 1 × 10 − 7                                                           i f       6 α l V g r i d π 3 < d s ∗ (13) where α l is the liquid volume fraction and V g r i d is the volume of the grid.

The k − ω SST (shear-stress transport) model (Menter, 1994) was used as the turbulence model in the present study. This model, which combines the best characteristics of the k − ω and k − ε turbulence models, behaves like a k − ω -type model near the wall, otherwise it behaves as k − ε -type model, which avoids strong freestream sensitivity.

The simulation of the liquid film off-take near the broken cold leg should consider the interface turbulence damping so that the estimated interfacial friction between the gas and liquid phases leads to reasonable predictions of the liquid film behavior and bypass phenomenon. One of the most widely used ways of damping the interface turbulence is to apply the Egorov model (Egorov et al., 2004). The Egorov model is only available with the k − ω model, however, other models such as the k − ε model may also benefit from an Egorov-like turbulence damping term near large-scale interface (Frederix et al., 2018). In this model, a source term is added to the ω equation of the k − ω model, which enhances the specific turbulence dissipation rate as follows: S i = α A i · Δ n β ρ i ( ω i ) 2 = α A i · 36 B 2 μ i 2 β ρ i Δ n 3 (14) ω i = B 6 μ i β ρ i Δ n 2 (15) where B is a damping coefficient and ∆n is the cell height across the interface. This source term is only activated in the interfacial region by introducing an indicator A i which is the interfacial area density. As the value of the damping coefficient B increases, the specific dissipation rate increases, which makes the eddy viscosity decrease.

Thus, determining an appropriate value for B is crucial to predicting how the free surface behavior is affected by the interfacial friction. However, there has been no general guideline provided for the selection of B in previous studies. In most cases, B has been tuned to match the experimental results well (Hohne et al., 2002; Egorov et al., 2004; Lo and Tomasello, 2010; Gada et al., 2017; Fan et al., 2019). Therefore, in the present study, a sensitivity analysis was performed to investigate the effect of B on the simulation results, and the value of B that yielded reasonable results for each simulation case was determined.

Because the effect of interface turbulence damping is also dependent on the cell height, as shown in Eq. 15, Frederix et al. (2018) introduced a mesh-independent length scale δ , which was incorporated into the Egorov approach to obtain consistent results regardless of the grid cell size, as shown below: δ 2 = Δ n 2 6 B (16) then the ω i term in Eq. 15 can be expressed independent of the cell size by adopting δ as follows: ω i = μ i β ρ i δ 2 (17) Using this approach, Frederix simulated a co-current stratified flow in a parametric study of δ, and found that the maximum value of δ was 10^–4 m. In the present study, the effect of interface turbulence damping was first investigated in terms of B. Then, an appropriate value of B that yielded reasonable results was converted to δ following the approach used by Frederix.

The test section in the SNU experiment (Choi and Cho, 2019) was modeled as shown in Figure 2A. In the experiment, the test section describes 1/10th scale of APR1400 reactor vessel downcomer. There were two intact cold legs for air intake and a broken cold leg for air–water outtake. The DVI nozzle was placed above the broken cold leg, and water was injected through this nozzle. The liquid film sensor installed at the inner wall of annular test section measured the film thickness in the region between the DVI nozzle and the cold legs. To reduce the computational cost, the lengths of the cold legs and the DVI nozzle were reduced, and the velocity profiles obtained from the fully developed flow were used as inlet boundary conditions. The bottom of the domain was set as a water outlet that controlled the water out flow rate to keep the water level constant.

Figure 2B shows the meshing configuration. Trimmed meshes were used, and 13 prism layers were generated near the wall. Considering the requirements of the turbulence model used, the k − ω SST model, wall y + values were maintained at approximately 1. The total number of cells was 6,223,381. Mesh convergence tests were conducted by changing the number of near-wall meshes. The results are presented in Simulation Under No-Airflow Conditions Section, along with the liquid film thickness prediction results.

The flow conditions for the simulations are summarized in Table 1.

For the case of no airflow, a transient simulation with a time step of 0.2 ms was conducted, as shown in Figure 3 (W089A00). The time step was determined based on the local Courant number, which should be less than 1.0 for applying the High-Resolution Interface Capturing (HRIC) scheme suited for tracking sharp interfaces. The film interface was expressed by the iso-surface of the void fraction at 0.5. In the simulation, the injected water from the DVI nozzle impinged on the wall and spread in the form of a liquid film. Because the size of a computational cell in CFD is about 70 times smaller than that of a measurement area on the liquid film sensor, it was necessary to spatially average the simulation results for quantitative comparison with the experimental results. Therefore, the calculated film thickness in the areas corresponding to the measurement point on the sensor was extracted and spatially averaged. Then, the averaged thickness was compared with the experimental result and used to validate the CFD simulation.

The front view of the averaged film thickness is compared with the experiment results in Figure 4. It was found that the CFD predicted the position of the thick film boundary and the film distribution reasonably. A quantitative comparison of the liquid film thicknesses is shown in Figure 5A. The film thicknesses were comparable overall, but there were some points at which the liquid film thickness was greatly underestimated (see the area marked in gray). The points at which significant errors occurred are shown in Figure 5B. In the figure, x = 0 corresponds to the center of the broken cold leg, and the peak of the graph appears on the film boundary region. Large discrepancies in film thickness could be seen at positions where the change in the film thickness was large in the film boundary region. Nevertheless, the CFD model predicted the peak value of the film thickness and the position at which the value appears well, which is important for the prediction of the ECC bypass induced by entrainment.

To confirm the validity of using the current mesh configuration near the wall where the liquid film is formed, the effect of the mesh size on the film thickness prediction was investigated as shown in Figure 6. In the case of a coarse mesh (n = 3), the film thickness was predicted to be thinner than in other cases. When the number of near-wall meshes was greater than 10, the film thickness converged to a specific value. Consequently, taking into consideration the computational cost, the current mesh configuration (n = 13) was judged to be adequate for describing the film behavior.

To simulate the film off-take phenomenon with airflow, the boundary conditions for the two intact cold legs were set as the velocity inlets. Figure 7 shows the streamline of the airflow in the computational domain. The injected air first impinges against the inner wall of the downcomer, and radial airflow is formed near each intact cold leg. Then, the air flows to the broken cold leg, which narrows the spreading width of the liquid film and leads to film off-take.

The results obtained for the shape of the liquid film were compared with visual observations from the SNU experiment as shown in Figure 8. The visual observations (Figure 8A) show that the flow rate of the liquid film affects how much the film boundary is dragged for a given air flow rate. When the film flow rate is higher, the film boundary is dragged less by the airflow because of the greater momentum of the film flow. As shown in Figure 8B, the CFD simulation also successfully reproduced this film behavior.

Based on the simulation results, two mechanisms were found to contribute to ECC bypass, as shown in Figure 9. The first mechanism is the entrainment phenomenon that occurred at the thick film boundary. This mechanism involved the roll wave formed at the thick film being sheared off by the airflow and droplets then being generated from the roll wave. These generated droplets flowed into the broken cold leg with the airflow. If the liquid film boundary spreads farther from the broken cold leg, the bypass rate could be inferred to decrease because the position of the film boundary is related to the distance the droplets must travel toward the break. The second mechanism is the film off-take phenomenon that occurred near the broken cold leg. Under conditions of high airflow, the thick film region appeared near the broken cold leg as the gravitational and interfacial friction forces were balanced, creating a hanging liquid film. This thick film could easily be stretched by the airflow, and wisps in the form of large liquid lumps were generated. Then, the wisps that flowed to the broken cold leg retained their shape or broke up into droplets. These principal mechanisms of ECC bypass detected in the CFD simulation results were also confirmed in a previous experimental study (Choi and Cho, 2002). Thus, it was qualitatively demonstrated that the proposed CFD model can accurately predict the film off-take phenomenon, including the entrainment at the thick film boundary.

As described above, the prediction accuracy of the ECC bypass rate is determined by the behavior of the liquid film and the detached liquid (wisps and droplets). Therefore, further investigations of the parameters that affect the liquid phase behavior by means of quantitative assessment of simulation results are needed.

When airflow toward the broken cold leg occurred, the spreading width of the liquid film became narrower, and entrainment occurred from the liquid film boundary. Under the flow conditions for the simulation, the entrainment started to occur when j g , o reached 20 m/s regardless two different liquid flow rates, which was also confirmed in the experiment. This indicates that the sensitivity of the gas inertia related to the interfacial shear stress is dominant for determining the onset of entrainment. The generated droplets flowed to the broken cold leg with the airflow and contributed to the ECC bypass. Among these droplets, those that the mesh could not resolve because of their small size were treated with the mixture approach. To determine the slip velocity for the unresolved droplets, the droplet diameter had to be known, but it was not measured in the SNU experiment. It can however be estimated that the mean droplet diameter would be in the range of approximately 100–200 μm, based on existing empirical correlations (Tatterson et al., 1977; Kataoka et al., 1983) and the results of experiments (Zaidi et al., 1998; Hurlburt and Hanratty, 2002; Westende, 2008) in which the droplet diameter was measured under conditions similar to those of the SNU experiment. Therefore, in this study, simulations were performed for droplet diameters in the range of 100–200 μm, and a reasonable diameter was determined based on the calculated ECC bypass fraction. The simulations for determining the droplet diameter were carried out under relatively low airflow conditions ( j g , o ≤ 22   m / s ) to minimize the effect of wisps generated by strong airflow on the bypass fraction, so that the droplet size effect could be independently confirmed.

Figure 10 shows the variation in the ECC bypass fraction with the droplet diameter. The bypass fraction was obtained from the water injection rate and the water flow rate to the broken cold leg, as follows. ( B y p a s s   f r a c t i o n ) = m ˙ f , b r e a k m ˙ f , i n (18) As the droplet diameter increased, the drag coefficient and the phase slip velocity increased according to Eqs 6–10. The increase in the slip velocity indicates the decrease in the velocity of the unresolved droplets. Accordingly, the ECC bypass fraction decreased as the droplet diameter increased. It was confirmed that the bypass fraction was nearly saturated with at droplet diameter of 50 μm. When the VOF model was applied, the bypass fraction was overestimated by more than two times because phase slip was not considered for the dispersed phase flow. Based on fact that a 150 μm droplet diameter in the VOF–slip model yielded reasonable bypass fraction prediction results for this condition, a droplet diameter of 150 μm was used for all simulation cases.

Under high-airflow conditions ( j g , o > 22   m / s ), the spreading width of the liquid film became much narrower, and a thick film was formed near the broken cold leg. A strong airflow toward the broken cold leg made the thickened film stretch, generating wisps. This implies that not only the unresolved droplet behavior but also the liquid film and wisp behavior are major factors that affect the ECC bypass fraction under high-airflow conditions. This free-surface behavior is closely related to interfacial friction, which may have to be reduced using the Egorov damping model. Therefore, the effect of the damping coefficient B on the simulation results was investigated.

Figure 11A shows snapshots of the film spreading width for different values of B (W063A24). As B increases, the spreading width of the liquid film near the broken cold leg increases because of reduced interfacial friction. As shown in Figure 11B, the CFD model was able to predict the film boundary most comparable to the experiment results for B = 10. In the figure, the dotted white lines indicate the peak position of the film boundary confirmed in the SNU experiment.

The effect of B on the bypass fraction is shown in Figure 12A. When the turbulence was not damped at the film interface, the overestimated interfacial friction made the spreading width of liquid film excessively narrow, which caused a large error in predicting the bypass fraction. It should be noted that the calculated bypass fraction for B = 10 was the most comparable to the experiment results. This means that accurate simulation of the free-surface behavior achieved by adopting a suitable B value led to reliable prediction of the bypass fraction as well. The prediction of the film spreading width is directly related to the position at which droplets are generated, as well as the formation of wisp flow around the broken cold leg, as shown in Figure 12B. If an excessive damping effect occurs at B = 30, wisp flow is not generated around the broken cold leg because of the excessively wide liquid film width, and as a result, the bypass rate is reduced.

Figure 13 illustrates the calculated bypass fraction for B values of 0, 10, and 30 for each simulation case. The figure shows that the appropriate value of B depends on the flow conditions, which is consistent with the results of previous studies described in Interface Turbulence Damping Section . Under low-airflow conditions, for which the j g , o is less than 22 m/s, B has little effect on the bypass fraction, which can be predicted well even without the damping model. This means that there is no need to dampen the interface turbulence when the velocity difference between the phases at the interface is small, which can no longer be regarded as wall-like behavior. When j g , o exceeds 22 m/s, the appropriate value of the interface turbulence damping coefficient ranges from 10 to 30 in terms of B and from 2.7 × 10^–5 m to 5.4 × 10^–5 m in terms of δ. It is worth noting that the range of δ identified in the present study is consistent with the results reported by Frederix et al. (2018) (in which the maximum value of δ was 10^–4 m). However, finding an appropriate value for B or δ for each flow condition by means of multiple simulations is a very time-consuming process; hence, modeling of interface turbulence damping for various flow conditions is required in future research.

The predicted liquid film width and thickness obtained using a suitable value of B was compared with the experiment results, as shown in Figure 14. Although the spreading width of the liquid film was under-predicted in the CFD simulation, the peak position of the film boundary was comparable. Because entrainment mostly occurs at the thick film, it can be deduced that the accurate prediction of these peak positions results in reliable prediction of the ECC bypass phenomenon by the CFD simulation.