Boiling and Condensation in Two-Phase System Transients with Water Hammer

Abstract

Water hammer in two-phase systems, induced by direct steam condensation on subcooled water or by separation of subcooled water column, results in the most intensive pipeline pressure surges. Amplitudes of pressure spikes along the course of these dangerous transients strongly depend on the condensation and evaporation rates. The present paper provides a literature overview of thermal-hydraulic models for the prediction of water hammer phenomenon in two-phase systems, together with an original mechanistic approach for the prediction of phase transition rates, based on the shape and size of vapor-liquid interfacial area and the phase transition potential expressed through vapor and liquid phase temperature difference. Available water hammer experimental conditions were numerically simulated with the new modeling approach. Driving parameters of boiling and condensation rates at the steam-water interfaces are evaluated, and a good agreement is shown between numerical results and experimental data of bulk two-phase flow parameters during water hammer transients.

1. Introduction

The water hammer is recognized as a very dangerous phenomenon, and therefore, the prediction of its occurrence is necessary in order to prevent accidents. Depending on the mechanisms of their origin, we can divide water hammers into two types [1]. One type of water hammer, in single-phase flow, is caused by a quick valve opening or closure or when a pump suddenly stops. During the course of a transient with rapid flow rate changes, a column separation might occur, which is characterized with the rapid gaseous and vaporous cavitation and the formation of a two-phase system [2]. The other type of water hammer is caused by rapid condensation of steam in direct contact with subcooled liquid in a pipe or vessel under pressure. This second type is called condensation induced water hammer (CIWH). The CIWH and water hammer with column separation take place in two-phase system with rapid boiling and condensation. An overview of investigation of these two mechanisms of water hammers in two-phase systems follows.

1.1 Water hammer with column separation

The safety of various hydraulic systems depends on the accuracy of the prediction of water hammer with gaseous and vaporous cavitation. The first pressure surge in the single-phase system, caused by a sudden valve closure, causes a maximum peak pressure, close to the pressure rise in liquid obtained by the Joukowsky relation,$\mathrm{\Delta }p=-\mathit{\rho c}\mathrm{\Delta }u$, where ρ is the fluid density, c is the sonic velocity in the fluid-pipe system, and Δu is the fluid velocity change. The intensities of the subsequent pressure peaks during the water hammer transient are greatly affected by the presence of dissolved air in the liquid and vapor generation due to the evaporation of liquid if pressure is lower or equal to the saturation pressure. In all systems where water comes into contact with air, the air dissolves in the water so that water contains entrained air microbubbles. The air content in untreated tap water is 1.13 × 10⁹ microbubbles per m³. The most probable air microbubble diameter is 6 × 10⁻⁶ m. In degassed water, the air content is reduced to 0.911 × 10⁹ microbubbles per m³ with the same most probable microbubble diameter [3]. A detailed review of solubility of air and solubility of other gases in water can be found in [4]. Before the occurrence of the water hammer, entrained air microbubbles have the same velocity as a water volume and practically no influence on the thermophysical parameters of air-water mixture. If during the pressure transient pressure drops and water is degassed, a significant influence on the hydrodynamics of the mixture of water and air is observed. In cases with pressure drops below the saturation pressure along the course of the transient, the liquid becomes superheated, and adiabatic evaporation occurs, i.e., bubbles are generated although no heat is added to the liquid. The main difference between the two-phase mixture patterns with adiabatic evaporation and diabatic wall boiling is in the locations of bubbles nucleation and rise and corresponding void distribution. In case of adiabatic evaporation, the nucleation of bubbles occurs both on the heated microscopically rough wall surface and on the impurities within the bulk of liquid phase, i.e., a rapid bubbling or so-called flashing occurs within the whole liquid volume, while in case of wall boiling, the liquid is superheated within the thin liquid layer on the heated wall, and the bubbles are mainly generated on the wall surface and hydrodynamically transferred to the bulk of liquid volume [5]. Therefore, in cases of wall boiling, the void fraction is mainly higher in the vicinity of wall surface, which especially holds in cases of subcooled boiling with bubbles condensation within the bulk of subcooled liquid.

Bergant and Simpson [6] compared results of several numerical models with the data measured within an experimental test of the water hammer caused by the rapid valve closure. The system is described with a set of one-dimensional equations. This set includes water hammer equations for single liquid flow, two-phase flow equations for a distributed vaporous cavitation region, shock equations for condensation of liquid-vapor mixture back to the liquid and equations for a discrete vapor cavity separating a liquid and a vaporous cavitation region. The occurrence of cavitation during the transient was simulated using three models: the discrete vapor cavitation model (DVCM), the discrete gas cavitation model (DGCM) and the generalized interface vapor cavitation model (GIVCM). In the DVCM, it is assumed that the location of the gas phase is at the grid nodes, the liquid phase fills the space between nodes, the speed of wave propagation between adjacent nodes is equal to the sonic velocity in the liquid phase, the minimum pressure during the transient is determined by the pressure at which the first vapor bubble forms (liquid/vapor saturation pressure at the liquid temperature), and there is no wave propagation in the two-phase flow [7]. The model is applied to discrete cavities in which vapor bubbles as well as vapor cavities are placed. The DVCM is the most widely used for modeling vapor cavitation during transient in hydraulic systems, but its accuracy depends on the ratio of the cavitation volume and the volume between the numerical nodes, which is its main drawback [8]. The DGCM model is simple and gives good results for a wide range of input parameters if the gas fraction in the working fluid is small [9]. It is fully specified with its characteristic equations, the continuity equation for the ith numerical node and the ideal gas law. DGCM has been successfully used to model both vapor and gas cavitation. The GIVCM explicitly describes discrete cavities and vapor cavitation regions. As a basis for the development of the cavitation vapor interface, the algorithm of the DVCM was used, which allows the cavities to be formed into a branched network, which is calculated using the method of characteristics (MOC). The weakness of a model, in comparison with DVCM, is a long computational time. The main difference between the three models used is in the physical interpretation and description of the formation of the vapor cavitation region that spreads from the middle cavity towards the closed valve.

The limitations of DVCM, DGCM and GIVCM cavitation models, which are a consequence of adopted assumptions, are overcome in the presented mechanistic modeling approach, based on non-equilibrium gas release and absorption as well as on vaporous cavitation. Previous research has shown the lack of data on the air content in the liquid during the pressure transient and that it is necessary to explain the gas release in more detail. In this book chapter, the homogeneous gas-liquid two-phase flow model is applied to the simulation of water hammer with gaseous and vaporous cavitation, which includes non-equilibrium gas-liquid and vapor-liquid mass transfer at the interface.

1.2 Condensation induced water hammer

Direct contact of steam and subcooled water leads to CIWH. Since the specific volume of steam is significantly greater than the specific volume of liquid, there is a pressure drop in a part of the pipe occupied by the steam. The pressure difference in the parts of the pipe occupied by the liquid phase and the parts occupied by the vapor causes the liquid column movement and acceleration towards the area occupied by the vapor. Consequently, the vapor condensation continues towards the propagating liquid and vapor interface. A liquid column accelerates and eventually hits the obstacle, such as a valve, the closed end of a pipe or another liquid column, and reflects from the obstacle. A sudden pressure pulse, generated at the moment of impact of a liquid column onto an obstacle and the consequent propagation of pressure waves, can cause severe mechanical damage to the equipment, such as damage to pipe walls, fittings and hangers or pressure vessels, and might endanger the safety and cause serious injury to operating personnel. The dangerous effects of CIWH were shown by Milivojevic et al. [10] with simulations of destructive pressure peaks greater than 10 MPa in systems that were initially at a low pressure close to the atmospheric.

