A study of rewetting and conjugate heat transfer inﬂuence on dryout and post-dryout phenomena with a multi-domain coupled CFD approach

In the present study, a multi-domain coupled computational ﬂuid dynamics (CFD) approach is developed for the modeling of dryout and post-dryout heat transfer. For the ﬂuid part, the thin ﬁlm and gas core are modeled by the liquid ﬁlm model and two-ﬂuid model, respectively. For the solid part, the heating process is modeled by either using a constant heat source or solving the Joule heating problem. The ﬂuid-solid conjugate heat transfer is calculated by using carefully designed coupling schemes which can automatically determine the operation mode for pre-and post-dryout regions. Unlike standalone simulations where only the inner wall temperature is predicted, coupled simulations are able to predict the outer wall temperature, allowing a direct comparison with experiments. Simulations were carried out for a wide range of ﬂow conditions and validated against the corresponding steady state experiments. By newly introducing a ﬁlm rewetting model, the current CFD code is capable of simulating the transient behavior of dryout. With the rewetting model, the coupled code successfully predicted the dryout hysteresis.


Introduction
Flow boiling is widely used in high heat flux systems, e.g., electronic components and boiling water reactors (BWRs), since vaporization allows the coolant to efficiently remove heat without a significant temperature increase.A main challenge for such systems is the so-called boiling crisis phenomenon, where the efficient heat transfer deteriorates catastrophically due to changes in the underlying heat transfer mechanisms.This can happen to annular flow which is a two-phase flow regime characterized by the coexsitence of a thin liquid film and a fast-moving gas core.When the liquid film is depleted, the efficient film-wall heat transfer is replaced by the much less efficient vapor-wall convection, causing the wall temperature to increase dramatically to an extent that could damage the heating wall.This type of boiling crisis is referred to as dryout and should be accurately predicted due to its potentially dangerous consequences.
The challenge of modeling dryout stems from the fact that it involves multiple complex phenomena.In the pre-dryout region, there are three main heat and mass transfer mechanisms between the liquid film and the gas core: (a) droplet entrainment from the film to the gas core; (b) droplet deposition from the gas core to the film; (c) evaporation from the film to the gas core.In the post dryout region, the vapor phase becomes over-heated, causing those dispersed droplets to evaporate.Such processes should be properly modeled in the prediction of dryout.
The traditional approach for simulating dryout-related phenomena is using 1D empirical correlations where macroscopic quantities, e.g., superficial velocities and cross-sectional void fraction, are used [1][2][3][4] .Since there is no need to simulate flow details, this approach requires negligible computing resources.In addition, such correlations can be very accurate provided that there are adequate experimental data.However, there are two main drawbacks of this approach.On the one hand, empirical correlations are usually constructed based on certain apparatus and flow conditions.Therefore, their application to other flow configurations should be always conducted with caution.On the other hand, such correlations are usually constructed based on steady state data.Therefore, their validity in transient simulations is questionable.
A promising approach to overcome such limitations is to use the computational fluid dynamics (CFD) method, where localized quantities are obtained by numerically solving the corresponding governing equations.Given good enough grids, interface capturing methods, e.g., the volume of fluid method and the level-set method, are able to simulate the complex interface structure in annular flows [5][6][7][8][9] .However, such a fine grid requirement is too demanding to achieve for practical dryout predictions due to the presence of tiny droplets and thin films in long flow channels.Therefore, detailed interface structures are usually not modeled in practical CFD simulations, which can be categorized into two groups.In the first group, a two-fluid multi-field framework is used to describe the annular flow, where a separate volumetric fraction [10,11] or the total liquid phase volumetric fraction [12,13] is used to model the thin film.The occurrence of dryout is determined by the local vapor quality or void fraction.In the second group, the thin film is modeled by the liquid film model [14] and the gas core is modeled by either the Eulerian-Eulerian or Eulerian-Lagrangian methods [15][16][17][18][19][20] .This group of multi-domain simulations can be further divided into three categories: (a) only pre-dryout annular flow is modeled [15,16,19] ; (b) dryout is predicted but no postdryout heat transfer is modeled [20] ; (c) dryout is predicted with the post-dryout heat transfer being modeled [17,18] .In the last two categories, the dryout occurrence is determined by the calculated film thickness.It should be noted that even though all these practical CFD simulations were carried out in a transient manner, only the (quasi) steady state results were compared to the corresponding experiments.Therefore, further investigations are required to make such models capable of simulating the transient behavior of dryout.
Another limitation of previous studies using both empirical correlations and CFD is that most simulations were carried out on the fluid side with the heating structure being modeled as a prescribed heat flux.Therefore, the quality of such heat fluxes plays an important role in the accurate prediction of the wall temperature, making this treatment obviously unsuitable for complex heat structures [21] .Actually, this approach should be revisited even for simple heating structures, e.g., a round tube, due to several reasons.First, around the dryout location, there is a large temperature gradient which leads to an axial heat flux inside the heating structure.Second, in experiments using Joule heating [22][23][24] , temperature changes cause the heating power to redistribute due to the temperature dependency of the solid electrical conductivity.This effect might be strong for the post-dryout heat transfer.Third, during power transients, wall heat flux does not change synchronously with the heating power due to the thermal inertia of the heating structure.There is an additional issue with such a heat flux treatment when coming to the model validation.In dryout experiments, there are two types of wall temperatures, namely, inner wall and outer wall temperatures.Inner wall temperature refers to the temperature at the fluid-solid interface.While outer wall temperature is the temperature at the outer surface of the solid.Only outer wall temperatures are accessible in the experiments, which are then used to deduce the inner wall temperatures by using certain assumptions.Consequently, such experimental inner wall temperatures may themselves be subject to large uncertainties, making them unfavorable for model validations.Therefore, it is necessary and crucial to accurately model the heating structure in dryout simulations.For the multi-field approach, there has been some recent development to model the solid heat transfer [25,26] , where steady state results have been presented.However, such conjugate heat transfer modeling method is yet to be developed for the multi-domain CFD approach.
In the present study, we further develop the multi-domain CFD approach for the modeling of dryout and post-dryout heat transfer.New models are introduced such that the code is capable of simulating dryout transients as well as the heat transfer inside the heating structure.

