Abstract
The subcooled decompression under temperature gradient experiment performed by Takeda and Toda in 1979 has been reproduced using the in-house code WAHA version 3. The sudden blowdown of a pressurized water pipe under temperature gradient generates a travelling pressure wave that changes from decompression to compression, and vice versa, every time it reaches the two-phase region near the orifice break. The pressure wave amplitude and frequency are obtained at different locations of the pipe's length. The value of the wave period during the first 20 ms of the experiment seems to be correct but the pressure amplitude is overpredicted. The main three parameters that contribute to the pressure wave behavior are: the break orifice (critical flow model), the ambient pressure at the outlet, and the number of volumes used for the calculation. Recent studies using RELAP5 code have reproduced the early pressure wave (transient) of the same experiment reducing the discharge coefficient and the bubble diameter. In the present paper, the long-term pipe pressure, that is, 2 seconds after rupture, is used to estimate the break orifice that originates the pressure wave. The numerical stability of the WAHA code is clearly proven with the results using different Courant numbers.
1. Introduction
Propagation of pressure waves through the piping systems of nuclear power plants is of key interest in the field of reactor safety analyses. Pressure surges are important also in other fields of process engineering [1]. During the loss-of-coolant-accident in a water-cooled reactor, the subcooled decompression results in the propagation of a rarefaction wave through the reactor system and can even result in further structural damage and failure to maintain core geometry and core cooling [2].
One of the key experiments historically well known as a benchmark case for thermal hydraulics system code RELAP5 [3] was the one of Edwards and O’Brien [4]. They ruptured a horizontal pipe filled with water at pressures and temperatures characteristic for water-cooled reactors. Later, Takeda and Toda [5] experimentally observed the pressure behavior in a vertical pipe ruptured at the top end with a temperature gradient (higher temperature at the top end) and analyzed it with simple equilibrium and nonequilibrium wave propagation models. Takeda and Toda showed that for a pipe with a temperature gradient, flashing in the hotter section results from the passing of a rarefaction wave with a large enough amplitude, whereas the colder section, with a lower vapor pressure, remains subcooled liquid. A large discontinuity in the sound speed arises between the two-phase and single-phase regions serving as a reflective surface. The data from Takeda and Toda has been compared with a RELAP5 simulation by Lafferty et al. [6].
The RELAP5 computer code, widely used in nuclear reactor transient analysis, is based on a two-fluid model consisting of the mass, momentum, and energy conservation equations for each phase solved using a semi-implicit finite-difference technique [3, 7]. In the present paper, computer code WAHA [8], developed for simulations of acoustic phenomena in piping systems, is used to simulate the experiment of Takeda and Toda. WAHA code uses a similar mathematical model as RELAP5, but it is based on a characteristic upwind finite difference scheme. This allows more accurate tracking of pressure waves with less numerical diffusion and modeling of the critical two-phase flow with the basic two-fluid model without any special models.
The paper is divided in two main groups. Description of the experimental facility and introduction to the WAHA computer program are given in Section 2, and the results of the simulations are presented in Section 3 and 4. Section 3 starts with an explanation of the WAHA piping model and it is followed by a descriptive presentation of the simulations’ results testing different numerical parameters. The influence of the parameters in the results gives an insight in the WAHA code numerical implementation. The influence of the Courant number, order of accuracy of the calculation and the number of volumes, used to divide the pipe’s length, are tested initially. In Section 4 the sensitivity of different physical models is studied, first testing the pipe elasticity model and wall thickness influence, Section 4.1. Using the best of these parameters, the simulations are improved further studying the consequences of a smaller effective orifice diameter and higher values of the tank pressure in Sections 4.2 and 4.3 respectively. At the end, the conclusions are drawn.
2. Description of the Experiment and the Waha Code
The subcooled blowdown under temperature gradient experimental facility [5] comprised a test section of stainless steel pipe of 3.2 m length and 53.5 mm internal diameter. An orifice of 15 mm diameter was located at the top of the pipe resembling the pipe break and it was followed by an exhaust tube minimizing back pressure effects after the orifice. The orifice was covered with 0.03~0.75 mm thick mylar paper. By rupturing the mylar paper, the sudden blowdown experiments started. The water filling the pipe section was heated by electrical heaters until the desired linear temperature distribution which was measured with thirteen temperature measuring terminals was attained. Five different experiments were performed with different temperature distributions. The temperature at the bottom of the pipe was fixed at room temperature (RT) for all five experiments. The temperatures at the top, close to the orifice, were RT, 95°C, 120°C, 140°C, and 160°C. Six different pressure transducers (PT) were mounted along the pipe length section to capture the fluid pressure changes. The positions of the transducers can be seen in Table 1. After heating up the water to the specific test temperature distribution, the pressure was set at 0.855 MPa at the position PT1 before the beginning of the experiments. Figure 1 shows a schematic view of the pipe.
The experiment with fluid temperature at the top of the pipe equal to 160°C has been modeled using the in-house code WAHA version 3 [8]. The WAHA computer program is a one-dimensional six-equation two-fluid solver developed for simulations of two-phase flow water hammer phenomena. The continuity, momentum and energy equations for liquid and gas phases are used in a nonconservative form. The continuity and energy equations are derived to account for the pipe elasticity and variable pipe cross-section effects. The pipe cross section () can change due to variable geometry, which depends on the position, and the material elasticity, (1):
The pressure pulse () produces a variation of the pipe diameter () following the linear relation [9] where stands for the pipe wall thickness and E the modulus of elasticity of the material.
The WAHA code uses its own steam tables and contains correlations for mass, momentum, and heat transfer between phases and for wall friction. The correlations are flow regime dependent. The system of WAHA equations can be written in the form: where is the vector of independent variables (pressure, vapor volume fraction, velocities, and specific internal energy for the liquid and gas phases, resp.), and are the matrices of the system, and is a vector which contains the nondifferential terms in the equations. The numerical scheme of the WAHA computer program uses the splitting operator with implementation consisting of two steps explained below. First, the source term vector is split in two parts. The “non-relaxation” sources which contain the wall friction, volumetric forces, sources due to the variable pipe cross-section and sources due to the pipe motion (not implemented in version 3). The second part of the source term vector includes the relaxation sources that tend to establish thermal and mechanical equilibrium between liquid and vapor phases. The relaxation source terms are often stiff, meaning that their characteristic times are much shorter than other phenomena in two-phase flow. In the first step of the splitting operator method, the convection term with the nonrelaxation sources is treated with shock-capturing numerical scheme with low numerical diffusion: In the second step, the system of ordinary differential equations—relaxation source terms—is solved with a variable time step: where now contains the temperature of both phases instead of the specific internal energy. Matrices , , and are fully developed in [8].
3. Waha Model—Sensivity to Numerical Parameters
The experiment with water temperature at the top of the pipe equal to 160°C has been reproduced using the WAHA computer program.
3.1. Model and Initial Conditions
The WAHA code uses ASCII input data files that contain all geometrical, physical, computational and temporal characteristics of a modeled piping system. Figure 1 shows the two modeled components to simulate the experimental facility. The first component is the pipe which consists of N volumes dividing its length. The lower part of the pipe is a closed end ensuring fluid velocity equal to 0. The last volume (Figure 1) at the top of the pipe has a cross section that is equivalent to the break/orifice diameter. This volume is connected to the second component, a tank. The tank component is a constant pressure boundary condition, and it is set to ambient pressure. The fluid initial conditions of temperature, pressure, and vapor pressure can be seen in Figure 2. The temperature at the top of the pipe is 160°C. The total simulated time is two seconds (2 s) after the blowdown with special attention to the first 20 ms of the transient. The simulations’ outputs are the pressure values of those volumes corresponding to the pressure transducers, locations (Table 1).
(a)
(b)
(c)
3.2. Results of the Simulations
In order to assess the capability of the WAHA computer code to simulate the two-phase flow decompression blowdown, the influence that different parameters have on the results has been tested. In the next subsections, the numerical parameters analyzed are the Courant number, order of accuracy of the calculation and the number of volumes dividing the pipe’s length. In the next section, sensitivity to various physical models is shown: role of the pipe elasticity model and wall thickness, and the break flow modeling. Other parameters (results not shown in the present paper) which have no influence on the results were also tested. Those parameters are the wall friction model, material roughness or friction factor. The results presented in this work are compared with the experimental results from [5].
3.2.1. Influence of Courant Number
The convective part of the system of equations in WAHA code, (4) is derived following a finite difference scheme [8]. The stability domain for its integration is limited by the CFL (Courant-Friedrichs-Levy) condition: where is the Courant number, the spatial increment, the time increment, and the eigenvalues of the diagonal matrix , which come from the diagonalization of the Jacobian matrix in (4). The eigenvalues provide a good approximation for the sum of the speed of sound in the fluid media () and the fluid velocity (). After the blowdown, a small region with low void fraction appears close to the orifice where the speed of sound is one order of magnitude lower, that is, from approximately 1400 m/s in the one-phase liquid region to 50–500 m/s in the two-phase region. The fluid velocity is almost zero in the one-phase region and close to 30 m/s in the two-phase region. Therefore, the most limiting value for resides in the one-phase region where the time increment follows:
Figure 3 shows the 1 ms time window when the rarefraction wave reaches the two-phase region and then it is reflected as a compression wave. The moment it reaches the two-phase region the pressure drop increases the void fraction near the break [10]. Figure 4 shows the results of the simulations performed with WAHA code using 100 volumes and pipe wall thickness of 3 mm, for different Courant numbers. They are compared with the experimental results and the numerical simulations of the same experiment using the RELAP5 code, published recently by Lafferty et al. [6]. The RELAP5 simulations were performed with the default thermal nonequilibrium model and 99 nodes dividing the pipe section.
(a)
(b)
A general observation that can be extracted from Figure 4 is that the behavior of the simulated pressure wave is correct when compared with the experimental results from [5]. However, the results in Figure 4 show important differences between the WAHA and RELAP5 codes. While the WAHA simulations are not affected by the Courant number, the RELAP5 results show different pressure drops and subsequent average pressure levels for the two Courant numbers shown here. The first pressure drop value (≈0.6 MPa) simulated with WAHA code follows the experimental results closer. Furthermore, WAHA predicts higher pressure amplitudes than RELAP5.
3.2.2. Order of Accuracy of the Calculation
The WAHA code has implemented a combination of first-order upwind discretization and second-order discretization of (3). The problems arising from the pure second-order discretization, like oscillations near the vicinity of non-smooth solutions can be solved using a combination of the first- and the second-order accurate discretization, that is, with the use of limiters. Limiters account for the smoothness of the solutions. If the solutions are smooth, a larger part of the second-order discretization is used, otherwise larger part of the first-order discretization is used [8, 11]. Figure 5 shows the simulation results for pure first-order upwind discretization and the combination of first- and second-order discretization using the limiters described in [8, 11]. In the transient solved here, the pure second-order implementation gives unstable results and they are not shown. The simulations have been performed with 100 volumes, Courant number equal to 0.8 and pipe wall thickness equal to 3 mm.
The period of the pressure oscillations does not change significantly between the first-order and the combined-order with different limiters. First-order upwind discretization, which is supposed to give the least accurate solution, overpredicts excessively the first pressure drop at the beginning of the simulation. On the other hand, the pressure amplitudes are better estimated. The three limiter options give reasonable results for the first pressure drop, the Superbee limiter the one which follows the experimental results more accurately at the beginning of the simulation. Superbee limiter is supposed to give better results for steepest waves while the most smeared waves are better reproduced with MINMOD. Van-Leer’s limiter solution should lie between the solutions obtained with the MINMOD and the Superbee limiters. This can be clearly observed in Figure 5. Superbee limiter gives the highest pressure amplitudes differing from the experimental results. As well, Superbee and Van-Leer solutions give non-physical oscillations after steep pressure changes. Pressure changes with MINMOD limiter are smooth and the pressure amplitudes, although overestimated, are not as high as with Superbee or Van-leer. MINMOD will be used in further sections.
3.2.3. Influence of Number of Volumes
The influence of the number of volumes used to divide the pipe’s length can be seen in Figure 6. Three possibilities with N equal to 100, 1000, and 2000 have been tested, despite the general recommendation for 1D models, according to which, the use of volume length lower than the volume diameter often does not make too much sense. Almost no differences can be seen between 1000 and 2000 volumes except the much higher computing time needed to perform the simulations. The required time using 1000 volumes is two hours while the time needed for 2000 volumes goes up to ten hours. The period of the oscillations is practically the same for 100 or 1000 volumes. A small change in the wave behavior can be seen at the lowest pressure value when the rarefaction wave passes by the position PT3, at time 0.005 s. The two main differences that can be observed in Figure 6 are the lower pressure level reached with 1000 volumes after the first pressure drop and the higher pressure amplitudes using 100 volumes. While the first pressure drop is better reproduced using 100 volumes, the pressure amplitude is excessively overpredicted. The number of volumes used to divide the pipe’s length gives the height of an individual volume. Consequently the last volume, which has a cross section equivalent to the break orifice, has a height of 3.2 cm for 100 volumes and of 3.2 mm for 1000 volumes. In the 100 volumes case, the longer pipe wall at the orifice region produces extra friction effects on the fluid, which finds more difficulties to leave the pipe after the blowdown, decreasing the first pressure drop value. Further on, the pressure in the pipe is kept at higher levels and slightly higher pressure amplitudes are generated after each wave reflection.
No data has been found showing the dimensions of the wall thickness at the top of the pipe where the orifice was machined. It seems more likely that the wall thickness should be similar to the pipe’s wall thickness, that is, 3 mm (see Section 4.1), and so the number of volumes to divide the pipe should be 1000.
After the results shown so far, there seems to be no parameter available through the input file, which could further adjust the simulation results to the experimental ones. Further adjustments could be made by changing the physical models inside the code, however, this was avoided in the present work. Slightly lower pressure level is reached after the first blowdown and further pressure amplitudes are clearly overpredicted. Everything seems to point towards extra friction effects, that influenced the experimental results but that aren’t considered so far in our simulations, or 3D flow effects at the region close to the break which cannot be modeled with one-dimensional codes. Parameters such as material roughness, the wall friction model, or multiplying friction factor [8] have been tested as well but they did not produce any difference in the results. The fluid velocity in the one-phase region is practically 0. Consequently the influence of any friction-related parameter is negligible in this region. The fluid velocity in the two-phase region obtained during the simulations is in the order of 30 m/s. The use of reasonable values for the material roughness and friction factor parameter has no influence on the results in this region either. Furthemore, the pressure levels after the two second simulations were compared with the experimental values. The break diameter and downstream pressure (tank pressure) were modified to adjust the long term pressure.
4. Sensitivity of Physical Models
Various physical models were tested to check the sensitivity of the results. Pipe elasticity and break modeling turned out to be relevant and are described below. As mentioned previously, the wall friction model, material roughness or friction factor were tested too, but with no major influence on the results.
4.1. Elasticity Model and Wall Thickness
The pipe elasticity can be taken into account in the WAHA code equations [1, 2]. Pipe elasticity reduces the effective sonic velocity in the pipe [9]: where as in (2), is the effective sound velocity and is the sonic velocity in the pure fluid.
In Figure 7 one can observe the influence that the material elasticity has on the results. Three pipe wall thicknesses (2, 3 and 4 mm) and the results without elasticity (or stiff pipe) are compared, because the wall thickness is not available in the experiment’s description. The simulations were performed using 100 volumes dividing the pipe’s length, Courant number equal to 0.8 and pipe’s material with elastic modulus of 2.0 1011 Pa. The first pressure drop at the beginning of the experiment is slightly overpredicted but the average pressure follows the experimental data quite well. The pressure wave amplitude is clearly overpredicted at PT3 location but it is less at PT4 and PT5. The period of the simulated wave observed at PT4 and PT5 positions agrees with the experimental results as well. The influence of the pipe elasticity is clear observing the results from the figure. When a completely stiff pipe is considered, meaning that material elastic modulus (or the pipe wall thickness) is infinity, the period at PT4 and PT5 closely matches the experimental data. On the other hand, the pressure drop generated by the rarefraction wave is too fast. This can be seen in PT3 location at 0.005 s. If the pipe elasticity is considered, the pressure wave is slower. Thinner the pipe wall, higher is the wave period. The results for pipe thickness equal to 3 mm match the time of the experimental lower pressure value at 0.005 s. Furthermore, for a 53.5 mm internal diameter stainless steel pipe, a thickness of 3 mm seems to be the best assumption [12].
(a)
(b)
4.2. The Break Diameter
The results of the fluid pressure after two seconds transient obtained with the WAHA code are half the value of the experimental results (Figure 8). One of the reasons to explain this behavior could be that the effective orifice diameter was smaller than 15 mm due to failure of the Mylar paper to break completely generating unaccounted pressure drops. The pressure value obtained in [5] after two seconds blowdown was approximately 0.3 MPa, and it can be seen as a horizontal line in Figure 8(a). We can also observe in the figure the pressure values at the six pressure transducers locations using different simulated break diameters. The simulations were performed using 1000 volumes, MINMOD limiter, courant number 0.8 and pipe wall thickness 3 mm. It can be seen that with a 9 mm break diameter the final experimental pressure can be reached while with the 15 mm diameter, the pressure is exactly half the experimental value. Diameters from 6 to 9 mm are modeled slightly different using the WAHA code. Recommendation for simulations with smooth area change model in WAHA code manual (Section in [8]) is to avoid cross-section ratios, between neighboring volumes, exceeding a factor of roughly 2. Thus, due to the big difference between the dimensions of the pipe and break diameters, the top three volumes have decreasing cross sections, modeling a kind of “bottle neck” at the top of the pipe. For smaller than approximately 10 mm, the ratio with the pipe diameter is . Figure 8(b) shows the first 20 ms of the transient modeled using 15 mm and 9 mm break diameter and they are compared with the experimental results. Figure 9 shows the same comparison but for the 2 s time simulations.
(a)
(b)
The short-term effects using 9 mm break diameter are a slightly lower initial pressure drop, but higher pressure amplitudes, Figure 8(b). Like with the number of volumes, the results are not conclusive. On the other hand, looking at the 2 s results (Figure 9), at PT4 and PT6 locations the experimental pressure level is followed quite well by the simulations using mm. It has to be pointed out that the simulated pressure oscillations at those locations are kept for a longer time and higher in magnitude compared to the experimental data. The final pressure level at PT1 is the same as in the experiment. WAHA code simulated a smooth pressure decrease during the 2 s transient while the experiment showed slightly lower pressure levels between 0.5 and 0.75 s and then the pressure went up to 0.3 MPa. The initial pressure decrease tendency at PT1 during first 0.5 s seems to be followed more accurately modeling the break using 15 mm diameter.
4.3. Downstream—Tank Pressure Influence
Another reason for the 2s experimental pressure values to be higher than the simulated ones could be a back pressure downstream of the break orifice higher than 0.1 MPa. The experimental facility had mounted a ball valve and an exhaust pipe downstream of the break orifice. The exhaust pipe diameter was not much bigger than the pipe diameter and not far from the break orifice the exhaust pipe had an elbow. Most likely the pressure right after the break was not exactly constant to 0.1 MPa, as it was initially modeled in WAHA by the use of the tank boundary condition. No pressure readings were taken at locations after the break orifice in [5]. The pipe downstream of the break has not been modeled with WAHA but the tank pressure has been set at higher values to match the experimental results, as can be seen in Figure 10(b). The break diameter used in the simulations is 15 mm. It can be seen that with a tank pressure of 0.3 MPa, the pressure at PT1 and PT6 locations follows the experimental results closely and at PT4 the pressure is a bit underpredicted. Furthermore, at PT1 the experimental results follow initially the simulations using a tank pressure of 0.1 MPa and then move towards those with values using 0.3 MPa. Despite the discrepancies, the results look quite promising.
(a)
(b)
Looking at Figure 10(a), one can observe the pressure at PT3 for 20 ms transient with tank pressure of 1 and 3 atm. They are compared with the experiment and the RELAP5 results [6]. In this case RELAP5 simulations were performed with a reduced discharge coefficient of 0.5 and different bubble diameter by setting a Laplace length factor of 0.1. RELAP5 approaches slightly better the first pressure drop value from the experiment but the rarefraction wave passes PT3 location too soon, at t = 0.005 s. The WAHA results instead do not seem to be affected excessively during the 20 ms transient due to an external pressure of 3 atm. We can see that the pressure amplitude is a bit lower with a tank pressure of 3 atm which is the desired behavior. The time for the rarefraction wave is better predicted compared to the RELAP5. Also the average pressure level seems to be better reproduced with the WAHA code.
Due to these last results we can conclude that a smaller effective break diameter together with downstream unexpected back pressure and 3D fluid flow effects before the pipe break are probably the reasons for the overestimated pressure amplitudes and the lower value of the first pressure drop obtained with the simulations.
5. Conclusions
In the present paper, the in-house computer program WAHA version 3 is used to reproduce the subcooled decompression under temperature gradient experiment performed by Takeda and Toda in 1979. The WAHA code is a one-dimensional six-equation two-fluid solver developed for simulations of the two-phase flow water hammer phenomena. The equations for the two-phase flow are implemented in non-conservative form and the pipe elasticity effects are taken into account.
Different WAHA parameters have been tested such as the pipe elasticity model and wall thickness, the Courant number, order of accuracy of the calculation and the number of volumes that divide the pipe. The WAHA code can’t reproduce perfectly the experimental results. The main problem resides in the overpredicted pressure amplitudes. The first pressure drop after the blowdown is also slightly overpredicted. The period of the wave though follows very well the experimental one. Extra friction effects that had influenced the experimental results but that are not considered so far in our simulations or 3D flow effects at the region close to the break which can’t be modeled with one-dimensional codes seem to be the reasons for the disagreement. Parameters such as material roughness, the wall friction model or multiplying friction factor have been tested as well but did not produce any significant difference in the results.
In order to improve the results, the influence of different effective break diameters is checked and it is found that the 2 s transient experimental pressure value can be obtained with a break diameter of 9 mm. Further on, the influence of the downstream pressure after the break is also tested. With a boundary condition of 0.3 MPa tank pressure after the break, the final experimental pressure is also obtained. The first pressure drop is the same and the pressure amplitudes, although still overpredicted, are slightly lower. The comparison of the WAHA results with recent RELAP5 simulations of the same experiment shows that WAHA code follows the experimental data more accurately although the pressure amplitudes are highly overpredicted. In order to further improve the WAHA model of the experiment, pressure readings downstream of the break orifice would be of a big advantage.
Acknowledgment
The authors gratefully acknowledge the financial contribution of the Slovenian Research Agency through the Research Project J2-4078 and Research Program P2-0026.