© F. Világi et al., published by EDP Sciences 2019

1 Introduction

The safety of nuclear power plants is the most important aspect during harnessing the energy from the atom. While a lot of people may think that the safety of nuclear energy is sufficient, the nuclear accidents during the last few decades indicates the contrary. Especially the Fukushima disaster on 11 March 2011, showed that there is a need for a form of passive safety for nuclear power plants. The new main objective is to minimize the dependence on electrical supply. The new reactor solutions will also rely on passive systems. The use of passive safety systems eliminates the costs associated with installation, maintenance and operation of active safety systems that require multiple pumps with independent and redundant power supplies [1]. Another incentive to use passive security systems is the potential for increased security through increased system reliability, because passive systems don't have any moving parts which could fail during operation. Keller [3] first introduced the idea of using a natural circulation loop to transmit the heat between hot and cold piping branch. Keller also provided mathematical description of such a device backed by experimental measurements. The stability of the process had yet to be presented, where the works of Welander [4] showed instabilities in the meaning of oscillating flow-rate and pressure in the device. The research of NCL was further developed by theoretical work of Zvirin [5]. All the measurements conducted were only at laboratory scale and the mathematical models had to be adapted to the industrial scale. In order to accomplish this task, proper scaling criteria had to be established by Zuber and Heisler [6,7]. Scaling laws, proposed by Zuber, are also known as the philosophy of power/volume ratio. This philosophy has some inherent flaws, which surpasses some specific phenomenons in the NCL, such as flow instability [8]. This scaling contained a lot of dimensionless factors, and the conversion to an NCL with another device layout (heater, cooler) was troublesome. The problems had to be solved and therefore new scaling laws were introduced by Vijayan [9], who proposed scaling laws, which included the instability nature of NCLs. This method is currently widely used as the state of the art.

2 Description of the experimental facility

The main reason for building the experimental facility is to evaluate the capability of a natural circulation loop in removing the decay heat produced from a nuclear reactor in case of a failure of the main cooling system. This safety mechanism is supposed to be implemented into nuclear reactors of IV generation. The most reasonable medium used is helium, due to its inert properties regarding chemical reactions, also as radiological processes [2]. Helium is supposed to be heated directly from the nuclear reactor and then to be cooled down in the heat exchanger (water). The flow of helium should be ensured only by the difference of densities in hot and cold piping branch. The driving power is thus obtained directly from the heat supplied [2]. The experimental loop is 10 m height and 3 m wide construction, the flow process is single-phase and the device orientations are vertical heater and horizontal cooler. The experimental facility allows us to observe the helium flow at various conditions. The facility stands in the area of Energomont s.r.o. in Trnava. It consists of two main elements, the GFR (reactor substitution) and the DHR (heat sink), along with two piping branches (hot and cold). The GFR is considered to be used for energy input in form of heat, it consists of electrical heating rods with maximum installed power of 500 kW, with adjustable power range of 220 kW. The projected heat output of the DHR is also 220 kW with possible water flow regulation. The pressure difference on the elements, helium temperatures, surface temperatures, water flow, inlet power and velocities are measured during the experiment at several points along the loop. Helium can achieve the following operating values based on the settings of GFR and DHR devices:

• Helium temperature at the inlet to the GFR: 400–520 °C,

• Helium temperature at the outlet from the GFR: 150–250 °C,

• Operating pressure: 3–7 MPa.

The experimental facility provides an opportunity to verify and research the thermodynamic and hydraulic properties of self-driven, highly compressed, helium flow. The facility enables modelling of changes in input and output power of elements, as well as changes in operating pressure. It is used for the research of steady and unsteady helium flow in a given range of pressures and temperatures [10].

Figure 1 represents the diagram of the experimental facility. Besides the geometry it also shows the positions of measuring elements for temperature, pressure and velocity. The facility encompasses an advanced system of automatized data collecting and evaluating. As the heat conduction through pipes to the surroundings is generally neglected, the most important values from the measurement are temperatures and pressures from the control points 1–4, which are marked with red colour. The measurement of velocities happen in points marked by blue colour. Te small difference between the velocities in red points and the velocities in blue points was neglected.

