Experimental Investigation and Computational Validation of Thermal Stratification in Piping Systems of PWR Reactors

1. Introduction

One phase thermally stratified flow occurs in horizontal piping where two different layers of the same liquid flow separately without appreciable mixing due to the low velocities and difference in density (and temperature). This condition results in a varying temperature distribution in the pipe wall and in an excessive differential expansion between the upper and lower parts of the pipe walls. This phenomenon can induce thermal fatigue in the piping system threatening its integrity. In some safety related piping systems of pressurized water reactors (PWR) plants, temperature differences of about 200 ^oC can be found in a narrow band around the hot and cold water interface. To assess potential piping damage due to thermal stratification, it is necessary to determine the transient temperature distributions in the pipe wall (Häfner, 2004) (Schuler and Herter, 2004).

Aiming to improve the knowledge on thermally stratified flow and increase life management and safety programs in PWR nuclear reactors, experimental and numerical programs have been set up at Nuclear Technology Development Center, a researcher institute of the Brazilian Nuclear Energy Commission (CDTN/CNEN) (Rezende, 2012), (Rezende et al. 2012). The Thermal Stratification Experimental Facility (ITET) was built to allow the study of the phenomenon as broadly as possible. The first test section was designed to simulate the steam generator injection nozzle and has the objective of studying the flow configurations and understanding the evolution of the thermal stratification process. The driving parameter considered to characterize flow under stratified regime due to difference in specific masses is the Froude number. Different Froude numbers, from 0.019 to 0.436, were obtained in different testes by setting injection cold water flow rates and hot water initial temperatures.

The use of Computational Fluid Dynamics (CFD) in nuclear reactor safety analyses is growing due to considerable advancements made in software and hardware technology. However, it is still necessary to establish quality and trust in the predictive capabilities of CFD methodologies. A validation work requires comparisons of CFD results against experimental measurements with high resolution in space and time. Recently, some research laboratories have been implementing experimental programs aiming to assist this demand.

The organization of the XVII ENFIR – Seventeenth Meeting on Nuclear Reactor Physics and Thermal Hydraulics propose the Special Theme on Thermal Hydraulics for CFD Codes as a contribution for the validation of CFD methodologies (Rezende el al. 2011a). The experimental results of thermal stratification developed at the Thermal Hydraulics Laboratory of CDTN/CNEN were used for comparisons with CFD results. This Chapter shows the results of the validation done by the Brazilian researcher (Rezende, 2012). The purposes of this special theme are: CDF simulations of a transient with a coolant thermally stratified single phase flow in the steam generation injection nozzle simulating experimental facility; performing comparisons of different CFD simulations and comparisons of CFD simulations with experimental results. Two sets of experimental data are proposed for the numerical simulation.

Numerical simulation was performed with the commercial finite volume Computational Fluid Dynamic code CFX. A vertical symmetry plane along the pipe was adopted to reduce the geometry in one half, reducing mesh size and minimizing processing time. The RANS two equations RNG k- turbulence model with scalable wall function and the full buoyancy model were used in the simulation. In order to properly evaluate the numerical model a Verification and Validation (V&V) process was performed according to an ASME standard. Numerical uncertainties due to mesh and time step were evaluated. The performed validation process showed the importance of proper quantitative evaluation of numerical results. In past studies a qualitative evaluation of the results would be considered sufficient and the present model would be (as it has been) considered very good for the prediction and study of thermal stratification. However, with the present V&V study it was possible to identify objectively the strengths and weaknesses of the model.

Results show the influence of Froude number on the hot and cold water interface position, temperature gradients and thermal striping occurrence. Results are presented in terms of wall temperature, internal temperature, vertical probe temperature, temperature contours and velocity fields.

2. The thermal stratification experimental facility at CDTN

The Thermal Stratification Experimental Facility (ITET) wasbuilt in the Thermal-hydraulic Laboratory at Nuclear Technology Development Center (CDTN) (Fig. 1) to allow the study of the phenomenon as broadly as possible. The first test section was designed to simulate the steam generator injection nozzle. Figure2shows a drawing of this test section that consists of a stainless steel tube (AISI 304 L), 141.3 mm in outside diameter and9.5 mm thick.

