Numerical Investigation of Water Film Evaporation with the Countercurrent Air in the Asymmetric Heating Rectangular Channel for Passive Containment Cooling System

Passive containment cooling system (PCCS) is an important passive safety facility in the large advanced pressurized water reactor. Using the physical laws, such as gravity and buoyancy, the water film/air countercurrent flow is formed in the external annular channel to keep inside temperature and pressure below the maximum design values. Due to the large curvature radius of the annular channel, one of the short arc segments is taken out, as a rectangular channel, to analyze the main water film evaporation heat transfer characteristics. Two numerical methods are used to predict the water film evaporative mass flow rate during the heat transfer process in the large-scale rectangular channel with asymmetric heating when the water film temperature is not saturated. At the same time, these numerical simulation results are validated by the experiment which is set up to study water film/air countercurrent flow heat transfer on a vertical back heating plate with 5m in length and 1.2m in width. It is shown that the maximum deviation between numerical simulation and experiment is 30%. In addition, the influences on these parameters, such as heat flux, evaporative mass flow rate, and water film thickness, are evaluated under the different tilted angles of the rectangular channel and horizontal plane, water/air inlet flow rates, water/air inlet temperatures, heating surface temperatures, and air inlet relative humidities. All these results can provide a good guidance for the design of PCCS in the future.


