Entropy Assessment on Direct Contact Condensation of Subsonic Steam Jets in a Water Tank through Numerical Investigation

1 School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China; jiyu1994joe11@163.com (Y.J.); zhangmdc@163.com (Y.-N.Z.) 2 Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China 3 Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China; tongjf@ihep.ac.cn 4 State Nuclear Power Technology R & D Center, Beijing 102209, China; wangxuwei@snptc.com.cn (X.-W.W.); wanghan@snptc.com.cn (H.W.) * Correspondence: zhc5@vip.163.com; Tel.: +86-451-8641-2328


Introduction
Direct contact condensation (DCC) phenomena exists in many industrial facilities, such as the suppression pools of boiling water reactors, the emergency core cooling systems (ECCSs) for nuclear reactors, direct contact heat exchangers, and thus have attracted great interest and received much attention in the past several decades.Generally, steam is discharged into subcooled water at high speed, then contacts with water and condenses at the interface directly in a DCC process.Knowledge including the jet shape, plume length, temperature and pressure fluctuation is primarily required for engineering design, and a reliable prediction of the process is essential, because the rapid condensation could induce accidents under certain conditions during the DCC process.A wide variety of theoretical and experimental studies have been performed to investigate the DCC phenomenon, in attempting to provide technical support for equipment design and operation.Kerney et al. [1] and Weimer et al. [2] investigated the submerged sonic steam jet in quiescent sub-cooled water tanks by experimental and theoretical approaches, and a semi-empirical correlation to calculate the steam penetration length was developed.Aya and Nariai [3] proposed the correlation of heat transfer coefficients under low mass flux conditions through the experimental investigation of steam jets injected into cold water.Chan and Lee [4] proposed the regime maps at low steam mass flux, and the condensation mode was categorized into several classes.On top of the experimental investigations, Computational Fluid Dynamics (CFD) technique is also a powerful tool to analyze two phase flow processes.With the rapid development of CFD, several numerical studies have been performed to investigate the DCC phenomenon.Gulawani et al. [5,6] analyzed the temperature profile, flow pattern and heat transfer characteristics by using a thermal phase change model of commercial CFD code under a three dimensional simulation of DCC process, and the corresponding experiment was carried out to validate the results by CFD prediction.Dahikar et al. [7] investigated the flow pattern and temperature distribution of DCC using PIV (Particle Image Velocimetry), PLIF (Planar Laser Induced Fluorescence) and CFD, and during the simulation, two turbulence models, i.e., k-ε and Large Eddy Simulation (LES) were employed, through which the flow field, temperature field and phase change information can be monitored instantaneously.Shah et al. [8] simulated the DCC of supersonic steam jet under the Euler-Euler framework of FLUENT 6.3, and the results agreed fairly well with the published experimental data.Recently, Li et al. [9] performed the simulation using the Volume of Fluid (VOF) model for multiphase flow and LES turbulent model in FLUENT, through which, good consistance with the available experiment data was obtained.
As an effective method to estimate transport phenomena, the thermodynamic approach has been widely applied to viscous flow, heat and mass transfer processes [10][11][12][13][14].The traditional way in two phase flow emphasizes the essential thermal-hydraulic parameters, while the second law of thermodynamics provides another powerful tool to investigate phase changes [15][16][17].The condensation process is accompanied by thermodynamic irreversibility resulting from fluid friction, heat transfer over finite temperature gradients and inner phase change in fluid flow domain; therefore, a relationship exists between entropy generation and the energy devaluation in the process [18,19].Naterer and Adeyinka [20] proposed a correlation for entropy generation and exergy in film condensation; and in this way, entropy generation can be used as an important parameter in optimizing a two-phase system.Thiel [21] examined the mechanism of entropy generation in a condenser with high fractions of non-condensable gases using scaling and boundary layer technique, aiming at providing a criterion for minimum entropy generation, which is valuable in engineering analysis.Revellin et al. [22] applied the local entropy generation rate formulation catering to the saturated two phase flow to analyze the performance of heat exchangers in the various conditions.Furthermore, Herwig [23] and Wenterodt [24] introduced the concept of entropic potential of energy, which consider the energy transfer within a certain process from a broader perspective.
From the literature review above, it should be highlighted that, although numerous investigations for DCC have been performed, essentially, it is still not well understood.In the current work, a transient simulation for DCC process of steam is performed through the mixture model in the commercial CFD code, ANSYS-FLUENT, and a phase change model derived from the kinetic theory of gas is implemented to the solver through User Defined Functions (UDF).Moreover, the entropy generation analysis approach is used to estimate the irreversibility of DCC, which is expected to be a powerful tool for designing and optimizing the equipment involved in the phenomenon.

Geometry Model
DCC processes appear in various types of facilities, but it can be generally considered as a fundamental transport processes in a specified geometry for the same mechanism.In this work, the simulation is performed on a simplified hexagonal tank of 560 mm ˆ560 mm ˆ870 mm dimensions.
The tank is full of subcooled water, while saturated steam is ejected from a vertical downward pipe with a diameter of 5 mm.The pipe is centrally located in the top boundary of water tank and it is submerged vertically 250 mm in the water.A schematic view of the geometry model is shown in Figure 1.
Entropy 2016, 18, 21 3 with a diameter of 5 mm.The pipe is centrally located in the top boundary of water tank and it is submerged vertically 250 mm in the water.A schematic view of the geometry model is shown in Figure 1.

Mathematical Model
In the current work, the mixture model concerning the slip velocity is employed, through which the intensive shear interaction can be simulated in the non-homogenous two phase flow.Liquid water is treated as the primary phase, while the vapor is secondary phase.In addition, both phases are regarded as inter-penetrating continua in the simulation.The condensation and evaporation process occurring at the interface is described by the phase change model based on the kinetic theory of gas.For the turbulence model, the standard k-ε model with standard wall functions treatment is adopted.Moreover, an entropy generation assessment approach is proposed to analyze the irreversibility during the DCC process.