It was made of two pieces of this tube connected each other by a 90^o curve, a vertical and a horizontal piece respectively 500 mm and 2000 mm length. A flanged extension of the tube was placed inside a pressure vessel, which simulates the steam generator. Thermocouples were placed in four Measuring Stations along the length of the test section tube. Measuring Stations I, II and III, located in the horizontal length of the tube were instrumented with thermocouples, measuring both fluid and wall temperature at several positions of each Measuring Station. Measuring Station A, positioned in the vertical length of the tube, was instrumented with three thermocouples just to determine the moment when the injected cold water reaches its position.

Figure 1.

Thermal-hydraulic Laboratory at Nuclear Technology Development Center (CDTN)

Experiment	Flow rate [kg/s]	Pressure [Pa]	Initial system temperature [°C]	Cold water injection temperature [°C]
2146	0.76	2.11 x 10^6	32	220
2212	1.12	2.14 x 10^6	28	221

Table 1.

Input data for the proposed experiments

Before the beginning of each test the whole system is filled with cold water. Then it is pressurized and heated by steam supplied by a boiler. A temperature equalization pump ensures that the entire system is heated in a homogeneous way. After the heating process, the equalization pump is turned off and both the steam supply and equalization lines are isolated by closing valves V3, V5 and V6. The test itself begins then by injecting cold water from the lower end of the vertical tube after opening valve V4. The cold water flow rate was previously adjusted at a value planned in the test matrix. This flow rate and the system pressure are maintained stable through a set of safe (V1) and relieve (V2) valves at the upper side of the pressure vessel, which controls upstream pressure. The water flows from the injection nozzle simulator pipe to the steam generator simulator vessel through 11 holes at the upper side of the extension tube placed inside the vessel. These holes are 12 mm in diameter and they are displaced 42 mm from each other. The center of the first hole is 20 mm from the end of the tube.

Figure 2.

Position of the Measuring Stations A, I, II and III in the steam generator injection nozzle simulating test section

2.1. The instrumentation

Measurement Stations I, II and III, positioned along the longitudinal length of the tube simulating the steam generator injection nozzle, as shown in Figure 2, were used for temperature measurements. Figures 3, 4 and 5 show the thermocouples distribution in Measuring Stations I, II and III, respectively. To measure fluid temperature on Measurement Station I a set of 12 thermocouples was angularly distributed along the tube’s internal wall (3 mm from the wall), shown in Figure 3 by circle symbols. These internal thermocouples were named clockwise starting from the highest vertical position as T1I01, T1I02, …, T1I11 and T1I12. To measure the tube’s wall temperature another set of 12 thermocouples was brazed on the outside wall at the same angular position as the internal thermocouples, displayed by triangle symbols in Figure 3. These external thermocouples were named clockwise starting from the highest vertical position as T1E01, T1E02, …, T1E11 and T1E12. Finally, a removable probe was placed along the cross section’s vertical diameter, containing a set of 9 fluid thermocouples placed at the same vertical position of each of the internal thermocouples, shown by square symbols in Figure 3. These probe thermocouples were named from the highest to the lowest vertical position as T1S01, T1S02,.....T1S08 and T1S09.

Figure 3.

Positions of the thermocouples at Measurement Station I

Figure 4 shows the thermocouple distribution on Measurement Station II. A set of 19 thermocouples was angularly distributed along the tube’s internal wall (3 mm from the wall) to measure fluid temperature, shown in Fig. 4 by circle symbols. Close to the angular position of 90° a set of 5 internal thermocouples was positioned in close proximity, displaced 2 mm from each other, to capture fluctuations of the cold-hot water interface. In the opposite side 2 internal thermocouples were positioned in the same manner to capture asymmetrical behaviors of the interface. These internal thermocouples were named clockwise starting from the highest vertical position as T2I01, T2I02, …, T2I18 and T2I19. Another set of 14 thermocouples was brazed on the outside wall at the same angular position as the internal thermocouples (only 1 external thermocouple was positioned at the angular positions of 90° and 270°), shown in Fig. 4 by triangle symbols. These external thermocouples were named clockwise starting from the highest vertical position as T2E01, T2E02, …, T2E13 and T2E14. Finally, a removable probe was placed along the cross section’s vertical diameter, containing a set of 10 fluid thermocouples placed at the same vertical position of each of the internal thermocouples, as shown in Fig. 4 by square symbols. These probe thermocouples were named from the highest to the lowest vertical position as T2S01, T2S02,...., T2S09 and T2S10.

Figure 4.

Positions of the thermocouples at Measurement Station II