Introduction
In recent years, more and more passive safety facilities have been applied in the large advanced pressurized water reactor. For example, the passive containment cooling system (PCCS) is used in AP600 [1], AP1000 [2], and CAP1400 [3], mainly to remove the decay heat from the postulated accident scenario and prevent the release of the radioactive materials.
e main characteristic of PCCS is that the gravity-driven water film flows downward along the outer surface of containment wall with the countercurrent buoyancy-driven air flowing above the water film in the annular channel. e water film evenly covering the heated shell wall cools down the containment together with naturally circulating air, and the water film evaporation will dominate the heat transfer process [4]. Because PCCS does not rely on electric power or active mechanics and the number of equipment simplified, it is adopted in the larger advanced power reactor. In order to determine the actual heat removal capacity of PCCS, it is very important to understand the relevant influencing factors by establishing the experiment facility with extended test conditions and carrying out corresponding numerical research studies.
For the experimental research, Forgione et al. [5] experimentally investigated water film evaporation on the heating plate with 2.0 m long and 0.6 m wide with countercurrent air flow. ey mainly analyzed the influence of channel depth on the heat transfer and fitted a heat transfer correlation considering the short duct entrance effect (L/De). Later, Ambrosini et al. [1] also conducted a series of evaporation tests in the channel with 2.0 m long, 0.6 m wide, and 0.1 m high. ey confirmed that heat and mass analogy method could accurately predict the heat flux. Similar to the study of Ambrosini et al. [1], Kang and Park [6] set up a test facility to study falling water film evaporation with the countercurrent air flow in the channel which had a 2.0 × 0.5 m 2 heating plate. e facility had an adjustable gap height ranging from 5 cm to 20 cm. Based on their experiment results, they developed a correlation considering vapor concentration on the film surface. ey also found that water film wave which depended on the film temperature and air velocity affected the heat transfer characteristics. Huang et al. [7] conducted the experiment about the influence of water film wave on the process of water film/air countercurrent flow heat transfer. ey believed that both the water film Reynolds number and the air Reynolds number affected the water film wave, but the air Reynolds number had a greater effect. Recently, Hu et al. [4] carried out the water film evaporation tests in a large-scale rectangular channel considering the countercurrent air flow. e channel was 5.0 m long, 1.2 m wide, and 0.3 m high. e bottom surface temperature of the channel was changed by adjusting the input electric power of the conducting oil heater. Based on their experimental results, they fitted an evaporation mass transfer relationship that mainly considered the influence of the channel gap height on mass transfer.
For the numerical research, both lumped parameters (LP) codes and computational fluid dynamics (CFD) codes are used to investigate water film evaporation characteristics. Ren and Wan [8], respectively, utilized the 1-D LP model and 2-D CFD model to study the water film evaporation model with the upward moist air in a vertical rectangular channel by analyzing the average Nusselt number (Nu) and Sherwood number (Sh). ey found that 1D model could be utilized to retrieve the 2D model results. Besides, they also investigated the influence of different wall boundary conditions [9] on evaporation heat transfer characteristics. Based on the assumption that the water film thickness was very thin and evaporation only occurred on the liquid film surface, Du et al. [10] proposed a 1D simplified PCCS steady heat transfer model coupled with the steam condensation process on the inside surface of the containment wall. Comparing with the experimental data of Kang and Park [6], this 1D model was reliable for the heat and mass transfer process. ere had been a lot of analysis codes that used to analyze the heat removal capacity of PCCS, such as MEL-COR [11], COCOSYS [12], WGOTHIC [13], MAAP4 [14], COMMIX [15], and CAST3M [16]. However, these current codes were lack of the experiment validation and did not consider the gradient distribution of parameters in the direction of the channel gap height, so their application was limited. Recently, in the 3-D GASFLOW-MPI CFD code, a dynamic water film model was developed by Xiao et al. [17] through adding essential interphase exchange source terms in the governing equations of mass, momentum, and energy. ey believed that this film model in GASFLOW-MPI code could efficiently analyze the heat removal capacity of PCCS, from the excellent comparisons of the calculation results between CFD code and exact numerical solution (ENS). In addition, Hong et al. [18] pointed out that water film evaporation could be simulated by the Lee model [19] of FLUENT. However, the Lee model could not accurately predict the water evaporation heat transfer characteristics when the temperature of water film was not saturated. Ambrosini et al. [20,21] thought the water film evaporation on the heated wall was a transpiration problem, and the water vapor concentration of this heated wall was saturated. erefore, they set up a 2D CFD numerical model adding the equivalent mass and energy sources in the first cell layer close to the wall, but they did not model the water film. Wang et al. [2] also used FLUENT code to analyze AP1000 PCCS outside cooling, utilizing the Eulerian wall film (EWF) model simultaneously coupled with the Eulerian multiphase flow model and species transport model. eir numerical results were validated by the experimental data from EFFE [22] and showed that the EWF model could well simulate the falling film evaporation behavior. However, in this study, the water film mass flow rate was very small, and the sensitivity of related parameter was not discussed.
Although there are lots of experimental studies on the external cooling capacity of PCCS currently, the channel gap height in these studies, except for Hu et al. [4], is relatively small. Due to the high cost and limited measurement points of experimental research, it is quite necessary to carry out numerical validation work based on Hu's experiment. Besides, as for the water film/air countercurrent flow and heat transfer process, there are few public 3D CFD numerical investigations which are performed under varied conditions of water/air inlet flow rates, water/air inlet temperatures, heating surface temperatures, air inlet relative humidities, and the tilted angles between the rectangular channel and horizontal plane. erefore, in this study, a 3D rectangular channel geometric model is set up. Two numerical methods are used to calculate the mass and energy exchanged during the water film evaporation heat transfer process using FLUENT. ese numerical results are validated by the experimental data from the study of Hu et al. [4]. In addition, the sensitivity of related parameters is analyzed. Figure 1, the rectangular channel is formed by a bottom steel plate (red), a top glass cover plate (transparent), and two side steel walls (gray). e channel has a length of L (5.0 m), width of W (1.2 m), gap height of H (0.3 m), and tilted angle of θ. In order to simulate the water film flow and heat transfer characteristics in the zones of the arc-shaped dome of PCCS, θ ranges from 0 to 90 deg. Because the vertical side area of the containment is much larger than the area of the dome, θ is usually equal to 90 deg. e water film enters the calculation domain along the negative direction of the z-axis under the action of gravity, and its initial thickness is δ (0.01 m). Meanwhile, the air is above the water film, forming a countercurrent flow with the water. e inlet cross-section area is the same as the outlet, both for the air and water film. In addition, the inlet width of the water film is W so that the water film can cover the entire bottom steel plate. Because there is a certain amount of the water vapor in the dry air, the capability of the dry air to absorb the water vapor is reflected by the air relative humidity φ or the water vapor mass fraction yi. e bottom steel plate has the uniform temperature T w , and the top glass cover plate and two side steel walls are insulated. e velocity inlet is assigned to the boundary of air/water film inlet, where the velocity direction is perpendicular to the cross-section and the profiles of velocity and temperature are uniform. e parameters of the air inlet velocity and temperature are V g,in and T g,in , respectively. In addition, V l,in and T l,in , respectively, refers to the water film velocity and temperature. V l,in is derived from the water film inlet mass flow rate _ m l,in . e air inlet relative humidity is φ. Because the ends of the rectangular channel are open to the ambient, the outlet boundary of air/water film is the pressure outlet. e different parameter values (θ, T w , T l,in , _ m l,in , T g,in , V g,in , and φ) in each research target are shown in Table 1.

Mathematical Model.
In this study, it is very critical to accurately calculate the mass exchange during the process of the water film evaporation heat transfer with the countercurrent air flow. On the one hand, the water film is converted into the water vapor; on the other hand, the generated water vapor will be diffused in the air. erefore, the multiphase flow models must be coupled with the species transportation model. Next, two calculation methods will be introduced separately. at is, the first method is the VOF multiphase flow model with the user-defined functions about the mass and energy source (VOF_UDF), and the second method is the Eulerian multiphase flow model coupled with the Eulerian wall film model (Eulerian_EWF).