Methods
The present CFD model solves different governing equations in the fluid and solid regions.

Multi-domain modeling of fluids
The thin film region is modeled by the liquid film model [14] and the gas core is modeled using the Eulerian-Eulerian approach.

Liquid film
For any quantity φ, we introduce the following depth-averaged quantity: where δ is the film thickness and n denotes the wall-normal direction.The depth-averaged governing equations are as follows: where ∇ s is surface nabla operator; S ρδ , S ρU δ , and S ρh δ represent corresponding mass, momentum and energy source terms, respectively.

Gas core
The following governing equations are used for the Eulerian-Eulerian modeling of the gas core based on the two-fluid model [27] ; where k denotes phases; k is the total mass source; M k is the total momentum source; Q k is the total energy source.
For simplicity, we introduce the following notations: subscripts dep, ent and eva denote deposition, entrainment and evaporation, respectively; in the superscripts, f, d and v mean film, droplet and vapor, respectively, and " → " is used to denote the direction of the transfer process.Therefore, we can write out source terms as follows: where M v → d and M d → v are momentum sources that are not related to mass sources and satisfy the following relation In the present study, only drag and turbulent dispersion forces are considered for the droplets since it has been shown that lift and wall lubrication forces are only important for bubbly flows and can be neglected for pre-and post-dryout simulations [13] .The virtual mass is also neglected since droplets are much heavier than the vapor phase.The standard k − model is used for the vapor phase, while no turbulence model is used for the droplet phase.Source terms related to the coupling between the gas core and the film are discussed in Appendix A .The heat transfer between droplets and the vapor phase is modeled as follows where T d is the droplet temperature and is kept at saturation temperature; a i is the interfacial area density which is obtained by solving a simple interfacial area transport equation [28][29][30] ; h eva is the heat transfer coefficient, which is calculated as follows [31] : where κ v is the thermal conductivity of the vapor phase; d 32 is the Sauter mean diameter calculated based on a i ; U r is the relative velocity between droplets and vapor.

