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

Authors: Kashuai Du; Po Hu; Zhen Hu
Journal: Science and Technology of Nuclear Installations
 (2020)

---

## Abstract

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 5 m in length and 1.2 m 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.


---

## Body

## 1. 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. The 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. The 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. They 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. They 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 m2 heating plate. The 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. They 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. They 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. The channel was 5.0 m long, 1.2 m wide, and 0.3 m high. The 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). They 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. There had been a lot of analysis codes that used to analyze the heat removal capacity of PCCS, such as MELCOR [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. They 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. Therefore, 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. Their 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. Therefore, 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. These 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.


## 2. CFD Model

### 2.1. Problem Description

As shown in Figure1, the rectangular channel is formed by a bottom steel plate (red), a top glass cover plate (transparent), and two side steel walls (gray). The 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. The 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. The 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.Figure 1
Schematic of 3D computational domain for the rectangular channel.The bottom steel plate has the uniform temperatureTw, and the top glass cover plate and two side steel walls are insulated. The 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. The parameters of the air inlet velocity and temperature are Vg,in and Tg,in, respectively. In addition, Vl,in and Tl,in, respectively, refers to the water film velocity and temperature. Vl,in is derived from the water film inlet mass flow rate m˙l,in. The 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. The different parameter values (θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, and φ) in each research target are shown in Table 1. The based case parameter conditions are θ = 90 deg, Tw=85°C, Tl,in = 70°C, m˙l,in=0.18kg/s, Tg,in=45°C, Vg,in=6m/s, and φ = 27.5%.Table 1
Numerical simulation plan.
Research goalTilted angle,θ (deg)Temperature of bottom wall,Tw (°C)Temperature of water film inlet,Tl,in (°C)Temperature of air inlet,Tg,in (°C)Velocity of air inlet,Vg,in (m/s)Mass flow rate of water film inlet,m˙l,in (kg/s)Relative humidity of air inlet,φ (%)Grid analysis9084.727044.5920.1836.0Evaluation aboutθ15, 30, 45, 60, 75, 9085704560.1827.5Evaluation aboutTw9075, 80, 85, 90, 95, 98704560.1827.5Evaluation aboutTl,in908520, 30, 40, 50, 60, 70, 804560.1827.5Evaluation aboutm˙l,in9085704560.09, 0.18, 0.36, 0.72, 1.4427.5Evaluation aboutVg,in908570452, 4, 6, 80.1827.5Evaluation aboutTg,in9085705, 25, 45, 65, 8560.1827.5Evaluation aboutφ9085704560.1810.0, 27.5, 50.0, 70.0, 90.0


### 2.2. 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. Therefore, the multiphase flow models must be coupled with the species transportation model. Next, two calculation methods will be introduced separately. That 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).


#### 2.2.1. The 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. The water film (liquid) is the secondary phase.The convection and diffusion behavior of the water vapor in air can be described by(1)∂∂tαgρgyi+∇⋅αgρguyi=−∇⋅J+Sm,J=−αgρgD∇yi,D=D0TT01.5,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; Sm is the mass exchange source term; D is the water vapor mass diffusion coefficient in air; T is the grid cell node temperature; D0 and T0 are constants (D0 = 2.56 × 10−5 m2/s and T0 = 298 K); and ∇ is the gradient operator.In (1), the water vapor mass fraction of the moist air yi can be calculated as follows:(2)yiT=MvpvTMvpvT+Map−pvT,pvT=φp,where Mv and Ma are the molar mass of the water vapor and air, respectively; φ is the relative humidity of the moist air;pvT is the water vapor partial pressure of the moist air at the temperature of T; and p is the total gas phase pressure. The connection between yi and φ can be established through the humidity ratio d of the moist air [23] as shown in the following equation:(3)yiT=d1+d,d=0.622φpsatp−φpsat,where psat is the saturated water vapor partial pressure of the moist air.Besides,ρg in (1) can be given by [24](4)ρg=10−3p0RT1−yiTMa+yiTMv−1,yiT≤yiTsat,10−3p0RT1−yiTMa+yiTsatMv+yiT−yiTsatρlTsat−1,yiT>yiTsat,where p0 is the standard atmospheric pressure of 0.1013 MPa; R is the universal gas constant; Tsat is the water vapor saturated temperature; and ρl (Tsat) is the water film density under the temperature of Tsat.The energy exchange sourceSΦ during the water evaporation process should be given by(5)SΦ=Smγ,γ=hvTsat−hlTsat,where γ is the latent heat of water vapor and hvTsat and hlTsat 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,Sm 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 t0. As the time increases, the thickness of the water film is gradually reduced to δ1 at time t1 due to the water film evaporation. The entire water film surface is divided into many small enough cells. The area and volume of each cell are Ai and Vi, respectively. In addition, the 1-1 surface is the water film surface closest to the phase interface, whose temperature is Tl,i. Meanwhile, the 0-0 surface is the moist air surface closest to the phase interface, where the water vapor partial pressure is saturated. Thus, the water vapor mass fraction of the 0-0 surface is yi (Tsat,i), which can be obtained from (2). The water vapor mass flow of the 0-0 surface equals to the sum of the convective mass flow Sc,i and the diffusion mass flow Sd,i. Because Ai is sufficiently small, the water vapor can only exchange mass in the normal direction of the 0-0 surface, i.e., the vector n. The evaporative water film mass flow of the 1-1 surface Si equals to the water vapor mass flow of the 0-0 surface. Therefore, Si can be calculated as follows:(6)Si=Sc,i+Sd,i=ρgyiTi,ivg−ρgD∂yiTi,i∂nAiVi,where Ti,i is the temperature of the grid cell at the gas-liquid phase interface; vg is the 0-0 surface normal velocity of the moist air; yi (Ti,i) is the water vapor mass fraction of the moist air on the 0-0 surface; and ∂yiTi,i/∂n is the gradient of yi (Ti,i) on the 0-0 surface.Figure 2
Schematic diagram of mass transfer for water film evaporation on the heated wall.Moreover, because the net mass increase is zero for the dry air on the 0-0 surface,vg in (6) should also satisfy(7)ρg1−yiTi,ivg−ρgD∂1−yiTi,i∂n=0,Then, substituting (7) into (6), Si can be obtained as follows:(8)Si=−ρgD1−yiTi,i∂yiTi,i∂nAiVi,Therefore, the total mass exchangeSm during the water film evaporation process can be obtained by accumulating all the Si of the phase interface cells as follows:(9)Sm=∑i=1NSi,where N is the total number of cells in the gas-liquid phase interface.The continuity equations for the gas phase and the liquid phase of this VOF multiphase model are shown as (10) and (11):(10)∂αgρg∂t+∇⋅αgρg⋅u=Sm,(11)∂αlρl∂t+∇⋅αlρl⋅u=−Sm,where αl is the liquid phase volume fraction. The sum of αl and αg is 1, i.e., αg+αl=1.The momentum equation of this VOF multiphase model can be written asfollows:(12)∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇uT−23μ∇⋅u⋅I+ρg+f,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:(13)f=2σρkl∇αlρg+ρl,kl=∇⋅n,n=∇αl∇αl,n=nwcosθw+twsinθw,where σ is the surface tension coefficient between liquid phase and gas phase; kl is the curvature of the interface; nw and tw are the unit normal vector and unit tangent vector of the wall, respectively; and θw is the contact angle of the wall.The energy equation of this VOF multiphase model can be written as follows:(14)∂∂tρE+∇⋅uρE+p=∇⋅λ+Cpμtσt∇T,μt=ρCμk2ε,where E and λ are the average total energy and thermal conductivity of the control volume weighted by the volume fraction, respectively, with E=αgρgEg+αlρlEl/αgρg+αlρl and λ=αgλg+αlλl; Cp 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; Cμ is a constant; and k and ε are the turbulent kinetic energy and turbulent dissipation rate, respectively.