Figure 5 shows the thermocouple distribution on Measurement Station III. Close to the angular position of 90° a set of 4 internal thermocouples, named from the highest to the lowest vertical position as T3I01, T3I02, T3I03, and T3I04, was positioned 3 mm from the internal wall and displaced 2 mm from each other to measure fluid temperature. A fifth internal thermocouple, named T3I05, was placed at the angular position of 180°, shown in Fig. 5 by circle symbols. Two thermocouples, named T3E01 and T3E02, were brazed on the outside wall of the tube at the angular positions of 90° and 180° respectively, shown by triangle symbol in Fig. 5. Finally, a removable probe was placed along the cross section’s vertical diameter containing a set of 6 fluid thermocouples, shown as square symbols in Fig. 5. These probe thermocouples were named T3S01, T3S02, T3S03, T3S04, T3S05 and T3S06 from top to bottom. They were placed respectively at the same vertical positions of thermocouples T2S03, T2S04, T2S05, T2S07, T2S08 and T2S10.

A set of three thermocouples was positioned at Measuring Station A to detect the instant when the injected cold water reaches its position. The thermocouples were placed inside the tube 3 mm from the wall, at the center of the cross section by a probe and at the external wall

Figure 5.

Positions of the thermocouples at Measurement Station III

Figure 6 shows a photograph of the test section pipe after the brazing of the thermocouples. Figure 6 shows in detail the outside of Measuring Station I. The external thermocouples were brazed directly to the pipe and the internal thermocouples were brazed through special stainless steel injection needles. Some aluminum brackets for the thermocouples are seen in the back, which were only used during the assembly of the experimental facility. Figure 7 shows the Measuring Station I internal thermocouples. Figure 7 shows a photograph of the ITET, including the horizontal tube of the injection nozzle, the pressure vessel simulating the steam generator and the cold water tank.

Other measurements performed were:

• injection flow rate of cold water, using a set of orifice plate and differential pressure transmitter;

• water temperature in the cold water tank, using an isolated type K thermocouple of 1 mm in diameter;

• water temperature in the cold water injection pipe, both close to the orifice plate and also close to the point of injection to the nozzle simulation tube, using two isolated type K thermocouples of 1 mm in diameter;

• temperature inside the steam generator simulation vessel, using an isolated type K thermocouple of 1 mm in diameter;

• pressure inside the steam generator simulation vessel, using a gauge pressure transducer;

• pressure in water injection line, using a gauge pressure transducer;

• water level in the cold water tank using a differential pressure transmitter.

Figure 6.

The test section’s horizontal pipe after the thermocouples brazing, and detail of the measuring station

Figure 7.

The internal thermocouples in the Measuring Station I, and the Thermal Stratification Experimental Facility (ITET) during assembly

2.3. The measuring uncertainty

The measuring uncertainties for the main parameters, obtained according to ISO (1993), were:

• 2.4°C for the temperature measurements;

• 2.4 % of the measured value for the flow rate measurements; and,

• 1.5 % for the gauge pressure measurements.

3. Simulation results

The experimental results of thermal stratification developed at the Thermal Hydraulics Laboratory of CDTN/CNEN were used for comparisons with CFD results. In recent theoretical evaluations, CFD (Computational Fluid Dynamic) analysis using three dimensional Reynolds Averaged Navier Stokes (RANS) has been used, which is due to several reasons, from the ease of use of commercial codes and development of low costs computational systems of reasonable processing capacity, to the speed at which results are obtained.

However, before CFD can be considered as a reliable tool for the analysis of thermal stratification there is a need to establish the credibility of the numerical results. Procedures must be defined to evaluate the error and uncertainty due to aspects such as mesh refinement, time step, turbulence model, wall treatment and appropriate definition of boundary conditions. These procedures are referred to as Verification and Validation (V&V) processes (Roache, 2010). In 2009 a standard was published by the American Society of Mechanical Engineers (ASME) establishing detailed procedures for V&V of CFD simulations (ASME, 2009).

According to the Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer – V&V 20 (ASME, 2009), the objective of validation is to estimate the modeling error within an uncertainty range. This is accomplished by comparing the result of a simulation (S) and an experiment (D) at a particular validation point. The discrepancy between these two values, called comparison error (E), can be defined by Equation 1 as the combination of the errors of the simulation (${\mathrm{\delta }}_{\mathrm{s}}=\mathrm{S}-\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{ }\mathrm{V}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}$) and experiment (${\mathrm{\delta }}_{\mathrm{e}\mathrm{x}\mathrm{p}}=\mathrm{D}-\mathrm{T}\mathrm{r}\mathrm{u}\mathrm{e}\mathrm{ }\mathrm{V}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}$) to an unknown True Value.