Heat transfer in solid
For solid regions, the governing equation for heat conduction reads where ρ s , c p,s and κ s are the solid density, heat capacity and thermal conductivity, respectively; Q s is the volumetric heat source in the solid region.By assuming Q s is uniformly distributed, its value could be calculated based on the wall heat flux and the solid geometry.However, this assumption could be weak in the postdryout related experiments, where Joule heating is used to heat up the solid.The uniform heat source treatment is equivalent to assume that the electrical conductivity, σ e , is constant, which is actually temperature dependent.In order to model the Joule heating effect accurately and to evaluate the uniform heat source assumption, Q s is also modeled in the following way: where V e is the electric potential field obtained by solving the following equation: with a constraint on the voltage difference, V e , across the heated length.However, V e is unfortunately not documented in [22] .Therefore, an iterative procedure is adopted to determine the value of V e such that the total heating power matches the experiment.

Film-gas core coupling for pre-dryout regions
The coupling between the liquid film and the gas core is achieved by modeling different phenomena with closures.The mass source term consists of source terms for droplet deposition, droplet entrainment and film evaporation as follows The liquid film momentum source term is split into two components, namely, shear-stress-based source terms and pressure based source terms.Shear-stress-based momentum source terms consider momentum transfer in the interface-tangential direction as follows The pressure-based momentum sources handle momentum transfer in the interface-normal direction as follows: Energy source terms for the liquid film are: Details of the individual source terms are given in Appendix A

Dryout criteria
In the simulation, when the local film thickness is smaller than a critical film thickness, δ critical , it is assumed that dryout occurs at this spot.δ critical could be calculated using the following correlation [32] : Even though this criterion has been used in previous simulations [17,18,20] , it is not used in the present study since it can be substituted with a much simpler one.
A minimum film thickness, δ min , is often used in the liquid film simulations to maintain the numerical stability.In such a treatment, only patches with δ > δ min are considered wet.Then we can directly use δ critical = δ min for dryout detection.This criterion will be further discussed in Section 4.1 .

Dry patch treatment
The identified dry patches should be excluded from the liquid film calculation, i.e., the corresponding film thickness should be set to zero.This is handled on the matrix level [33] .It should be noted that this dry patch treatment implicitly assumes that once a patch turns dry, it will remain dry, since a zero-thickness film always satisfies two aforementioned dryout criteria.In order to avoid this issue, a rewetting model is introduced to allow certain dry patches to become wet again, as shown in Fig. 1 .The idea is that a dry patch is allowed to become wettable if it has at least one wet neighbor with an excessive film thickness.Therefore, for a given dry patch P , we set it to wettable if there is a neighboring patch > γ , where γ is set to 1.5 in the present study.Then, based on the governing balance equations, the film thickness at P can increase via droplet deposition and film advection.
We note that handling dry patches on the matrix level has a side effect.When dryout occurs, film boundary conditions set at the outlet become redundant since their contribution to the matrix is neglected.In such cases, the real boundary is the wet-dry  borders as shown in Fig. 1 , which are moving during the simulation.Therefore, proper conditions must be applied to these moving boundaries.In the present study, this is achieved by using the first-order upwind scheme for advection terms in the film region.As a result, the wet-dry border will always take values from the wet patch, which is equivalent to setting a no gradient condition to a fixed boundary.

Film-gas core coupling for post-dryout regions
In the present study, the transition boiling region is not considered due to our limited understanding of this complex process.Therefore, for the whole post-dryout region, we assume that the vapor is in direct contact with the tube wall and absorbs all the heat flux.