#### 2.2.2. The 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. The water film is the secondary phase.The 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.The governing equations of continuity, momentum, and energy for all the control volumes are shown as equations (15)–(17):(15)∂∂tαqρq+∇⋅αqρquq=m˙pq−m˙qp+Sq,(16)∂∂tαqρquq+∇⋅αqρququq=−αq∇p+tq+αqρqg+rpq+m˙pqupq−m˙qpuqp+ftotal+sq,(17)∂∂tαqρqhq+∇⋅αqρquqhq=tq∇uq−∇⋅qq+Φpq+m˙pqhpq−m˙qphqp+Φs,where m˙pq represents the mass exchanged from the p phase to the q phase, with m˙pq=−m˙qp; uq is the q phase velocity; upq is the phase slip velocity (when m˙pq>0,upq=up, when m˙pq<0,upq=up); tq is the q phase stress tensor, satisfying tq=∇⋅μ∇u+∇uT−2/3μ∇⋅u⋅I; rpq is the internal force between phases; ftotal is the sum of other forces, including the external volume force fq, the lift force flift,q, the wall lubrication force fwl,q, the virtual mass force fvm,q, and the turbulent diffusion force ftd,q; hq is the q phase specific enthalpy; qq is the heat flux vector, which represents the thermal conduction diffusion term; Φpq is the energy exchanged between phases, with Φpq=−Φqp; and Sq, sq, and Φs are the source term of mass, force, and energy, respectively.Therefore, 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. The 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):(18)∂δ∂t+∇⋅δul=m˙sρl,m˙s=αlρlul,n,(19)∂δul∂t+∇⋅δulul=−δ∇pLρl+gδ+32ρlt−3νlulδ+qρl,pL=pgas+ph+pσ,ph=−ρδn⋅g,pσ=−σ∇⋅∇δ,(20)∂∂tδTl,ave+∇⋅ulδTl,ave=1ρCp2λlTi+Twδ−2Tl,aveδ+q˙imp+m˙phaseγTi,where δ is the water film thickness; ul is the water film velocity; m˙s represents the change of the water film mass flow due to the formation, separation, and splashing; ul,n is the normal component of ul; pL is the total surface force, mainly includes the pressure gradient loss of the air pgas, the gravity component perpendicular to the wall surface ph 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˙sul; Tl,ave is the water film average temperature; Ti is the phase interface temperature; Tw 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 γ (Ti) represents the latent heat of the water vapor at the temperature Ti.In addition,m˙phase in equation (20) can be calculated as [2, 19](21)m˙phase=ρgD/δρgD/δ+CphaseCphaseyiTsat−yiTi,Cphase=Ccondensation,yiTsat>yiTi,10−4Cevaporationδ,yiTsat≤yiTi,where Cphase, Ccondensation, and Cevaporation are constants, with Ccondensation=Cevaporation=1010, and Tsat is the water vapor saturated temperature, which depends on the water vapor saturated partial pressure psat.The expression ofpsat is shown as follows:(22)log10psat=−2.1794+0.02953T−273.17−9.1837×10−5T−273.172+1.4454×10−7T−273.173,where T is the thermodynamic temperature of the water film surface.


### 2.3. 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 correspondingTl,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. The 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).


#### 2.3.2. Physical Parameter

These 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. The 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.


#### 2.3.3. Discrete Formats, Algorithms, and Residuals

As for the VOF_UDF method, the numerical calculation is based on the finite volume method. The transient pressure-based solver is used. In order to improve the convergence of the solution, the implicit volume force option is turned on. The pressure implicit with splitting of operators (PISO) is utilized as the pressure velocity coupling method. The equations of momentum, energy, turbulent transportation, and the water vapor mass fraction are discretized with the second-order upwind scheme. The time term is adopted by the first-order implicit scheme. In addition, the phase interface is captured with the geometrical reconstruction method. The 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. The 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. The 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. The governing equations of momentum, energy, turbulence, and water vapor mass fraction are selected by the second-order upwind scheme. The 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. The time term of equations (18)∼(20) is discretized with the first-order explicit scheme. The water film governing equations of mass, momentum, and energy are discretized with the first-order upwind scheme. The 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.


### 2.4. Grid Independent Analysis