With a more precise insight into the mechanism of the occurrence of water hammer, it is possible to improve the protection systems, implement necessary safety measures and thus prevent its consequences. Some of the facilities where this undesirable thermal–hydraulic phenomenon can occur are steam power plant units [11], nuclear power plants (NPPs) [12, 13, 14], district heating systems [15] and the ammonia refrigeration system [16]. In the previous research, whether experimental or numerical, the goal was mostly to record the highest peak pressure value due to safety reasons. An experimental and analytical investigation was performed to estimate the impulse generated during the large steam bubble collapse in a vertical pipe between the lower stagnant hot water column and the upper downward accelerating column of cold water [17]. The liquid column was observed as a rigid body, and a simple mechanical model was derived for the prediction of the water column velocity and impulse at the moment of impact and the resulting pressure peak. The CIWH in a vertical pipe filled with steam and closed at the top was reported in [18]. In this study, the pipe filled with steam is immersed in the reservoir of subcooled water. The direct contact of steam and subcooled water is caused by fast opening valve at the bottom of the pipe. The experiments showed cases with great pressure pulses, from initial several bars to approximately 16 MPa. These experimental conditions were simulated with the thermal-hydraulic nuclear reactor safety code TUF [12]. The CIWH simulations required the TUF code upgrading with the model for the steam-water interfacial area concentration at the water column head. It was assumed that the water column and steam interface consist of certain bubbles and droplets, which are formed during the interface movement. This interfacial area was kept constant during the water column movement.

The safety of NPP steam generators and feedwater systems is of great importance. Consequently, CIWH caused by countercurrent flow of steam and subcooled water in the horizontal pipe, or in pipes with small horizontal inclination, is most common, and also the most experimentally investigated in the latest decade [1, 2, 19, 20, 21, 22], as well as earlier [23, 24, 25]. In such systems, stratified flow occurs, followed by steam bubbles formation and the rapid condensation of steam bubble entrapped by the subcooled water. Some guidelines are given for the prevention of CIWH in a two-phase flow. Barna and Ezsol [25] performed the water hammer experiment in the countercurrent flow of steam and water in the horizontal pipe in the Hungarian PMK-2 facility and also in the Japanese Rig-of-safety Assessment (ROSA) facility. In addition, they performed numerical simulation by WAHA3 code for transients in single- and two-phase flow. Large system codes like RELAP5, Trac and CATHARE are used for the safety analysis of transients in NPP, and they model two-phase flow transients. Numerical simulations of CIWH are performed by using large system codes as ATHLET, WAHA3 [25], RELAP5 [20], Trac, CATHARE, in-house codes [26] and OpenFOAM, open source CFD software [21]. A type of the CIWH known as the water cannon phenomenon can occur during the discharge of steam into a large volume of cold water. In case when the exhaust valve is closed, steam can be found trapped in a pipeline. The rapid condensation of steam on cold water is followed by subcooled water suction into the pipe. The formed water slug eventually hits the valve and causes a large pressure pulse. Yeung et al. [27] performed a simulation of the water cannon phenomenon by nuclear reactor safety code RELAP5/MOD3. Dynamics of CIWH depends on complex thermal and hydraulic effects, such as the steam condensation on subcooled water, the disintegration of the water column head and the droplets entrainment from the water column head to the steam space and the transient friction of the water column. Experimental investigations have shown that practically the same experimental conditions can result in significant scattering of the experimental test data, as it was reported in [12, 22, 28]. This scattering of the results is caused by the stochastic nature of the water column head disintegration. The steam–cold water contact area affects the way the condensation at the water column head develops which further depend on the acceleration of the water column and the peak pressure at the impact during CIWH. The entrainment of water droplets from the accelerated water column head and the condensation of steam onto the droplets and the disturbed water column head were investigated experimentally in [29]. The Taylor instability of the accelerated water column head was observed. The relative velocity of the entrained droplets and the water column head was correlated with the column head acceleration.

Hou et al. [26] experimentally and numerically investigated flashing instability induced water hammer (FIIWH), which belongs to the class of CIWH, and occurs in open natural circulation systems (NCS). The occurrence of FIIWH is more likely in NCS with long horizontal pipes. NCS are widely used in NPPs as the components of the passive safety systems. Sun et al. [19] experimentally investigated CIWH in NCS, regarding flow parameters and fluid temperature. Since they used tap water as a medium, in which there are dissolved non-condensable gases, when a CIWH occurred, these gases were released. The volume of non-condensable gases did not change during the experiment, which shows that all the vapor in the bubbles, which occurred after CIWH, completely evaporated. The presence of non-condensable gases in the entrapped bubble reduces the heat transfer coefficient and lowers condensation rate, so that the speed of the water columns is lower, which results in smaller pressure peaks when they collide. A prerequisite for the formation of CIWH is the existence of an entrapped bubble, and the degree of subcooling of the water determines whether or not a CIWH will occur. The steam in the entrapped bubble will condense in contact with surrounding subcooled water, the pressure in the bubble will suddenly drop, then the water at the bottom of the bubble will instantly evaporate, and steam condensation induced flashing (SCIF) will occur. The SCIF occurs when water on one side of the entrapped bubble is saturated and on the other is subcooled. The SCIF diminishes the pressure drop in the bubble due to condensation and slows down the occurrence of CIWH. After CIWH in NCS water column moves in the opposite direction of the initial flow. Sun et al. [19] conducted 67 CIWH events, and they identified three types of CIWH in NCS. Type I CIWH occurs when there is a large difference in velocity of water and steam, slip ratio. The velocity change results in pressure change which increases the disturbances at the interface. When pressure change prevails over gravity and surface tension, Kelvin-Helmholz instability occurs. Type II CIWH is caused by the solitary wave formation and reverse flow of subcooled water into the pipe. Interaction of solitary and interface wave is also important. Type III CIWH is induced by the pressure wave generated after the occurrence of CIWH. The most probable type is Type II CIWH, and the least probable is Type III.

The transient friction strongly affects water column acceleration during CIWH, especially in the vicinity of the pressure wave front, because it affects the propagation speed of the wave front and the evolution of its shape. By modeling vapor cavitation in transient fluid flow, Shu [30] concluded that by applying the unsteady friction model, better damping is achieved, in cases where the cavitation is not too strong and when the assumption of laminar single-phase flow with the appropriate weighting function is applicable. An unsteady turbulent skin friction model for one-dimensional smooth pipes, where the wall shear stress is proportional to the fluid instantaneous acceleration, is developed by Vardy and Brown [31]. Errors resulting from the numerical simulation of friction in unsteady flows in small pipe networks were evaluated in [32] in the case of applying the MOC with a fixed node distance in the numerical grid and without interpolation. The application of the friction relaxation model to the prediction of the water hammer in two-phase flow is presented in [33].

The dynamics of deaeration and the fluid structure interaction (FSI) after the occurrence of the first pressure pulse are also effects that are distinguished as very influential on CIWH. Neuhaus and Dudlik in [34] experimentally investigated the air release process. The FSI model was implemented into fluid dynamic equations. Taking into account the air release as well as the FSI model predictions are improved.

In order to improve the process safety and to prevent equipment damage and accidents, it is very important to single out the thermal–hydraulic conditions that can lead to CIWH, as well as the maximum pressure peaks and periods of generated pressure pulses during CIWH. The dynamics of generated pressure pulses during CIWH event is predominantly governed by the intensity of the condensation rate. Direct contact of subcooled liquid and steam leads to intensive condensation. The condensation rate is highly influenced by the interfacial area concentration and condensation heat transfer coefficient. The interfacial area has a very irregular shape. During CIWH transient liquid jets are formed, and liquid column is disintegrated so that droplets from the head are entrained in steam volume. This droplet entrainment increases interfacial area and condensation rate. In the available literature, there is a lack of information on the determination of the condensation rate. Liu presented in [12] model of condensation during CIWH event, but without using values for interfacial area concentrations of steam and liquid droplets and also of liquid and steam bubbles in close to the water column head. In studies by [27, 28, 33] the rate of condensation was only mentioned without providing further details. Yeung et al. [27] referred that a combination of mechanistic models and experimental correlations was used in the calculation using the RELAP code. Barna et al. [28] mentioned the quick condensation model in WAHA code, but no detailed information was provided. Kucienska et al. [33] outlined the use of the heat and mass transfer model for dispersed flow. The model is derived from the homogeneous relaxation approach, assuming a large heat transfer coefficient at the vapor side of phase interface. No further details about the model application are given.

The next chapter presents our own numerical model for the prediction of CIWH and for the prediction of water hammer with gaseous and vaporous cavitation.

2. Modeling approach

2.1 Governing equations