Coupling schemes
All the aforementioned submodels are coupled together via carefully designed schemes.For instance, the coupling scheme for fluid-only simulations is illustrated in Fig. 2 .Such simulations are conducted within a single fluid solver implemented in OpenFOAM, which handles both calculation and data exchanging.Therefore, these simulations are referred to as standalone simulations in the present study.
The coupling scheme for fluid-solid simulations is more sophisticated due to the fact that film patches need to exchange information with the other two regions simultaneously.This is quite challenging to achieve in a single OpenFOAM solver.Therefore, we adopt a partitioned coupling strategy for the fluid-solid coupling, as illustrated in Fig. 3 .The solid solver solves Eq. ( 17) .The fluid solver is essentially the solver developed for standalone simulations but operates in a different mode.These two OpenFOAM solvers are coupled together via the preCICE library [34] .The fluidsolid simulations are referred to as coupled simulations in the present study.
The coupled simulation can be performed in two modes, namely, with and without thermal inertia of the solid.For the latter, we just need to set ρ s c p,s to zero to deactivate the temporal term in Eq. ( 17) .This is very useful when we are only interested in the steady state results.For this type of coupling, we only need to replace the heat flux in Fig. 2 with the heat flux calculated by the solid solver, as shown in Fig. 3 a.While a different coupling scheme is used for coupled simulations with thermal inertia, as depicted in Fig. 3 b.This is based on the fact that heat flux decreases dramatically right after a wet film patch becomes dry, which cannot be captured by the former coupling scheme.