e First Method (VOF_UDF).
In this method, the moist air (the mixed gas) which consists of the air and the water vapor is the first phase of the VOF multiphase flow model. e water film (liquid) is the secondary phase. e convection and diffusion behavior of the water vapor in air can be described by where ρ g and yi are the density and the water vapor mass fraction of the moist air, respectively; u is the velocity vector; α g is the gas phase volume fraction; J is the diffusion flux; S m is the mass exchange source term; D is the water vapor mass diffusion coefficient in air; T is the grid cell node temperature; D 0 and T 0 are constants (D 0 � 2.56 × 10 − 5 m 2 /s and T 0 � 298 K); and ∇ is the gradient operator. In (1), the water vapor mass fraction of the moist air yi can be calculated as follows: where M v and M a are the molar mass of the water vapor and air, respectively; φ is the relative humidity of the moist air; p v (T) is the water vapor partial pressure of the moist air at the temperature of T; and p is the total gas phase pressure. e connection between yi and φ can be established through the humidity ratio d of the moist air [23] as shown in the following equation: where p sat is the saturated water vapor partial pressure of the moist air. Besides, ρ g in (1) can be given by [24]  Science and Technology of Nuclear Installations where p 0 is the standard atmospheric pressure of 0.1013 MPa; R is the universal gas constant; T sat is the water vapor saturated temperature; and ρ l (T sat ) is the water film density under the temperature of T sat . e energy exchange source S Φ during the water evaporation process should be given by where c is the latent heat of water vapor and h v (T sat ) and h l (T sat ) are the enthalpy of saturated water vapor and the water film, respectively, which can be determined by the Prakash formula [25].
Since the evaporation of water film and the condensation of water vapor are opposite physical process, S m in (1) and (5) can be calculated according to the formation and growth of droplets on the surface of the heat exchanger analyzed by Zhuang et al. [26]. As shown in Figure 2, it is assumed that the water film thickness on the bottom heated wall is δ 0 at the time t 0 . As the time increases, the thickness of the water film is gradually reduced to δ 1 at time t 1 due to the water film evaporation. e entire water film surface is divided into many small enough cells. e area and volume of each cell are A i and V i , respectively. In addition, the 1-1 surface is the water film surface closest to the phase interface, whose temperature is T l,i . Meanwhile, the 0-0 surface is the moist air surface closest to the phase interface, where the water vapor partial pressure is saturated. us, the water vapor mass fraction of the 0-0 surface is yi (T sat,i ), which can be obtained from (2). e water vapor mass flow of the 0-0 surface equals to the sum of the convective mass flow S c,i and the diffusion mass flow S d,i . Because A i is sufficiently small, the water vapor can only exchange mass in the normal direction of the 0-0 surface, i.e., the vector n. e evaporative water film mass flow of the 1-1 surface S i equals to the water vapor mass flow of the 0-0 surface. erefore, S i can be calculated as follows: where T i,i is the temperature of the grid cell at the gas-liquid phase interface; v g is the 0-0 surface normal velocity of the moist air; yi (T i,i ) is the water vapor mass fraction of the moist air on the 0-0 surface; and (z[yi(T i,i )]/zn) is the gradient of yi (T i,i ) on the 0-0 surface. Moreover, because the net mass increase is zero for the dry air on the 0-0 surface, v g in (6) should also satisfy en, substituting (7) into (6), S i can be obtained as follows: erefore, the total mass exchange S m during the water film evaporation process can be obtained by accumulating all the S i of the phase interface cells as follows: where N is the total number of cells in the gas-liquid phase interface. e continuity equations for the gas phase and the liquid phase of this VOF multiphase model are shown as (10) and (11): where α l is the liquid phase volume fraction. e sum of α l and α g is 1, i.e., α g + α l � 1. e momentum equation of this VOF multiphase model can be written as follows: where I is the unit stress tensor; ρ and μ are the average density and dynamic viscosity of the control volume weighted by the volume fraction, respectively, with ρ � α g ρ g + α l ρ l and μ � α g μ g + α l μ l ; g is the gravity acceleration vector; and f is the volume force source term, which mainly refers to the surface tension. Because the surface tension acts on the phase interface, f can be calculated by the continuous surface force (CSF) model [27] as follows: where σ is the surface tension coefficient between liquid phase and gas phase; k l is the curvature of the interface; n w and t w are the unit normal vector and unit tangent vector of the wall, respectively; and θ w is the contact angle of the wall. e energy equation of this VOF multiphase model can be written as follows: where E and λ are the average total energy and thermal conductivity of the control volume weighted by the volume fraction, respectively, with E � (α g ρ g E g + α l ρ l E l )/(α g ρ g + α l ρ l ) and λ � α g λ g + α l λ l ; C p is the specific heat at constant pressure; σ t is the turbulent Prandtl number for energy that equals to 0.85; μ t is the turbulent eddy viscosity coefficient;

Moist air
Water film Phase interface Science and Technology of Nuclear Installations C μ is a constant; and k and ε are the turbulent kinetic energy and turbulent dissipation rate, respectively.