Mixture Model
The mixture model [25] is used to simulate separated phase flow, but it assumes local equilibrium over short spatial length scales.Compared to the full Eulerian model in [6][7][8], the model solves smaller number of variables; and therefore, it is assumed to be a good substitute for Eulerian model in several cases, including the simulation in the present work.

Continuity Equation
The continuity equation for mixture describing the mass flux into and out of a control volume boundary as well as the internal viration of mass is: where ρm and vm are the density and mass-averaged velocity for mixture, which are defined respectively as follows: where αq is the volume fraction of phase q.

Momentum Equation
The momentum equation for mixture solved throughout the calculating domain is given by:

Mathematical Model
In the current work, the mixture model concerning the slip velocity is employed, through which the intensive shear interaction can be simulated in the non-homogenous two phase flow.Liquid water is treated as the primary phase, while the vapor is secondary phase.In addition, both phases are regarded as inter-penetrating continua in the simulation.The condensation and evaporation process occurring at the interface is described by the phase change model based on the kinetic theory of gas.For the turbulence model, the standard k-ε model with standard wall functions treatment is adopted.Moreover, an entropy generation assessment approach is proposed to analyze the irreversibility during the DCC process.

Mixture Model
The mixture model [25] is used to simulate separated phase flow, but it assumes local equilibrium over short spatial length scales.Compared to the full Eulerian model in [6][7][8], the model solves smaller number of variables; and therefore, it is assumed to be a good substitute for Eulerian model in several cases, including the simulation in the present work.

Continuity Equation
The continuity equation for mixture describing the mass flux into and out of a control volume boundary as well as the internal viration of mass is: where ρ m and v m are the density and mass-averaged velocity for mixture, which are defined respectively as follows: where α q is the volume fraction of phase q.

Momentum Equation
The momentum equation for mixture solved throughout the calculating domain is given by: The terms on the left side of equation are the inner variation of momentum with time and convective momentum flux crossing the boundaries of control volumes, respectively.The terms on the right side are overall pressure, stresses, gravitational force, and interaction forces between phases due to slip velocity, in proper order.The interaction force is defined as follows: where µ m is the viscosity of mixture, defined as: q"1 α q µ q (5)

Energy Equation
The energy equation for the mixture is: where E q is the total energy of phase q; κ e f f is the effective thermal conductivity, which can be denoted as: in which, κ q is the thermal conductivity of phase q; κ t is the turbulent thermal conductivity, determined by the turbulence closure model; S Eq is the volumetric heat source for phase q.

Turbulence Model
The standard k-ε model [25] is chosen for turbulence model.In the formulation, the conservation equations for turbulent kinetic energy k and the turbulent dissipation ε are expressed as: where µ t,m is the turbulent viscosity of mixture, defined as: G k,m is the production of turbulent kinetic energy, and it is computed from: The value for the constant above are: σ k " 1.0, σ ε " 1.3, C µ " 0.09, C 1ε =1.44, C 2ε =1.92.

Phase Change Model
According to the kinetic theory of gas, the molecules emitted from both vapor and water will be absorbed totally when they move to the equilibrium interface, and the mass transfer rate between the phases can be calculated in Hertz-Knudsen relation; however, the interface is usually in non-equilibrium state, and the reflection exists there, as is shown in Figure 2.
Entropy 2016, 18, 21 the phases can be calculated in Hertz-Knudsen relation; however, the interface is usually in non-equilibrium state, and the reflection exists there, as is shown in Figure 2. Then the simplified Hertz-Knudsen-Schrage Equation gives the following form for a smooth phase interface [26,27]: In Equation (12), the evaporation coefficient and condensation coefficient are defined respectively by Knudsen and Prüger [28,29] as: The pressure and temperature can be correlated in the vicinity of the saturation condition using Clapeyron-Clausius equation, which is defined as: in which, L is the latent heat.In this view, the condensation and evaporation rate can be rewritten as: where T is the real temperature, Tsat (pg) is the saturated temperature corresponding to the pressure pg, and Tsat (pl) is the saturated temperature corresponding to the pressure pl.Then the simplified Hertz-Knudsen-Schrage Equation gives the following form for a smooth phase interface [26,27]: In Equation (12), the evaporation coefficient and condensation coefficient are defined respectively by Knudsen and Prüger [28,29] as: number of molecules transferred to the vapor phase number of molecules emitted from the liquid phase γ c " number of molecules absorbed by the liquid phase number of molecules impinging on the liquid phase The pressure and temperature can be correlated in the vicinity of the saturation condition using Clapeyron-Clausius equation, which is defined as: in which, L is the latent heat.In this view, the condensation and evaporation rate can be rewritten as: where T is the real temperature, T sat (p g ) is the saturated temperature corresponding to the pressure p g , and T sat (p l ) is the saturated temperature corresponding to the pressure p l .
It is also assumed that vapor bubble is in dispersed phase when it condenses, and liquid droplet is in dispersed phase when it evaporates at the interface.Furthermore, the dispersed phases contain Entropy 2016, 18, 21 6 of 23 spherical particles of the same diameter, and then the volumetric interfacial area can be derived as the following form: in which, d p is the diameter of dispersed phase, and α p is volume fraction of dispersed phase.Therefore, the volumetric condensation rate and evaporation rate are given by: In addition, volumetric phase change rate is:

Entropy Generation Analysis Model
The irreversibility of DCC process consists of dissipations due to viscous and turbulence fluctuation, finite temperature difference in heat transfer as well as the phase change within the domain.From the view of irreversible thermodynamics, the transport equation for entropy in a Cartesian coordinate system can be written as tensor form [30]: where J i s is the entropy flux in i th direction, and .