As shown in Figure3, hexahedral grids are used in the numerical simulation. The 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.Figure 3
Grid detail for the numerical simulation.These parameters, such as the average water film temperatureTl,aveTl,ave=Tl,in+Tl,ave,out/2, average air temperature TfTf=Tg,in+Tg,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. Thus, the expressions of the average heat flux of heated wall qw,ave,numerical, water film evaporative ratio ωnumerical, and mass transfer coefficient hm are shown as equations (23)–(25):(23)qw,ave,numerical=Φw,numericalA,(24)ωnumerical=m˙evaporation,numericalm˙l,in,(25)hm=m˙evaporation,numericalAρi−ρ∞,where A is the area of the bottom wall and ρi and ρ∞ are the density of the water vapor at the temperature of Tl,ave and Tf, respectively.Therefore, 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, Tl,ave, Tf, and hm. The specific parameter values of θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,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 Eulerian_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.Table 2
Grid independence analysis for the VOF_ UDF method.
Total number of grid cells1,800,0002,700,0004,000,0006,000,000Height of the 1st layer cell near the heated wall (mm)0.2500.1000.0400.016Φw,numerical (kW)72.4385.2285.4686.34m˙evaporation,numerical, (kg/s)0.027680.033140.033250.03364hm (m/s)0.033760.035960.036040.03621Table 3
Grid independence analysis for the Eulerian_ EWF method.
Total number of grid cells640,000960,0001,440,0002,160,000Φw,numerical (kW)102.9091.4090.9790.77m˙evaporation,numerical, (kg/s)0.040700.035730.035620.03553Tl,ave (K)348.98347.21347.34347.31Tf (K)322.40321.02321.79321.77


## 2.1. Problem Description

As shown in Figure1, the rectangular channel is formed by a bottom steel plate (red), a top glass cover plate (transparent), and two side steel walls (gray). The 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. The 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. The 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.Figure 1
Schematic of 3D computational domain for the rectangular channel.The bottom steel plate has the uniform temperatureTw, and the top glass cover plate and two side steel walls are insulated. The 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. The parameters of the air inlet velocity and temperature are Vg,in and Tg,in, respectively. In addition, Vl,in and Tl,in, respectively, refers to the water film velocity and temperature. Vl,in is derived from the water film inlet mass flow rate m˙l,in. The 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. The different parameter values (θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, and φ) in each research target are shown in Table 1. The based case parameter conditions are θ = 90 deg, Tw=85°C, Tl,in = 70°C, m˙l,in=0.18kg/s, Tg,in=45°C, Vg,in=6m/s, and φ = 27.5%.Table 1
Numerical simulation plan.
Research goalTilted angle,θ (deg)Temperature of bottom wall,Tw (°C)Temperature of water film inlet,Tl,in (°C)Temperature of air inlet,Tg,in (°C)Velocity of air inlet,Vg,in (m/s)Mass flow rate of water film inlet,m˙l,in (kg/s)Relative humidity of air inlet,φ (%)Grid analysis9084.727044.5920.1836.0Evaluation aboutθ15, 30, 45, 60, 75, 9085704560.1827.5Evaluation aboutTw9075, 80, 85, 90, 95, 98704560.1827.5Evaluation aboutTl,in908520, 30, 40, 50, 60, 70, 804560.1827.5Evaluation aboutm˙l,in9085704560.09, 0.18, 0.36, 0.72, 1.4427.5Evaluation aboutVg,in908570452, 4, 6, 80.1827.5Evaluation aboutTg,in9085705, 25, 45, 65, 8560.1827.5Evaluation aboutφ9085704560.1810.0, 27.5, 50.0, 70.0, 90.0


## 2.2. 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. Therefore, the multiphase flow models must be coupled with the species transportation model. Next, two calculation methods will be introduced separately. That 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).


### 2.2.1. The 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. The water film (liquid) is the secondary phase.The convection and diffusion behavior of the water vapor in air can be described by(1)∂∂tαgρgyi+∇⋅αgρguyi=−∇⋅J+Sm,J=−αgρgD∇yi,D=D0TT01.5,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; Sm is the mass exchange source term; D is the water vapor mass diffusion coefficient in air; T is the grid cell node temperature; D0 and T0 are constants (D0 = 2.56 × 10−5 m2/s and T0 = 298 K); and ∇ is the gradient operator.In (1), the water vapor mass fraction of the moist air yi can be calculated as follows:(2)yiT=MvpvTMvpvT+Map−pvT,pvT=φp,where Mv and Ma are the molar mass of the water vapor and air, respectively; φ is the relative humidity of the moist air;pvT is the water vapor partial pressure of the moist air at the temperature of T; and p is the total gas phase pressure. The connection between yi and φ can be established through the humidity ratio d of the moist air [23] as shown in the following equation:(3)yiT=d1+d,d=0.622φpsatp−φpsat,where psat is the saturated water vapor partial pressure of the moist air.Besides,ρg in (1) can be given by [24](4)ρg=10−3p0RT1−yiTMa+yiTMv−1,yiT≤yiTsat,10−3p0RT1−yiTMa+yiTsatMv+yiT−yiTsatρlTsat−1,yiT>yiTsat,where p0 is the standard atmospheric pressure of 0.1013 MPa; R is the universal gas constant; Tsat is the water vapor saturated temperature; and ρl (Tsat) is the water film density under the temperature of Tsat.The energy exchange sourceSΦ during the water evaporation process should be given by(5)SΦ=Smγ,γ=hvTsat−hlTsat,where γ is the latent heat of water vapor and hvTsat and hlTsat 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,Sm 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 t0. As the time increases, the thickness of the water film is gradually reduced to δ1 at time t1 due to the water film evaporation. The entire water film surface is divided into many small enough cells. The area and volume of each cell are Ai and Vi, respectively. In addition, the 1-1 surface is the water film surface closest to the phase interface, whose temperature is Tl,i. Meanwhile, the 0-0 surface is the moist air surface closest to the phase interface, where the water vapor partial pressure is saturated. Thus, the water vapor mass fraction of the 0-0 surface is yi (Tsat,i), which can be obtained from (2). The water vapor mass flow of the 0-0 surface equals to the sum of the convective mass flow Sc,i and the diffusion mass flow Sd,i. Because Ai is sufficiently small, the water vapor can only exchange mass in the normal direction of the 0-0 surface, i.e., the vector n. The evaporative water film mass flow of the 1-1 surface Si equals to the water vapor mass flow of the 0-0 surface. Therefore, Si can be calculated as follows:(6)Si=Sc,i+Sd,i=ρgyiTi,ivg−ρgD∂yiTi,i∂nAiVi,where Ti,i is the temperature of the grid cell at the gas-liquid phase interface; vg is the 0-0 surface normal velocity of the moist air; yi (Ti,i) is the water vapor mass fraction of the moist air on the 0-0 surface; and ∂yiTi,i/∂n is the gradient of yi (Ti,i) on the 0-0 surface.Figure 2
Schematic diagram of mass transfer for water film evaporation on the heated wall.Moreover, because the net mass increase is zero for the dry air on the 0-0 surface,vg in (6) should also satisfy(7)ρg1−yiTi,ivg−ρgD∂1−yiTi,i∂n=0,Then, substituting (7) into (6), Si can be obtained as follows:(8)Si=−ρgD1−yiTi,i∂yiTi,i∂nAiVi,Therefore, the total mass exchangeSm during the water film evaporation process can be obtained by accumulating all the Si of the phase interface cells as follows:(9)Sm=∑i=1NSi,where N is the total number of cells in the gas-liquid phase interface.The continuity equations for the gas phase and the liquid phase of this VOF multiphase model are shown as (10) and (11):(10)∂αgρg∂t+∇⋅αgρg⋅u=Sm,(11)∂αlρl∂t+∇⋅αlρl⋅u=−Sm,where αl is the liquid phase volume fraction. The sum of αl and αg is 1, i.e., αg+αl=1.The momentum equation of this VOF multiphase model can be written asfollows:(12)∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇uT−23μ∇⋅u⋅I+ρg+f,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:(13)f=2σρkl∇αlρg+ρl,kl=∇⋅n,n=∇αl∇αl,n=nwcosθw+twsinθw,where σ is the surface tension coefficient between liquid phase and gas phase; kl is the curvature of the interface; nw and tw are the unit normal vector and unit tangent vector of the wall, respectively; and θw is the contact angle of the wall.The energy equation of this VOF multiphase model can be written as follows:(14)∂∂tρE+∇⋅uρE+p=∇⋅λ+Cpμtσt∇T,μt=ρCμk2ε,where E and λ are the average total energy and thermal conductivity of the control volume weighted by the volume fraction, respectively, with E=αgρgEg+αlρlEl/αgρg+αlρl and λ=αgλg+αlλl; Cp 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; Cμ is a constant; and k and ε are the turbulent kinetic energy and turbulent dissipation rate, respectively.