$$\mathrm{E}=\mathrm{S}-\mathrm{D}={\mathrm{\delta }}_{\mathrm{s}}-{\mathrm{\delta }}_{\mathrm{e}\mathrm{x}\mathrm{p}}$$E1

The simulation error can be decomposed in input error (${\mathrm{\delta }}_{\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}$) that is due to geometrical and physical parameters, numerical error (${\mathrm{\delta }}_{\mathrm{n}\mathrm{u}\mathrm{m}}$) that is due to the numerical solution of the equations and modeling error (${\mathrm{\delta }}_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}$) that is due to assumptions and approximations. Splitting the simulation error in its three components and expanding Equation 1 to isolate the modeling error gives Equation 2.

$${\mathrm{\delta }}_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}=\mathrm{E}-({\mathrm{\delta }}_{\mathrm{n}\mathrm{u}\mathrm{m}}+{\mathrm{\delta }}_{\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}-{\mathrm{\delta }}_{\mathrm{e}\mathrm{x}\mathrm{p}})$$E2

The standard applies then to this analysis the same concepts of error and uncertainty used in experimental data analysis, defining a validation standard uncertainty, u_val as an estimate of the standard deviation of the parent population of the combination of the errors in brackets in Equation 2, in such a way that the modeling error falls within the range$\left[\mathrm{E}+{\mathrm{u}}_{\mathrm{v}\mathrm{a}\mathrm{l}},\mathrm{E}-{\mathrm{u}}_{\mathrm{v}\mathrm{a}\mathrm{l}}\right]$, or using a more common notation:

$${\mathrm{\delta }}_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}=\mathrm{E}\pm {\mathrm{u}}_{\mathrm{v}\mathrm{a}\mathrm{l}}$$E3

Supposing that the errors are independent, u_val can be defined as Equation 4.

$${\mathrm{u}}_{\mathrm{v}\mathrm{a}\mathrm{l}}=\sqrt{{\mathrm{u}}_{\mathrm{n}\mathrm{u}\mathrm{m}}^2+{\mathrm{u}}_{\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}}^2+{\mathrm{u}}_{\mathrm{e}\mathrm{x}\mathrm{p}}^2}$$E4

The estimation of these uncertainties is at the core of the process of validation. The experimental uncertainty can be estimated by well established techniques (ISO, 2003). Input uncertainty is usually determined by any propagation techniques or analytically (ASME, 2009). The numerical uncertainty, on the other hand, poses greater difficulties to access.

The estimation of the numerical uncertainty is called verification and is usually split into two categories: code and solution verification. Code verification evaluates the mathematical correctness of the code and is accomplished by simulating a problem that has an exact solution and verifying if that solution is obtained. This activity requires extensive programming access to the core of the code which is not available in commercial codes, due to this it is common practice to take commercial codes as verified by the supplier.

Solution verification is the process of estimating the numerical uncertainty for a particular solution of a problem of interest. The two main sources of errors here are the discretization and iteration processes. The discretization error is the difference between the result of a simulation using a finite grid in time and space and that obtained with an infinitely refined one. The methods developed to evaluate it are based on a systematic grid refinement study where the solution is expected to asymptotically approximate the exact value as the grid is refined, at a rate proportional to the discretization order of the solution. The iteration error is present in codes that use iterative solvers, where the result must converge to the exact value as the iterations develop. It is usually estimated using the residual root mean square (RMS) between subsequent iterations of a variable over all the volumes of the domain.

The numerical simulation of the Experiment 1 shown in Table 2 was performed by Resende et al (2011b and 2011c) using CFX 13.0 (ANSYS, 2010) code in a simplified geometry. The geometry in Fig. 2 was simulated with the omission of the flanges and most of the lower inlet geometry, as shown in Fig. 8. These simplifications have no significant influence on the results. A second flow condition showed in Table 2 (Experiment 2) was also simulated to further evaluate the numerical methodology.

Experiment	Flow rate [kg/s]	Pgauge [bar]	Thot [^oC]	Tcold[^oC]
1	0.76	21.1	219.2	31.7
2	1.12	21.4	217.7	28.7
^(1)	0.03	0.5	2.4	2.4