S 3
gen is the entropy generation rate.The governing entropy equation for transport process including single-mixture flow, Fourier heat conduction, and incompressible fluid is considered to deduce the entropy generation model over the calculation domains in turbulent multiphase flows, as is in [31]: In the equation above, the last three terms on the right side are the entropy generation owing to viscous, phase change and the heat transfer with finite temperature difference, in proper order.Besides, the instantaneous variable can be divided into time-averaged part and fluctuating part in the RANS approach, such as: The terms in the Equation ( 19) are substituted by time-averaged variables, and the main processes are shown as follows.

Convective Terms
The time-averaged convective term for entropy on the left side of Equation ( 19) is written as: Entropy 2016, 18, 21 7 of 23

Entropy Generation by Dissipation
Time-averaged entropy generation due to viscous dissipation contains two groups of terms, i.e., averaged flow part and fluctuating part, which read as: - in which, the T 1 m in the denominator is neglected [31].

Entropy Generation by Heat Transfer
In the time averaging process of entropy generation by heat transfer, the T 1 m is neglected for its existence in higher order terms.Then, the time-averaged entropy generation terms with respect to finite temperature gradients can be expressed as: Here, the first term on the right side is the entropy generation due to heat transfer with time-averaged temperature gradients, while the second term on the right side is due to fluctuating temperature gradients.

Entropy Generation by Inner Phase Change
Considering the fluctuation of volumetric phase change rate, the time-averaged entropy generation by inner phase change is given by the following form: where L is the latent heat of water.

Time-Averaged Transport Equation for Entropy
Substituting the time-averaged terms into the Equation ( 19), the transport equation for the mean entropy in turbulence conditions leads to: Obviously, the term (II) and the term (IV) concerning the velocity and temperature fluctuations in different entropy generation sources above are still unclosed.Then, the information contained in the k-ε turbulence model can be used to replace the two terms in a closed form, which will be shown as follows: (1) The exact turbulent dissipation approximately equals to the production of density ρ m and the turbulent dissipation rate ε, therefore, the entropy generation rate due to turbulent dissipation reads as: .
(2) Use the Boussinesque-like approach [31], and then the entropy generation due to fluctuating temperature gradients is: .
with a thermal diffusivity and a t the turbulent thermal diffusivity.Assume that the turbulent thermal diffusivity is related as thermal diffusivity through: Then, the entropy generation due to mean temperature gradients and entropy generation due to fluctuating temperature gradients can be combined as the entropy generation due to heat transfer: .
From the derivation above, the local volumetric entropy generation is concluded as follows: . looooooooooooooooooooooooooomooooooooooooooooooooooooooon

Simulation Details
In the simulation, the SIMPLE algorithm [32] is used for pressure-velocity coupling.In addition, the second order upwind discretization scheme is used for solving the momentum equations and energy equation, while a first order upwind discretization scheme is used for solving turbulent kinetic energy equations, turbulent dissipation rate equations and volume fraction equations.Considering the instability resulting from two phase flow and condensation, the under-relaxation factors are initially set to be 0.1 for all variables at the first 0.1 s, and they subsequently change to be default value for the remainder of simulation time, while the time step is 0.001 s based on the previous work conducted by other researchers [5][6][7]9].All the simulations were performed on a workstation with four Intel computer cores and 8 GB RAM.
As for the solutions residuals, two criteria are selected to determine whether they are convergent: (1) the residuals for all solutions except for energy solution should be below 1 ˆ10 ´3, while the residual for solution of energy should be below 1 ˆ10 ´6; (2) the entropy generation rate for all two phase domains hardly change with iterations in one time steps.
The boundary conditions are set as follows: the mass flux with value of 60 kg/m 2 s for vapor is enforced at the inlet section, the temperature and pressure for vapor is 374.15K and 3 kPa.The pressure outlet with value 0 Pa (gauge pressure) is employed for mixture at the outlet section.The rest boundaries are set to be adiabatic walls with no velocity slip.
The initial conditions are set as follows: inside the steam pipe is patched with vapor, the temperature and gauge pressure of which are 374.15K and 3 kPa separately, whilst the temperature and gauge pressure of accumulated water inside the tank are 300 K and 0 Pa, and the operating Entropy 2016, 18, 21 9 of 23 pressure is 101.325kPa, i.e., 1 atm.For simplify, the phase change temperature is set to be constant, 373.15 K.
Steam in the fluid domain is treated as an incompressible ideal gas, which means that the density and specific heat capacity would change with temperature only, while other properties such as viscosity and thermal conductivity are constant value, which are taken from the FLUENT material database.By contrast, liquid water is regarded as incompressible fluid; and its properties are all specified as constant value, as is shown in Table 1.In the simulation, the values of evaporation/condensation coefficient are treated as constant, 0.01, derived from the [33,34].
) are the default coefficients in ANSYS FLUENT.