e Second Method (Eulerian_EWF).
In this method, the air and the water vapor are also considered the mixed gas (the moist air) which is the first phase of the Eulerian multiphase flow model. e water film is the secondary phase.
e species transport equation that describes the diffusion of the water vapor in the air can also be given by (1), and the amount of the water vapor in the moist air is weighted by the water vapor mass fraction yi. e governing equations of continuity, momentum, and energy for all the control volumes are shown as equations (15)- (17): where _ m pq represents the mass exchanged from the p phase to the q phase, with _ m pq � − _ m qp ; u q is the q phase velocity; u pq is the phase slip velocity (when _ m pq > 0, u pq � u p , when _ m pq < 0,u pq � u p ); t q is the q phase stress tensor, satisfying ; r pq is the internal force between phases; f total is the sum of other forces, including the external volume force f q , the lift force f lift,q , the wall lubrication force f wl,q , the virtual mass force f vm,q , and the turbulent diffusion force f td,q ; h q is the q phase specific enthalpy; q q is the heat flux vector, which represents the thermal conduction diffusion term; Φ pq is the energy exchanged between phases, with Φ pq � − Φ qp ; and S q , s q , and Φ s are the source term of mass, force, and energy, respectively. erefore, the mass and energy exchanged during the water film evaporation process can be determined by solving these equations (15)-(17) together with the models of surface tension, turbulent, and species transportation. However, this solution method is relatively complex and difficult. In this research, the water film thickness is very thin, so it is suitable to use the EWF model of the FLUENT to predict the creation and flow of liquid film on the wall surface [2]. To simplify equations (15)-(17), some assumptions are made in the EWF model. e water film is parallel to the wall surface; the internal velocity of the water film is a parabolic distribution; the internal temperature of the water film satisfies the piecewise linear distribution, and the turning point is located at the half thickness of the water film; the numerical condition is steady; and the radiation heat transfer is negligible.
Hence, the water film governing equations of mass, momentum, and energy are shown as equations (18)-(20): where δ is the water film thickness; u l is the water film velocity; _ m s represents the change of the water film mass flow due to the formation, separation, and splashing; u l,n is the normal component of u l ; p L is the total surface force, mainly includes the pressure gradient loss of the air p gas , the gravity component perpendicular to the wall surface p h and the surface tension p σ ; t is the shear stress tensor; ] l is the water film kinematic viscosity; q indicates the force related to the collection or separation of the water film, satisfying q � _ m s u l ; T l,ave is the water film average temperature; T i is the phase interface temperature; T w is the wall temperature; _ q imp represents the heat released by the water film when impinging the wall; _ m phase is the phase change mass flow caused by the evaporation of the water film or the condensation of the water vapor; and c (T i ) represents the latent heat of the water vapor at the temperature T i .
In addition, _ m phase in equation (20) can be calculated as where C phase , C condensation , and C evaporation are constants, with C condensation � C evaporation � 10 10 , and T sat is the water vapor saturated temperature, which depends on the water vapor saturated partial pressure p sat . e expression of p sat is shown as follows: where T is the thermodynamic temperature of the water film surface.

Numerical Solution Setting
2.3.1. Boundary and Initial Condition. As for the setting of the boundary condition, the water vapor mass fraction of the wall is zero. Because the water film flows close to the bottom surface of the channel, the water film temperature only needs to be set on this wall according to the corresponding T l,in in Table 1 when the EWF model is used. Besides, the turbulence parameters of the inlet and outlet boundary for the air and water film are set by turbulence intensity and hydraulic equivalent diameter for each case of Table 1.
As for the setting of the initial condition, the standard initialization calculation is performed from the water film inlet. e entire computational domain is filled with the moist air at a certain temperature, and the water film covers the entire wall with a certain thin thickness. In addition, the governing equations of the EWF model should also be initialized due to the time term of equations (18)∼(20).

Physical Parameter.
ese parameters, such as the density, the dynamic viscosity, the thermal conductivity, and the specific heat, are calculated as the function of the temperature for the air, water vapor, and water film. e density of the moist air is determined by the volume-weighted mixing law of the air and water vapor. In addition, the mass diffusion coefficient of water vapor in the moist air, which is implemented into the solver by means of UDF, is calculated by (1). Besides, the standard state enthalpy of water film and water vapor are equal to zero and latent heat, respectively.