Simulated cases
Simulations were carried out based on post-dryout experiments by Becker et al. [22] , where the outer wall temperature was measured for a wide range of flow conditions in different electrically heated round tubes.Uncertainties for heating power, flow rate, temperature measurements were ± 1%, ± 1%, ± 3.6 K, respec-  tively.The following factors have been considered when selecting the six flow conditions shown in Table 1 .
Regarding the operating pressure, values around 7 MPa have been almost always used in previous studies [12,13,[15][16][17][18][19][20]26] since it is the operating pressure for BWRs.In this study, we would like to test our models with a wide range of operating pressures.For each selected operating pressure, different flow rates and wall heat fluxes have been chosen such that we can get both saturated vapor and superheated vapor (if possible).Also, two different test sections were selected.The one with smaller inner and outer diameters is denoted as Tube-1, and the other as Tube-2.The thermal conductivities for Tube-1 and Tube-2 are calculated as and respectively [22] .The corresponding electrical conductivities are calculated as follows [22,35] σ ( Such a difference in electrical conductivities allows us to investigate the role of Joule heating effect modeling.

Computational domains
A proper computational domain should be firstly decided since the present model is only capable of modeling annular and mist flows.The straightforward option would be using a computational domain starting from the onset of annular flow.For example, the onset steam quality of an annular flow, x tr , can be determined as follows [36,37] : The required length to get such a quality, L tr , is calculated by: where h in is the water enthalpy at the test-section inlet in the experiment.Then the simulations can start at z = L tr .In the present study, the inlet position is evaluated as follows z = max (L tr , 3 m ) .(31) This allows us to avoid using excessively long computational domains where relatively thick liquid films may exist, as will be discussed in Section 4.3 .The inlet locations are given in Table 1 .

Inlet conditions
Proper boundary conditions must be specified at the computational domain inlets, which is achieved by using energy balance and empirical correlations.The detailed procedure is provided in Appendix B .

Discretization schemes
All the equations are solved numerically using OpenFOAM, which is based on the finite volume (FV) method.A brief description on how to solve the liquid film model using the FV approach is given in Appendix C .Linear interpolations are used to calculate face values based on values at cell centroids.The least squares method is used to calculate (full) gradient terms.Facenormal gradients are calculated by the central-differencing scheme with an explicit non-orthogonal correction.These face-normal gradients are then used for the calculation of Laplacian terms using the Gauss divergence theorem.For the gas core region, divergence terms in the volumetric fraction equation are discretized using the vanLeer scheme.Divergence terms for other scalar and vector variables are discretized using the linearUpwind and linearUpwindV schemes, respectively.These divergence schemes are all second-order accurate.Temporal terms are discretized using the implicit Euler scheme, which is first-order accurate.

Verification
Our verification tests focus on the mesh size dependency study based on the standalone simulations.Each flow condition is simulated with three successively refined hexahedral O-grids (numbered from A to C), as shown in Table 2 .For instance, the crosssectional view of mesh-C for Run-251 is given in Fig. 4 .y * = ρv u τ y μv is the dimensionless boundary distance for a near-boundary cell where boundary refers to film surface and solid wall for the preand post-dryout regions, respectively.y is the corresponding physical distance between the cell centroid and the boundary.u τ is the friction velocity which is calculated based on k .
All simulations were run in transient mode and finally reached steady states.The axial wall temperature and film thickness profiles are plotted in Fig. 5 .It should be noted that, in order to compare results in a compact manner, such profiles are averaged along  the circumferential direction.There are two important conclusions that we can draw from these simulations.On the one hand, the dryout location is insensitive to mesh refinement in all cases.On the other hand, a mesh size dependency is spotted for the post-dryout heat transfer that finer grids tend to predict lower wall temperatures.This discrepancy disappears when droplets are totally evaporated, e.g., regions near outlets for Run-277 and Run-334.Actually, this issue is caused by the modeling of the evaporation of droplets.As shown in Eq. ( 15) , there is a positive correlation between Q evq and T v .For a near-wall cell in the post-dryout region, T v is larger for finer grids since the cell centroid is closer to the heating surface.As a result, there is a higher evaporation of droplets, meaning that less energy will be taken away by the less efficient vapor-wall convective heat transfer.Therefore, the wall temperature will be lower for finer grids.Actually, the enhanced evaporation also increases the convective heat transfer coefficient by increasing the vapor flow rate, which also contributes.This issue will be further investigated in future studies.In the remaining parts of the present study, we will use mesh-B for further discussions unless otherwise stated.

Simulations without thermal inertia
Standalone simulations use the heat flux boundary condition to mimic the solid heating effect.This heat flux is usually calculated based on the steady state heating power which does not consider the thermal inertia effect of the solid.The coupled simulations were first carried out without thermal inertia and compared with the corresponding standalone simulations, as shown in Fig. 6 .For solid regions, there are always ten layers of cells in the thickness direction.In addition, two heating approaches, i.e., uniform heating and Joule heating, were both used for the solid region.
The most obvious advantage of the coupled simulation is that the outer wall temperature is available, allowing a direct comparison with the experiment.
Two test sections respond differently to the solid heating approaches.For Tube-1, almost the same results are predicted by the two heating approaches.The reason is that, even though Tube-1 has a temperature-dependent electrical conductivity, as shown in Eq. ( 27) , σ 1 e varies in a tiny range, which is essentially equivalent to the uniform heating approach.However, for Tube-2, quite dif-ferent results are predicted.The reason is that σ 2 e varies in a much larger range and decreases with temperature.As a result, in comparison with the uniform heating approach, post-dryout regions have larger heating powers in the Joule heating approach, resulting in higher wall temperatures.This is proved by inspecting the heat flux at the inner wall, as shown in Fig. 7 .For the uniform heating approach, the heat flux only shows different values around the dryout location, which is caused by the internal heat flux in the solid region that goes from high-temperature dry regions to lowtemperature wet regions.However, for the Joule heating approach, in addition to this abrupt change around the dryout location, other regions also have different heat fluxes.Pre-dryout regions have rather uniform heat fluxes since the solid temperature is uniform along the flow direction.The heat flux for the post-dryout region decreases in the flow direction since electrical resistance drops in the same direction.The heat flux difference between the pre-and post-dryout regions is as high as 15%.We note that Tube-2 was made of stainless steel, which is widely used in experiments and industrial applications.Such findings challenge the standalone approach since the uniform heat flux assumption could be poor.
Since the transition boiling is not considered in the present study, the gradual temperature increase in this region cannot be captured in the simulations.Instead, wall temperatures increase sharply at the dryout point.For the fully developed post-dryout region, qualitatively correct temperature trends have been predicted for all the cases.However, only simulations for Run-277 can quantitatively match the experiment.Relatively large errors are spotted for Run-497 and Run-510 where the operating pressure was 15 MPa.In our opinion, the reason might be that deposition and entrainment closures described in Appendix A become inaccurate for high operating pressures.

Simulations with thermal inertia
A postulated power transient based on Run-277 is used to assess the code performance when the thermal inertia effect is considered.For the standalone simulation, there is no thermal inertia.Consequently, all the power is taken away by the fluid region.Therefore, we can use the wall heat flux in the standalone simulation to show how the input power changes during the transient, as depicted in Fig. 8 .The transient starts from a fully developed temperature and flow fields with a low heating power (20% of the heating power in Run-277).Then it undergoes two step increases and two step decreases, where the time interval between power changes is set to 20 s.This heat flux transient is used to calculate the corresponding uniform heating source in the coupled simulation, where ρ s = 8370 kg/m 3 and c p,s = 461 J/kg/K are used in Eq. (17) .
As shown in Fig. 8 , in comparison with the standalone simulation, the wall heat flux at the fluid-solid interface has an obvious delayed response to both power increases and decreases, which is caused by the thermal inertia of the solid region.In addition, there are more interesting phenomena in Fig. 8 that are better explained by inspecting Fig. 9 , where the wall film coverage is used as an indicator for the dryout location.
First of all, we could clearly see the effect of the rewetting model described in Section 2.5 , without which the film coverage would not increase during the power decrease transients.
Another finding is that, during the power increase transients, there are two slight heat flux drops around 5 s and 23 s in Fig. 8 , which are caused by newly occurring dry patches that worsen the heat transfer at the wall.After the film coverage stabilizes, the heat flux starts to increase again.
The third interesting phenomenon is the so-called dryout hysteresis.Dry patches form rather fast in the power increase transients.However, the disappearance of dry patches is much slower     when power decreases, as clearly seen in Fig. 9 .The solid thermal inertia plays an important role in this process, since it takes time for the wall heat flux to decrease, as shown in Fig. 8 .

Sensitivity studies
Sensitivity studies have been carried out for several parameters to help us better understand the performance of present models.A sensitivity study on this parameter has been carried out using Run-251, as shown in Fig. 10 .As expected, dryout occurred earlier for some large δ min values.However, the same results have be obtained for δ min ≤ 5 × 10 −6 mm.Therefore, one can choose any value in this range.
δ min also plays a role in the rewetting modeling, as shown in Fig. 1 .For a given γ , a smaller δ min means that the rewetting event is easier to be triggered.The postulated transient has been carried out to check this effect with different δ min for Run-277 with mesh-A, as shown in Fig. 11 .It can be seen in the zoomed-in view that rewetting occurs a little earlier for smaller δ min values.However, this difference is barely visible when inspecting the film coverage   evolutions for the whole power transient.Therefore, the value of δ min has a rather limited effect on dryout transients.

γ
In terms of the rewetting modeling, γ plays a role similar to δ min that a smaller value makes it easier for rewetting to occur, as confirmed by the zoomed-in view in Fig. 12 .Nevertheless, this difference is negligible for the whole rewetting process.

Inlet elevation
One important reason for using Eq. ( 31) to determine the elevation of the inlet in simulations is that the film might be relatively thick at the beginning of annular flow regions.In such cases, the liquid film model described by Eqs. ( 2) -( 4) may become erroneous.A sensitivity study on z has been carried out for Run-334, where L tr ≈ 0.76 m.As shown in Fig. 13 , for the z = 1 m case, the film thickness could exceed 0.8 mm, which is rather thick, considering the fact the tube inner diameter is only 14.9 mm.In comparison with higher elevations, different results have been predicted for the dryout location and post-dryout wall temperature.The consistency within high-elevation results also indicates that the inlet conditions described in Appendix B are rather robust.

Conclusions and outlook
In the present study, the multi-domain coupled CFD approach is further developed for the modeling of dryout and post-dryout heat transfer, which is achieved by introducing several key submodels.A rewetting model has been developed which enables the code to simulate both the occurrence and disappearance of dryout.Advanced coupling schemes have been developed such that the gas core-film-solid coupling can be carried out efficiently in different modes.In com parison with standalone simulations, coupled simulations could predict the outer wall temperatures of the test sections, allowing a direction comparison with the experiment.The coupled simulation also shows its capability of modeling more complex heat transfer phenomena, e.g., power redistribution caused by the Joule heating effect and the dryout hysteresis.In addition, sensitivity studies show that the newly introduced models perform quite robust despite the fact that few tunable coefficients are used.
More developments are still required in the field.One direction is trying to develop a better droplet evaporation model to eliminate the grid size dependency issue for the post-dryout heat transfer.The other is modeling the gas core with the Eulerian-Lagrangian approach such that there is no need to introduce closures for droplet deposition.It will also be useful to develop models for the transition boiling region to get a better agreement with the experiment.In addition, this would allow us to further introduce wall temperature into the rewetting modeling.Nevertheless, more localized mechanistic closures should be developed for the coupling between the film and the gas core, especially for the droplet entrainment process, such that better predictions could be obtained for high-pressure flow conditions.and the velocity at the interface between the film and the gas core is 28) This is used as the boundary conditions for the gas-core velocities.The film-vapor shear stress term is calculated as where C f,i is the interfacial friction factor which is calculated as follows [36] :