Grid Independent Verification
In the current work, hexahedral elements are used for meshing, and grids with good quality are ensured using ICEM, a mesh generation tool.For the cylindrical pipe in the geometry, the block with O-Grid is adopted.The values of Y + is required to locate in the interval (11.225~200) for the adoption of standard wall functions, then the mesh near wall region is adapted to ensure 30 < Y + < 50.Totally, three computational grids consisted with 38,6972, 633,452, and 914,597 nodes are constructed.Before the simulation, the effect of grid density on the results is investigated.Figure 3a,b show the transverse distribution of longitudinal velocity at selected locations (z 1 = 265 mm, z 2 = 280 mm) with the three different grids above.All three grids can predict similarly the velocity profile but the medium size mesh, i.e., 633,452 nodes (« 80 ˆ80 ˆ100) is adopted considering the balance of accuracy and efficiency, as is shown in Figure 4. 9 enforced at the inlet section, the temperature and pressure for vapor is 374.15K and 3 kPa.The pressure outlet with value 0 Pa (gauge pressure) is employed for mixture at the outlet section.The rest boundaries are set to be adiabatic walls with no velocity slip.
The initial conditions are set as follows: inside the steam pipe is patched with vapor, the temperature and gauge pressure of which are 374.15K and 3 kPa separately, whilst the temperature and gauge pressure of accumulated water inside the tank are 300 K and 0 Pa, and the operating pressure is 101.325kPa, i.e., 1 atm.For simplify, the phase change temperature is set to be constant, 373.15 K.
Steam in the fluid domain is treated as an incompressible ideal gas, which means that the density and specific heat capacity would change with temperature only, while other properties such as viscosity and thermal conductivity are constant value, which are taken from the FLUENT material database.By contrast, liquid water is regarded as incompressible fluid; and its properties are all specified as constant value, as is shown in Table 1.In the simulation, the values of evaporation/condensation coefficient are treated as constant, 0.01, derived from the [33,34].

Grid Independent Verification
In the current work, hexahedral elements are used for meshing, and grids with good quality are ensured using ICEM, a mesh generation tool.For the cylindrical pipe in the geometry, the block with O-Grid is adopted.The values of Y + is required to locate in the interval (11.225~200) for the adoption of standard wall functions, then the mesh near wall region is adapted to ensure 30 < Y + < 50.Totally, three computational grids consisted with 38,6972, 633,452, and 914,597 nodes are constructed.Before the simulation, the effect of grid density on the results is investigated.Figure 3a,b show the transverse distribution of longitudinal velocity at selected locations (z1 = 265 mm, z2 = 280 mm) with the three different grids above.All three grids can predict similarly the velocity profile but the medium size mesh, i.e., 633,452 nodes (≈ 80 × 80 × 100) is adopted considering the balance of accuracy and efficiency, as is shown in Figure 4.

Verification & Validation
The method in the current study is based on the mixture model and the kinetic theory of gas, which is different from the previous work.To validate the current models, a simulation case with the same geometry as the model in the study by Takase et al. [35] has been performed.In this simulation, a velocity boundary which has a value of 25 m/s is adopted, the initial temperature of steam and accumulated water are 111 °C and 20 °C, respectively, and the gauge pressure inside the tank is 2.3 kPa.In addition, the steam changes to liquid water at 110.5 °C.
Figure 5 shows the comparisons of experimental results derived by Takase et al. [35] and the results by numerical predictions in the present work.Here, Figure 5a is the initial view, i.e., at time t = 0.0 s, while Figure 5b is the result at time t = 0.5 s.From Figure 5b, it can be seen that the length, spreading angle and the width of plume jet are quite similar to those in the experimental observation, but the experiment shows a wavy and turbulent flow structure, while the shape of plume jet by present CFD prediction is smooth.The most possible reason is the adoption of RANS model, in which the turbulent kinetic energy formulation considers that all the normal components of stresses are isotropic, i.e., the equal spreading of the plume jet in all directions.In general, the simulation steam plume shape agrees qualitatively well with the experimental observations.

Verification & Validation
The method in the current study is based on the mixture model and the kinetic theory of gas, which is different from the previous work.To validate the current models, a simulation case with the same geometry as the model in the study by Takase et al. [35] has been performed.In this simulation, a velocity boundary which has a value of 25 m/s is adopted, the initial temperature of steam and accumulated water are 111 ˝C and 20 ˝C, respectively, and the gauge pressure inside the tank is 2.3 kPa.In addition, the steam changes to liquid water at 110.5 ˝C.
Figure 5 shows the comparisons of experimental results derived by Takase et al. [35] and the results by numerical predictions in the present work.Here, Figure 5a is the initial view, i.e., at time t = 0.0 s, while Figure 5b is the result at time t = 0.5 s.From Figure 5b, it can be seen that the length, spreading angle and the width of plume jet are quite similar to those in the experimental observation, but the experiment shows a wavy and turbulent flow structure, while the shape of plume jet by present CFD prediction is smooth.The most possible reason is the adoption of RANS model, in which the turbulent kinetic energy formulation considers that all the normal components of stresses are isotropic, i.e., the equal spreading of the plume jet in all directions.In general, the simulation steam plume shape agrees qualitatively well with the experimental observations.

Verification & Validation
The method in the current study is based on the mixture model and the kinetic theory of gas, which is different from the previous work.To validate the current models, a simulation case with the same geometry as the model in the study by Takase et al. [35] has been performed.In this simulation, a velocity boundary which has a value of 25 m/s is adopted, the initial temperature of steam and accumulated water are 111 °C and 20 °C, respectively, and the gauge pressure inside the tank is 2.3 kPa.In addition, the steam changes to liquid water at 110.5 °C.
Figure 5 shows the comparisons of experimental results derived by Takase et al. [35] and the results by numerical predictions in the present work.Here, Figure 5a is the initial view, i.e., at time t = 0.0 s, while Figure 5b is the result at time t = 0.5 s.From Figure 5b, it can be seen that the length, spreading angle and the width of plume jet are quite similar to those in the experimental observation, but the experiment shows a wavy and turbulent flow structure, while the shape of plume jet by present CFD prediction is smooth.The most possible reason is the adoption of RANS model, in which the turbulent kinetic energy formulation considers that all the normal components of stresses are isotropic, i.e., the equal spreading of the plume jet in all directions.In general, the simulation steam plume shape agrees qualitatively well with the experimental observations.Figure 6 shows the comparisons of transverse temperature distribution at selected longitudinal location between the present numerical predictions and the published work by Takase et al. [35], where X + is the dimensionless transverse location, and Z + is the dimensionless distance departure from the steam pipe exit (X + = x/d 0 , Z + = z/d 0 , d 0 is the inner diameter of steam pipe).Both the results show steep value in the vicinity of the jet, and smooth temperature change in the outward zone.Besides, the CFD prediction agrees quantitatively well with the numerical results derived by Takase et al. [35] (within 5%) in the central location, whilst it follows the same trend of the results by Takase et al. [35] with marginal deviation within 14% in the outer location.
Entropy 2016, 18,21 Figure 6 shows the comparisons of transverse temperature distribution at selected longitudinal location between the present numerical predictions and the published work by Takase et al. [35], where X + is the dimensionless transverse location, and Z + is the dimensionless distance departure from the steam pipe exit (X + = x/d0, Z + = z/d0, d0 is the inner diameter of steam pipe).Both the results show steep value in the vicinity of the jet, and smooth temperature change in the outward zone.Besides, the CFD prediction agrees quantitatively well with the numerical results derived by Takase et al. [35] (within 5%) in the central location, whilst it follows the same trend of the results by Takase et al. [35] with marginal deviation within 14% in the outer location.