The experiment starts by pressurising the loop with helium at room temperature. The heater, as well as the cooler power output is then gradually adjusted, based on a predefined time course. By this procedure, the steady-state was achieved within around 4–5 hours depending on the experiment scenario. The steady state was kept for a relatively short period of time (30–60 minutes). After the steady state period, the shutdown of the system was initiated. The shutdown began with the cool-down of the whole system to 30 °C. Then, the helium was transferred to the low-pressure vessel and then to the pressure cylinders.

Fig. 1 Diagram of the experimental facility.

3 Measurement of velocity

There are four velocity sensors installed on the natural circulation loop (BS1-BS4), but for the most measurements only one of the sensors gave valid data (BS1). The method of velocity measurement is based on the application of the Bernoulli equation − the velocity is measured through the determination of the dynamic portion of pressure. The dynamic portion of pressure p_d is determined as the difference between the total pressure p_t and the static pressure p_s . The static pressure is measured by the sampling of static pressure at the wall of the pipe. A pitot tube, connected to the capacitive pressure difference sensor was used for the measurement of total pressure (Fig. 2).

On the basis of theoretical derivation a relation between the total pressure (pressure of the braked fluid) and other properties of the flow, for example the velocity and the static pressure, can be determined. Obviously, the measured total pressure value p_c differs from the local pressure of the braked fluid p_ci , which was theoretically determined to stop the isoentropic flow to zero. Assuming of an ideal fluid, the relation between the properties of flow can be expressed as:$$\frac{p_{ci}}{p_s}={\left(1+\frac{k-1}{2}{M}_{\infty }^2\right)}^{\frac{k}{k-1}}\text{.}$$(1)

For a real fluid flowing at low Mach numbers the equation can be adjusted to:$$u=C_p\sqrt{\frac{2\left(p_{ci-p_s}\right)}{\rho }}\text{,}$$(2)where C_p is a calibration coefficient, experimentally estimated in the calibration tunnel. According to the described process the instantaneous point velocity u was measured. For the comparison of measurement with mathematical model, we need an relation between the measured point velocity u and the cross-sectional velocity v. Assuming that the velocity u is measured at the pipe axis, this point velocity takes the maximum value of u _max. Then we can express the cross-sectional velocity as:$$v=\beta u_{\text{max}}$$(3) where β is the coefficient of fullness of velocity profile. The value of coefficient β depends on the regime of the flow. During a laminar flow, the velocity profile in a round pipe is parabolic. The ratio of (cross-sectional) velocity to the maximal point velocity at the pipe axis is expressed as the coefficient of fullness of velocity profile by equation:$${\beta }_L=\frac{v}{u_{\max }}=0.5\text{.}$$(4)

The flow in a round pipe is considered as fully developed turbulent flow at velocities corresponding to Reynolds number greater than 4 × 10³. The velocity profile can be approximated by the power law in the form (Tab. 1):$$\frac{v}{u_{\max }}\cong {\left(\frac{y}{R}\right)}^{\frac{1}{n}}\text{.}$$(5)

However, in the binding sub-layer (near the channel wall), the velocity increases linearly with increasing distance from the wall. The exponent n in the power law of the velocity profile increases slightly with increasing Reynolds number.

It is obvious that during the flow of helium in experimental natural circulation loop, there is a period of time when the flow is not laminar anymore, but also can't be considered as fully developed turbulent flow. For the calculation of velocity in this area, it is necessary to consider the coefficient value of β as a function of Re. This dependence is determined by the measurement in the calibration tunnel. By this, the resulting velocity is determined with real parameters of helium flow taken into account, by an iterative process.

Fig. 2 Helium velocity measurement.

4 Representation of data from measurement number X