A2.2. Gravity terms
The film-normal term is calculated as and the film-tangential contribution is calculated as (A.32)

A2.4. Gas core pressure
The gas core pressure p c is taken from the adjacent gas core cell.

B1. Primary quantities
For a given distance, z , to the inlet in the experiment, the corresponding steam quality is calculated based on the energy balance: At the inlet of the computational domain, the total flux is split into vapor flux G v , droplet flux G d and film flux G f : where j v = G v ρv , and ρ c is calculated according to In the above entrained fraction model, an assumption is made as follows: 11) In order to calculate the film thickness, the void fraction is introduced, which is defined as: Then, the void fraction α is calculated using the following correlation [43] : With all such information, α d , U v , U d , U f and δ f could all be calculated for the inlet of the computational domain.

B2. Auxiliary quantities
There are several auxiliary quantities to be calculated at the inlet.One is d 32 which is calculated based on the correlation by Caraghiaur and Anglart [44] :

Appendix C. FV discretization of the liquid film model
In order to solve Equations ( 2) -(4) using the FV method, the investigated surface is first extruded in its normal direction with a thickness of δ 0 to form a 3D domain, which is then discretized into a finite number of non-overlapping and convex control volumes, as shown in Fig. 15 .It should be emphasized that we could only have one layer of cells in the surface-normal direction.
The next step is to integrate the governing equations over each discretized cell, P .The integrals of temporal and common source terms can be easily calculated.We only focus on the calculation of the integrated advction term V P ∇ s • φ d V, where φ = ρ f δU f φ.
One basic assumption for the liquid film model is that the filmnormal advection is negligible:  Therefore, the integrated advection term is calculated as follows using the Gauss divergence theorem: where φ i is the value of φ at the centroid of face i .
In order to enforce the condition stated by Eq. (C.1) , the surfacenormal fluxes, i.e., φ so • S so and φ ex • S ex , should be set to zero.This is achieved by setting a no gradient boundary condition for the velocity field.Actually, this boundary condition is used for all the involved variables such that we can use ∇ and ∇ 2 to calculate ∇ s and ∇ 2 s , respectively.The extruded mesh not only introduces the extruded surface, but also shifts the position of the investigated surface, as illustrated in Fig. 16 .In the FV framework, discretized integral-form equations are solved for cell centroids ( C FV in Fig. 16 ).In fact, the locations of such cell centroids determine the FV surface.Therefore, there is always a δ 0 2 gap between the FV surface and the source surface.For curved surfaces, such gaps lead to curvature differences between the FV surface and source surface.Therefore, it is necessary to take this gap into consideration when generating the mesh.
For film-only simulations, there is an obvious solution that we could explicitly specify the locations of the source and extruded surfaces such that the FV surface is exactly located at where it should be.While when the liquid film model is coupled with the gas core, the source surface is always selected as the physical boundary (tube wall in this case), meaning that we cannot move the source surface as we do in the film-only cases.Since it is no longer possible to eliminate the gap effect, the idea is to reduce this effect.Then, the solution is straightforward that we can use a small value for δ 0 such that the FV surface is close to the physical boundary.
In order to evaluate the influence of δ 0 , we use the following  17 .Qualitatively, the result is just as expected that a smaller δ 0 leads to a better prediction since the FV surface is closer to the source surface.Quantitatively, the result is very encouraging since the difference is rather small.If δ 0 R is set to 0.1, the largest difference is around -0.2%, which is a very decent result for such simplified modelings.In the present study, δ 0 = 0 . 1 mm has been used, resulting in a maximum δ 0 R ≈ 0 .013 .Therefore, the corresponding gap effect is believed to be negligible.