Discrete Formats, Algorithms, and Residuals.
As for the VOF_UDF method, the numerical calculation is based on the finite volume method. e transient pressure-based solver is used. In order to improve the convergence of the solution, the implicit volume force option is turned on. e pressure implicit with splitting of operators (PISO) is utilized as the pressure velocity coupling method. e equations of momentum, energy, turbulent transportation, and the water vapor mass fraction are discretized with the second-order upwind scheme. e time term is adopted by the first-order implicit scheme. In addition, the phase interface is captured with the geometrical reconstruction method. e iterative residuals of the governing equations, such as mass, momentum, turbulent transportation, and water vapor mass fraction, are set to 10 − 4 , and the residual of the energy governing equation is 10 − 6 . e iterative relaxation factors of these governing equations are kept as the default values. Besides, the transient time step is 10 − 4 s.
As for the Eulerian_EWF method, the numerical calculation is also based on the finite volume method. e steady pressure-based solver is used. In order to obtain the stability of the solution, the implicit volume fraction is selected. In the pressure and velocity coupling method, the phase coupled SIMPLE algorithm is utilized and the least square cell-based gradient term is adopted. e governing equations of momentum, energy, turbulence, and water vapor mass fraction are selected by the second-order upwind scheme. e volume fraction equation is used by the QUICK scheme. In addition, the transient iteration mode is used in the EWF model with the time step 0.001 s. e time term of equations (18)∼(20) is discretized with the first-order explicit scheme. e water film governing equations of mass, momentum, and energy are discretized with the first-order upwind scheme. e residuals of all the governing equations are set to 10 − 4 . Besides, the maximum and minimum water film thickness of the EWF model are set to 0.01 and 5.0 × 10 − 7 m, respectively. Figure 3, hexahedral grids are used in the numerical simulation. e grid nodes are refined at the near-wall regions, entrance, and exit of the rectangular channel. In addition, for the VOF_UDF method, the grid nodes are refined in the regions, where y is less than 15 mm, because the VOF multiphase flow model focuses on the phase interface position.

Grid Independent Analysis. As shown in
ese parameters, such as the average water film temperature T l,ave (T l,ave � (T l,in + T l,ave,out )/2), average air temperature T f (T f � (T g,in + T g,ave,out )/2), water film evaporative mass flow _ m evaporation,numerical , and total heat transfer rate of heated wall Φ w,numerical , can be directly obtained in the FLUENT code. us, the expressions of the average heat flux of heated wall q w,ave,numerical , water film Science and Technology of Nuclear Installations evaporative ratio ω numerical , and mass transfer coefficient h m are shown as equations (23)-(25): where A is the area of the bottom wall and ρ i and ρ ∞ are the density of the water vapor at the temperature of T l,ave and T f , respectively. erefore, the grid independent analysis is conducted by comparing the influence of the total grid number or the grid node distribution on these parameters, such as Φ w,numerical , _ m evaporation,numerical , T l,ave , T f , and h m . e specific parameter values of θ, T w , T l,in , _ m l,in , T g,in , V g,in , and φ are shown in Table 1. For the VOF_UDF method, as shown in Table 2, when the first layer cell height of the heated wall is 0.1 mm high (total number of grid cells is 2.70 million, and y + � 0.4), the numerical results will not be affected. For the Euler-ian_EWF method, as shown in Table 3, when the total number of grid cells is 0.96 million (first layer cell height of the heated wall is 0.5 mm, and y + � 1.0), the numerical results will not be affected. In addition, the maximum cell size in both the length and width direction of the channel is 50 mm for these two numerical methods, i.e., VOF_UDF and Eulerian_EWF. Moreover, the calculation results of Φ w,numerical and _ m evaporation,numerical from the VOF_UDF method agrees with those from the Eulerian_EWF method, and the differences of Φ w,numerical and _ m evaporation,numerical between these two methods are 7.25% and 7.82%, respectively.