The results of one measurement (nr. X) are presented below. The time course of temperature, pressure and velocity was processed into charts. Other properties of flow such as mass flow and energy balances were derived from the measured properties. Figure 3 shows the time course of selected physical properties of the helium (temperatures in hot and cold branch, absolute pressure).

The energy balance was established form the measured data. The chart below shows the time course of thermal power accumulated into the helium (GFR) and thermal power drained from the coolant in (DHR). The third line then represents the electrical power input to the GFR (Figs. 4–7).

Fig. 3 Time course of temperature and pressure of the coolant.

Fig. 4 Time course of velocity in the hot branch (the most accurate).

Fig. 5 Time course of the outside insulation temperature.

Fig. 6 Time course of mass flow.

Fig. 7 Time course of thermal power gained and given to/from the coolant inside GFR and DHR.

5 Hydraulic model of flow

The mathematical model described in this article has been named as “Hydraulic” model, because it doesn't calculated the heat transfer between the elements and it needs the temperature distribution as input [2,11]. For a closed natural circulation loop the flow depends on the difference of densities in the cold and hot branch:$$\left({\rho }_{CB}-{\rho }_{HB}\right)gh={\displaystyle \sum _{i=1}^n}{p}_{li}\text{.}$$(6)where the left side of the equation represents the potential energy available due to different densities in the cold ρ_CB and the hot ρ_HB branch. During steady flow, the equation must be valid as the flow can not continuously accelerate or slow down, so the potential energy gained must be spent on the pressure losses p_l. The sum of pressure change across the loop (i.e., 1-2-3-4-1) must be zero (steady flow):$${\displaystyle \sum _{i=1}^3}\left(p_i-p_{i+1}\right)+\left(p_4-p_1\right)=0\text{.}$$(7)

The method uses a relation for calculation of pressure change between two points during steady flow [11].$$p_{i+1}=\sqrt{\frac{p_i^2}{e^{\alpha }}-\lambda \frac{L}{d}\frac{zr{T}_m}{A^2}{Q}_m^2\left(\frac{e^{\alpha }-1}{\alpha e^{\alpha }}\right)}$$(8) $$\alpha =\frac{2g\left(h_{i+1}-h_i\right)}{zrT}\text{.}$$(9)

This method is supposed to observe and evaluate the effect of pressure loss on the mass flow during steady state flow. The equation contains two unknowns that are dependent on each other, namely the mass flow Q_m and the pressure at the next point p _i+1. Thus, the equation needs to be calculated by iterative method [11].

We calculate the pressure p_{i +1} with equation (8) from the pressure at the previous point. Then we can sequentially calculate the pressure at each next control point on the loop as illustrated in (10). To start the calculation we need an initial pressure. The total operating pressure − p _1,est, which is either chosen as the condition at which we want to calculate the steady state velocity, or in our case, this is the measured pressure at point 1. The calculated pressure p _1,calc is the estimated pressure which should be in point 1 according the equation (8) and can be obtained as follows:$$p_{1,est}{\displaystyle \stackrel{}{}}{p}_2{\displaystyle \stackrel{}{}}{p}_3{p}_4{\displaystyle \stackrel{}{}}{p}_{1,calc}\text{.}$$(10)

When calculating the pressure at points 2, 3, 4, 1, the relationship (8) is used. At places where the local losses are significant (GFR and DHR) the friction factor λ is expressed in terms of local losses coefficient ξ as described in section 5.2 equation (17). During steady flow, the pressure difference must be zero, so we can write$$f\left(Q_m\right)=p_{1,calc}-p_{1,est}=0\text{.}$$(11)