### 2.2.2. The 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. The water film is the secondary phase.The 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.The governing equations of continuity, momentum, and energy for all the control volumes are shown as equations (15)–(17):(15)∂∂tαqρq+∇⋅αqρquq=m˙pq−m˙qp+Sq,(16)∂∂tαqρquq+∇⋅αqρququq=−αq∇p+tq+αqρqg+rpq+m˙pqupq−m˙qpuqp+ftotal+sq,(17)∂∂tαqρqhq+∇⋅αqρquqhq=tq∇uq−∇⋅qq+Φpq+m˙pqhpq−m˙qphqp+Φs,where m˙pq represents the mass exchanged from the p phase to the q phase, with m˙pq=−m˙qp; uq is the q phase velocity; upq is the phase slip velocity (when m˙pq>0,upq=up, when m˙pq<0,upq=up); tq is the q phase stress tensor, satisfying tq=∇⋅μ∇u+∇uT−2/3μ∇⋅u⋅I; rpq is the internal force between phases; ftotal is the sum of other forces, including the external volume force fq, the lift force flift,q, the wall lubrication force fwl,q, the virtual mass force fvm,q, and the turbulent diffusion force ftd,q; hq is the q phase specific enthalpy; qq is the heat flux vector, which represents the thermal conduction diffusion term; Φpq is the energy exchanged between phases, with Φpq=−Φqp; and Sq, sq, and Φs are the source term of mass, force, and energy, respectively.Therefore, 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. The 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):(18)∂δ∂t+∇⋅δul=m˙sρl,m˙s=αlρlul,n,(19)∂δul∂t+∇⋅δulul=−δ∇pLρl+gδ+32ρlt−3νlulδ+qρl,pL=pgas+ph+pσ,ph=−ρδn⋅g,pσ=−σ∇⋅∇δ,(20)∂∂tδTl,ave+∇⋅ulδTl,ave=1ρCp2λlTi+Twδ−2Tl,aveδ+q˙imp+m˙phaseγTi,where δ is the water film thickness; ul is the water film velocity; m˙s represents the change of the water film mass flow due to the formation, separation, and splashing; ul,n is the normal component of ul; pL is the total surface force, mainly includes the pressure gradient loss of the air pgas, the gravity component perpendicular to the wall surface ph 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˙sul; Tl,ave is the water film average temperature; Ti is the phase interface temperature; Tw 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 γ (Ti) represents the latent heat of the water vapor at the temperature Ti.In addition,m˙phase in equation (20) can be calculated as [2, 19](21)m˙phase=ρgD/δρgD/δ+CphaseCphaseyiTsat−yiTi,Cphase=Ccondensation,yiTsat>yiTi,10−4Cevaporationδ,yiTsat≤yiTi,where Cphase, Ccondensation, and Cevaporation are constants, with Ccondensation=Cevaporation=1010, and Tsat is the water vapor saturated temperature, which depends on the water vapor saturated partial pressure psat.The expression ofpsat is shown as follows:(22)log10psat=−2.1794+0.02953T−273.17−9.1837×10−5T−273.172+1.4454×10−7T−273.173,where T is the thermodynamic temperature of the water film surface.