Table 2.

Setup parameters for the experiments and simulations

(1) Global uncertainty

Figure 8.

Computational model domains and boundary conditions.

The computational model was generated with two domains: one solid, corresponding to the pipes, and one fluid for the water in its interior. A vertical symmetry plane along the pipe was adopted to reduce the mesh size in one half, minimizing processing time. The walls in the vessel region were considered adiabatic as the external tube walls. Mass flow inlet and outlet conditions were defined at the bottom end of the pipe and high end of the vessel, respectively. Figure 2 shows the computational model’s details.

The initial conditions shown in Table 2 were used in the simulations. Water properties like density, viscosity and thermal expansivity were adjusted by regression as function of temperature with data extracted from Table IAPWS-IF97, in the simulation range (25 ^oC to 221 ^oC). The RANS - Reynolds Averaging Navier-Stokes equations, the two equations of the RNG k- turbulence model, with scalable wall functions, the full buoyancy model and the total energy heat transfer model with the viscous work term were solved. The simulations were performed using parallel processing with up to six workstations with two 4 core processor and 24 GB of RAM. All simulations were performed using the high resolution numerical scheme (formally second order) for the discretization of the conservation and RNG k- turbulence model equations terms and second order backward Euler scheme for the transient terms. A root mean square (RMS) residual target value of 10^-6 was defined as the convergence criteria for the simulations in double precision. By using this RMS target the interactive error is minimized and can be neglected in the uncertainty evaluation as its contribution are usually many orders lower that of other sources like discretization (Roache, 2010).

A mesh and time step study described in the following section were performed according to ASME V&V 20 standard to assess the numerical uncertainty (ASME, 2009).

A solution verification study was performed according to ASME CFD Verification and Validation standard to evaluate mesh and time step uncertainties (ASME, 2009).

Three gradually refined non-structured tetrahedral meshes with prismatic near wall elements (inflated) were generated for the model presented in Fig. 2 to evaluate mesh related uncertainty. Progressive grid refinements were applied to edge sizing of the piping elements. The ratio between the height of the last prismatic layer and the first tetrahedral was kept equal to 0.5 for all meshes. Three layers of prismatic structured volumes were built close to the surfaces in the solid and fluid domains. The growth factor between prismatic layers was maintained constant with a value of 1.2. A localized mesh edge sizing of 5 mm was applied at the inlet nozzle of the vertical pipe and vessel outlet nozzle for all meshes. At the outlet holes of the horizontal pipe an edge sizing of 2 mm was also used for all meshes. Element sizing in the vessel was set to expand freely with a growth factor of 1.2.

The characteristics of the generate meshes are shown in Table 3. The table includes the resulting grid refinement ratio (r_i) and representative grid edge size (h_i) defined by Equations 5 and 6, respectively. Figure 9 shows some details of the generated meshes.

$${\mathrm{r}}_{\mathrm{i}}={\mathrm{h}}_{\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{ }\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{e}\mathrm{ }\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{ }\mathrm{i}+1}/{\mathrm{h}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{ }\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{ }\mathrm{i}}$$E5

$${\mathrm{h}}_{\mathrm{i}}={\left(\mathrm{M}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{ }\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}/\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{ }\mathrm{o}\mathrm{f}\mathrm{ }\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{ }\mathrm{o}\mathrm{f}\mathrm{ }\mathrm{i}\mathrm{ }\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}\right)}^{1/3}$$E6

Mesh i	hi [mm]	No. of elements / nodes	ri	Element Edge Length [mm]
1	2.84	2,809,114 / 13,533,642	1.83	2.5
2	5.22	583,012 / 2,191,174	1.67	5.0
3	8.70	198,152 / 472,909	-	10.0

Table 3.

Meshes characteristics

Figure 9.

Mesh details.

To evaluate time step related uncertainty, three gradually refined time steps shown in Tab. 4 were used for the simulation of the model with mesh 2 presented in Fig. 3. Table 4 includes the resulting time step refinement ratio (r_j) defined by Equation 7.

$$r_j=t_{lastcoarsetimesyepj+1}/t_{presenttimestep}$$E7

Time j	tj [s]	rj
1	0.075	1.51
2	0.113	1.50
3	0.169	-

Table 4.

Time steps characteristics