Single-phase vapor or liquid and two-phase vapor-liquid flow in a pipe is considered as the one-dimensional, transient and compressible fluid flow of homogeneous fluid. The velocity and thermal equilibrium are assumed between the vapor and liquid phase in the two-phase flow. The following mass, momentum and energy conservation equations are applied:

• Mass conservation

  $$\frac{\mathit{D\rho }}{\mathit{Dt}}+\rho \frac{\partial u}{\partial x}=0$$E1

• Momentum conservation

  $$\frac{\mathit{Du}}{\mathit{Dt}}+\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mathit{fu}\left|u\right|}{2d_{\mathrm{H}}}+f_{\mathrm{u}}sign\left(u\right)\left|\frac{\partial u}{\partial t}\right|+gsin\theta =0$$E2

• Energy conservation

  $$\frac{\mathit{Dh}}{\mathit{Dt}}-\frac{1}{\rho }\frac{\mathrm{D}p}{\mathrm{D}t}-\frac{{\mathit{fu}}^2\left|u\right|}{2d_{\mathrm{H}}}-f_{\mathrm{u}}usign\left(u\right)\left|\frac{\partial u}{\partial t}\right|+\underset{x_{\mathrm{i}}-\epsilon }{\overset{x_{\mathrm{i}}+\epsilon }{\int }}\mathit{sign}\left(h-h^{\text{'}}\right)\frac{{\mathrm{\Gamma }}_c{h}^{"}}{\rho }\mathit{dx}=0$$E3

where the dependent variables are velocity u, pressure p and specific enthalpy h, and independent variables are time t and spatial coordinate x. The hydraulic diameter is denoted with $d_{\mathrm{H}}$, the pipe inclination angle from the horizontal axis with $\theta$, gravity acceleration with g, the Darcy friction coefficient with f and the unsteady friction coefficient with f_u. The third and fourth terms on the left-hand side of momentum and energy conservation equations, Eqs. (2) and (3), are related to pressure drop due to steady-state and transient friction, respectively. The last term on the left-hand side of the energy Eq. (3) determines interfacial heat transfer per unit mass of fluid due to condensation.

2.2 Closure laws for condensation induced water hammer

The pressure drop due to friction is determined using Darcy friction coefficient f. The unsteady friction coefficient is determined according to [31], as

$$f_{\mathrm{u}}=2\sqrt{\frac{12.86}{{\mathbf{Re}}^{\left(1.1844-0.0567{log}_{10}\mathbf{Re}\right)}}}$$E4

where Re is the Reynolds number.

The condensation rate in Eq. (3) is determined as

$${\mathrm{\Gamma }}_c=\frac{q_c{a}_i}{r}$$E5

where r is the latent heat of condensation, and the condensation heat flux $q_c$ is determined as the product of the condensation heat transfer coefficient $h_c$ and the difference between the saturation temperature $T_{\mathit{sat}}$ and the subcooled liquid temperature $T_1$

$$q_c=h_c\left(T_{\mathit{sat}}-T_1\right).$$E6

The condensation takes place from vapor or vapor-liquid two-phase mixture to the liquid column at the interface, whose position is denoted with x_i (Figure 1). Parameter ε represents infinitesimal distance from the interface to the subcooled liquid and from the interface to the two-phase mixture or vapor. The transfer of thermal energy of condensation through the interface from the two-phase mixture or vapor to the liquid column is determined by the function $\mathit{sign}\left(h-h^{\text{'}}\right)$.

Figure 1.

Vapor-liquid interface position and the direction of the condensation process.

The prediction of interfacial area concentration $a_i$ is the main task when determining the condensation rate, besides determining the condensation heat transfer coefficient. Throughout the CIWH event transient liquid droplets or jets separate from the liquid column head and entrain in the vapor volume, so that the interfacial area increases, as well as the condensation rate. The developed model assumes that disintegration process of liquid column head, the entrainment of droplets into vapor and the intensity of the heat transfer at the interface depend on the liquid column head acceleration. The product of the condensation heat transfer coefficient and the interfacial area concentration is divided to two additions, first related to liquid column head (LCH) and the other to entrained droplets (ED)

$$h_c{a}_i={\left(h_c{a}_i\right)}_{\mathit{LCH}}+{\left(h_c{a}_i\right)}_{\mathit{ED}}.$$E7

Bloemeling [35] used correlations based on surface renewal theory to predict the condensation of steam at the turbulent liquid column head, but greater uncertainty occurs in predicting turbulent characteristic length and velocity. These turbulent scales are correlated with turbulent kinetic energy and dissipation rate, and accuracy of their calculations requires application of 3D models of steam–water interface behavior in transient conditions. In this chapter, the condensation heat transfer coefficient at the LCH is calculated using Dittus and Boelter (1930) correlation, i.e., turbulent heat transfer at the LCH is considered as convective heat transfer at the pipe wall

$${\left(h_c\right)}_{\mathit{LCH}}=\left({\lambda }_1/d_H\right)0.023{Re}_1^{0.8}{Pr}_1^{0.4}$$E8

where index 1 denotes liquid, $\lambda$ is the thermal conductivity, and Re and Pr are Reynolds and Prandtl number, respectively. The interfacial area concentration at the LCH is calculated with the assumption that it is equal to the cross section of the pipe

$${\left(a_i\right)}_{\mathit{LCH}}=\frac{A_i}{A_i\mathrm{\Delta }x}=\frac{1}{\mathrm{\Delta }x}.$$E9

The condensation heat transfer coefficient on the ED is determined from the Nusselt number, whose constant value is theoretically predicted and experimentally confirmed [36, 37, 38] as

$$\mathrm{Nu}=\frac{{\left(h_c\right)}_{\mathit{ED}}{d}_D}{{\lambda }_1}=C.$$E10

The adopted constant value for C is 17.9 according to Kronig and Brink [36]. The entrained droplet diameter depends on parameter Y

$$Y=\frac{\sigma {\text{We}}_{\mathrm{cr}}}{{\rho }_2{\left(u_2-u_1\right)}^2}$$E11

where $\sigma$ is the surface tension, index 2 denotes superheated vapor, and the value of the critical Weber number We_cr = 0.799 is taken from [39]. The entrained droplet diameter is calculated as

$$d_D=\left\{\begin{array}{l}{10}^{-4}\mathrm{m},Y\le {10}^{-4}\mathrm{m}\\ Y,{10}^{-4}\mathrm{m}<Y<5\cdot {10}^{-4}\\ 5\cdot {10}^{-4}\mathrm{m},Y\ge 5\cdot {10}^{-4}\mathrm{m}\end{array}\right.\mathrm{m}.$$E12

The major challenge here exists in determination of the interfacial area concentration between entrained droplets from the liquid column head and the vapor. It is assumed that the interfacial area concentration depends on the acceleration of the liquid column head as follows

$${\left(a_i\right)}_{\mathit{ED}}=a_{D,0}+4\cdot {10}^{-3}{\left(\mathit{Du}/\mathit{Dt}\right)}^3$$E13

Where a_D,0 takes the values between 0 and 40 m²/m³. The determination of the parameters on the right-hand side of Eq. (13) is performed by comparing the results of numerical simulations of CIWH with experimental data available in [12, 18] and data in [27]. The term a_D,0 is determined by dynamics of propagation of the liquid column towards the vapor volume, and it depends on the initial conditions under which CIWH occurs, e.g., the fast opening of the valve or removing the obstacle between superheated liquid and vapor. The formation of the interfacial area between liquid and steam during the opening of the valve is stochastic in nature. Interfacial area is very irregular with entrained droplets in vapor and formation of liquid jets, and therefore, a_D,0 is in the range between 0 and 40 m²/m³. If its value is 0, that means that there are no entrained droplets which significantly increase the interfacial area. The second term in Eq. (13) takes into account the influence of the acceleration of the liquid column head on the disintegration of the column head and the droplet entrainment in the vapor volume.

The thermodynamic quality is used to determine the phase state of the fluid

$$x_t=\frac{h-h^{\text{'}}}{r}$$E14

where for $x_t\le 0$ liquid phase takes place, for $x_t\ge 0$ vapor phase, and for $0<x_t<1$ a two-phase mixture. The fluid density is calculated as a reciprocal value of the specific volume determined by applying the equations of state for subcooled liquid $v_1\left(p,h\right)$, superheated vapor $v_2\left(p,h\right)$, saturated liquid $v^{\text{'}}\left(p\right)$ and saturated vapor $v^{\text{'}\text{'}}\left(p\right)$ (data from the steam tables [40] are used)

$$\rho =\left\{\begin{array}{l}1/v_1\left(p,h\right),\text{if}{x}_t\le 0\\ 1/\left(v^{\text{'}}\left(p\right)+x_t\left(v^{\text{'}\text{'}}\left(p\right)-v^{\text{'}}\left(p\right)\right)\right),\text{if}0<x_t<1\\ 1/v_2\left(p,h\right),\text{if}{x}_t\ge 1\end{array}\right..$$E15

2.3 Closure laws for water hammer with gaseous and vaporous cavitation

Liquid single-phase and gas/vapor-liquid two-phase flow in a pipe is observed as the one-dimensional, transient and compressible fluid flow of homogeneous fluid. The velocity equilibrium is assumed between the gas/vapor and liquid phase in the two-phase flow. The flow is isenthalpic. Mass and momentum conservation equations Eqs. (1) and (2) are applied, where the unsteady friction coefficient is calculated according to Eq. (4) and the fluid density is evaluated by Eq. (15).

The reasonable assumption is that the liquid density is constant and determined by the initial liquid temperature and the initial average pressure in the flow channel. The gas/vapor density is calculated with the ideal gas law

$${\rho }_2=\frac{p_2}{R_g{T}_2}$$E16

where $T_2$ is the gas/vapor temperature approximated with the initial liquid temperature $T_{1,0}$, while $R_g$ is the gas constant. In the case of bubbly flow, the gas/vapor pressure is calculated by taking into account the surface tension σ, according to the Laplace eq. [38]

$$p_2=p_1+\frac{2\sigma }{r_b}.$$E17

Assuming that there is no heat exchange between the flow channel and the surroundings and that the single-phase flow and homogeneous two-phase flow are isenthalpic, gas/vapor density only depends on pressure

$${\rho }_2={\rho }_2\left(p_2\right).$$E18

The infinitesimal change in gas/vapor quality in the case of water hammer with gaseous or vaporous cavitation is determined as [41, 42, 43].

$$\frac{\mathit{dx}}{\mathit{dt}}=\frac{\mathrm{\Gamma }}{\rho }$$E19

where $\mathrm{\Gamma }$ represents the rate of interfacial mass transfer of gas/vapor (due to evaporation and condensation in case of vaporous cavitation and due to absorption and desorption in case of non-condensable gaseous cavitation) per unit volume and per unit of time, and $\rho$ is the gas/vapor-liquid two-phase mixture density. The rate of interfacial mass transfer of gas/vapor is determined as

$$\mathrm{\Gamma }=j_i{a}_i$$E20

where $a_i$ is the gas/vapor-liquid interfacial area concentration, and $j_i$ is the interfacial mass transfer flux of vapor or non-condensable gas.

The relation between the gas/vapor void fraction and quality is [38].

$${\alpha }_2=\frac{1}{1+\frac{1-x}{x}\frac{{\rho }_2}{{\rho }_1}}.$$E21

The calculation of gas/vapor-liquid interfacial area concentration depends on gas/vapor quality and the two-phase flow pattern. It is assumed that for voids lower than or equal to 0.74 the flow pattern is bubbly [44], while for voids greater than 0.74 the annular pattern is assumed. The interfacial area concentration in the bubbly flow is calculated as [37].

$$a_i=\frac{6{\alpha }_2}{d_b}$$E22

with the assumption that the bubbles are spheres with uniform diameter $d_b$. The bubble diameter is calculated from the relation

$${\alpha }_2=\frac{n_b\pi d_b^3}{6}$$E23

where the gas/vapor void fraction is related to the bubble concentration $n_b$. The interfacial area concentration for the annular flow pattern is calculated from the following expression by assuming that the liquid phase wets the tube wall as liquid film [37].

$$a_i=\frac{4}{d}\sqrt{{\alpha }_2}$$E24

where $d$ is the tube inner diameter.

The interfacial mass transfer flux of vapor on the surface of bubble, in case of vaporous cavitation, is

$$j_i={\rho }_2\frac{{\mathit{dr}}_b}{\mathit{dt}}.$$E25

where ${\rho }_2$ denotes vapor density, and $r_b$ is the bubble radius. The time change of the bubble radius is described by the correlation [45].

$$\frac{{\mathit{dr}}_b}{\mathit{dt}}=\frac{{\lambda }_1}{{\rho }_{2}r}\left(\frac{1}{{\left(a_{1}t\right)}^{1/2}}+\frac{1}{r_b}\right)\left(T_1-T_{\mathit{sat}}\right).$$E26

where ${\lambda }_1$ denotes thermal conductivity of liquid, $r$ is latent heat, $a_1$ is the thermal diffusivity of liquid, $T_1$ and $T_{\mathit{sat}}$ are, respectively, liquid temperature and saturation temperature. During the bubble growth, the bubble radius $r_b$ is much greater than value of ${\left(a_{1}t\right)}^{1/2}$so that the equation for time change of the bubble radius can be reduced to

$$\frac{{\mathit{dr}}_b}{\mathit{dt}}=\frac{{\lambda }_1}{{\rho }_{2}r{\left(a_{1}t\right)}^{1/2}}\left(T_1-T_{\mathit{sat}}\right).$$E27

The interfacial mass transfer rate of vapor in bubbly flow is determined by introducing Eq. (25) and Eq. (22) into Eq. (20) and then introducing the bubble diameter from Eq. (23) and the time change of the bubble radius from Eq. (27) into obtained relation. The following expression is derived

$$\mathrm{\Gamma }=4.835\frac{{\alpha }_2^{2/3}{n}_b^{1/3}{\lambda }_1}{r{\left(a_{1}t\right)}^{1/2}}\left(T_1-T_{\mathit{sat}}\right).$$E28

In the case of annular flow, the evaporation and condensation interfacial mass flux at the liquid film surface is determined from the surface heat flux equation for heat conduction in semi-infinite solid with the initial temperature $T_1$, that is suddenly lowered and maintained at a temperature $T_{\mathit{sat}}$ as

$$q=\frac{{\lambda }_1}{{\left(\pi a_{1}t\right)}^{1/2}}\left(T_1-T_{\mathit{sat}}\right).$$E29

By dividing Eq. (29) with latent heat, the evaporation and condensation interfacial mass flux is obtained

$$j_i=\frac{{\lambda }_1}{r{\left(\pi a_{1}t\right)}^{1/2}}\left(T_1-T_{\mathit{sat}}\right).$$E30

The interfacial mass transfer rate of vapor in annular flow is determined by introducing Eq. (24) and Eq. (30) into Eq. (20) as

$$\mathrm{\Gamma }=\frac{4{\lambda }_1{\alpha }_2^{1/2}}{\mathit{rd}{\left(\pi a_{1}t\right)}^{1/2}}\left(T_1-T_{\mathit{sat}}\right).$$E31

The absorption and desorption rate of non-condensable gas in liquid, in case of gaseous cavitation, is determined in the following manner. If non-condensable gas is in contact with liquid for a long time period, the liquid becomes saturated with gas, and the equilibrium condition is reached. The molar fraction of absorbed non-condensable gas in liquid ${\tilde{x}}_{\mathit{eq}}$ is then determined by Henry’s law [4].

$${\tilde{x}}_{\mathit{eq}}=\frac{p_2}{H_c}$$E32

where $p_2$ is non-condensable gas pressure, and $H_c$ is Henry’s constant, which depends on the type of gas, as well as the pressure and temperature of the liquid.

Under non-equilibrium conditions, when the mass fraction of the dissolved gas in the liquid is different from the equilibrium gas saturation in the liquid, the gas transfer occurs at the interface between the liquid and gas phase. Assuming that the relative velocity between the gas and liquid phase is negligible in the case of water hammer, the convective mechanisms of gas transfer are neglected, and interface gas transfer is determined by the gas diffusion on the liquid side of the interface. In bubbly flow, which occurs during water hammer with gas cavitation, small bubbles are dispersed in the liquid, gas–liquid relative velocity can be neglected. The mass balance equation of the non-condensable gas in liquid in the coordinate system connected with the moving boundary, for isothermal absorption or desorption [46], is

$$\frac{\partial C}{\partial t}+u_i\frac{\partial C}{\partial x}=D_1\frac{{\partial }^{2}C}{\partial x^2}$$E33

where $C$ is the mass fraction of dissolved non-condensable gas in liquid, $D_1$ is the liquid diffusivity, $x$ is the coordinate, and $u_i$ is the boundary displacement velocity defined by the condition of interface impenetrability for the liquid phase

$$j_{i,1}=\left(1-C_i\right){\rho }_1{u}_i-\rho D_1\frac{\partial \left(1-C\right)}{\partial x}=0$$E34

as

$$u_i=-\frac{D_1}{1-C_i}\frac{\partial C}{\partial x}\left(\mathrm{at}\mathrm{x}=0\right)$$E35

and $j_{i,1}$ in Eq. (34) is the mass flux of the liquid phase which contains convection and diffusion components. The solution to Eq. (33), for boundary displacement velocity defined by Eq. (35), with following boundary conditions: if $x=0$, $C=C_i$ and when $x\to \infty$, $C=C_1$ is

$${\left.\frac{\partial C}{\partial x}\right|}_{x=0}=-\frac{\left(C_i-C_1\right)}{\sqrt{\pi D_{1}t}}.$$E36

The interfacial mass transfer flux of non-condensable gas

$$j_i=-\frac{{\rho }_1\sqrt{D_1}\left(C_i-C_1\right)}{\sqrt{\mathit{\pi t}}}$$E37

is determined by the difference between the gas mass fraction in the liquid $C_1$ and the gas mass fraction at the interface $C_i$. The interface gas mass fraction, if the liquid temperature is constant, is determined as a function of pressure in following form

$$C_i=k_1+k_{2}p$$E38

where coefficients $k_1$ and $k_2$ depend on water temperature. Eq. (38) is applicable for air absorption in water and for pressures from low vacuum to several bars. According to Eq. (38), the pressure of air saturation in water is calculated as

$$p_{\mathit{sat}}=\frac{\left(C_1-k_1\right)}{k_2}.$$E39

The value of the liquid pressure determines whether absorption or desorption takes place, if it is lower than the pressure of gas saturation in liquid, desorption occurs, and if it is higher, then absorption occurs.

By including Eqs. (22), (23) and (37) in Eq. (20), the interface mass transfer rate of non-condensable gas in bubbly flow can be calculated as

$$\mathrm{\Gamma }=-\frac{6^{2/3}{\rho }_1{D}_1^{1/2}{n}_b^{1/3}{\alpha }_2^{2/3}}{{\pi }^{1/3}{t}^{1/2}}\left(C_i-C_1\right)$$E40

and in case of annular flow, by including Eqs. (24) and (37) in Eq. (20), follows

$$\mathrm{\Gamma }=-\frac{4{\rho }_1{D}_1^{1/2}{\alpha }_2^{1/2}}{{\left(\mathit{\pi t}\right)}^{1/2}d}\left(C_i-C_1\right).$$E41

Taking into account that the change in the non-condensable gas mass fraction in the liquid is equal to the change in non-dissolved gas quality in the two-phase gas-liquid mixture, ${\mathit{dC}}_1=-\mathit{dx}$ from Eq. (19) follows

$$\frac{{\mathit{dC}}_1}{\mathit{dt}}=-\frac{\mathrm{\Gamma }}{\rho }.$$E42

Since the bubbly flow in air-water mixture occurs during water hammer with gaseous cavitation [6, 47], the Eq. (40) is used for the estimation of the maximum value of air mass transfer during gaseous cavitation in [48].

In technical systems with degassed water, there is a certain amount of non-condensable gas dissolved in the liquid, while a small amount of gas is undissolved and is dispersed in liquid in the form of microbubbles. When, during the transient, the rarefaction wave propagates, the volume of microbubbles increases due to the pressure drop. After the passage of the rarefaction wave, the gas microbubbles can interact with each other and they may merge. The dynamic change of the bubble radius is expressed by the Rayleigh-Plesset equation, obtained in [48],

$$\frac{p_b-p_1}{{\rho }_1}=r_b\frac{d^2{r}_b}{{\mathit{dt}}^2}+\frac{3}{2}{\left(\frac{{\mathit{dr}}_b}{\mathit{dt}}\right)}^2+\frac{4{\mu }_1}{{\rho }_1{r}_b}\frac{{\mathit{dr}}_b}{\mathit{dt}}+\frac{2\sigma }{{\rho }_1{r}_b}$$E43

where $p_b$ is pressure in the gas bubble, $p_1$ is the water pressure, $r_b$ is the bubble radius, $\sigma$ is the surface tension, ${\mu }_1$ is dynamic viscosity of water, and ${\rho }_1$ is density of water. A detailed numerical solution, for air bubble growth in water at room temperature under a sudden pressure drop, is presented in [48], as well as the analysis of the influence of the initial bubble radius and the pressure drop value on the bubble growth dynamics. It is concluded that the rate of mass transfer is significantly higher in the initial period of gaseous cavitation than during the rest of the transient. This period is approximately equal to the time steps of the numerical integration of governing differential equations. For instance, in the case of the MOC, the time step of integration is determined by the Courant criterion, see Eq. (55) below. The time step of integration of the governing balance equations is no smaller than 10⁻⁴ s, which is of the same order of magnitude as the time period of inertial bubble growth. With regard to these effects, for the first step of integration of the balance equation when the liquid pressure drops below the pressure of gas saturation in liquid, the rate of mass transfer is calculated with the following empirical equation

$$\mathit{dx}=-{\mathit{dC}}_1=\frac{\mathrm{\Gamma }}{\rho }\mathit{dt}=-k\left(C_i-C_1\right)$$E44

where the value of coefficient k = 0.7735 is determined in [48], by comparing the measured and calculated pressure changes during the transient with gaseous cavitation. After this first step of integration, it is assumed that the finite number of air bubbles is formed, which results in the bubble number density $n_b\approx {10}^8$, while Eq. (40) is later applied for the mass transfer rate calculation.

The pipes wall elasticity is taken into account in determination of the speed of pressure wave propagation

$$c=a{\left(1+\frac{\rho a^2}{E}\frac{d}{\delta }\right)}^{-1/2}$$E45

where E is Youngs’s modulus of elasticity of the pipe, $\delta$ is the pipe wall thickness, and a is the sonic velocity. In case of the single-phase water flow,

the sonic velocity [38] is calculated as

$$a={\sqrt{\left(\frac{\mathit{dp}}{\mathit{d\rho }}\right)}}_s$$E46

and in the homogeneous two-phase flow without evaporation or condensation, the sonic velocity is calculated from the so-called frozen sonic velocity expression [49]

$$a={\left[\rho \left(\frac{{\alpha }_2}{{\rho }_2{c}_2^2}+\frac{1-{\alpha }_2}{{\rho }_1{c}_1^2}\right)\right]}^{-\frac{1}{2}}$$E47

where $c_1$ and $c_2$ are the sonic velocities in liquid and vapor, respectively, and they are calculated using Eq. (46), for known values of thermodynamic parameters of liquid and vapor.

Presented closure laws for boiling and condensation are incorporated into the hydrodynamic model presented in Section 2.1, and the results obtained with the numerical solutions show an interplay of hydrodynamic effects of pressure waves generation and propagation and interfacial mass and energy transfer by phase transition and diffusion. It is interesting to note that an interaction of the thermal and hydrodynamic effects is observed even at the nanoscale level by the molecular dynamic simulations of boiling [50].

2.4 Numerical solution

Governing equations Eqs. (1)–(3) are transformed, including application of material derivative, introducing sonic velocity in one-phase fluid, and thus in homogeneous two-fluid model, as a function of pressure and enthalpy

$$c={\left({\left(\frac{\partial \rho }{\partial p}\right)}_h+\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial h}\right)}_p\right)}^{-1/2}$$E48

and grouping of all partial differentials over time and coordinate on the left-hand side of equations, so that the system of quasi-linear hyperbolic partial differential equations is obtained in following form

$$\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}+c^2\rho \frac{\partial u}{\partial x}=X$$E49

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=Y$$E50

$$\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}-\frac{1}{\rho }\left(\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}\right)=Z$$E51

where

$$X=-c^2{\left(\frac{\partial \rho }{\partial h}\right)}_p\left(\frac{{\mathit{fu}}^2\left|u\right|}{2d_{\mathrm{H}}}+f_{\mathrm{u}}usign\left(u\right)\left|\frac{\partial u}{\partial t}\right|\right.\left.-\mathit{gu}sin\theta -\underset{x_{\mathrm{i}}-\epsilon }{\overset{x_{\mathrm{i}}+\epsilon }{\int }}\mathit{sign}\left(h-h^{\text{'}}\right)\frac{{\mathrm{\Gamma }}_c{h}^{"}}{\rho }\mathit{dx}\right)$$E52

$$Y=-\frac{\mathit{fu}\left|u\right|}{2d_{\mathrm{H}}}-f_{\mathrm{u}}sign\left(u\right)\left|\frac{\partial u}{\partial t}\right|-gsin\theta$$E53

$$Z=\frac{{\mathit{fu}}^2\left|u\right|}{2d_{\mathrm{H}}}+f_{\mathrm{u}}usign\left(u\right)\left|\frac{\partial u}{\partial t}\right|-\underset{x_{\mathrm{i}}-\epsilon }{\overset{x_{\mathrm{i}}+\epsilon }{\int }}\mathit{sign}\left(h-h^{\text{'}}\right)\frac{{\mathrm{\Gamma }}_c{h}^{"}}{\rho }\mathit{dx}.$$E54

The system of equations, Eq. (49)–(51), is solved for the appropriate initial and boundary conditions by the MOC. Initial conditions are defined by thermal–hydraulic characteristics of the fluid at the initial time, before the disturbance happens. The boundary conditions are determined by the state of the fluid at the beginning and end of the sections, as well as at the boundaries of the observed system. Three characteristic paths are used, where two correspond to the pressure wave propagation (C⁺ and C⁻, Figure 2) and the third to the propagation of the fluid particle enthalpy front (C^P, Figure 2). The spatial step of integration $\mathrm{\Delta }x$, that is the distance between two adjoining nodes, is constant within one pipe segment. The time step of integration is determined according to the Courant criterion:

Figure 2.

Spatial coordinate (x)–time (t) plane and characteristic paths (in case of the fluid particle velocity u_P > 0).

$$\mathrm{\Delta }t\le min\left(\frac{\mathrm{\Delta }x}{c_{i,j}+\left|u_{i,j}\right|}\right),i=1,2,\dots N\left(j\right),j=1,2,\dots ,M$$E55

where index i denotes the node within the pipe j.

Conservation equations Eq. (49)–(51) are transformed into the system of ordinary differential equations along three characteristic paths.

$$\mathit{dp}+\mathit{\text{ρcdu}}=\left(X+\mathit{\rho cY}\right)\mathit{dt}\text{along}{C}^{+}:\frac{\mathit{dt}}{\mathit{dx}}=\frac{1}{u+c}$$E56

$$\mathit{dp}-\mathit{\text{ρcdu}}=\left(X-\mathit{\rho cY}\right)\mathit{dt}\text{along}{C}^{-}:\frac{\mathit{dt}}{\mathit{dx}}=\frac{1}{u-c}$$E57

$$\mathit{dh}-\frac{1}{\rho }\mathit{dp}=\mathit{Zdt}\text{along}{C}^P:\frac{\mathit{dt}}{\mathit{dx}}=\frac{1}{u}.$$E58

Approximation of total derivatives by finite differences along the characteristic paths transforms differential equations, Eqs. (56)–(58), into difference equations.

A, B and C in Figure 2 denote three successive nodes in the flow channel, which are used for calculation of initial values of the dependent variables at time level t. The pressure and velocity values in points R and L are calculated by the linear interpolation of their initial values in A, B and C. The point D denotes the intersection point of all characteristic paths C⁺, C⁻ and C^P and the x-axis at time level t + Δt. Hence, D denotes the node where disturbance arrives in the next time level, t + Δt. Coordinate x_P is determined using the slope of the characteristic path C^P and the linear interpolation of the velocity between nodes A and B for positive flow direction and between nodes B and C for negative flow direction, Figure 3.

Figure 3.

Determination of nodes for the Lagrange’s interpolation polynomial of the third degree, for positive (u_P > 0) and negative (u_P < 0) flow direction.

The Lagrange’s interpolation polynomial (LIP) of the third degree is used for determination of the initial enthalpy value in point P with the aim of reducing the numerical diffusion of the enthalpy front propagation. For the purpose of derivation of LIP of the third degree it is necessary to use the enthalpy values at four nodes. The choice of nodes depends on the flow direction. It is always necessary to use the enthalpy value in one node downstream, in two nodes upstream, as well as in the observed node (designated with i in Figure 3). The following calculation algorithm, suitable for the computer programming of LIP, is applied:

$$L_m\left(x\right)={\mathrm{\Pi }}_{m+1}\left(x\right)\sum _{j=0}^m\frac{h_j}{D_j}$$E59

where

$${\mathrm{\Pi }}_{m+1}\left(x\right)=\left(x-x_0\right)\dots \left(x-x_m\right)$$E60

$$D_j=\left(x_j-x_0\right)\left(x_j-x_1\right)\dots \left(x_j-x_{j-1}\right)\left(x-x_j\right)\left(x_j-x_{j+1}\right)\dots \left(x_j-x_m\right),\mathrm{j}=0,1,\dots ,\mathrm{m}.$$E61

The error of the interpolation with the LIP is [51].

$$\left|h\left(x\right)-L_m\left(x\right)\right|\le \frac{M_{m+1}}{\left(m+1\right)!}\left|{\mathrm{\Pi }}_{m+1}\left(x\right)\right|$$E62

and

$$M_{m+1}=\underset{a\le x\le b}{max}\left|f^{\left(m+1\right)}\left(x\right)\right|.$$E63

The use of LIP of the third degree to determine the enthalpy at the point P gives the truncation error of the fourth order O[(Δx)⁴] for the numerical discretization of the enthalpy along the x coordinate. The integration of the energy equation Eq. (58) with respect to time is performed along the characteristic path C^P with the Euler explicit method, which gives truncation error of the first order O(Δt).

3. Results and discussion

The developed model and computer program of transient compressible fluid flow is applied for calculation of thermal-hydraulic parameters during the CIWH event. The validation of the developed program was performed by comparing the obtained numerical results with the available results of experimental measurements from the literature [27] and comparing them with the numerical results obtained using commercial programs TUF [12] and RELAP5/MOD3 [27].

The numerical simulation of the water hammer with air cavitation caused by the rapid valve closure and the liquid column separation in case of low-velocity flow is carried out for experimental conditions of water hammer test in [6]. The saturation pressure of water is not reached, and evaporation does not occur. Measured and calculated values of pressure change at different distances from the tube entrance are compared.

The upstream type of vaporous cavitation is simulated with the developed model and presented numerical method. Calculated values of pressure head and fluid velocity, at different distances from the isolating valve, are compared with measured values in experimental installation [52].

3.1 Condensation induced water hammer in a test facility

Numerical simulation of the steam–water interface propagation and the CIWH caused by direct steam condensation was performed for conditions of a simple experimental apparatus consisting of a tank, a horizontal and a vertical pipe, of the same diameter (Figure 4). Zaltsgendler et al. [18] experimentally investigated CIWH in the test facility. The horizontal pipe is at one end connected to the tank. The tank and the horizontal pipe are filled with subcooled water, at initial pressure 0.551 MPa and temperature 22°C. The vertical pipe is filled with saturated steam at the pressure of 0.382 MPa. At one end, the pipe is connected by a fast-acting ball valve to the horizontal pipe, while at the other end is closed.

Figure 4.

Schematic view of the experimental apparatus for CIWH testing [18].

By opening the fast-acting ball valve at the moment t = 0 s, the subcooled water and saturated steam are brought into direct contact, and conditions are met for CIWH event. The intensive condensation of steam onto the head of the water column occurs, and the water column starts moving towards the space which was filled with steam before condensation occurred. Due to condensation, the pressure in the steam drops sharply, and the water column accelerates towards the closed end of the vertical pipe. At the moment when all the steam in the vertical tube is condensed, the head of the water column splashes the closed end of the pipe at 1.34 s and causes a pressure increase of 8.53 MPa (Figure 5a). Afterwards, the pressure wave propagates towards the tank, where it is reflected from the water mass, and its amplitude attenuates due to the friction on the pipe walls. This wave with attenuated amplitude moves towards the closed end of the vertical pipe and hits it, but this time with less intensity. This process is repeated periodically. The numerically obtained results are in acceptable agreement with the measured values presented in Figure 5b. In Figure 5a, it is shown that the calculated pressure pulses diminish at about 3.8 s, while the measured pressure pulses diminish after 4 s, Figure 5b. The cause of the enlarged difference between measured and calculated pressure pulses, after the first peak occurrence, may be in the method of determination of the evaporation and condensation rate and the sonic velocity in the two-phase mixture after the pressure wave reflection from the closed pipe end. The first pressure pulse is the most dangerous from the point of view of equipment safety because it has the largest amplitude. By applying the developed model, a good prediction of the time of occurrence and the amplitude of the first and most dangerous pressure pulse was achieved. The presented calculation results were obtained with spatial step of 0.1 m.

Figure 5.

Pressure change near the closed end of the vertical pipe in the experimental apparatus for CIWH testing. Results obtained by applying the developed model (up) and the results of experimental measurements [18] (down).

3.2 Water cannon test

The CIWH in the vertical pipe for steam discharge into a pool with subcooled water, known as water cannon, was experimentally investigated in [27]. The experimental apparatus, shown in Figure 6, consists of a vertical metal pipe 0.7112 m long with an inner diameter of 0.0381 m and a large tank filled with subcooled water in which the pipe is immersed with its bottom to a depth of several centimeters. Pressure gauges are mounted on the top of the pipe. Saturated steam at a pressure of 0.1023 MPa is introduced into the vertical pipe at a constant speed through a small pipe, located at its top, and the lower end of the vertical pipe is immersed in a tank with water at a temperature of 49°C, at the same pressure of 0.1023 MPa. The problem is observed from the moment when the vertical pipe is completely filled with steam and the valve at its top is closed. The initial velocity in all parts of the system is equal to 0 m/s.

Figure 6.

Schematic view of the experimental apparatus for simulating a water cannon.

This situation can occur at the steam turbine exit, during the discharge of steam into a larger volume of cold water. Direct contact of steam and subcooled water leads to intensive steam condensation. The steam pressure drops, and a liquid column in the form of a plug enters the pipe. As the water level in the pipe rises, due to the difference in pressure, the condensation process is carried further towards its top, where the water column eventually hits and causes a large pressure pulse on the valve. The maximal pressure pulse occurs at the first impact.

Pressure change over time near the closed end of the vertical pipe obtained using the developed program (Figure 7a) is compared with pressure change obtained by the RELAP code [27] (Figure 7b). The intensity of the first pressure peak, which represents the risk of causing mechanical damage to the pipeline, matches the results from the literature with which they were compared. Other peaks predicted by the RELAP code (Figure 7b) are caused by the acoustic propagation of pressure waves in the vertical pipe between the tank and the closed end of the pipe, while the developed program (Figure 7a) predicts only four pulses caused by the movement of the water column in the pipe.

Figure 7.

Pressure change near the closed end of the vertical pipe in the water cannon test by Yeung et al. [27]. The result obtained by applying the developed model (up) and result of RELAP5/MOD3 simulation (down).

3.3 Water hammer with non-equilibrium gas release

The numerical simulation of the water hammer caused by the rapid valve closure and the liquid column separation was carried out for experimental conditions of water hammer test [6]. The experimental installation, shown in Figure 8, consists of two reservoirs under different pressures that are connected by a straight copper pipeline with a constant inclination to the horizontal of 3.2° and a total length of 37.23 m. The inner diameter of the tube is 0.0221 m. Demineralized water is used as the working fluid. Each performed experiment consisted of two phases. First, a stationary flow in the tube is established, the initial conditions are determined, and then, a transient is induced by rapidly closing the valve. The water hammer test is performed for the case where the fast-closing valve is upstream and water flows upwards. The initial velocity in the installation is equal to 0.3 m/s (i.e., low-velocity flow). The pressure in the reservoir 2 is 0.32 MPa, and the assumed initial water temperature is 20°C. The fast-closing valve with spring, located in front of the reservoir 1, was closed in 0.009 s in the experiment. Since the effective flow reduction in the experiment is 0.004 s, the reduction of the flow velocity in the calculation starts at 0.005 s, and at 0.009 s, the value of the velocity drops to zero. The initial value of the volume fraction of the gas phase (i.e., air micro-bubbles) is 10⁻⁷, which is the value assumed in [6], and the measured sonic velocity before closing the valve is 1319 m/s. The pressure is measured at four locations: in front of the fast-closing valve (p_v1 in Figure 8), at three quarters of the pipe length from the Reservoir 2 (p_q1 in Figure 8), at the half-length of the tube (p_mp in Figure 8) and at one quarter of the pipe length from the reservoir 2 (p_q2 in Figure 8).

Figure 8.

Schematic view of the experimental installation for water hammer with air cavitation [6].

The presented case of low-velocity flow leads to water hammer with air cavitation, where the saturation pressure of water has not been reached, and the evaporation does not occur. Measured and calculated values of pressure change at different distances from the tube entrance are shown in Figure 9. The calculated pressure values are in very good agreement with the measured values in both amplitudes and oscillation periods during the entire duration of the recorded transient process of 1.5 s. Both the measured and calculated data show a sudden pressure jump due to the rapid closing of the valve at the beginning of the transient (Figure 9a). The amplitude of the pressure increase is about 0.44 MPa, which corresponds to the value determined according to the Joukowsky equation. The compression wave generated in front of the closed valve propagates towards the reservoir 2 and is reflected as rarefaction wave that travels back towards the closed valve. At 0.066 s, the rarefaction wave is reflected at the closed valve as a wave of the same sign, and the pressure drops to a low value of 0.006 MPa, which is lower than the saturation pressure of air in water, but still higher than the saturation pressure of water and steam, which is 0.0023 MPa at 20°C. The liquid column separation occurs when the gap in the water flow appears due to the presence of air, resulting from gaseous cavitation, when pressure drops below the saturation pressure of air in water. Since the conditions for water evaporation have not been met, there is no vaporous cavitation.

Figure 9.

Comparison of measured and calculated pressure: (a) in front of the closing valve, (b) at three quarters of the pipe length from the inlet, (c) at half of the pipe length, (d) at one quarter of the pipe length from the inlet.

Here, presented results are obtained with a numerical discretization of the tube with 300 nodes. In the presented calculation, the initial mass fraction of dissolved air in water is estimated at 8 × 10⁻⁶, which corresponds to the saturation of air in water at a pressure of 0.035 MPa and a temperature of 20°C. This value was chosen so as to obtain a good agreement between the calculated peak pressures and the measured data. The appropriate prediction of the friction effect is important for a more precise determination of the pressure change. It is not enough to take into account only steady friction [48].

3.4 Water hammer with vaporous cavitation

Sanada et al. [52] have experimentally investigated water hammer with vaporous cavitation in the experimental installation shown in Figure 10, which consists of the horizontal acrylic pipeline of 200 m length with the inner diameter of 0.0152 m connected with two reservoirs. The pipeline is filled with water. The initial water parameters are given in [30]. The focus here will be on the upstream type of vaporous cavitation.

Figure 10.

Schematic view of the experimental installation for column separation investigation—upstream type [52].

The fast closure of upstream valve causes the pressure drop behind the valve. The pressure drops to and below the saturation pressure, the pressure at which the intensive evaporation occurs. After certain time, reflected pressure wave leads to collapse of vapor bubbles, formed in the evaporation process, and it causes pressure surge at the valve. Vapor bubbles form and disappear in cycles until the minimal steady pressure, higher than the saturation pressure, sets at the valve. This is column separation which occurs downstream, behind the obstacle which stops fluid flow.

The upstream vaporous cavitation is simulated with the model and numerical solution method presented in Section 2. Calculated values of absolute pressure head in meters and fluid velocity, at different distances from the isolating valve, are compared with measured values in experimental installation [52] in Figures 11 and 12. Figure 13 shows the comparison between measured pressure values in experimental installation [52] and numerical predictions from the various models for upstream cavitation obtained by [30]. For this severe case of vaporous cavitation, all predictions of Shu show similar characteristics. The best agreement is shown in case of frequency-dependent friction model, Figure 13d. Comparison of measured and calculated values of the pressure change in Figure 11 and of the change in velocity in Figure 12 shows satisfactory agreement for the amplitudes of the pressure wave and the periods of their oscillation. Compared to the results by Shu in [30] using a DVCM, a homogeneous equilibrium model and an improved transient friction model (Figure 13), the developed model provides a better prediction of the attenuation of the pressure wave amplitudes in period after 6 s.

Figure 11.

Comparison of measured and calculated absolute pressure head during column separation at different distances from the isolating valve: (a) x = 0 m, (b) x = 40 m and (c) x = 120 m. (d) Calculated absolute pressure head during column separation at the distance x = 200 m from the isolating valve.

Figure 12.

(a) Calculated fluid velocity during column separation at the location of the isolating valve. Comparison of measured and calculated fluid velocity during column separation at different distances from the isolating valve: (b) x = 40 m, (c) x = 120 m and (d) x = 200 m.

Figure 13.

Upstream vaporous cavitation. The absolute pressure changes at the isolating valve: (a) experimental results [52], (b) column separation model, (c) two-phase homogeneous equilibrium vaporous cavitation model and (d) frequency-dependent friction model [30].

4. Conclusions

The one-fluid thermal-hydraulic model is developed. The same governing equations are used for simulation and analysis of CIWH and water hammer with gaseous and vaporous cavitation. The model is based on one-dimensional, transient and compressible fluid flow of homogeneous fluid. Closure laws are developed with certain specificities for both phenomena. The model is solved by the application of MOC.

In the case of CIWH simulation, the energy equation is solved with the application of the LIP of the third degree for the calculation of initial enthalpy values at the characteristic paths of the fluid particles propagation. The prediction of the pressure pulses in CIWH events is validated by comparing with the available experimental measurements. The importance of the liquid column head tracking and the transient friction calculation for CIWH prediction is demonstrated. Previous experimental observations showed that the scatter of the test data with the same test conditions exists. This scattering is obtained also by numerical simulations in this study, and it is taken into account by modeling of the product of the condensation heat transfer coefficient and the interfacial area concentration which are correlated with the acceleration of the liquid column head and vapor interface. From the point of view of plant safety, the most important outcome of the numerical simulation of CIWH is the prediction of the first pressure pulse caused by the liquid column splashing.

A new model for water hammer with gaseous and vaporous cavitation is presented. The difference in modeling using the presented model compared to standard models is reflected in the fact that in this model the two-phase mixture exists anywhere along the pipe length, and consequently, the speed of wave propagation is equal to the sonic velocity in the two-phase mixture. Also, the model has no restrictions on the minimum pressure value during the transient, which is determined by the dynamics of the propagation of the pressure waves and the intensity of rarefaction waves. The new modeling approach is validated by comparing the obtained numerical results with the measurement results of an experimental test of the water hammer caused by the rapid valve closure on the experimental installation [6]. The present model is applicable to water hammer with gaseous and vaporous cavitation. The closure laws consist of the non-equilibrium model of gas release and absorption and desorption in water in case of non-condensable gaseous cavitation, as well as evaporation and condensation model in cases of vaporous cavitation. This modeling approach is sensitive to the spatial integration step, as a consequence of the large nonlinearity of the change of two-phase flow parameter, like the speed of pressure wave propagation and the two-phase mixture density, with changes of gas volume fraction. One of the advantages of this method is a simpler algorithm, than in standard methods, that can be easily implemented in computer programs. The gas release rate has uppermost value during the propagation of the first rarefaction wave, which causes the pressure to drop below the value of the saturation pressure of the gas in the water. In the sequel of the transient, when the gas is already released and the gas bubbles or pockets are formed, the rates of gas degassing or absorption are of a smaller order of magnitude. The difference in the gas generation rate could be the consequence of disturbed conditions in water, with dissolved gas which undergoes sudden pressure drop caused by rarefaction wave propagation during the transient. The developed model of gas release during water hammer with gaseous cavitation is justified by the comparison of the obtained numerical results with the experimental data for one research case.

Acknowledgments

The authors are grateful to Dr. Anton Bergant for providing the experimentally measured data.

This research was supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia (grant 451-03-68/2022-14/200105 and grant 451-03-68/2022-14/200213 of 4.2.2022.).

Nomenclature

A	area, m^2
a	interfacial area concentration, 1/m
a	sonic velocity, m/s
a	thermal diffusivity, m^2/s
C	constant in Eq. (10)
C	mass fraction
c	speed of pressure wave propagation, m/s
D	diffusion coefficient, m^2/s
d	inside diameter, m
E	Young’s modulus of elasticity, N/m^2
f	friction coefficient
g	gravity, m/s^2
H	pressure head, m
Hc	Henry’s constant, Pa
h	specific enthalpy, J/kg
h	heat transfer coefficient, W/m^2 K
j	mass flux, kg/m^2 s
k	coefficient of gas interface mass fraction in Eq. (38), k1 and k2, 1/Pa
k	coefficient of mass transfer in Eq. (44)
L	characteristic length, m
n	concentration, 1/m^3
p	pressure, Pa
q	surface heat flux, W/m^2
Rg	gas constant, J/kg K
r	latent heat of evaporation/condensation, J/kg
r	radius, m
T	temperature, K
t	time, s
t	renewal time in Eqs. (40) and (41), s
u	velocity, m/s
v	specific volume, m^3/kg
x	spatial coordinate, m
x	mass fraction
x	quality
$\tilde{x}$	molar fraction in Eq. (32)

Greek letters

$\alpha$	void fraction
$\mathrm{\Gamma }$	mass transfer rate, kg/m^3 s
$\delta$	thickness, m
$\epsilon$	infinitesimal distance, m
$\theta$	angle of inclination, °
$\lambda$	thermal conductivity, W/m K
$\mu$	dynamic viscosity, Pa s
$\nu$	kinematic viscosity, m^2/s
$\rho$	density, kg/m^3
$\sigma$	surface tension, N/m

Non-dimensional numbers

$\mathrm{Nu}=\frac{\mathit{hL}}{\lambda }$	Nusselt number
$Pr=\frac{\nu }{a}$	Prandtl number
$Re=\frac{\mathit{ud}}{\nu }$	Reynolds number
$\text{We}=\frac{{\rho }_1{\left(u_2-u_1\right)}^2{d}_b}{\sigma }$	Weber number

Subscripts or superscripts

${}^{\text{'}}$	saturated liquid state
${}^{\text{'}\text{'}}$	saturated vapor state
0	initial condition
1	liquid
2	gas/vapor
b	bubble
c	condensation
D	droplet
ED	Entrained Droplet
eq	equilibrium condition
H	hydraulic
h	isenthalpic process
i	interface
LCH	Liquid Column Head
p	isobaric process
s	isentropic process
sat	saturation
t	thermodynamic
u	unsteady

Acronyms and abbreviations

CFD	Computational Fluid Dynamics
CIWH	Condensation Induced Water Hammer
DGCM	Discrete Gas Cavitation Model
DVCM	Discrete Vapor Cavitation Model
ED	Entrained Droplets
FIIWH	Flashing Instability Induced Water Hammer
FSI	Fluid Structure Interaction
GIVCM	Generalized Interface Vapor Cavitation Model
LCH	Liquid Column Head
LIP	Lagrange’s Interpolation Polynomial
MOC	Method of Characteristics
NCS	Natural Circulation System
NPP	Nuclear Power Plant
ROSA	Rig-Of-Safety Assessment (ROSA)
SCIF	Steam Condensation Induced Flashing