## 2.2.1. The 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. The water film (liquid) is the secondary phase.The convection and diffusion behavior of the water vapor in air can be described by(1)∂∂tαgρgyi+∇⋅αgρguyi=−∇⋅J+Sm,J=−αgρgD∇yi,D=D0TT01.5,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; Sm is the mass exchange source term; D is the water vapor mass diffusion coefficient in air; T is the grid cell node temperature; D0 and T0 are constants (D0 = 2.56 × 10−5 m2/s and T0 = 298 K); and ∇ is the gradient operator.In (1), the water vapor mass fraction of the moist air yi can be calculated as follows:(2)yiT=MvpvTMvpvT+Map−pvT,pvT=φp,where Mv and Ma are the molar mass of the water vapor and air, respectively; φ is the relative humidity of the moist air;pvT is the water vapor partial pressure of the moist air at the temperature of T; and p is the total gas phase pressure. The connection between yi and φ can be established through the humidity ratio d of the moist air [23] as shown in the following equation:(3)yiT=d1+d,d=0.622φpsatp−φpsat,where psat is the saturated water vapor partial pressure of the moist air.Besides,ρg in (1) can be given by [24](4)ρg=10−3p0RT1−yiTMa+yiTMv−1,yiT≤yiTsat,10−3p0RT1−yiTMa+yiTsatMv+yiT−yiTsatρlTsat−1,yiT>yiTsat,where p0 is the standard atmospheric pressure of 0.1013 MPa; R is the universal gas constant; Tsat is the water vapor saturated temperature; and ρl (Tsat) is the water film density under the temperature of Tsat.The energy exchange sourceSΦ during the water evaporation process should be given by(5)SΦ=Smγ,γ=hvTsat−hlTsat,where γ is the latent heat of water vapor and hvTsat and hlTsat 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,Sm 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 t0. As the time increases, the thickness of the water film is gradually reduced to δ1 at time t1 due to the water film evaporation. The entire water film surface is divided into many small enough cells. The area and volume of each cell are Ai and Vi, respectively. In addition, the 1-1 surface is the water film surface closest to the phase interface, whose temperature is Tl,i. Meanwhile, the 0-0 surface is the moist air surface closest to the phase interface, where the water vapor partial pressure is saturated. Thus, the water vapor mass fraction of the 0-0 surface is yi (Tsat,i), which can be obtained from (2). The water vapor mass flow of the 0-0 surface equals to the sum of the convective mass flow Sc,i and the diffusion mass flow Sd,i. Because Ai is sufficiently small, the water vapor can only exchange mass in the normal direction of the 0-0 surface, i.e., the vector n. The evaporative water film mass flow of the 1-1 surface Si equals to the water vapor mass flow of the 0-0 surface. Therefore, Si can be calculated as follows:(6)Si=Sc,i+Sd,i=ρgyiTi,ivg−ρgD∂yiTi,i∂nAiVi,where Ti,i is the temperature of the grid cell at the gas-liquid phase interface; vg is the 0-0 surface normal velocity of the moist air; yi (Ti,i) is the water vapor mass fraction of the moist air on the 0-0 surface; and ∂yiTi,i/∂n is the gradient of yi (Ti,i) on the 0-0 surface.Figure 2
Schematic diagram of mass transfer for water film evaporation on the heated wall.Moreover, because the net mass increase is zero for the dry air on the 0-0 surface,vg in (6) should also satisfy(7)ρg1−yiTi,ivg−ρgD∂1−yiTi,i∂n=0,Then, substituting (7) into (6), Si can be obtained as follows:(8)Si=−ρgD1−yiTi,i∂yiTi,i∂nAiVi,Therefore, the total mass exchangeSm during the water film evaporation process can be obtained by accumulating all the Si of the phase interface cells as follows:(9)Sm=∑i=1NSi,where N is the total number of cells in the gas-liquid phase interface.The continuity equations for the gas phase and the liquid phase of this VOF multiphase model are shown as (10) and (11):(10)∂αgρg∂t+∇⋅αgρg⋅u=Sm,(11)∂αlρl∂t+∇⋅αlρl⋅u=−Sm,where αl is the liquid phase volume fraction. The sum of αl and αg is 1, i.e., αg+αl=1.The momentum equation of this VOF multiphase model can be written asfollows:(12)∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇uT−23μ∇⋅u⋅I+ρg+f,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:(13)f=2σρkl∇αlρg+ρl,kl=∇⋅n,n=∇αl∇αl,n=nwcosθw+twsinθw,where σ is the surface tension coefficient between liquid phase and gas phase; kl is the curvature of the interface; nw and tw are the unit normal vector and unit tangent vector of the wall, respectively; and θw is the contact angle of the wall.The energy equation of this VOF multiphase model can be written as follows:(14)∂∂tρE+∇⋅uρE+p=∇⋅λ+Cpμtσt∇T,μt=ρCμk2ε,where E and λ are the average total energy and thermal conductivity of the control volume weighted by the volume fraction, respectively, with E=αgρgEg+αlρlEl/αgρg+αlρl and λ=αgλg+αlλl; Cp 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; Cμ is a constant; and k and ε are the turbulent kinetic energy and turbulent dissipation rate, respectively.


## 2.2.2. The 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. The water film is the secondary phase.The 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.The governing equations of continuity, momentum, and energy for all the control volumes are shown as equations (15)–(17):(15)∂∂tαqρq+∇⋅αqρquq=m˙pq−m˙qp+Sq,(16)∂∂tαqρquq+∇⋅αqρququq=−αq∇p+tq+αqρqg+rpq+m˙pqupq−m˙qpuqp+ftotal+sq,(17)∂∂tαqρqhq+∇⋅αqρquqhq=tq∇uq−∇⋅qq+Φpq+m˙pqhpq−m˙qphqp+Φs,where m˙pq represents the mass exchanged from the p phase to the q phase, with m˙pq=−m˙qp; uq is the q phase velocity; upq is the phase slip velocity (when m˙pq>0,upq=up, when m˙pq<0,upq=up); tq is the q phase stress tensor, satisfying tq=∇⋅μ∇u+∇uT−2/3μ∇⋅u⋅I; rpq is the internal force between phases; ftotal is the sum of other forces, including the external volume force fq, the lift force flift,q, the wall lubrication force fwl,q, the virtual mass force fvm,q, and the turbulent diffusion force ftd,q; hq is the q phase specific enthalpy; qq is the heat flux vector, which represents the thermal conduction diffusion term; Φpq is the energy exchanged between phases, with Φpq=−Φqp; and Sq, sq, and Φs are the source term of mass, force, and energy, respectively.Therefore, 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. The 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):(18)∂δ∂t+∇⋅δul=m˙sρl,m˙s=αlρlul,n,(19)∂δul∂t+∇⋅δulul=−δ∇pLρl+gδ+32ρlt−3νlulδ+qρl,pL=pgas+ph+pσ,ph=−ρδn⋅g,pσ=−σ∇⋅∇δ,(20)∂∂tδTl,ave+∇⋅ulδTl,ave=1ρCp2λlTi+Twδ−2Tl,aveδ+q˙imp+m˙phaseγTi,where δ is the water film thickness; ul is the water film velocity; m˙s represents the change of the water film mass flow due to the formation, separation, and splashing; ul,n is the normal component of ul; pL is the total surface force, mainly includes the pressure gradient loss of the air pgas, the gravity component perpendicular to the wall surface ph 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˙sul; Tl,ave is the water film average temperature; Ti is the phase interface temperature; Tw 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 γ (Ti) represents the latent heat of the water vapor at the temperature Ti.In addition,m˙phase in equation (20) can be calculated as [2, 19](21)m˙phase=ρgD/δρgD/δ+CphaseCphaseyiTsat−yiTi,Cphase=Ccondensation,yiTsat>yiTi,10−4Cevaporationδ,yiTsat≤yiTi,where Cphase, Ccondensation, and Cevaporation are constants, with Ccondensation=Cevaporation=1010, and Tsat is the water vapor saturated temperature, which depends on the water vapor saturated partial pressure psat.The expression ofpsat is shown as follows:(22)log10psat=−2.1794+0.02953T−273.17−9.1837×10−5T−273.172+1.4454×10−7T−273.173,where T is the thermodynamic temperature of the water film surface.