Velocity Profile
Figure 7 shows the velocity streamline in the x = 0 plane at different time, t = 4 ms, t = 8 ms, t = 12 ms, t = 16 ms, t = 44 ms, and t = 120 ms.From the figures, it can be seen that, once the steam discharging into the accumulated water, the interaction between steam and water due to momentum exchange decays the velocity of steam to 0 immediately in certain distance, which drives the surrounding water flow slightly.At the jet interface, surrounding water is entrained into the plume region due to buoyancy effect, and then a recirculation zone as shown in figures is formed.In this way, intense circulation in the tank enhance the thermal mixing effect during DCC processes.

Velocity Profile
Figure 7 shows the velocity streamline in the x = 0 plane at different time, t = 4 ms, t = 8 ms, t = 12 ms, t = 16 ms, t = 44 ms, and t = 120 ms.From the figures, it can be seen that, once the steam discharging into the accumulated water, the interaction between steam and water due to momentum exchange decays the velocity of steam to 0 immediately in certain distance, which drives the surrounding water flow slightly.At the jet interface, surrounding water is entrained into the plume region due to buoyancy effect, and then a recirculation zone as shown in figures is formed.In this way, intense circulation in the tank enhance the thermal mixing effect during DCC processes.
Transverse profile of V z for selected longitudinal positions (z = 0.26 m and z = 0.30 m) at different time are displayed in Figure 8, which indicate that V z of steam in the central location near pipe exit decreases rapidly after it reaches its maximum value.The observations correspond to the results obtained by Dahikar et al. [7] in their experiments and CFD simulation.Besides, it can be observed that the velocity in the centerline is not at the maximum all the time, and sometimes it increases in the jet zone deviated from the center line.Furthermore, the velocity changes sharply with the variation of transverse location near the region of pipe exit, while it changes smoothly in the outward zone, due to finite spreading width of jets in the pool water.obtained by Dahikar et al. [7] in their experiments and CFD simulation.Besides, it can be observed that the velocity in the centerline is not at the maximum all the time, and sometimes it increases in the jet zone deviated from the center line.Furthermore, the velocity changes sharply with the variation of transverse location near the region of pipe exit, while it changes smoothly in the outward zone, due to finite spreading width of jets in the pool water.

Temperature Field
The heat transfer process is modeled by associating the energy source with phase change in control volume, and transient simulation is carried out to determine the temperature field.Predictions of temperature in the x = 0 plane at different time and transverse profile of temperature for selected longitudinal position at different time are shown in Figures 9 and 10, respectively.From the Figure 9, it can be observed that steam temperature in the pipe is almost a constant, 374.15 K, and decreases dramatically when the steam is discharged into the water, due to the energy exchange between steam and pool water.Large temperature variation during the DCC process is characterized by Figure 9a-d

Temperature Field
The heat transfer process is modeled by associating the energy source with phase change in control volume, and transient simulation is carried out to determine the temperature field.Predictions of temperature in the x = 0 plane at different time and transverse profile of temperature for selected longitudinal position at different time are shown in Figures 9 and 10 respectively.From the Figure 9, it can be observed that steam temperature in the pipe is almost a constant, 374.15 K, and decreases dramatically when the steam is discharged into the water, due to the energy exchange between steam and pool water.Large temperature variation during the DCC process is characterized by Figure 9a-d Entropy 2016, 18, 21 13 obtained by Dahikar et al. [7] in their experiments and CFD simulation.Besides, it can be observed that the velocity in the centerline is not at the maximum all the time, and sometimes it increases in the jet zone deviated from the center line.Furthermore, the velocity changes sharply with the variation of transverse location near the region of pipe exit, while it changes smoothly in the outward zone, due to finite spreading width of jets in the pool water.

Temperature Field
The heat transfer process is modeled by associating the energy source with phase change in control volume, and transient simulation is carried out to determine the temperature field.Predictions of temperature in the x = 0 plane at different time and transverse profile of temperature for selected longitudinal position at different time are shown in Figures 9 and 10, respectively.From the Figure 9, it can be observed that steam temperature in the pipe is almost a constant, 374.15 K, and decreases dramatically when the steam is discharged into the water, due to the energy exchange between steam and pool water.Large temperature variation during the DCC process is characterized by Figure 9a-d In Figure 10, a similar transverse profile of temperature can be observed, which indicates that the temperature field (heat transfer) is relevant to the velocity profile.In the central location of steam jets, the entrainment of subcooled water into the plume generates small perturbation, resulting in lower temperature in the plume, e.g., z = 0.26 m at t = 4 ms, z = 0.30 m at t = 8 ms.Moreover, a sharp In Figure 10, a similar transverse profile of temperature can be observed, which indicates that the temperature field (heat transfer) is relevant to the velocity profile.In the central location of steam jets, the entrainment of subcooled water into the plume generates small perturbation, resulting in lower temperature in the plume, e.g., z = 0.26 m at t = 4 ms, z = 0.30 m at t = 8 ms.Moreover, a sharp In Figure 10, a similar transverse profile of temperature can be observed, which indicates that the temperature field (heat transfer) is relevant to the velocity profile.In the central location of steam jets, the entrainment of subcooled water into the plume generates small perturbation, resulting in lower temperature in the plume, e.g., z = 0.26 m at t = 4 ms, z = 0.30 m at t = 8 ms.Moreover, a sharp decrease of temperature occurs at the vapor-liquid interface, while smooth decrease of temperature within 5 K is observed far away from pipe exit, e.g., liquid region.The results of numerical predictions show the similar profile as the work by Gulawani [5,6].