Solution verification was performed using the three generated meshes and three simulated time steps based on the Grid Convergence Index method (GCI) of the ASME V&V 20 standard (ASME, 2009). The theoretical basis of the method is to assume that the results are asymptotically converging towards the exact solution of the equation system as the discretization is refined with an apparent order of convergence (p) that is in theory proportional to the order of the discretization scheme. The objective of the method is to determine p utilizing three systematically refined discretizations and determine relative to the finest discretization result a 95% confidence interval (±U_{num 95%} = ±GCI) where the exact solution is. In other word, the objective is to determine the expanded uncertainty interval due to the discretization.

Considering the representative grid edge sizes ${\mathrm{h}}_{\mathrm{i}-1}<{\mathrm{h}}_{\mathrm{i}}<{\mathrm{h}}_{\mathrm{i}+1}$ and grid refinement ratios${\mathrm{r}}_{\mathrm{i}}={\mathrm{h}}_{\mathrm{i}+1}/{\mathrm{h}}_{\mathrm{i}}$, the apparent order of convergence p can be determined by Equations 8, 9 and 10. In an analogous manner similar equations can be obtained for time discretization, however these will be omitted for brevity.

$${\mathrm{p}}_{\mathrm{i}}=\frac{1}{\mathrm{l}\mathrm{n}\left({\mathrm{r}}_{\mathrm{i}}\right)}\left|\mathrm{l}\mathrm{n}\left|{\mathrm{\epsilon }}_{\mathrm{i}+1}/{\mathrm{\epsilon }}_{\mathrm{i}}\right|+\mathrm{q}\left({\mathrm{p}}_{\mathrm{i}}\right)\right|$$E8

$$\mathrm{q}\left({\mathrm{p}}_{\mathrm{i}}\right)=\mathrm{l}\mathrm{n}\left(\frac{{\mathrm{r}}_{\mathrm{i}}^{{\mathrm{p}}_{\mathrm{i}}}-\mathrm{s}}{{\mathrm{r}}_{\mathrm{i}+1}^{{\mathrm{p}}_{\mathrm{i}}}-\mathrm{s}}\right)$$E9

$$\mathrm{s}=1\bullet \mathrm{s}\mathrm{g}\mathrm{n}\left({\mathrm{\epsilon }}_{\mathrm{i}+1}/{\mathrm{\epsilon }}_{\mathrm{i}}\right)$$E10

where${\mathrm{\epsilon }}_{i+1}={\mathrm{\varphi }}_{\mathrm{i}+2}-{\mathrm{\varphi }}_{i+1}$, ${\mathrm{\epsilon }}_i={\mathrm{\varphi }}_{i+1}-{\mathrm{\varphi }}_i$, ${\mathrm{\varphi }}_{\mathrm{k}}$denotes the variable solution on the k^th grid and sgn is the signal function (sgn(x) = -1 for x < 0; 0 for x = 0 and 1 for x > 0).

It is recommended by the standard ASME (2009) that the obtained value of p be limited to the maximum theoretical value, which for the used high resolution and Euler discretization scheme is 2. Also the value of p can be limited to a minimum of 1 to avoid exaggerations of the predicted uncertainty, however when limited it is recommended that the obtained value is presented for comparison.

With the value of p the expanded uncertainty GCI can be calculated using Equation 11 using an empirical Factor of Safety (Fs), equal to 1.25, that is recommended for studies with more than three meshes (ASME, 2009).

$${\mathrm{G}\mathrm{C}\mathrm{I}}_{\mathrm{i}}=\frac{\mathrm{F}\mathrm{s}\bullet {\mathrm{\epsilon }}_i}{{\mathrm{r}}_{\mathrm{i}}^{{\mathrm{p}}_{\mathrm{i}}}-1}$$E11

When the presented procedure is applied to obtain the GCI for local variables, such as a temperature profile, an average value of p should be used as to represent a global order of accuracy.

Mesh and time step uncertainties are considered independent in this study and the total numerical expanded uncertainty is calculated through Equation 12.

$${\mathrm{U}}_{\mathrm{n}\mathrm{u}\mathrm{m}}=\sqrt{{\mathrm{G}\mathrm{C}\mathrm{I}}_{\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}}^2+{\mathrm{G}\mathrm{C}\mathrm{I}}_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{ }\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}}^2}$$E12

In this study the temperature profiles along time in several positions of the test section were evaluated. Figure 10 displays the analyzed positions that are equivalent to the thermocouple positions of the experiments.

Figure 10.

Thermocouples positions