## 2.3. 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 correspondingTl,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. The 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).


### 2.3.2. Physical Parameter

These 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. The 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.


### 2.3.3. Discrete Formats, Algorithms, and Residuals

As for the VOF_UDF method, the numerical calculation is based on the finite volume method. The transient pressure-based solver is used. In order to improve the convergence of the solution, the implicit volume force option is turned on. The pressure implicit with splitting of operators (PISO) is utilized as the pressure velocity coupling method. The equations of momentum, energy, turbulent transportation, and the water vapor mass fraction are discretized with the second-order upwind scheme. The time term is adopted by the first-order implicit scheme. In addition, the phase interface is captured with the geometrical reconstruction method. The 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. The 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. The 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. The governing equations of momentum, energy, turbulence, and water vapor mass fraction are selected by the second-order upwind scheme. The 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. The time term of equations (18)∼(20) is discretized with the first-order explicit scheme. The water film governing equations of mass, momentum, and energy are discretized with the first-order upwind scheme. The 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.


## 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 correspondingTl,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. The 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).


## 2.3.2. Physical Parameter

These 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. The 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.


## 2.3.3. Discrete Formats, Algorithms, and Residuals

As for the VOF_UDF method, the numerical calculation is based on the finite volume method. The transient pressure-based solver is used. In order to improve the convergence of the solution, the implicit volume force option is turned on. The pressure implicit with splitting of operators (PISO) is utilized as the pressure velocity coupling method. The equations of momentum, energy, turbulent transportation, and the water vapor mass fraction are discretized with the second-order upwind scheme. The time term is adopted by the first-order implicit scheme. In addition, the phase interface is captured with the geometrical reconstruction method. The 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. The 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. The 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. The governing equations of momentum, energy, turbulence, and water vapor mass fraction are selected by the second-order upwind scheme. The 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. The time term of equations (18)∼(20) is discretized with the first-order explicit scheme. The water film governing equations of mass, momentum, and energy are discretized with the first-order upwind scheme. The 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.


## 2.4. Grid Independent Analysis

As shown in Figure3, hexahedral grids are used in the numerical simulation. The 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.Figure 3
Grid detail for the numerical simulation.These parameters, such as the average water film temperatureTl,aveTl,ave=Tl,in+Tl,ave,out/2, average air temperature TfTf=Tg,in+Tg,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. Thus, the expressions of the average heat flux of heated wall qw,ave,numerical, water film evaporative ratio ωnumerical, and mass transfer coefficient hm are shown as equations (23)–(25):(23)qw,ave,numerical=Φw,numericalA,(24)ωnumerical=m˙evaporation,numericalm˙l,in,(25)hm=m˙evaporation,numericalAρi−ρ∞,where A is the area of the bottom wall and ρi and ρ∞ are the density of the water vapor at the temperature of Tl,ave and Tf, respectively.Therefore, 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, Tl,ave, Tf, and hm. The specific parameter values of θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,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 Eulerian_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.Table 2
Grid independence analysis for the VOF_ UDF method.
Total number of grid cells1,800,0002,700,0004,000,0006,000,000Height of the 1st layer cell near the heated wall (mm)0.2500.1000.0400.016Φw,numerical (kW)72.4385.2285.4686.34m˙evaporation,numerical, (kg/s)0.027680.033140.033250.03364hm (m/s)0.033760.035960.036040.03621Table 3
Grid independence analysis for the Eulerian_ EWF method.
Total number of grid cells640,000960,0001,440,0002,160,000Φw,numerical (kW)102.9091.4090.9790.77m˙evaporation,numerical, (kg/s)0.040700.035730.035620.03553Tl,ave (K)348.98347.21347.34347.31Tf (K)322.40321.02321.79321.77


## 3. 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. The schematic of this experimental facility is shown in Figure 4. The 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. The 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. The 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. The 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. The 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. The 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.Figure 4
Schematic of the experimental facility.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 convection, the average heat flux of the bottom plate surfaceqload can be calculated by(26)qload=ΦloadStotal,Φload=Φair+Φl,sensible+Φl,evaporation,Stotal=LW,Φair=Cp,airTfm˙gTg,ave,out−Tg,in,m˙g=ρairTg,ave,outVg,ave,outLW,Φl,sensible=Cp,waterTl,avem˙l,outTl,ave,out−Tl,in,m˙l,out=m2−m1t2−t1,Φl,evaporation=m˙evaporationhv−hl,m˙evaporation=m˙l,in−m˙l,out,where Stotal 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; Cp,air (Tf) is the air specific heat at the temperature Tf; m˙g is the air mass flow; Tg,ave,out is the air outlet average temperature; ρair (Tg,ave,out) is the air density at the temperature Tg,ave,out; Vg,ave,out is the air outlet average velocity; Cp,water (Tl,ave) is the water specific heat at the temperature Tl,ave; m˙l,out is the water film recovery mass flow; Tl,ave,out is the water film outlet average temperature; t1 and t2 are the start and end time for the water film weighing, respectively; and m1 and m2 are the readings of weighing instrument, respectively, corresponding to the t1 and t2.There are five cases listed in Table4 to validate these two numerical methods. The results of the average bottom plate surface heat flux qw,ave and the water film evaporation mass flow ratio ω obtained from the experiment and the numerical simulation are shown in Figure 5. The 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%. The reason is that the heat flux measured from the experiment is relatively small due to the lower Tw 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 qw,ave and ω from these two numerical methods are relatively small, within ±12%. Therefore, 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.Table 4
Thermodynamic parameters of the analyzed test.
Test no.Tilted angle,θ (deg)Temperature of bottom wall,Tw (°C)Temperature of water film inlet,Tl,in (°C)Temperature of air inlet,Tg,in (°C)Velocity of air inlet,Vg,in(m/s)Mass flow of water film inlet,m˙l,in(kg/s)Relative humidity of air inlet,φ (%)19060504540.1833.62857020.4268.03857060.4230.84857040.1833.85857060.1827.5Figure 5
Comparison of the average bottom plate surface heat flux and the water film evaporation mass flow ratio obtained from the experiment and numerical simulation.
(a)(b)Moreover, the VOF_UDF method uses the transient iteration to accurately capture the water film and moist air phase interface. The time step size cannot be larger, so this calculation process is time consuming. The 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. Thus, 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 asqw,ave, ω, and δ¯. Therefore, it is wise to select the Eulerian_EWF method in the parameter sensitivity analysis of θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, and φ. This numerical method not only improves the calculation speed but also reduces the number of grids.