Plume Shape
A vapor cavity called plume jet can be observed in the region near pipe exit after the steam ejecting from the nozzle.Figure 11 shows the contours plot of vapor void fraction in the x = 0 plane at different time.According to the plume shapes shown in the figure, three stages can be discriminated, i.e., the initial stage (Figure 11a,b), developing stage (Figure 11c-e) and the oscillatory stage (Figure 11f).In the initial stage, the steam spreads immediately after discharged into the subcooled water, and the plume shows no fixed shape; then in the developing stage, the plume begin to change its form from "pear" shape to elliptical boundary, and the size of plume grows gradually; in the oscillatory stage, the small perturbation makes the plume shape become ellipsoidal shape with wavy boundary.
decrease of temperature occurs at the vapor-liquid interface, while smooth decrease of temperature within 5 K is observed far away from pipe exit, e.g., liquid region.The results of numerical predictions show the similar profile as the work by Gulawani [5,6].

Plume Shape
A vapor cavity called plume jet can be observed in the region near pipe exit after the steam ejecting from the nozzle.Figure 11 shows the contours plot of vapor void fraction in the x = 0 plane at different time.According to the plume shapes shown in the figure, three stages can be discriminated, i.e., the initial stage (Figure 11a,b), developing stage (Figure 11c-e) and the oscillatory stage (Figure 11f).In the initial stage, the steam spreads immediately after discharged into the subcooled water, and the plume shows no fixed shape; then in the developing stage, the plume begin to change its form from "pear" shape to elliptical boundary, and the size of plume grows gradually; in the oscillatory stage, the small perturbation makes the plume shape become ellipsoidal shape with wavy boundary.There are some possible reasons to explain the phenomenon: (1) steam is discharged into the stagnant water at high speed, and strong fluctuation occurs in the initial stage, which prevents the plume from forming a fixed shape; (2) with the proceeding of steam condensation, the surrounding water is heated up and the degree of water subcooling decreases, which results in the growth of plume; (3) simultaneous decrease of fluctuation and subcooling make the plume grow faster, while instantaneous condensation gives rise to the formation of oscillation.

Mass Transfer
Figure 12 presents the instantaneous condensation rate in the DCC process of steam.If the boundary with vapor void fraction 0.5 is defined as vapor-liquid interface, it is found that steam mainly condenses at the interface compared with the contour plot of volume fraction of vapor, as is shown in Figure 11.The condensation rate in the numerical simulation is rather small, and the possible reason may be that the low steam mass flux generates less turbulent fluctuation, which obviously weakens the effect of thermal mixing associating with the fluctuation.There are some possible reasons to explain the phenomenon: (1) steam is discharged into the stagnant water at high speed, and strong fluctuation occurs in the initial stage, which prevents the plume from forming a fixed shape; (2) with the proceeding of steam condensation, the surrounding water is heated up and the degree of water subcooling decreases, which results in the growth of plume; (3) simultaneous decrease of fluctuation and subcooling make the plume grow faster, while instantaneous condensation gives rise to the formation of oscillation.

Mass Transfer
Figure 12 presents the instantaneous condensation rate in the DCC process of steam.If the boundary with vapor void fraction 0.5 is defined as vapor-liquid interface, it is found that steam mainly condenses at the interface compared with the contour plot of volume fraction of vapor, as is shown in Figure 11.The condensation rate in the numerical simulation is rather small, and the possible reason may be that the low steam mass flux generates less turbulent fluctuation, which obviously weakens the effect of thermal mixing associating with the fluctuation.There are some possible reasons to explain the phenomenon: (1) steam is discharged into the stagnant water at high speed, and strong fluctuation occurs in the initial stage, which prevents the plume from forming a fixed shape; (2) with the proceeding of steam condensation, the surrounding water is heated up and the degree of water subcooling decreases, which results in the growth of plume; (3) simultaneous decrease of fluctuation and subcooling make the plume grow faster, while instantaneous condensation gives rise to the formation of oscillation.

Mass Transfer
Figure 12 presents the instantaneous condensation rate in the DCC process of steam.If the boundary with vapor void fraction 0.5 is defined as vapor-liquid interface, it is found that steam mainly condenses at the interface compared with the contour plot of volume fraction of vapor, as is shown in Figure 11.The condensation rate in the numerical simulation is rather small, and the possible reason may be that the low steam mass flux generates less turbulent fluctuation, which obviously weakens the effect of thermal mixing associating with the fluctuation.In addition, the condensation rate in the initial stage is larger than the condensation rate in other stages, due to gradual decreasing subcooling of pool water.