Table 5 shows the some of the obtained results of the performed verification process. Average values for p and GCI are presented as the maximum GCI of the entire profile. These maximums were all located in regions of steep temperature gradients, which explain the very high observed values.

Position in the pipe	Mesh			Time step:
pm*	GCIm* [^oC]	Maximum GCIm [^oC]	pt*	GCIt* [^oC]	Maximum GCIt[^oC]
Internal
T1I01	1.58	14.012	41.608	1.00	0.056	0.173
T1I05	1.88	1.174	35.020	1.32	0.415	17.257
T1I10	1.52	0.496	75.830	1.27	0.578	45.339
T2I01	1.32	7.394	22.829	1.31	0.030	0.218
T2I06	1.87	1.377	51.166	1.22	0.594	21.577
T2I15	1.48	1.489	99.554	1.20	0.748	70.006
T3I02	1.80	1.099	64.103	1.21	0.544	13.988
T3I05	1.47	1.220	122.122	1.23	0.584	51.247
Probe
T1S05	1.64	2.034	60.014	1.23	0.546	14.964
T2S04	1.65	1.766	58.885	1.16	0.902	15.776
T3S02	1.44	2.381	45.358	1.32	0.796	20.076
External
T1E02	1.61	2.198	6.347	1.05	0.010	0.095
T1E10	1.56	0.199	0.601	1.21	0.087	0.537
T2E01	1.34	0.164	1.115	1.17	0.003	0.015
T2E11	1.16	1.599	3.302	1.47	0.0185	0.621

Table 5.

Verification process results for several thermocouple positions.

* Time averaged values.

It can be observed in Table 5 that uncertainties due to the mesh are in average greater than those due to the time step. One reason for these values could be attributed to the course mesh used in the study that could lead to overestimation of the total uncertainty of the refined mesh. In average only thermocouples T1I01 and T2I01 displayed uncertainties above the experimental one of 2.4 ^oC, both located in the upper region of the vertical tube which indicates that this region is the most affected by the mesh refinement.

An example of the obtained results from the verification process is shown in Fig. 11 that displays the temperature profiles at thermocouple position T2S04 obtained by the simulated meshes and time steps. The figure also shows the uncertainties obtained along the simulated time. It is observed that the mesh contribution to uncertainty is much greater than that of the time step. The grater values of uncertainty were obtained the abrupt temperature drop region and at the subsequent temperature oscillation period.

Figure 11.

Numerical uncertainty evaluation due to the mesh and time step

The obtained uncertainty prediction through the solution verification process proposed by ASME V&V 20 standard ASME (2009) showed very variable and sometimes incoherent results for the uncertainty prediction. The method takes in account three discretizations for the estimate of GCI and requires that the results between these discretizations be “well behaved” to produce coherent uncertainty estimates. Convergence must be “well behaved” due to the core assumption made by the method that the solution is converging asymptotically as the mesh is refined. This is a very strong assumption as it has been concluded in recent studies that it is safest to assume that the numerical data are not within the asymptotic regime (Eça et al., 2009). It is in fact questionable that even the finest meshes in use today can produce solutions that are in this regime (Lockard, 2010). However, the obtained method gives a good insight as to how results are varying as the mesh is refined and the estimated uncertainty may not be accurate but is a quantification on how “well behaved” is the solution.

Following the solution verification, a validation process was performed comparing the numerical results with experimental data. To determine the validation expanded uncertainty, U_val (Eq. 4), only the estimated numerical and experimental uncertainties were considered neglecting the input contribution. Although the input uncertainty is in fact non-neglectable, its evaluation is beyond the purpose of this study as it is extremely complex requiring hundreds of simulations taking in account fabrication tolerances and uncertainties in all measurable variables.

Figures 12 shows a comparison between numerical and experimental results as well as the validation error (E = S – D) and validation uncertainty for the upper thermocouples. Very high validation uncertainty after the beginning of the temperature drop can be observed. This high uncertainty is attributed to the mesh that influences greatly the results in this region. Validation become poor after 150 s of simulation, however before that time numerical and experimental results agree well.

Figure 12.

Validation results for the upper thermocouples

Figure 13.

Validation results for the probe thermocouples

Figure 14.

Validation results for the lower thermocouples