Experimental Validation of CFD Model
In this section, the numerical results are validated by the experimental data obtained from the WAFT facility [4] which was conducted to study the water film evaporation heat transfer in a large vertical plate with the countercurrent air flow. e schematic of this experimental facility is shown in Figure 4. e facility is made up of the inlet wind passage, the water supply and recovery system, the main rectangular channel test section, the oil circulation and heating system, the hydraulic rotating frame, and the measuring system. e rectangular channel consists of a 5.0 m long and 1.2 m wide oil heating bottom steel plate, a top glass cover plate, and two side insulated walls. e gap height between the steel plate and the glass cover plate is 0.3 m. Water film can cover the entire bottom wall with the gravity from the top water film distribution tank. e countercurrent airflow is driven by the centrifugal fan. Several heat flux sensors are directly attached to the bottom wall surface to measure the heat flux. e evaporative water film mass flow equals to the water film mass flow difference between the supply and recovery system. Several heaters are mounted on the pipes and devices to adjust the temperatures of the heat transfer oil, air, and water. e power of these heaters is measured by the power analyzer (HIOKI-3390). In addition, in order to study the heat transfer characteristics of PCCS in the dome zones of the containment, the channel can be rotated from the horizontal position to the vertical position with the hydraulic rotating frame.
During this test, the heat is transferred through the bottom plate surface by the ways of radiation and water flow convection, while the heat is transferred through the water film surface by the ways of radiation, evaporation, and convection. Because the radiation heat transfer rate is much smaller than the heat transfer rate of evaporation and Φ air � C p,air T f _ m g T g,ave,out − T g,in , _ m g � ρ air T g,ave,out V g,ave,out LW, Φ l,sensible � C p,water T l,ave _ m l,out T l,ave,out − T l,in , where S total is the area of the bottom plate; Φ load is the total heat transfer rate of the bottom plate; Φ air is the air convection heat transfer rate; Φ l,sensible is the water film convection heat transfer rate; Φ l,evaporation is the water film evaporation heat transfer rate; C p,air (T f ) is the air specific heat at the temperature T f ; _ m g is the air mass flow; T g,ave,out is the air outlet average temperature; ρ air (T g,ave,out ) is the air density at the temperature T g,ave,out ; V g,ave,out is the air outlet average velocity; C p,water (T l,ave ) is the water specific heat at the temperature T l,ave ; _ m l,out is the water film recovery mass flow; T l,ave,out is the water film outlet average temperature; t 1 and t 2 are the start and end time for the water film weighing, respectively; and m 1 and m 2 are the readings of weighing instrument, respectively, corresponding to the t 1 and t 2 .
ere are five cases listed in Table 4 to validate these two numerical methods. e results of the average bottom plate surface heat flux q w,ave and the water film evaporation mass flow ratio ω obtained from the experiment and the numerical simulation are shown in Figure 5. e computational time that is needed to simulate these experiments is 300 s for the VOF_UDF numerical method. As can be seen from Figure 5, these numerical results well agree with the corresponding experimental measurements, and the maximum deviations are within 30%. In fact, in only two cases (test nos. 1 and 3 of Table 4), the deviation between experiment and simulation result exceeds 20%. e reason is that the heat flux measured from the experiment is relatively small due to the lower T w or the larger _ m l,in , so the water film evaporation is not obvious. In addition, the deviation of experiment and numerical calculation in the research literature of Ambrosini et al. [20] is also more than 30% when the water film evaporative mass flow is small. Besides, for all the cases of Table 4, the discrepancies of q w,ave and ω from these two numerical methods are relatively small, within ±12%. erefore, these current numerical methods are feasible and can well predict the heat and mass transfer characteristics of the water film evaporation process with the countercurrent air.
Moreover, the VOF_UDF method uses the transient iteration to accurately capture the water film and moist air phase interface. e time step size cannot be larger, so this calculation process is time consuming. e Eulerian_EWF method could not capture the phase interface, but it uses the steady iteration and the calculation efficiency is higher. Due to the larger size of experiment facility as well as the PCCS channel, especially the channel height, the water film thickness of heated plate surface is too small. us, it is not very important to accurately capture the phase interface. In the meanwhile, parameter sensitivity studies generally focus on the physical quantities that reflect the overall heat transfer charateristics of PCCS, such as q w,ave , ω, and δ. erefore, it is wise to select the Eulerian_EWF method in the parameter sensitivity analysis of θ, T w , T l,in , _ m l,in , T g,in , V g,in , and φ. is numerical method not only improves the calculation speed but also reduces the number of grids.

Sensitivity Analysis of Operation Conditions
In order to quantify the influence of operation conditions (listed in Table 1) on the heat and mass characteristics, the calculated values of the average bottom plate surface heat flux q w,ave , the water film evaporation mass flow ratio ω, and the average water film thickness on the bottom plate surface δ are deduced based on the simulation results from the Eulerian_EWF method q w,ave and ω: where δ i is the local value of the water film thickness for every grid node on the bottom plate surface and n is the total number of grid nodes on the bottom plate surface. Figure 6 shows the influence of the tilted angle θ on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, q w,ave and ω increase with the increase of θ. e rising slope of q w,ave and ω are very small when θ is larger than 60 deg, while δ decreases with the increase of θ, and its change is also very small for θ > 60 deg. e reason is that θ mainly affects the component of the gravity acceleration along the direction of water film falling flow (g · sin θ). When θ becomes larger, the water film surface speed will increase and its trend is the same as sin θ. us, the water film thickness will get thinner. In this way, the heat transfer resistance of the water film will also get smaller, and the bottom plate surface heat flux will increase with the increase of θ. In addition, the bottom plate surface heat flux is mainly absorbed by the water film. erefore, _ m evaporation and ω will increase with the increase of θ, while the water film inlet mass flow rate remains the same. Figure 7 shows the influence of the water film inlet mass flow rate _ m l,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, δ increases significantly with the increase of _ m l,in because the increase of _ m l,in will cause the water film velocity to increase simultaneously in both the length and width direction of the channel. In addition, q w,ave first increases with the increase of _ m l,in and then decreases, and it reaches its maximum value when _ m l,in is 0.18 kg/s. ere are two reasons for this phenomenon. On the one hand, the water film evaporation heat transfer plays a major role in the total bottom plate surface heat flux because δ is less than 1.0 mm for all these cases. On the other hand, the convection heat transfer increases due to the increase of turbulent intensity when the range of _ m l,in is 0.09-0.18 kg/s, while it reduces due to the increase of the water film thermal resistance when the range of _ m l,in is 0.18-1.44 kg/s. Moreover, ω decreases with the increase of _ m l,in because _ m l,in rapidly increases while _ m evaporation slightly changes. Figure 8 shows the influence of the water film inlet temperature T l,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, q w,ave first decreases with the increase of T l,in and then increases, and it reaches its minimum value when T l,in is 70°C. e reason is that the heat removed by the water film evaporation is more than that by the water film convection. When T l,in is less than 70°C, the water film convection heat transfer is more important, and the convection heat transfer difference decreases with the increase of T l,in . However, when T l,in is more than 70°C, the water film evaporation heat transfer begins to dominate, and the water vapor partial pressure at the water film surface increases with the increase of T l,in , which will promote the water film evaporation. erefore, when T l,in increases, _ m evaporation increases. As shown in Figure 8, under the same _ m l,in , ω increases with T l,in increasing, and δ reduces with the increase of T l,in . However, the effects of T l,in on ω and δ are relatively weak due to the small change of the water vapor partial pressure of the phase interface.