Fig. 1 .
Fig. 1.Rewetting model for the liquid film region: dry patch P is set to wettable if it has at least one wet neighbor with an excessive film thickness.

Fig. 5 .
Fig. 5. Film thickness and wall temperature profiles calculated by different grids.

Fig. 8 .
Fig. 8. Mean normalized wall heat flux at the fluid-solid interface during the power transient.

Fig. 9 .
Fig. 9. Film coverage at the wall during the power transient.

. 5 )
where the entrained fraction e = G d G d + G f is calculated based on the following correlation by Cioncolini and Thome [42] : e = 1 + 279 .6 W e −0 .

Fig. 15 .
Fig. 15.Finite volume discretization using one-layer polyhedral control volumes.Each bounding surface of cell P has a corresponding vector, whose direction and magnitude denote the orientation and area of the surface, respectively.S so and S ex are the vectors for the source and extruded surfaces, respectively.Other vectors are denoted by S i .

Fig. 16 .
Fig. 16.Side view of the FV mesh generated from a source surface.

Fig. 17 .
Fig. 17.Evaluating the influence of the extrusion thickness in the finite volume approach based on the analytical solution for laminar cylindrical falling film.

4 ) 2 ,
analytical solution for U for a free falling film on a vertical cylindrical surface R is the radius of the cylindrical surface.We further define the following liquid film Reynolds number:Using this analytical solution, we could evaluate the influence of the extrusion thickness, δ 0 , as shown in Fig.16.For a given combination of R and δ, U and Re are calculated according to Eq. (C.3) and Eq.(C.4) , respectively.In the finite volume approach, the radius of the FV surface is R + δ 0 as shown in Fig.16.Therefore, we substitute R with R + δ 0 2 in Eq. (C.3) to get U δ 0 , which is the depth-averaged velocity corresponding to the FV surface.The relative difference η = U δ 0 −U U is calculated for different R and Re values, as shown in Fig.

Table 1
Key parameters for simulated experiments.
Then we can calculate the inlet condition for a i based on d 32 and α d .Other quantities are related to the turbulence modeling.Inlet values for k are calculated based a turbulent intensity of 5%.The inlet turbulent length scale is estimated to be 0.07 D h during calculation.