Figures 13 and 14 show a comparison between numerical and experimental results with the validation uncertainty and the validation error (E = S – D) for the probe and lower thermocouples, respectively. For both regions results show a high validation error and uncertainty for the beginning of the temperature drop and subsequent oscillations. It is observed that at the cold water front reaches the center of the pipe (probe thermocouples) before the experiment and that the temperature drop in the lower region of the pipe is quicker in the simulation. Although considerable validation error is observed the qualitative agreement between experiment and simulation can be considered good as most of the behavior observed was reproduced.

Figure 15 shows the evolution of the temperature differences between the average temperatures on the highest and lowest positions of the horizontal tube calculated through the Equation 13 for the internal and external thermocouples.

$$DT={\left[\left(T_1^u+T_2^u\right)-\left(T_1^l+T_2^l\right)\right]}^{}/2$$E13

where the superscripts u and l are relative to the upper and the lower positions and subscripts 1 and 2 to the first and the second measuring station of the horizontal tube.

Figure 15.

Validation results for the temperature difference between upper and lower thermocouples positioned internally and externally.

Figure 15 shows that for the region of highest temperature difference, and therefore, most critical for the piping integrity, the validation error is relatively low and well predicted. It is also observed that the external temperature difference agreement between experimental and numerical results is very good during the evaluated time.

The performed validation process showed the importance of proper quantitative evaluation of numerical results. In past studies a qualitative evaluation of the results would be considered sufficient and the present model would be (as it has been) considered very good for the prediction and study of thermal stratification. However, with the present V&V study it was possible to identify objectively the strengths and weaknesses of the model.

Figure 16.

Comparison of numerical results obtained for two flow conditions

Figure 16 shows a comparison between numerical results obtained by the presented model for two flow conditions. By the results it is possible to conclude that thermal stratification occurs for both flow rates with similar intensity and temperature differences levels.

Figure 17 shows the cold water front evolution obtained numerically for the flow rate of 1.12kg/s. It can be observed that a cold water “head” is formed as the cold water front advances in the horizontal pipe. It can also be observed a change in the direction of the cold water front after reaching the end of the tube.

Figure 17.

Temperature contours along time for flow rate 1.12 kg/s.

Figure 18 shows details of the flow behavior and flow velocity evolution in the simulation of flow rate 1.12 kg/s. It can be observed that as the cold water front enters the horizontal pipe it accelerates due to stratification and that the front induces a recirculation flow of the hot water at the top of the pipe, as mass must be conserved.

Figure 18.

Flow velocity and behavior along time for flow rate 1.12 kg/s.

Figure 18 highlights the previously observed behavior, i.e., as the cold front reaches the end of the pipe it starts filing the pipe in the inversed direction eliminating almost all of the recirculating hot water. However, some hot water remains imprisoned at the top of the pipe as the injected cold water takes control of all water exits. This phenomenon is observed experimentally and causes the thermal stratification at the top of the pipe to persist for many minutes depending on the flow rate

4. Conclusion

The numerical simulation of one phase thermally stratified flow experiments in a pipe, similar to the steam generator injection nozzle at the secondary loop of a Pressurized Water Reactor (PWR), was proposed. The simulations were done using CFD codes (Rezende et al. 2011b).

A V&V evaluation of the numerical CFD methodology based on ASME (2009) standard was performed. Solution verification was also performed using three progressively refined meshes and time steps. Temperature profiles in several positions inside and outside the piping system were evaluated. In average the uncertainties due to the mesh were greater than those due to the time step. One reason for these values could be attributed to the course mesh used in the study that could lead to overestimation of the total uncertainty of the refined mesh. In average only thermocouples located in the upper region of the vertical tube displayed uncertainties above the experimental one (2.4 ^oC), which indicates that this region is the most affected by the mesh refinement.

The performed validation process showed the importance of proper quantitative evaluation of numerical results. In past studies a qualitative evaluation of the results would be considered sufficient and the present model would be (as it has been) considered satisfactory for thermal stratification prediction and study. However, with the present V&V study it was possible to identify objectively the strengths and weaknesses of the model.

Although considerable validation error was observed the qualitative agreement between experiment and simulation can be considered good as most of the behavior observed was reproduced. The performed validation process showed the importance of proper quantitative evaluation of numerical results.

Acknowledgment

This research project is supported by the following Brazilian institutions: Nuclear Technology Development Centre (CDTN), Brazilian Nuclear Energy Commission (CNEN), Research Support Foundation of the State of Minas Gerais (FAPEMIG), and Brazilian Council for Scientific and Technological Development (CNPq).