Entropy Generation
The volumetric entropy generation rate (entropy generation rate per unit volume) is determined in this section, based on basic parameters required.Figure 13 shows the contour plot of total entropy generation rate (EGR) per unit volume in the x = 0 plane at different time.As is shown in Figure 13, two strips with larger value of the volumetric EGR occur in the vicinity of steam pipe exit.Upon the steam discharging into the stagnant water, there is no certain boundary to restrict the steam flow.From the velocity streamline plot (Figure 7), it can be found that the large velocity gradients exist near the pipe exit.And moreover, to the best of our knowledge, the turbulence fluctuations will appear at the vapor-liquid interface due to the Kelvin-Helmholtz instabilities.These two terms lead to great dissipation initially.Compared with the other stages, the local value of volumetric EGR in the initial stage is much larger, while the region possessing the considerable EGR is smaller due to the shorter penetration length and finite spreading width of the plume jet, as is shown in Figure 11.In the developing stage, the values of volumetric EGR reduces shaply, and the zones with lagre volumetric EGR continues to grow, which corresponds to the behaviours of In addition, the condensation rate in the initial stage is larger than the condensation rate in other stages, due to gradual decreasing subcooling of pool water.

Entropy Generation
The volumetric entropy generation rate (entropy generation rate per unit volume) is determined in this section, based on basic parameters required.Figure 13 shows the contour plot of total entropy generation rate (EGR) per unit volume in the x = 0 plane at different time.As is shown in Figure 13, two strips with larger value of the volumetric EGR occur in the vicinity of steam pipe exit.Upon the steam discharging into the stagnant water, there is no certain boundary to restrict the steam flow.From the velocity streamline plot (Figure 7), it can be found that the large velocity gradients exist near the pipe exit.And moreover, to the best of our knowledge, the turbulence fluctuations will appear at the vapor-liquid interface due to the Kelvin-Helmholtz instabilities.These two terms lead to great dissipation initially.Compared with the other stages, the local value of volumetric EGR in the initial stage is much larger, while the region possessing the considerable EGR is smaller due to the shorter penetration length and finite spreading width of the plume jet, as is shown in Figure 11.In the developing stage, the values of volumetric EGR reduces shaply, and the zones with lagre volumetric EGR continues to grow, which corresponds to the behaviours of plume jet in this stage.In the oscillatory stage, the profile of local EGR becomes asymmetrical owing to the perturbation of flow and heat transfer stated in the previous sections.The volumetric EGR is integrated to the whole two phase flow domain V, which gives the entropy generation rate: The volumetric EGR is integrated to the whole two phase flow domain V, which gives the entropy generation rate: .S gen " ż V .S 3 gen dV (30) representing the total dissipation in the whole domain.Table 2 has listed out the EGR due to four kinds of irreversibility and the total EGR, which are used for describing the energy dissipation at different time during the DCC processes.From the table, it can be seen that the total EGR is higher initially for the violent turbulence, and declines quickly in the following time, accompanying some oscillation.In the process, the absolute value of EGR resulting from turbulence dissipation decreases apparently, while the EGR on account of heat transfer and inner phase change increase.The possible reason accounting for the phenomenon may be that the heat transfer and phase change can not response in time when the steam is discharged into the stagnant water.The DCC process performs turbulence dominant firstly.And subsequently, the heat transfer and phase change start gradually.Moreover, the value volumetric EGR due to viscous dissipation is relatively low during the whole process.The decrease of EGR proves that the process is becoming more and more energy-economical.Figure 14 shows the variation of contributions to total EGR with time resulting from the four kinds of irreversibility.In the figure, the EGR due to turbulence fluctuation occupies the largest proportion among the four irreversibility element in the initial stage, then decreases sharply.The effects of heat transfer irreversibility and inner phase change irreversibility on the EGR increase, which makes DCC process become heat dominant in the developing stage.In the oscillatory stage, the contributions of four kinds of irreversibility to total EGR hardly changes with time in the macroscopic view, and the process is still heat dominant.However, some oscillation still exists here if a more detailed look was taken.From begining to the end, the effect of viscous dissipation is minor.The reason is assumed that: the process is a phenomenon with violent perturbation, and mean flow element has lost its significance among all irreversibility, which is a common view for steam jet in many researches.The results here are accordance with the derivation of local volumetric EGR above.The analysis above can contribute to understand and quantify the major reversibility during the whole process.From the table and figures, the three stages can be clearly discriminated according to different EGR characteristics and their partial proportions with time, which may be difficultly distinguished from the thermal-hydraulic information, such as velocity streamline and temperature profile solely.Based on the EGR approach, the dissipation characteristics of the DCC is derived, which may contribute to better understand the process.Then, it may also contribute to optimize the design for facilities, including pipe diameter and pipe length, to minify the size of two strips of volumetric EGR.In addition, it may help optimization of configuration for operation conditions, including temperature of water and mass flux of steam, to reduce the value of EGR.Moreover, other factors such as the pipe with or without diffusor can be investigated using the EGR assessment approach for thermal mixing analysis concerning DCC phenomenon.The discussions will be further conducted in future works.

Conclusions
In the current work, the DCC process of steam jet in a subcooled water tank is investigated through numerical simulation and the entropy generation assessment model based on the local equilibrium hypothesis is developed to evaluate the energy dissipation characteristics.The main inferences can be summarized as follows: (1) A different mathematical model from the previous studies by other researchers is proved to be valid through the comparisons with the results obtained by Takase et al [35], Dahikar et al [7] and Gulawani [5,6].(2) Three distinct stages of DCC are discriminated clearly at the present conditions, i.e., initial stage, developing stage and oscillatory stage.In the initial stage, the plume shows no fixed shape.In the developing stage, the plume begins to act as an elliptical boundary, and the size of the plume grows quickly.In the oscillatory stage, the plume shape becomes ellipsoidal shape with disturbed structure.(3) The local volumetric EGR in the initial stage is much larger than those in other stages, but the region possessing considerable entropy generation rate is smaller than other stage.The decrease of EGR proves that the process conform to increasingly economical energy utilization.(4) The largest proportion in total EGR is occupied by turbulence fluctuation in the initial stage, and then it decreases apparently in the following time, meanwhile, the contributions of heat transfer irreversibility and inner phase change irreversibility to the local entropy generation increase, which makes DCC process become heat dominant in the developing and the From the table and figures, the three stages can be clearly discriminated according to different EGR characteristics and their partial proportions with time, which may be difficultly distinguished from the thermal-hydraulic information, such as velocity streamline and temperature profile solely.Based on the EGR approach, the dissipation characteristics of the DCC is derived, which may contribute to better understand the process.Then, it may also contribute to optimize the design for facilities, including pipe diameter and pipe length, to minify the size of two strips of volumetric EGR.In addition, it may help optimization of configuration for operation conditions, including temperature of water and mass flux of steam, to reduce the value of EGR.Moreover, other factors such as the pipe with or without diffusor can be investigated using the EGR assessment approach for thermal mixing analysis concerning DCC phenomenon.The discussions will be further conducted in future works.