Influence of the Bottom Plate
Surface Temperature. Figure 9 shows the influence of the bottom plate surface temperature T w on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, q w,ave increases with T w increasing because of the increase of the heat transfer difference between the water film and the bottom plate surface. us, the water film will absorb more heat from the bottom plate surface, which results in the increase of the water film average temperature. At this time, as T w increases, the water vapor diffusion velocity from the water film surface to the mainstream moist air will become larger and _ m evaporation increases. erefore, as shown in Figure 9, ω increases with the increase of T w under the same _ m l,in , but δ decreases. In addition, when T w changes from 75°C to 98°C, the ranges of q w,ave , ω, and δ are 60%, 50%, and 7%, respectively. It can be seen that T w has a smaller effect on δ, but has a larger effect on q w,ave and ω. Besides, the phenomenon of the small change of δ during the water film evaporation process also exists in the study of Wang et al. [2]. Figure 10 shows the influence of the air inlet velocity V g,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. When V g,in increases, q w,ave and ω increases, but δ decreases. ese trends of q w,ave , ω, and δ in Figure 10 are the same as Figure 9. Because the increase of V g,in makes the mainstream moist air refreshing faster, the higher water vapor concentration of the water film surface will be replaced and diluted quickly, so the saturated water vapor will become undersaturated again. erefore, _ m evaporation increases with the increase of V g,in , so does q w,ave . When _ m l,in keeps unchanged, ω becomes larger and δ gets smaller, with V g,in increasing. In addition, the decrease of δ can result in the decrease of water film heat transfer thermal resistance, which can also cause the increase of q w,ave . Figure 11 shows the influence of the air inlet temperature T g,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, the influences of T g,in on q w,ave , ω, and δ are negligible. When T g,in changes from 5°C to 98°C, the ranges of q w,ave , ω, and δ are 0.03%, 2.10%, and 0.18%, respectively.  e reason can be explained from two aspects. On the one hand, the heat transfer process of the water film evaporation with the countercurrent airflow is dominated by the water film evaporation, seen from equation (26), and the heat transfer rate of air convection is much smaller than that of water film evaporation. On the other hand, the increase of T g,in will cause slight deterioration of the air convective heat transfer. It is also shown that the slight increase of evaporation heat flux at low T g,in is counteracted by the decrease of the forced convection heat flux, even leading to a weak decreasing trend of total heat flux [10]. erefore, the dominant factor in this heat and mass transfer process is the V g,in rather than T g,in , comparing Figures 10 and 11. Figure 12 shows the influence of the air inlet relative humidity φ on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, q w,ave and ω decrease with the increase of φ, but δ increases with φ increasing; the ranges of q w,ave , ω, and δ are 12%, 8%, and 1.0%, respectively. Because the water vapor concentration (i.e., mass fraction yi) of the mainstream moist air in the rectangular channel increases with the increase of φ, the larger difference is formed between the water vapor partial pressure of the water film surface and the water vapor partial pressure of the mainstream moist air. erefore, the diffusion intensity of the water vapor is weakened, and _ m evaporation becomes smaller with the increase of φ. However, the water vapor partial pressure difference changes very slightly. In addition, there is a positive correlation between the transferred mass and the transferred heat during this water film evaporation process, so both q w,ave and ω decrease when φ increases. Besides, during this water film evaporation process, only the phase change mass loss exists. us, δ will become larger with the increase of φ. Table of ese Parameters. In order to clearly present the importance of these influencing factors (such as θ, T w , T l,in , _ m l,in , T g,in , V g,in , and φ) on the PCCS heat transfer process, the coefficient of sensitivity (COS) of the factor X is defined as the derivative of dimensionless q w,ave , ω, or δ, to the dimensionless variable of X itself, i.e.,

Influence
, where the superscript * stands for the dimensionless variable, the subscript o represents the based case conditions, and X can be referred to the θ, T w , T l,in , _ m l,in , T g,in , V g,in , or φ. Table 5 shows the relative strength of the COS for the different variables about the based case condition. e values of COS X,q w,ave , COS X,ω , and COS X,δ are obtained from equation (28), and the derivatives of dimensionless q w,ave , ω, or δ to the dimensionless variable X are derived from Figures 6-12. As seen from  Figure 10: Influence of the air inlet velocity on heat and mass transfer of the water film evaporation process. effect on q w,ave and ω, but they show lower effect on δ. e more important influence on δ are parameters of θ and _ m l,in . e parameter of T g,in has the little influence on q w,ave , ω, and δ. erefore, the efficient way to enhance the heat transfer capacity of PCCS depends on the control of V g,in to reduce the flow resistance.

