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.
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 Figure 1, the rectangular channel is formed by a bottom steel plate (red), a top glass cover plate (transparent), and two side steel walls (gray). The channel has a length of L (5.0 m), width of (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 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.
The bottom steel plate has the uniform temperature , 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 and , respectively. In addition, and Tl,in, respectively, refers to the water film velocity and temperature. is derived from the water film inlet mass flow rate . 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 (θ, , Tl,in, , , , and φ) in each research target are shown in Table 1. The based case parameter conditions are θ = 90 deg, , Tl,in = 70°C, , , , and φ = 27.5%.
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 bywhere and yi are the density and the water vapor mass fraction of the moist air, respectively; u is the velocity vector; 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:where and Ma are the molar mass of the water vapor and air, respectively; φ is the relative humidity of the moist air; 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:where psat is the saturated water vapor partial pressure of the moist air.
Besides, in (1) can be given by [24]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 source SΦ during the water evaporation process should be given bywhere γ is the latent heat of water vapor and and 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 , 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:where Ti,i is the temperature of the grid cell at the gas-liquid phase interface; 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 is the gradient of yi (Ti,i) on the 0-0 surface.
Moreover, because the net mass increase is zero for the dry air on the 0-0 surface, in (6) should also satisfy
Then, substituting (7) into (6), Si can be obtained as follows:
Therefore, the total mass exchange Sm during the water film evaporation process can be obtained by accumulating all the Si of the phase interface cells as follows: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):where αl is the liquid phase volume fraction. The sum of αl and is 1, i.e., .
The momentum equation of this VOF multiphase model can be written as follows:where I is the unit stress tensor; ρ and μ are the average density and dynamic viscosity of the control volume weighted by the volume fraction, respectively, with and ; g is the gravity acceleration vector; and f is the volume force source term, which mainly refers to the surface tension. Because the surface tension acts on the phase interface, f can be calculated by the continuous surface force (CSF) model [27] as follows:where σ is the surface tension coefficient between liquid phase and gas phase; kl is the curvature of the interface; and are the unit normal vector and unit tangent vector of the wall, respectively; and is the contact angle of the wall.
The energy equation of this VOF multiphase model can be written as follows:where E and λ are the average total energy and thermal conductivity of the control volume weighted by the volume fraction, respectively, with and ; 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):where represents the mass exchanged from the p phase to the q phase, with ; uq is the q phase velocity; upq is the phase slip velocity (when , when ); tq is the q phase stress tensor, satisfying ; 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 , the virtual mass force , 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; is the energy exchanged between phases, with ; and Sq, sq, and 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):where δ is the water film thickness; ul is the water film velocity; 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; is the water film kinematic viscosity; q indicates the force related to the collection or separation of the water film, satisfying ; Tl,ave is the water film average temperature; Ti is the phase interface temperature; is the wall temperature; represents the heat released by the water film when impinging the wall; 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, in equation (20) can be calculated as [2, 19]where Cphase, Ccondensation, and Cevaporation are constants, with , and Tsat is the water vapor saturated temperature, which depends on the water vapor saturated partial pressure psat.
The expression of psat is shown as follows: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 corresponding Tl,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 Figure 3, 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.
These parameters, such as the average water film temperature Tl,ave, average air temperature Tf, water film evaporative mass flow , and total heat transfer rate of heated wall , can be directly obtained in the FLUENT code. Thus, the expressions of the average heat flux of heated wall , water film evaporative ratio , and mass transfer coefficient hm are shown as equations (23)–(25):where A is the area of the bottom wall and ρi and ρ∞ are the density of the water vapor at the temperature of 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 , , Tl,ave, Tf, and hm. The specific parameter values of θ, , Tl,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 and from the VOF_UDF method agrees with those from the Eulerian_EWF method, and the differences of and between these two methods are 7.25% and 7.82%, respectively.
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.
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 surface qload can be calculated bywhere Stotal is the area of the bottom plate; is the total heat transfer rate of the bottom plate; is the air convection heat transfer rate; is the water film convection heat transfer rate; is the water film evaporation heat transfer rate; Cp,air (Tf) is the air specific heat at the temperature Tf; is the air mass flow; is the air outlet average temperature; ρair () is the air density at the temperature ; is the air outlet average velocity; Cp,water (Tl,ave) is the water specific heat at the temperature Tl,ave; 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 Table 4 to validate these two numerical methods. The results of the average bottom plate surface heat flux 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 or the larger , 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 and ω from these two numerical methods are relatively small, within %. 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.
(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 as , ω, and . Therefore, it is wise to select the Eulerian_EWF method in the parameter sensitivity analysis of θ, , Tl,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 Table 1) on the heat and mass characteristics, the calculated values of the average bottom plate surface heat flux , 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 :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
Figure 6 shows the influence of the tilted angle θ on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, and ω increase with the increase of θ. The rising slope of 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 (). 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, and ω will increase with the increase of θ, while the water film inlet mass flow rate remains the same.
4.2. Influence of the Water Inlet Mass Flow Rate
Figure 7 shows the influence of the water film inlet mass flow rate 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 because the increase of will cause the water film velocity to increase simultaneously in both the length and width direction of the channel. In addition, first increases with the increase of and then decreases, and it reaches its maximum value when 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 is 0.09–0.18 kg/s, while it reduces due to the increase of the water film thermal resistance when the range of is 0.18–1.44 kg/s. Moreover, ω decreases with the increase of because rapidly increases while slightly changes.
4.3. Influence of the Water Inlet Temperature
Figure 8 shows the influence of the water film inlet temperature 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, 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, increases. As shown in Figure 8, under the same , ω 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.
4.4. Influence of the Bottom Plate Surface Temperature
Figure 9 shows the influence of the bottom plate surface temperature 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 with 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 increases, the water vapor diffusion velocity from the water film surface to the mainstream moist air will become larger and increases. Therefore, as shown in Figure 9, ω increases with the increase of under the same , but decreases. In addition, when changes from 75°C to 98°C, the ranges of , ω, and are 60%, 50%, and 7%, respectively. It can be seen that has a smaller effect on , but has a larger effect on 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].
4.5. Influence of the Air Inlet Velocity
Figure 10 shows the influence of the air inlet velocity on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. When increases, and ω increases, but decreases. These trends of , ω, and in Figure 10 are the same as Figure 9. Because the increase of 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, increases with the increase of , so does . When keeps unchanged, ω becomes larger and gets smaller, with 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 .
4.6. Influence of the Air Inlet Temperature
Figure 11 shows the influence of the air inlet temperature 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 on , ω, and are negligible. When changes from 5°C to 98°C, the ranges of , ω, 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 will cause slight deterioration of the air convective heat transfer. It is also shown that the slight increase of evaporation heat flux at low 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 rather than , comparing Figures 10 and 11.
4.7. Influence of the Air Inlet Relative Humidity
Figure 12 shows the influence of the air inlet relative humidity φ on the heat and mass transfer of the water film evaporation process in the rectangular channel with the countercurrent air flow. As can be seen, and ω decrease with the increase of φ, but increases with φ increasing; the ranges of , ω, 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 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 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 φ.
4.8. Influence Table of These Parameters
In order to clearly present the importance of these influencing factors (such as θ, , Tl,in, , , , and φ) on the PCCS heat transfer process, the coefficient of sensitivity (COS) of the factor X is defined as the derivative of dimensionless , ω, or , to the dimensionless variable of X itself, i.e.,where the superscript ∗ stands for the dimensionless variable, the subscript o represents the based case conditions, and X can be referred to the θ, , Tl,in, , , , or φ.
Table 5 shows the relative strength of the COS for the different variables about the based case condition. The values of , , and are obtained from equation (28), and the derivatives of dimensionless , ω, or to the dimensionless variable X are derived from Figures 6–12. As seen from Table 5, and show higher effect on and ω, but they show lower effect on . The more important influence on are parameters of θ and . The parameter of has the little influence on , ω, and . Therefore, the efficient way to enhance the heat transfer capacity of PCCS depends on the control of to reduce the flow resistance.
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.
Nomenclature
| A: | Total water film coverage area on the heated wall (m2) |
| Ai: | Grid cell area of the phase interface (m2) |
| Cp: | Specific heat at constant pressure (J/kg/K) |
| Cphase: | Phase change constant |
| Cμ: | Constant |
| D: | Water vapor mass diffusion coefficient in the dry air (m2/s) |
| d: | Humidity ratio of the dry air (kg/kg) |
| E: | Total energy of the control volume (J/kg) |
| f: | Volume force source (N/m3) |
| ftotal: | Sum of other forces, such as the external volume force, lift force, wall lubrication force, virtual mass force, and turbulent diffusion force (N/m3) |
| g: | Gravity acceleration vector (m/s2) |
| h: | Specific enthalpy (J/kg) |
| H: | Height (m) |
| hm: | Mass transfer coefficient (m/s) |
| i: | Grid cell index of the phase interface |
| I: | Unit stress tensor |
| J: | Diffusion flux (kg/m2/s) |
| k: | Turbulent kinetic energy (m2/s2) |
| kl: | Curvature of the interface (m−1) |
| L: | Length (m) |
| m1: | Readings of weighing instrument corresponding to the time t1 (kg) |
| m2: | Readings of weighing instrument corresponding to the time t2 (kg) |
| Ma: | Molar mass of the dry air (g/mol) |
| : | Numerical calculated water film evaporative mass flow rate (kg/s) |
| : | Air mass flow rate (kg/s) |
| : | Water film inlet mass flow rate (kg/s) |
| : | Measured water film outlet mass flow rate (kg/s) |
| : | Phase exchanged mass flow caused by the evaporation of the water film or the condensation of the water vapor (kg/m2/s) |
| : | Mass exchanged between the phases (kg/m3/s) |
| : | Change of the water film mass flow due to the formation, separation, and splashing (kg/m2/s) |
| : | Molar mass of the water vapor (g/mol) |
| n: | Unit normal vector of the phase interface |
| N: | Total number of cells in the phase interface |
| : | Unit normal vector of the wall |
| p: | Total gas phase pressure (Pa) |
| pgas: | Pressure gradient loss of the air (Pa) |
| ph: | Gravity component perpendicular to the wall surface (Pa) |
| pL: | Total surface force (Pa) |
| psat: | Saturated water vapor partial pressure of the moist air (Pa) |
| : | Water vapor partial pressure of the moist air (Pa) |
| pσ: | Surface tension (Pa) |
| q: | Force related to the collection or separation of the water film (Pa) |
| : | Heat released by the water film when impinging the wall (W/m2) |
| qload: | Measured the average heat flux of the bottom plate surface (W/m2) |
| qq: | Conduction diffusion heat flux vector (W/m2) |
| : | Numerical calculated the average heat flux of the bottom plate surface (W/m2) |
| R: | Universal gas constant (J/mol/K) |
| rpq: | Internal force between phases (N/m3) |
| S: | Mass source term in equation (15) (kg/m3/s) |
| s: | Force source term in equation (16) (N/m3) |
| Sc,i: | Convective mass flow rate of the water vapor (kg/m3/s) |
| Sd,i: | Diffusion mass flow rate of the water vapor (kg/m3/s) |
| Si: | Exchanged mass flow rate of the grid cell of the phase interface (kg/m3/s) |
| Sm: | Total exchanged mass flow rate at the phase interface (kg/m3/s) |
| Stotal: | Total area of the heated wall (m2) |
| SФ: | Energy source term in equation (17) (W/m3) |
| T: | Temperature of grid cell (K) |
| t: | Time (s) |
| t: | Shear tensor (N/m2) |
| t1: | Start time of the weighing (s) |
| t2: | End time of the weighing (s) |
| Tf: | Average temperature of the air (K) |
| : | Air inlet temperature (K) |
| Ti: | Grid cell temperature of the phase interface (K) |
| Tl,ave: | Average temperature of the water film (K) |
| Tl,in: | Water film inlet temperature (K) |
| tq: | Shear tensor (N/m3) |
| Tsat: | Water vapor saturated temperature (K) |
| : | Unit tangent vector of the wall |
| : | Temperature of the heated wall (K) |
| u: | Velocity vector (m/s) |
| upq: | Slip velocity vector between the phases (m/s) |
| : | Normal moist air velocity scalar of the phase interface (m/s) |
| : | Air outlet average velocity (m/s) |
| : | Air inlet velocity (m/s) |
| : | Grid cell volume of the phase interface (m3) |
| : | Width (m) |
| X: | Influencing factor of PCCS heat and mass transfer process |
| x: | Cartesian coordinates of x direction, i.e., width direction of channel or the coverage width of water film on the wall (m) |
| y: | Cartesian coordinates of y direction, i.e., height direction of channel or the thickness of water film (m) |
| y+: | Dimensionless distance of the first cell to wall |
| yi: | Water vapor mass fraction of the moist air |
| z: | Cartesian coordinates of z direction, i.e., length direction of channel or the direction of airflow (m) |
| α: | Volume fraction of the phase |
| γ: | Latent heat of water vapor (J/kg) |
| δ: | Water film thickness (mm) |
| : | Average water film thickness of the heated wall (mm) |
| ε: | Turbulent dissipation rate (m2/s3) |
| : | Contact angle of the wall (deg) |
| λ: | Heat conductivity coefficient (W/m/K) |
| μ: | Dynamic viscosity (kg/m/s) |
| μt: | Turbulent eddy viscosity coefficient (kg/m/s) |
| ρ: | Density (kg/m3) |
| ρ∞: | Density of the water vapor at the average air temperature (kg/m3) |
| ρi: | Density of the water vapor at the average water film temperature (kg/m3) |
| σ: | Surface tension coefficient between liquid phase and gas phase (N/m) |
| σt: | Turbulent Prandtl constant |
| : | Water film kinematic viscosity (m2/s) |
| φ: | Relative humidity of the moist air (%) |
| : | Water film evaporative ratio (%) |
| : | Convection heat transfer rate of the air (W) |
| : | Evaporative heat transfer rate of the water film (W) |
| : | Convection heat transfer rate of the water film (W) |
| : | Measured total heat transfer rate of the bottom plate surface (W) |
| : | Exchanged energy between phases in equation (17) (W/m3) |
| : | Energy source term in equation (17) (W/m3) |
| : | Numerical calculated the total heat transfer rate of the bottom plate surface (W) |
| : | Gradient operator (m−1) |
| air: | Air |
| : | Gas phase |
| , ave, out: | Parameter average value of the air outlet section |
| l: | Liquid or water film |
| l, ave, out: | Parameter average value of the water film outlet section |
| max: | Maximum value of the variable X |
| min: | Minimum value of the variable X |
| o: | Based case |
| q: | qth phase of multiphase flow model |
| : | Vapor |
| T: | Transpose |
| : | Dimensionless. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The study was financially supported by the National Science and Technology Major Project (no. 2011ZX06002-005).