## 4. Sensitivity Analysis of Operation Conditions

In order to quantify the influence of operation conditions (listed in Table1) on the heat and mass characteristics, the calculated values of the average bottom plate surface heat flux qw,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 qw,aveandω:(27)δ¯=1n∑i=1nδi,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.


### 4.1. Influence of the Tilted Angle

Figure6 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, qw,ave and ω increase with the increase of θ. The rising slope of qw,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. The 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 θ. Thus, 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. Therefore, m˙evaporation and ω will increase with the increase of θ, while the water film inlet mass flow rate remains the same.Figure 6
Influence of the tilted angle on heat and mass transfer of the water film evaporation process.


### 4.2. Influence of the Water Inlet Mass Flow Rate

Figure7 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, qw,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. There 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 7
Influence of the water film inlet mass flow rate on heat and mass transfer of the water film evaporation process.


### 4.3. Influence of the Water Inlet Temperature

Figure8 shows the influence of the water film inlet temperature Tl,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, qw,ave first decreases with the increase of Tl,in and then increases, and it reaches its minimum value when Tl,in is 70°C. The reason is that the heat removed by the water film evaporation is more than that by the water film convection. When Tl,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 Tl,in. However, when Tl,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 Tl,in, which will promote the water film evaporation. Therefore, when Tl,in increases, m˙evaporation increases. As shown in Figure 8, under the same m˙l,in, ω increases with Tl,in increasing, and δ¯ reduces with the increase of Tl,in. However, the effects of Tl,in on ω and δ¯ are relatively weak due to the small change of the water vapor partial pressure of the phase interface.Figure 8
Influence of the water film inlet temperature on heat and mass transfer of the water film evaporation process.


### 4.4. Influence of the Bottom Plate Surface Temperature

Figure9 shows the influence of the bottom plate surface temperature Tw 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, qw,ave increases with Tw increasing because of the increase of the heat transfer difference between the water film and the bottom plate surface. Thus, 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 Tw increases, the water vapor diffusion velocity from the water film surface to the mainstream moist air will become larger and m˙evaporation increases. Therefore, as shown in Figure 9, ω increases with the increase of Tw under the same m˙l,in, but δ¯ decreases. In addition, when Tw changes from 75°C to 98°C, the ranges of qw,ave, ω, and δ¯ are 60%, 50%, and 7%, respectively. It can be seen that Tw has a smaller effect on δ¯, but has a larger effect on qw,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 9
Influence of the bottom plate surface temperature on heat and mass transfer of the water film evaporation process.


### 4.5. Influence of the Air Inlet Velocity

Figure10 shows the influence of the air inlet velocity Vg,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. When Vg,in increases, qw,ave and ω increases, but δ¯ decreases. These trends of qw,ave, ω, and δ¯ in Figure 10 are the same as Figure 9. Because the increase of Vg,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. Therefore, m˙evaporation increases with the increase of Vg,in, so does qw,ave. When m˙l,in keeps unchanged, ω becomes larger and δ¯ gets smaller, with Vg,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 qw,ave.Figure 10
Influence of the air inlet velocity on heat and mass transfer of the water film evaporation process.


### 4.6. Influence of the Air Inlet Temperature

Figure11 shows the influence of the air inlet temperature Tg,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 Tg,in on qw,ave, ω, and δ¯ are negligible. When Tg,in changes from 5°C to 98°C, the ranges of qw,ave, ω, and δ¯ are 0.03%, 2.10%, and 0.18%, respectively. The 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 Tg,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 Tg,in is counteracted by the decrease of the forced convection heat flux, even leading to a weak decreasing trend of total heat flux [10]. Therefore, the dominant factor in this heat and mass transfer process is the Vg,in rather than Tg,in, comparing Figures 10 and 11.Figure 11
Influence of the air inlet temperature on heat and mass transfer of the water film evaporation process.


### 4.7. Influence of the Air Inlet Relative Humidity

Figure12 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, qw,ave and ω decrease with the increase of φ, but δ¯ increases with φ increasing; the ranges of qw,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. Therefore, 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 qw,ave and ω decrease when φ increases. Besides, during this water film evaporation process, only the phase change mass loss exists. Thus, δ¯ will become larger with the increase of φ.Figure 12
Influence of the air inlet relative humidity on heat and mass transfer of the water film evaporation process.


### 4.8. Influence Table of These Parameters

In order to clearly present the importance of these influencing factors (such asθ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, and φ) on the PCCS heat transfer process, the coefficient of sensitivity (COS) of the factor X is defined as the derivative of dimensionless qw,ave, ω, or δ¯, to the dimensionless variable of X itself, i.e.,(28)COSX,qw,ave=dqw,ave∗dX∗=Xoqw,aveodqw,avedX,COSX,ω=dω∗dX∗=XoωodωdX,COSX,δ¯=dδ¯∗dX∗=Xoδ¯odδ¯dX,where the superscript ∗ stands for the dimensionless variable, the subscript o represents the based case conditions, and X can be referred to the θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, or φ.Table5 shows the relative strength of the COS for the different variables about the based case condition. The values of COSX,qw,ave, COSX,ω, and COSX,δ¯are obtained from equation (28), and the derivatives of dimensionless qw,ave, ω, or δ¯ to the dimensionless variable X are derived from Figures 6–12. As seen from Table 5, Tw and Vg,in show higher effect on qw,ave and ω, but they show lower effect on δ¯. The more important influence on δ¯ are parameters of θ and m˙l,in. The parameter of Tg,in has the little influence on qw,ave, ω, and δ¯. Therefore, the efficient way to enhance the heat transfer capacity of PCCS depends on the control of Vg,in to reduce the flow resistance.Table 5
Parameter influence table.
Variable,XBased case condition,XoMinimum value,XminMaximum value,XmaxCOSX,qw,aveCOSX,ωCOSX,δ¯θ (deg)9015900.14040.10330.7234Tw (°C)8575983.26571.58580.2075Tl,in (°C)7020800.20920.07400.0150m˙l,in (kg/s)0.180.091.440.01950.04880.1938Tg,in (°C)455850.00050.00600.0010Vg,in (m/s)6280.67800.55680.0295φ (%)27.510.090.00.03760.02590.0024