Conclusions
In the current work, the DCC process of steam jet in a subcooled water tank is investigated through numerical simulation and the entropy generation assessment model based on the local equilibrium hypothesis is developed to evaluate the energy dissipation characteristics.The main inferences can be summarized as follows: (1) A different mathematical model from the previous studies by other researchers is proved to be valid through the comparisons with the results obtained by Takase et al. [35], Dahikar et al. [7] and Gulawani [5,6].(2) Three distinct stages of DCC are discriminated clearly at the present conditions, i.e., initial stage, developing stage and oscillatory stage.In the initial stage, the plume shows no fixed shape.In the developing stage, the plume begins to act as an elliptical boundary, and the size of the plume grows quickly.In the oscillatory stage, the plume shape becomes ellipsoidal shape with disturbed structure.(3) The local volumetric EGR in the initial stage is much larger than those in other stages, but the region possessing considerable entropy generation rate is smaller than other stage.The decrease of EGR proves that the process conform to increasingly economical energy utilization.(4) The largest proportion in total EGR is occupied by turbulence fluctuation in the initial stage, and then it decreases apparently in the following time, meanwhile, the contributions of heat transfer irreversibility and inner phase change irreversibility to the local entropy generation increase, which makes DCC process become heat dominant in the developing and the oscillatory stage.The variation of EGR can be used to characterize the the dissipation and proceeding of DCC process.
In general, the current work presents a promising method to investigate the DCC process, in which an analytical solution to determine condensation and evaporation rate is introduced.And furthermore, it also proposes an effective approach based on the entropy generation for better understanding the DCC process, including the the thermal hydraulic and the energy dissipation information, which can contribute to design the thermodynamically efficient equipment involved in the DCC process in future works.

Figure 1 .
Figure 1.Schematic view of physical model used in the simulation of DCC process.

Figure 1 .
Figure 1.Schematic view of physical model used in the simulation of DCC process.

Figure 2 .
Figure 2. Molecular mechanisms of condensation and evaporation at vapor-liquid interface.
number of molecules transferred to the vapor phase number of molecules emitted from the liquid phase e   number of molecules absorbed by the liquid phase number of molecules impinging on the liquid phase c  

Figure 2 .
Figure 2. Molecular mechanisms of condensation and evaporation at vapor-liquid interface.

Figure 4 .
Figure 4. Schematic view of simulation mesh of the hexahedral tank.

Figure 5 .
Figure 5.The comparisons of simulation steam shape and the experimental observations at different time; (a) t = 0.0 s; (b) t = 0.5 s.

Figure 4 .
Figure 4. Schematic view of simulation mesh of the hexahedral tank.

Figure 5 .
Figure 5.The comparisons of simulation steam shape and the experimental observations at different time; (a) t = 0.0 s; (b) t = 0.5 s.

Figure 5 .
Figure 5.The comparisons of simulation steam shape and the experimental observations at different time; (a) t = 0.0 s; (b) t = 0.5 s.

Figure 6 .
Figure 6.The comparisons of the transverse temperature distribution at selected longitudinal locations between published data and present CFD predictions.

Figure 6 .
Figure 6.The comparisons of the transverse temperature distribution at selected longitudinal locations between published data and present CFD predictions.

Figure 7 .
Figure 7. Velocity streamline in the x = 0 plane at different time.(a) t = 4 ms; (b) t = 8 ms; (c) t = 12 ms; (d) t = 16 ms; (e) t = 44 ms; (f) t = 120 ms.Transverse profile of Vz for selected longitudinal positions (z = 0.26 m and z = 0.30 m) at different time are displayed in Figure 8, which indicate that Vz of steam in the central location near pipe exit decreases rapidly after it reaches its maximum value.The observations correspond to the results

Figure 8 .
Figure 8. Transverse profile of Vz for selected longitudinal position at different time; (a) z = 0.26 m; (b) z = 0.30 m.

Figure 8 .
Figure 8. Transverse profile of V z for selected longitudinal position at different time; (a) z = 0.26 m; (b) z = 0.30 m.
, while the predictions in Figure 9e-f indicate that temperature field varies little in the following time.The minor temperature change of pool water is due to low mass flux (60 kg/m 2 s) and low temperature (374.15K) of steam.

Figure 8 .
Figure 8. Transverse profile of Vz for selected longitudinal position at different time; (a) z = 0.26 m; (b) z = 0.30 m.

Figure 10 .
Figure 10.Transverse profile of temperature for selected longitudinal position at different time; (a) z = 0.26 m; (b) z = 0.30 m.

Figure 14 .
Figure 14.The variation of contributions of four kinds of irreversibility to total entropy generation with time.

Figure 14 .
Figure 14.The variation of contributions of four kinds of irreversibility to total entropy generation with time.

Table 1 .
The thermophysical properties of fluids.
* c p

Table 1 .
The thermophysical properties of fluids.