During the steady state, this condition must be kept, so we can iteratively calculate the mass flow as well as the flow velocity from measured temperatures and operating pressure. The mean temperature T_m that is in equation (8) is determined from the measurement as the average between two points. The calculation is performed by the iterative Newton method (12), which serves to find the root of the equation (11). The new mass flow rate after each iteration is determined as:$$Q_m^{j+1}=Q_m^j-\frac{f(Q_m)}{f\text{'}(Q_m)}\text{.}$$(12)

For the first iteration, the initial value of the mass flow is chosen in the order close to a predicted value. The derivation of the function is calculated by a numerical derivation:$$f\text{'}(Q_m)=\frac{f(Q_{m,\epsilon })-f(Q_m)}{\epsilon }\text{.}$$(13)

The member f(Q_m,ε ) expresses the difference between the selected pressure p _1,est and the pressure p _1,h. This pressure p _1,ε is calculated analogously to p _1,calc with increased mass flow by differential value: Q_m,h^j = Q_m ^j + ε where ε = 10⁻⁵.$$f(Q_{m,\epsilon })=p_{1,\epsilon }-p_{1,est}\text{.}$$(14)

The stop condition is set to Q_m ^j+1 –Q_m ^j < ε, where ε = 10⁻⁵. When this condition is reached, the iteration cycle ends and the velocity at any point on the loop is calculated from the continuity equation:$$v_i=\frac{Q_m}{{\rho }_i{A}_i}\text{.}$$(15)

5.1 Evaluation of pressure losses due to friction

The pressure losses due to friction were calculated by an existing method from literature [12]. The equation also known as the McKeon relation was used for calculation of friction factor:$$\frac{1}{\sqrt{\lambda }}=1.930\text{log}\left(\mathrm{Re}\sqrt{\lambda }\right)-0.537\text{.}$$(16)

The equation (16) is considered as one of the most accurate. The McKeon relation can be used only for smooth pipes and the literature [12] states the maximum error of 1.25% for low Reynolds values in range of 0 ⟨Re⟩ 310 × 10³, and error of 0.5% for range 310 × 10³ ⟨Re⟩ 18 × 10⁶. In our case, the Reynolds numbers are approximately in the range of 15 000 to 37 000 unless we count the start-up.

5.2 Evaluation of local pressure losses

The largest pressure losses occurs in the heater (GFR) and in the cooler (DHR). These local losses are introduced into the equation (8) by simple conversion of local loss friction factor to length loss friction factor.$$\lambda ={\xi }_{\text{DHR/GFR}}\frac{L}{d}\text{.}$$(17)

Equation (17) needs the local loss coefficient as input, which had to be calculated separately for both elements from measured values of velocities, temperatures and pressures. The calculation was by adjusting the well-known equation for pressure losses due to local resistance [2]:$$\sqrt{p_l}=\sqrt{\xi }\left(\sqrt{\frac{{\rho }_{ie}}{2}}{v}_{ie}\right)\text{.}$$(18)

As the flow in GFR and DHR is not isothermal, the density and the mean velocity change along the input to the output in both elements. Thus, the local loss coefficient ξ was related to the averaged values of the changing parameters v_ie and ρ_ie [2]. The Bell-Delaware method was used for calculation of pressure losses inside the elements for each measurement. This method was considered as sufficiently accurate. Subsequently, losses from sudden expansion and reduction at inputs and outputs were added. Function (18) is a straight line with$\sqrt{\xi }$ being its slope. By plotting the dependence into a graph and after data regression, the local loss coefficients ξ _DHR and ξ _GFR were determined.

Figure 8 shows nine plotted dependencies in which the pressure loss was calculated by the Bell-Delaware method and one dependency with directly measured pressure loss. The measurement is consistent with calculations. The final value of ξ _DHR was thereby determined at: 48.837 and the value for ξ _GFR at: 4.875.

Fig. 8 Determination of local loss coefficient for DHR.

6 Comparison between experiments and calculations

The most measurements started with helium being filled into the experimental facility at 2.5 MPa. After that, the GFR was started and as the heat input increased, the temperature in the hot branch rise resulting in decreased density, hence the flow was induced. By consequent helium cooling in the DHR exchanger, the temperature difference between two branches stayed constant at the steady state. In Figure 9, there are several pictures representing a comparison between measured values of velocities in the hot branch v_H and calculated values using the mathematical model. The temperature distribution, as well as operating pressure, is a necessary input to the calculation and was obtained from the measurement. The densities of helium were calculated by the Soave-Redlich-Kwong equation of state, which is considered as one of the most accurate [2,13]. Soave-Redlich-Kwong's (SRK) equation originated in 1972 when Soave [14], extended the Redlich-Kwong equation.

Figure 9 shows the temperatures progress along the length of the loop from measurements III, IV, VII and X during steady flow. The first temperature is the value from control point nr. 4 (bottom of the cold pipe), then the values from the temperature sensors in the interior of GFR are displayed, continuing with hot branch and DHR. The last displayed point is again the value from control point nr. 4.

Figure 10 shows a record of mean velocity v_H in the hot branch from four measurements (III, IV, VII and X). The frequency of measured velocity 1 Hz and the total duration of measurements varied from 22 000 to 28 000 seconds. The numerical model requires that a steady state is achieved, hence it was necessary to select a range where steady state was considered as achieved. For the selected measurements the ranges were: Measurement III. (12 500–18 000) with average operating pressure 6.156 MPa, measurement IV. (15 000–23 000) with average operating pressure 5.730 MPa, measurement VII. (10 000–17 300) with average operating pressure 4.992 MPa, and measurement X. (10 000–17 000) with average operating pressure 3.423 MPa. The measured velocity has a similar course as the calculated velocity. In the displayed steady-state range, the averaged deviation between the calculation and the measurements were as follows: Measurement III. − 4.964%, measurement IV. − 4.520%, measurement VII. − 6.311% and measurement X. − 4.456%.

Fig. 9 Distribution of temperatures along the loop for four selected measurements.

Fig. 10 Comparison of measured velocity (blue) with calculated values using numerical model (red).

7 The effect of leakages

Mathematical models, proposed by other authors, were tested on small scale models, where leakages are not prominent and the scaling factors didn't include possible leakages. One measurement (nr. X) was picked to represent the effect of leakage in this chapter. The steady state was achieved slowly, after approx. 11 000 seconds of operation.

The main reason of leakages is the small molecular volume of the helium, which easily leaks through any cracks on the loop. Having this in mind, the mass at every second of measurement was calculated simply as the product of mean density and the free volume at every section of the loop, using equation:$$m={\rho }_{m,GFR}{V}_{GFR}+{\rho }_{m,HB}{V}_{HB}+{\rho }_{m,DHR}{V}_{DHR}+{\rho }_{m,CB}{V}_{CB}\text{.}$$(19)

In the heater, hot and cold piping branch, the arithmetic mean of inlet and outlet temperatures was taken to calculate the density ρ_m . The density at the cooler was calculated from the logarithmic mean of inlet T ₁ and outlet T ₂ temperatures. Figure 11 shows the overall mass (violet) and leakage (black) time course during the measurement. The leakages Δm were calculated according to equation (20):$$\mathrm{\Delta }m=m_0-m\text{.}$$(20)

The device has openings, especially at the inserts of heating rods and cooling water pipes, which may not be sealed properly, and the coolant may leak out. This can be very dangerous, because if the coolant slowly leaks, it has to transfer more energy from the heater inside less mass, which results in its increased temperature. This effect of leakage can be observed in Figure 10, from approx. 12 000 s to 14 500 s.

The mass and leakages can be also represented by percentual values in every measured second as shown in Figure 12. The percentual mass value m_p and percentual leakage value Δm_p was calculated according to equations (21) and (22):$$m_p=\frac{m}{m_0}100\%$$(21) $$\mathrm{\Delta }{m}_p=\frac{\mathrm{\Delta }m}{m_0}100\%\text{.}$$(22)

As shown in Figures 11 and 12 the time course of leakage follows a seemingly linear character. The reason for this is yet uncertain, but the focus will be on unrevealing some exact relationship between the leakages and the state properties of coolant. The seemingly sudden rise/decline in mass during the measurement was caused by the quick change in the electrical input from GFR. This resulted in a short unsteady temperature distribution, where the temperature sensors at the GFR showed sudden temperature increase, but the temperature at the outlet of cooler had the same temperature.

Fig. 11 Mass and leakage time course during measurement X.

Fig. 12 Percentual mass and leakage time course during measurement X.

8 Conclusion

The article analyses a one-dimensional numerical model of flow in the helium loop with natural circulation. The flow is described as one-dimensional and steady. It turns out that a one-dimensional model can describe the flow in such a facility with sufficient accuracy. The difference between the calculated values and the selected measurements were about 4–7%. The deviations are apparently because the one-dimensional model does not take into account the flow to the y-axis and in particular, does not consider the backflows and swirls that almost certainly occur when helium enters into the DHR. We want to address the backflow problem in our further research, as well as the enhancement of the current numerical model into a thermodynamic model, capable of calculating the temperature distribution in the loop. The thermodynamic model is expected to have impaired accuracy due to deviations in temperature calculations. The thermodynamic model will be able to simulate the leakages, which appeared during measurements. This leakage will be incorporated only as an equation gained from a trendline, from measured data. Temperature and pressure have high influence on leakage and finding an equation for leakages is also in the scope of future work, but only after establishing the full thermodynamic model of flow in a natural circulation loop. No other researchers reported leakages on their own experimental devices, this may be due to the nature, how the heat input/output was solved in those devices. The heating and cooling of the coolant was made only by outside heaters/coolers around the pipes. As stated in the article, the presence of leakages means a high risk of overheating the coolant and the reason behind them have to be investigated more thoroughly.

Nomenclature

A : Cross-sectional surface [m²]

d : Diameter [m]

e : Euler's constant [–]

g : Constant of gravitational acceleration [m s⁻²]

h : Height [m]

L : Length [m]

m : Mass [kg]

p : Pressure [Pa]

Q_m : Mass flow [kg s⁻¹]

r : Specific gas constant [J kg⁻¹ K⁻¹]

T : Temperature [K]

V : Volume [m³]

v : Mean velocity [m s⁻¹]

z : Compressibility factor [–]

λ : Coefficient of length losses [–]

ρ : Density [kg m⁻³]

ξ : Coefficient of local losses [–]

Subscripts

0 : Initial value, value at 0 sec

Calc : Calculate value

CB : Reference to cold branch

DHR : Reference to DHR device

est : Estimated, or measured value as input

GFR : Reference to GFR device

H, HB : Reference to hot branch

i : Reference to spatial point on the loop

ie, m : Averaged value

j : Iteration cycle indicator

l : Reference to pressure losses

p : Percentual value

ε : Value calculated after small increment 10⁻⁵

Acknowledgements

This contribution was created on the basis of the project “Research centre ALLEGRO” (ITMS project code: 26220220198), supported by Operational Programme Research and Development funded by the European Regional Development Fund. The authors gratefully acknowledge the contribution of the Scientific Grant Agency of the Slovak Republic under the grant VEGA 1/0743/18.

References

All Tables

All Figures

	Fig. 1 Diagram of the experimental facility.
In the text

	Fig. 2 Helium velocity measurement.
In the text

	Fig. 3 Time course of temperature and pressure of the coolant.
In the text

	Fig. 4 Time course of velocity in the hot branch (the most accurate).
In the text

	Fig. 5 Time course of the outside insulation temperature.
In the text

	Fig. 6 Time course of mass flow.
In the text

	Fig. 7 Time course of thermal power gained and given to/from the coolant inside GFR and DHR.
In the text

	Fig. 8 Determination of local loss coefficient for DHR.
In the text

	Fig. 9 Distribution of temperatures along the loop for four selected measurements.
In the text

	Fig. 10 Comparison of measured velocity (blue) with calculated values using numerical model (red).
In the text

	Fig. 11 Mass and leakage time course during measurement X.
In the text

	Fig. 12 Percentual mass and leakage time course during measurement X.
In the text