## 4.1. Influence of the Tilted Angle

Figure6 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, qw,ave and ω increase with the increase of θ. The rising slope of qw,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. The 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 θ. Thus, 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. Therefore, m˙evaporation and ω will increase with the increase of θ, while the water film inlet mass flow rate remains the same.Figure 6
Influence of the tilted angle on heat and mass transfer of the water film evaporation process.


## 4.2. Influence of the Water Inlet Mass Flow Rate

Figure7 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, qw,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. There 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 7
Influence of the water film inlet mass flow rate on heat and mass transfer of the water film evaporation process.


## 4.3. Influence of the Water Inlet Temperature

Figure8 shows the influence of the water film inlet temperature Tl,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, qw,ave first decreases with the increase of Tl,in and then increases, and it reaches its minimum value when Tl,in is 70°C. The reason is that the heat removed by the water film evaporation is more than that by the water film convection. When Tl,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 Tl,in. However, when Tl,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 Tl,in, which will promote the water film evaporation. Therefore, when Tl,in increases, m˙evaporation increases. As shown in Figure 8, under the same m˙l,in, ω increases with Tl,in increasing, and δ¯ reduces with the increase of Tl,in. However, the effects of Tl,in on ω and δ¯ are relatively weak due to the small change of the water vapor partial pressure of the phase interface.Figure 8
Influence of the water film inlet temperature on heat and mass transfer of the water film evaporation process.


## 4.4. Influence of the Bottom Plate Surface Temperature

Figure9 shows the influence of the bottom plate surface temperature Tw 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, qw,ave increases with Tw increasing because of the increase of the heat transfer difference between the water film and the bottom plate surface. Thus, 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 Tw increases, the water vapor diffusion velocity from the water film surface to the mainstream moist air will become larger and m˙evaporation increases. Therefore, as shown in Figure 9, ω increases with the increase of Tw under the same m˙l,in, but δ¯ decreases. In addition, when Tw changes from 75°C to 98°C, the ranges of qw,ave, ω, and δ¯ are 60%, 50%, and 7%, respectively. It can be seen that Tw has a smaller effect on δ¯, but has a larger effect on qw,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 9
Influence of the bottom plate surface temperature on heat and mass transfer of the water film evaporation process.


## 4.5. Influence of the Air Inlet Velocity

Figure10 shows the influence of the air inlet velocity Vg,in on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. When Vg,in increases, qw,ave and ω increases, but δ¯ decreases. These trends of qw,ave, ω, and δ¯ in Figure 10 are the same as Figure 9. Because the increase of Vg,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. Therefore, m˙evaporation increases with the increase of Vg,in, so does qw,ave. When m˙l,in keeps unchanged, ω becomes larger and δ¯ gets smaller, with Vg,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 qw,ave.Figure 10
Influence of the air inlet velocity on heat and mass transfer of the water film evaporation process.


## 4.6. Influence of the Air Inlet Temperature

Figure11 shows the influence of the air inlet temperature Tg,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 Tg,in on qw,ave, ω, and δ¯ are negligible. When Tg,in changes from 5°C to 98°C, the ranges of qw,ave, ω, and δ¯ are 0.03%, 2.10%, and 0.18%, respectively. The 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 Tg,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 Tg,in is counteracted by the decrease of the forced convection heat flux, even leading to a weak decreasing trend of total heat flux [10]. Therefore, the dominant factor in this heat and mass transfer process is the Vg,in rather than Tg,in, comparing Figures 10 and 11.Figure 11
Influence of the air inlet temperature on heat and mass transfer of the water film evaporation process.


## 4.7. Influence of the Air Inlet Relative Humidity

Figure12 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, qw,ave and ω decrease with the increase of φ, but δ¯ increases with φ increasing; the ranges of qw,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. Therefore, 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 qw,ave and ω decrease when φ increases. Besides, during this water film evaporation process, only the phase change mass loss exists. Thus, δ¯ will become larger with the increase of φ.Figure 12
Influence of the air inlet relative humidity on heat and mass transfer of the water film evaporation process.


## 4.8. Influence Table of These Parameters

In order to clearly present the importance of these influencing factors (such asθ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, and φ) on the PCCS heat transfer process, the coefficient of sensitivity (COS) of the factor X is defined as the derivative of dimensionless qw,ave, ω, or δ¯, to the dimensionless variable of X itself, i.e.,(28)COSX,qw,ave=dqw,ave∗dX∗=Xoqw,aveodqw,avedX,COSX,ω=dω∗dX∗=XoωodωdX,COSX,δ¯=dδ¯∗dX∗=Xoδ¯odδ¯dX,where the superscript ∗ stands for the dimensionless variable, the subscript o represents the based case conditions, and X can be referred to the θ, Tw, Tl,in, m˙l,in, Tg,in, Vg,in, or φ.Table5 shows the relative strength of the COS for the different variables about the based case condition. The values of COSX,qw,ave, COSX,ω, and COSX,δ¯are obtained from equation (28), and the derivatives of dimensionless qw,ave, ω, or δ¯ to the dimensionless variable X are derived from Figures 6–12. As seen from Table 5, Tw and Vg,in show higher effect on qw,ave and ω, but they show lower effect on δ¯. The more important influence on δ¯ are parameters of θ and m˙l,in. The parameter of Tg,in has the little influence on qw,ave, ω, and δ¯. Therefore, the efficient way to enhance the heat transfer capacity of PCCS depends on the control of Vg,in to reduce the flow resistance.Table 5
Parameter influence table.
Variable,XBased case condition,XoMinimum value,XminMaximum value,XmaxCOSX,qw,aveCOSX,ωCOSX,δ¯θ (deg)9015900.14040.10330.7234Tw (°C)8575983.26571.58580.2075Tl,in (°C)7020800.20920.07400.0150m˙l,in (kg/s)0.180.091.440.01950.04880.1938Tg,in (°C)455850.00050.00600.0010Vg,in (m/s)6280.67800.55680.0295φ (%)27.510.090.00.03760.02590.0024


## 5. 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., Eulerian_EWF method). Meanwhile, these two numerical methods are validated by some experiments. The following conclusions could be drawn:(1)
The mass and heat transfer during the water film evaporation can be predicted by two numerical methods, i.e., VOF_UDF method and Eulerian_EWF method. The 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. The 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. Therefore, 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.



---

*Source: 5924325-2020-06-29.xml*