Conclusions
In this paper, the water film evaporation heat transfer characteristics in a large-scale rectangular channel with the countercurrent air flow are investigated mainly by the numerical simulation. Two numerical methods are utilized to calculate the mass and heat transfer when the water film temperature is not saturated. One is the VOF multiphase flow model loading the user-defined mass transfer model driven by the water vapor partial pressure difference and simultaneously coupled with the species transportation model (i.e., VOF_UDF method), the other is the Eulerian multiphase flow model coupled with the Eulerian wall film model and the species transportation model (i.e., Euler-ian_EWF method). Meanwhile, these two numerical methods are validated by some experiments. e following conclusions could be drawn: (1) e mass and heat transfer during the water film evaporation can be predicted by two numerical methods, i.e., VOF_UDF method and Euler-ian_EWF method. e maximum deviations between these obtained numerical results and the corresponding experimental results are 30%, so the accuracy of these numerical results can reach the engineering requirements. Besides, since the VOF_UDF method uses the transient iteration to capture the phase interface, the calculation efficiency of the VOF_UDF method is not as good as the Eulerian_EWF method with the steady iteration. For the heat and mass transfer problems that are not required to capture the phase interface, for example, the research of this paper, the Eulerian_EWF method is better. However, for the heat and mass transfer problems that required to capture the phase interface, for example, the research of the Ref. [25], only the VOF_UDF method can be used.
(2) Among these influencing factors, such as the tilted angle, the water film inlet mass flow rate, the water film inlet temperature, the bottom plate surface temperature, the air inlet velocity, the air inlet temperature, and the air inlet relative humidity, the water film evaporation heat transfer characteristics are not affected by the air inlet temperature, but are greatly affected by the air inlet velocity and the bottom plate surface temperature. e trend of the bottom plate surface average heat flux with these influencing factors is the same as that of the water film evaporation mass flow ratio, but is opposite to the trend of the water film average thickness on the bottom plate surface. Besides, during the water film evaporation process, the water film average thickness on the bottom plate surface changes a little with these influencing factors, except for the tilted angle and the water film inlet mass flow. (3) As for a wide range of PCCS thermal parameter, the numerical method presented in this paper can also accurately predict the heat and mass transfer rate as well as the correct trends of heat flux and water film thickness. In addition, the Eulerian wall film model can predict the flow and development transient process of the water film on the heated wall surface. erefore, an interesting study should be conducted in the future to simulate a full-scaled PCCS model applying the FLUENT during an accident evolution.
Overall, the numerical method in this paper can provide some guidance for the PCCS design in the future.

A:
Total water film coverage area on the heated wall (m 2 ) Ai: Grid cell area of the phase interface (m 2 ) Cp: Specific heat at constant pressure (J/kg/K) C phase : Phase change constant C μ : Constant D: Water vapor mass diffusion coefficient in the dry air (m 2 /s) d: Humidity ratio of the dry air (kg/kg) E: Total energy of the control volume (J/kg) f: Volume force source (N/m 3 ) f total : Sum of other forces, such as the external volume force, lift force, wall lubrication force, virtual mass force, and turbulent diffusion force (N/m 3 ) g: Gravity acceleration vector (m/s 2 ) Volume fraction of the phase c: Latent heat of water vapor (J/kg) δ: Water film thickness (mm) δ: Average water film thickness of the heated wall (mm) ε: Turbulent dissipation rate (m 2 /s 3 ) θ w : Contact angle of the wall (deg) λ: Heat conductivity coefficient (W/m/K) μ: Dynamic viscosity (kg/m/s) μ t : Turbulent eddy viscosity coefficient (kg/ m/s) ρ: Density (kg/m 3 ) ρ ∞ : Density of the water vapor at the average air temperature (kg/m 3 ) ρ i : Density of the water vapor at the average water film temperature (kg/m 3 ) σ: Surface tension coefficient between liquid phase and gas phase (N/m) σ t : Turbulent Prandtl constant ] l : Water film kinematic viscosity (m 2 /s) φ: Relative humidity of the moist air (%) ω numerical : Water film evaporative ratio (%) Φ air : Convection heat transfer rate of the air (W) Φ l,evaporation : Evaporative heat transfer rate of the water film (W) Φ l,sensor : Convection heat transfer rate of the water film (W) Φ load : Measured total heat transfer rate of the bottom plate surface (W) Φ pq : Exchanged energy between phases in equation (17)  Gradient operator (m − 1 ) air: Air g: Gas phase g, ave, out: Parameter average value of the air outlet section l: Liquid or water film l, ave, out: Parameter average value of the water film outlet section max: Maximum value of the variable X min: Minimum value of the variable X o: Based case q: qth phase of multiphase flow model v: Vapor T: Transpose * : Dimensionless.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.