Science and Technology of Nuclear Installations
Volume 2018 (2018), Article ID 7862847, 11 pages
https://doi.org/10.1155/2018/7862847

Research Article

Study of Fast Transient Pressure Drop in VVER-1000 Nuclear Reactor Using Acoustic Phenomenon

Soroush Heidari Sangestani,¹ Mohammad Rahgoshay,¹ Naser Vosoughi,² and Mitra Athari Allaf¹

¹Department of Nuclear Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran
²Department of Energy Engineering, Faculty of Engineering, Sharif University of Technology, Tehran, Iran

Correspondence should be addressed to Mohammad Rahgoshay; m.rahgoshay@gmail.com

Received 10 June 2017; Revised 3 September 2017; Accepted 10 September 2017; Published 31 January 2018

Academic Editor: Eugenijus Ušpuras

Copyright © 2018 Soroush Heidari Sangestani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This article aims to simulate the sudden and fast pressure drop of VVER-1000 reactor core coolant, regarding acoustic phenomenon. It is used to acquire a more accurate method in order to simulate the various accidents of reactor core. Neutronic equations should be solved concurrently by means of DRAGON 4 and DONJON 4 coupling codes. The results of the developed package are compared with WIMS/CITATION and final safety analysis report of Bushehr VVER-1000 reactor (FSAR). Afterwards, time dependent thermal-hydraulic equations are answered by employing Single Heated Channel by Sectionalized Compressible Fluid method. Then, the obtained results were validated by the same transient simulation in a pressurized water reactor core. Then, thermal-hydraulic and neutronic modules are coupled concurrently by use of producing group constants regarding the thermal feedback effect. Results were compared to the mentioned transient simulation in RELAP5 computer code, which show that mass flux drop is sensed at the end of channel in several milliseconds which causes heat flux drop too. The thermal feedback resulted in production of some perturbations in the changes of these parameters. The achieved results for this very fast pressure drop represent accurate calculations of thermoneutronic parameters fast changes.

1. Introduction

One of the most important aspects in design of different safety systems with sufficient preparation is simulation and analysis of transient states in reactor core. For these kinds of analysis basic equations in neutronic and thermal-hydraulic modules have to be coupled. Coupling of neutronic and thermal-hydraulic modules is different from each other, considering numerical solution methods and time and body meshing size. Thus, written codes for different transient states are mostly used for study of fuel and coolant temperature changes, power peak level, coolant pressure, stability time, and other parameters. Computer coding submits models with different degree of accuracy and validity. Most codes are not able to analyze coupled conditions of very fast transient (FT) states in very short time. This deficiency is associated with neutronic and thermal-hydraulic calculations or both. Therefore designing a code which is capable of analyzing FT conditions is highly needed.

One of the efficient ways in analyzing FT is using waveform method. Chan (1991) studied asymptotic waveform evaluation in analysis of time-dependent calculations [1]. Ooi et al. (2003) studied the finite element method (FEM) in thermal analysis that usually produces a formulation in the space/time domain. However, the sizes of the equations in FEM usually are large, and thus the conventional algorithms involve considerable computational time. The conventional methods have to take a very small time step size to avoid undesirable numerically induced oscillations or numerical instabilities. Thus, a new solution algorithm, named the asymptotic waveform evaluation scheme, was introduced by them to solve transient problems [2]. Liu et al. (2006) studied fast thermal simulation when power density increases by fast spectrum in frequency domain for computing steady state response. It indicates that the resulting thermal analysis algorithms lead to 10x–100x speedup over the traditional integration-based transient analysis with small accuracy loss. The studied parameters minimum time order is 15 ms [3]. Ham and Bathe (2012) used FEM to solve time-dependent two-dimensional wave propagation problems [4]. Ranaa et al. (2014) presented the well condition asymptotic waveform evaluation to solve heat conduction problem in the frequency domain. The method is presented for time-dependent problems [5]. In addition, Ishii et al. (2009) investigated the effect of acoustic phenomenon in steam dryer (in a Boiling water reactor) which indicates the process of pressure pulsations caused by hydroacoustic resonance propagation along the steam dryer [6]. Proskuryakov (2017) studied the effect of acoustic vibrations in the nuclear power plant (NPP) coolant such as VVER-1000 reactor and created new scientific direction “diagnosis, prognosis and prevention of vibration - acoustic resonances in the NPP equipment” shows that the developed methods can be used to predict and prevent the occurrence of vibration-acoustic resonances in the NPP equipment [7].

Various kinds of thermal-hydraulic and neutronic calculating models are put to use in transient calculating code development studies. Reducing costs and runtime and achieving required accuracy are three main purposes of them. Leung et al. (1981) studied acoustic impact techniques for increasing the accuracy of FT states modeling in CODA code [8]. However, regarding its very small meshing, acoustic methods are time-consuming and are not used any more.

In order to decrease complicated solving parameters in using compressible fluid method, geometry of every fuel assembly could be turned into a single heated channel (SHC). Four different methods are used to solve SHC transient equations. They are channel integral (CI), single mass velocity (SV), momentum integral (MI), and sectionalized compressible fluid (SC). SC model considers both sound impact and thermal expansion, while the other three models only consider thermal expansion [9]. The SHC method is used in costanza code (1967) in order to analyze the dynamics of liquid cooled nuclear reactor [10]. It is also used in DynCo code (2011) which is intended for complex neutron-physic and thermal-hydraulic dynamic calculation of the reactor core in 3D hex-Z geometry. DynCo code is designed to simulate the dynamic behavior of the Russian 3000-MWt pressurized water reactor (PWR) [11]. The coupling of CI and POWEX-K code (2011) simulates power excursion in reactor of Budapest University of Technology and Economics Reactor [12]. Hosseini et al. (2015) developed a coupled 3D neutronic with 1D thermal-hydraulic model in order to find reasonable power distribution. Neutronic module includes three-dimensional diffusion equation in two energy groups which was solved using analytic nodal method. Thermal-hydraulic module contains the conservation equations solved for 1D axial homogeneous downward flow through channel using SHC model (MI method). In conclusion, Xenon saturation transient analysis of a research reactor core was carried out [13].

The SHC compressible fluid method is one of the ways that make it possible to use wave propagation method and acoustic phenomenon. Shapiro (1953) considered the derivation and properties of the dynamics and thermodynamics of compressible fluid flow which is along with acoustic theory [14]. Meyer (1961) investigated that several approximations can be used to decouple the momentum and energy equations to facilitate solution of the transient problem. Additionally, the numerical solution of a transient problem would be particularly simplified if the compressibility of the fluid could be ignored [15]. Todreas and Kazimi (2001) investigated a rapid increase in the heat flux without change in the applied pressure drop in a PWR. The SC approach for a 10% heat flux step increase in the PWR channel shows that because the pressure begins to rise at internal channel points, a reduction in inlet flow rate and an accelerated exit flow rate occur [9].

The mechanism of compressible fluid method was published by Bar-Meir in 2007 [16]. Khola and Pandey (2013) studied the numerical simulation of transients in two-phase flow which is crucial to simulate accident-like conditions of nuclear reactors for safety analysis. In their work, a code for numerical computation of unsteady one-dimensional two-phase flow has been developed. The governing equations were obtained by the homogenous equilibrium mixture model and were decoupled and approximated using the SC and MI model [17].

Numerical considerations (i.e., the stability and/or accuracy) of the difference solution require that the time step of integration be less than the time interval for sonic wave propagation across the spatial grid points. Compared with the transport velocity, the fluid sonic velocity is large. It causes limitation of the time step in most numerical schemes to very small values. Acoustic phenomenon is used in accommodation of very short time and very small body meshing. This accommodation is determined by Courant’s criterion [8, 20]:

is time meshing period, is body meshing distance, is sound velocity in coolant, and is the mean transport velocity. The fluid sonic velocity is large comparing with the transport velocity. Therefore, limiting the time step in most numerical schemes to very small values leads to a computationally expensive solution of this problem. Regarding meshing size criterion, pressure wave propagation method mostly requires long-term computations for numeral stability. Therefore acoustic phenomenon features (despite the high accuracy) are not usually considered by researchers. However these features are very useful during core parameter vital changes.

FT pressure drop is a type of loss of coolant accident (LOCA). That is well known as the double ended guillotine break. When double ended guillotine break occurs (main coolant pipeline cold leg breaks at the reactor inlet), suddenly reactor coolant pressure decreases with leak coolant flow rate of 45000 Kg/s [18]. A pressure wave is also produced, while the large break LOCA occurs, which propagates across the channel. Therefore International Atomic Energy Agency (IAEA, 2003) studied the before-break vital moment (leak before break) [21].

The coolant fast pressure drop accident can lead the fluid to become two-phased and thermoneutronic parameters to change. Gonzalez-Santalo and Lahey Jr. (1972) investigated this matter by study of pressure drop transient in two-phase condition [22].

Calculation program of VVER-1000 reactor core FT pressure drop by means of SC method and acoustic phenomenon was developed in this investigation. In order to use the mentioned method every one of the 28 fuel assemblies should be considered as one SHC. Fuel assembly conversion into SHC and meshing method toward -axis are both shown in Figure 1.

Figure 1: Fuel assembly conversion into SHC and meshing method toward -axis.

2. Materials and Methods

2.1. An Overview of Neutronic and Thermal-Hydraulic Model

A computational program for simulation of VVER-1000 reactor core FT state of sudden pressure drop has been developed in this study. It includes the two thermal-hydraulic and neutronic basic models. Group constants of the fuel assemblies and reflectors are produced by DRAGON 4 code [23]. DRAGON 4 is a lattice code developed to solve Boltzmann transport equation in two and three dimensions, to apply self-shielding effects and to compute few-group macroscopic cross sections and diffusion coefficients. The mentioned program includes different models for simulation of fuel assembly behavior, such as interpolation of microscopic cross sections, resonance self-shielding calculation, different solvers for the Boltzmann transport equation with ability to take into account leakage effects, and calculation of condensed and homogenized parameters.

In this study, the integral transport equation is solved by the SYBILT module, self-shielding calculations are performed by the SHI module by means of generalized Stamm’ler method, and the CPO module is utilized for production of equivalent fuel assembly parameters in consistent format that can be used in forgoing calculation. After that, time-dependent multigroup neutron diffusion equations are solved by DONJON 4 models [24]. DONJON 4 is a multigroup diffusion solver. A three-dimensional simulation is performed by Thomas–Raviart–Schneider method using The TRIVAT module, the group constants of fuel assembly calculated by DRAGON 4 are recovered using the CRE module, and the FLUD module is to compute multiplication factor. The FEM is utilized for discretization of equation. Finally momentary amount of power distribution is achieved in 1/6 of reactor core, across every fuel assembly.

In thermal-hydraulic model, conservation of mass and momentum and energy dependent equations are solved by applying considered transients and use of channel axial division in MATLAB software. The channel is regarded as a typical coolant subchannel inside an assembly that only receives coolant through its bottom inlet. The fuel and clad heat transport equations are solved separately from coolant equations. One-dimensional transient transport equations of the coolant with radial heat input from the clad surfaces are resolved. The flow area is assumed axially uniform but pressure loss due to local area changes is still regarded. Any axial position of flow area could be considered as control area in order to derive radially averaged coolant flow equations. Mentioned equation is envisioned for SHC state and SC method of every 28 fuel assemblies (1/6 core). From the four SHC different methods only SC model considers both sound impact and thermal expansion, while the other three models considering thermal expansion [9]. Therefore it is able to study very FT states. Some of the thermodynamic features of water are also regarded [25].

Several approximations can be used to decouple the momentum and energy equations to facilitate solution of this transient. Additionally, the numerical solution of the discussed transient problem in this study is particularly simplified following the work of Meyer [15] and specific applications from Todreas and Kazimi [9]. Therefore the lateral variations of properties in the flow channels can be neglected. For this condition, (2), (3), and (4) are directly applicable:

where is enthalpy, is mass flux, is density, is fluid level passing, is of fuel assembly height, is time, is friction coefficient, is heat flux, is thermal-hydraulic diameter, is speed sound, is pressure, and is heated surface. Equations (2), (3), and (4) and appropriate constitutive relations provide the solutions for , , and for the initial and boundary conditions. The heat flux , which is a DONJON output, is specified as an input for the above equations. The boundary conditions for solving the momentum and energy equations are to be specified as and at the inlet and at the outlet. Furthermore, constitutive equations for and are required to complete definition of the problem. The equation of state for the density, assumed differentiable with respect to and , is specified as

The friction factor can be specified as

In SC model, numerical solution approach involves a set of difference equations representing the differential transport equations, arranged to consider , , and , and state variables at multiple points along the channel. Every one of the 28 fuel assemblies is considered as a single heated channel. The term sectionalized reflects the need to divide a channel into 360 segments to execute the numerical solution, regarding the height of every fuel assembly (3.6 m). Using (5) we get

where

From (2) and (7) we get

Equations (4) and (9) may be combined into two equations by eliminating and :

where we have defined as

Equations (3) and (10) are partial differential equations in , , and (i.e., the density does not appear as a differentiated variable). These equations were written in pointwise difference form (forward implicit method) and solved in , , and .

Sound velocity is different in various spaces; also the fluid is running in a specified direction with a specified speed. Fluid velocity is accordingly dependent on its other parameters. Regime pressure drop transient numerical considerations (i.e., the stability and/or accuracy) of the difference solution need that the time step of integration be less than the time interval for sonic wave propagation across the spatial grid points, that is, (1). The time discretization that is used in most codes varies from the semi-implicit scheme and multistep schemes or the fully implicit discretization.

Different codes have developed various methods for predicting theses states. CATHARE code [26] used critical flow rate correlation method; that way the geometry of the flow restriction is simplified and a coarse meshing is used. TRACE [27] and RELAP5 [28] codes use characteristic velocity; that way the geometry of the flow restriction is also simplified and a coarse meshing is used. A sound characteristic velocity is set to zero and simplified equations are used to predict flow evolution from upstream to the sonic section. The standard option in RELAP5 is a semi-implicit scheme. Only those governing parameters are evaluated at new time values which are needed to maintain stability and accuracy of the method or to avoid Courant time step size limitations based on pressure waves. CATHARE code [26] in another method assumes that the flow from upstream to sonic section is precisely calculated by 2-fluid equations using a very fine meshing in the vicinity of the throat. Such a method is only possible with an implicit numerical scheme which allows high velocity flow in small meshes without material Courant limitation.

All these codes, in which Courant criterion is not used, consider different ways in order to decrease the time of calculations but still they do not have the accuracy of the codes using Courant criterion, in time of critical flow conditions and very short times. Therefore Courant criterion was taken into account in the simulating codes of this study. Time mesh amount is obtained about 10 μs considering for solution stability and sound velocity of 900 m/s in coolant. DRAGON and DONJON/SHC developed code is regarded 1 μs for more accuracy.

2.2. Very FT Utility Analysis

Very FT coolant pressure drop in VVER-1000 reactor core is explained in

It is highly mentioned that axial fluid motion is considered steady in SC model which eliminates cross flow effect [9]. By considering core symmetry, one of the six parts of core is modeled, which includes 28 fuel assemblies’ channels. Therefore, despite the ones where core is assumed as one channel, radial accuracy of thermal-hydraulic and neutronic parameters is increased. The impact of thermal feedback on cross sections is utilized for coupling the two mentioned modules. Eventually standard deviation of predicted results for noted transient is estimated through accuracy evaluation of thermal-hydraulic and neutronic codes. Some advantages of SHC method by SC model are increasing the sufficient accuracy in sudden transient state simulations and more accurate analysis of different transient states. As mentioned before, acoustic phenomenon increases the number of time meshes, and therefore the time of calculation becomes longer. On that account, performing in FT states is the optimum condition for using noted code (considering user demands).

2.3. Coupling Procedure

Coupling algorithm of the two thermal-hydraulic and neutronic modules is as follows. Primarily every fuel assembly is modeled by DRAGON 4 cell calculating code according to Table 1. This code is used (considering fresh fuel) in order to determine the macroscopic cross sections sequences and other neutronic parameters of all reactor fuel assemblies. Instead of modeling the whole core, 28 fuel assemblies are modeled according to core symmetry. The next step is obtaining axial and radial power distributions by means of available triangular geometry module of DONJON 4 and numeral solving by means of FEM, also applying VVER-1000 reactor power. Considering nonpoisonous fresh fuel, reactor relative power radial distribution is measured in order to evaluate the calculation algorithm authenticity of the two codes’ correlation and validated by comparing to WIMS/CITATION [19] results and Bushehr VVER-1000 reactor FSAR [18].

Table 1: VVER-1000 geometry and operating conditions [18].

Operating condition
Value

Core height in the working state (cm)
355

Rod diameter (mm)
9.1

Pitch (mm)
12.75

Fuel pellet material
UO₂

Number of fuel rods in the fuel assembly
311

Average linear heat rate (W/cm)
166.7

Coolant flow rate (m³/h)
84800

Inlet pressure (MPa)
16

Outlet pressure (MPa)
15.70

Temperature at the inlet (°C)
291

In order to evaluate the authenticity of thermal-hydraulic module, the mentioned transient is modeled in PWR core fuel assembly (conditions are noted in Table 2) and its obtained results are compared to the Todreas and Kazimi [9] study results. Core is composed of 28 fuel assemblies; but only Hot Fuel Assembly (HFA) coupling results are noted in this assay. Thermal-hydraulic properties of coolant fast pressure drop accident in PWR core is simulated by SC method. The amounts of time periods which are caused by environmental sound velocity limitations are also considered for current simulation. After evaluating the authenticity of thermal-hydraulic module, coupling results of the two neutronic and thermal-hydraulic modules of coolant fast pressure drop accident in VVER-1000 reactor core could be studied for validation which is compared to RELAP5 (advanced-thermal-hydraulic code, coupled with point kinetic model to simulate neutronic behavior) [28]. The thermal-hydraulic and neutronic coupling module algorithm is as follows: the procedure starts with obtaining the axial power profile of neutronic module in steady state. The mentioned axial power profile and other core dependent parameters (sudden pressure drop transient, mass flux, pressure, temperature changes, etc.) are the entering data of thermal-hydraulic module program. Axial thermal changes affect the nuclear cross sections and the new power profile in neutronic module is produced. Then the new axial power profile is applied on thermal-hydraulic module. It is continued until the difference between axial thermal changes in every time step is less than a specified amount of convergence criteria. Neutronic and thermal-hydraulic coupling calculations algorithm are shown in Figure 2.

Table 2: PWR geometry and operating conditions [9].

Operating condition
Value

Channel length (m)
3.66

Rod diameter (mm)
9.70

Pitch (mm)
12.80

Flow area for rod (mm²)
90.00

Equivalent diameter (mm)
12.00

Linear heat (kW/m)
17.50

Mass flux (kg/m²s)
4125

Inlet pressure (MPa)
15.50

Outlet pressure (MPa)
15.42

Inlet enthalpy (kJ/kg)
1337.2

Figure 2: Neutronic and thermal-hydraulic coupling calculation algorithm.

As a result of time period selection limitations (short time periods), which leads to increase in the number of nodes, the rise in parameter numbers of defined study case and complicated equations or scales should not be considered. Therefore program solving time is lowered.

3. Results and Discussion

3.1. Results of Neutronic Calculations

Relative radial power distribution in 1/6 of VVER-1000 reactor core is shown in Table 3. Power Peaking Factor is compared to WIMS/CITATION [19] and Bushehr VVER-1000 reactor FSAR value [18], as observed. Obtained results show an unremarkable calculation error for the two neutronic coupled codes.

Table 3: Relative power radial distribution in 1/6 of VVER-1000 reactor core.

Fuel assembly number
Power peaking factor (FSAR) [18]
Power peaking factor
(DRAGON 4/DONJON 4)
Power peaking factor
(WIMS/CITATION) [19]

1
1.01
1.02
0.99

2
0.80
0.79
0.78

3
1.02
1.01
1.00

4
0.75
0.75
0.72

5
0.86
0.85
0.83

6
0.92
0.89
0.90

7
1.16
1.17
1.18

8
0.92
0.93
0.90

9
0.78
0.78
0.76

10
0.87
0.88
0.84

11
0.86
0.85
0.84

12
1.24
1.22
1.24

13
0.96
0.97
0.98

14
0.78
0.77
0.76

15
1.03
1.04
1.02

16
0.89
0.88
0.87

17
1.11
1.10
1.11

18
1.18
1.17
1.24

19
0.87
0.86
0.84

20
0.89
0.87
0.87

21
1.29
1.27
1.31

22
1.24
1.25
1.31

23
0.86
0.85
0.84

24
1.11
1.09
1.11

25
1.24
1.23
1.31

26
1.24
1.22
1.24

27
1.18
1.17
1.24

28
0.96
0.96
0.98

The maximum coupling calculation errors of DRAGON 4/DONJON 4 is 2.8% while it is 5.64% in WIMS/CITATION [19] codes. The obtained results of DRAGON 4 and DONJON 4 codes coupling shows higher accuracy than WIMS/CITATION codes as observed. Then relative axial power distribution in 1/6 of VVER-1000 reactor core is analyzed, which is shown in Figure 3. As observed, the peak of relative axial power distribution (in beginning of cycle) occurs lower than the middle height of core. Relative axial power distribution is more coordinated to Bushehr VVER-1000 reactor FSAR [18] when using DRAGON 4/DONJON 4 coupling rather than WIMS/CITATION coupling [19].

Figure 3: Relative axial power distribution in 1/6 of VVER-1000 reactor core.

3.2. Results of Thermal-Hydraulic Calculations

After the defined model is validated, thermal-hydraulic module preparation is started based on SHC by SC method. In order to evaluate the authenticity of newly made module, coolant fast pressure drop accident in PWR core fuel assembly is simulated by SC method (PWR core features are noted in Table 2). This simulation also considers different amounts of time periods which are caused by environmental limitations of sound velocity. Severe decreases in entrance pressure (also neutronic changes stability) cause mass flux reduction which is distributed in the channel during spent time [9, 29]. The outputs of thermal-hydraulic program, which are normalized mass flux changes during axial side of the channel, are shown in Figure 4. Sound velocity is calculated 900 m/s, causing channel transient state to last for 4.1 ms. The mass flux fast reduction of  ms would not reach the end of the channel; therefore upper sides will not be affected by pressure wave.

Figure 4: Short-term response of the PWR inlet pressure transient using the SC model.

For  ms, the reflected pressure wave affects the mass flux. It happens because of stable outer pressure theory, which means when opposite pole wave (with similar amplitude) moves in opposite direction, the diluted entrance wave is reflected like compressed wave in outer border. Pressure profile in  ms and longer times shows the effect of reflection wave progress all the way from exit to entrance. In that short amount of time, average mass flux rate only decreases from 1 to 0.98. Figure 4 shows the results of program performance and has an acceptable accommodation with reference diagram [9], as observed.

The written program accuracy in neutronic states is compared to Bushehr VVER-1000 reactor FSAR data [18]. Also it is compared to PWR data [9] in thermal-hydraulic states. After the authenticity evaluation of neutronic and thermal-hydraulic modules, correlated FT states and changes of various parameters during time are studied.

3.3. Results of Coupling Calculations

Changing diagrams of thermal-hydraulic parameters by using sound effect in pressure propagation, for very FT coolant pressure drop in VVER-1000 reactor HFA, are being studied. Pressure drop causes the mass flux reduction; therefore temperature increases and power drops in every node. The transfer of each node’s thermal data to DRAGON 4 module leads to the decrease of cross section moderation. Mentioned data and axial power profile are used as DONJON 4 module entrance. Finally lowered power distribution and pressure drop transient are applied as thermal-hydraulic module and thermal distribution of every one of 360 nodes is obtained again. The process continues until the convergence criteria of each node thermal difference are available. Considering the central pressure of each node and reactor core pressure changes, which is 0.3 MPa, alterations of pressure per time for the first node is shown in Figure 5; that reveals pressure drop of 0.0012 MPa per 1 ms.

Figure 5: Alterations of pressure per time for the first node.

Figure 6 shows mass flux changes during fuel assembly in three different times. Increasing the time from 10 μs to 2 ms causes a mass flux drop up to more than half of the fuel assembly. Mass flux drop in 10 μs, which is 0.35 Kg/m²s, occurs up to the beginning of fuel assembly. However, these changes are not sensed by the upper nodes. After 1 ms, mass flux drop arrives at the middle of fuel assembly, which increases from 0.35 Kg/m²s to 12 Kg/m²s for the first node and decreases to 25 Kg/m²s after 2 ms. The channel transient time of sonic wave is about 4 ms. The rapid decrease of the mass flux has not yet reached the end of the channel in this period of time; that means the upstream region is not yet affected by the pressure wave.

Figure 6: Mass flux changes during fuel assembly in different three times.

Mass flux changes in 5 ms indicate the wave reflux from end of the channel. From now on, for a few milliseconds, the reversible wave causes more mass flux drop in the end of channel than the beginning of it. Mass flux changes in 5 and 9 ms, which are calculated by RELAP5 [28], are shown in the diagram. RELAP5 does not sense pressure drop transient in times before 4.2 ms, and therefore it has some delay comparing to DRAGON 4 and DONJON 4/SHC codes and lesser time periods are not possible to be compared.

The observed difference between DRAGON 4 and DONJON 4/SHC diagram (in 5 ms) and RELAP5 diagram (in 5 ms) is because of the difference between neutronic and thermal-hydraulic modules accuracy. Also the effects of thermal feedback cause perturbation in mass flux alterations, in DRAGON 4 and DONJON 4/SHC code. RELAP5 [28] results only show the average mass flux drop caused by pressure drop. Average mass flux changes with steady slope have four reasons:

(1)Driving force drop caused by inlet pressure drop that reduces the channel average mass flux.(2)The fluid thermal expansion due to heating which causes a small slope change; that is, becomes less than .(3)The fluids being considered uncompressible (i.e., ).(4)Weak neutronic model and not using sufficient time meshes sizes for the mentioned transient.

The obtained results from mass flux changes simulation in primary, middle, and ending points of relative length of fuel assembly, in 9 ms, in DRAGON 4 and DONJON 4/SHC and RELAP5 code are shown in Figure 7. Mass flux in the first node initially faces a severe drop which continues for 2 ms. When this drop becomes more visible in upper nodes, diagrams slope turns more gradual along with average drop of 0.068 (Kg/m²s). Mass flux in central point of fuel assembly starts dropping after whereabouts 1.2 ms. However severe drop for ending points occurs after 4 ms. The profiles at  ms and more show the effect of the reflected wave progressing toward the channel exit from the inlet. During this short period of time, the mass flux only decreases 60 Kg/m²s. The net mass flux is the superposition of the forward wave and the reflected wave.

Figure 7: Mass flux changes during fuel assembly in different three nodes.

RELAP5 results show mass flux profile changes, with steady slope, in primary, middle, and ending points and does not present adequate description of the passing pressure wave.

HFA heat flux changing per time for beginning, middle, and ending nodes of fuel assembly, in “0–3 ms” period of time, is shown in Figure 8. The heat flux of first node drops whereabouts 40 KW/m² in that amount of time. However, mass flux of fuel assembly middle node drops whereabouts 30 KW/m² in 1.7 ms. It also occurs after 2.3 ms for ending node. Temperature changes per time for the first node are shown in Figure 9, which face drop of 0.01°C in 0.06 ms after the beginning of the accident. Mentioned drop is inevitable, based on constant axial power decrease.

Figure 8: HFA heat flux changing per time for beginning, middle, and ending nodes of fuel assembly.

Figure 9: Temperature changes per time for the first node.

4. Conclusions

Simulation of VVER-1000 reactor core FT pressure drop, by using acoustic phenomenon, was studied in this article. Coupling of the two DRAGON 4 and DONJON 4 neutronic modules, along with thermal-hydraulic module by means of SC method, was utilized in this simulation. It is the first time that neutronic and thermal-hydraulic modules are externally coupled, assuming that the coolant fluid is compressible, in order to obtain the best simulation for this transient.

In this transient, at the channel inlet, the mass flux begins to decrease because of the pressure reduction; it was observed in calculations that a perturbation was produced and propagates into the channel as time elapses. After pressure waves arrive at the end of the channel, reflected pressure waves affect the mass flux profile for short time. The incoming rarefaction wave will be reflected as a compression wave, due to the assumption of constant outlet pressure, at the exit boundary. Therefore a wave travels in the opposite direction with the same amplitude but opposite sign. The reverse wave causes more mass flux drop in upper nodes and as time passes, it decreases the difference between upper nodes mass flux drop and the ones in beginning of the channel. Thus it is concluded that, in such transients, mass flux drop rate does not follow stable conditions.

This transient (double ended guillotine break) is very fast, and therefore DRAGON 4 and DONJON 4 codes are used in neutronic module (comparing to WIMS/CITATION) considering their high speed and accuracy in analyzing the problem. Also SC method is utilized in thermal-hydraulic module for solutions of conservation equations by forward implicit method and Courant’s criterion in a SHC. Thus the very short time step (μs grades) was used in analyzing transient. The neutronic module was validated by WIMS/CITATION and FSAR. Thermal-hydraulic module was also validated by simulation of a PWR (comparing to Todreas and Kazimi [9] simulated one). Finally the simulated pressure drop transient (by DRAGON and DONJON/SHC code) was also simulated in RELAP5 code and results were compared. RELAP5 code was not capable of sensing the mentioned transient in times before 4.2 ms and, unlike DRAGON and DONJON/SHC code, coolant fluid is assumed uncompressible in it. So changes of density per pressure is considered zero. Of course this assumption is highly reasonable for classes of core transients that do not have the problem of losing a great amount of coolant. Therefore mass flux changes slope is almost steady along the channel and mass flux changes are only dependent on pressure changes in channel inlet and outlet. Thus it became completely obvious that, despite the fact that RELAP5 code uses “nearly implicit” scheme based on a “fractional step” approach for providing a sufficient degree of implicitness to eliminate the material Courant-type stability limit, but it does not give an adequate description of the mentioned studied transient.

The obtained results from DRAGON and DONJON/SHC for the times that a great volume of coolant is lost in a short time, assuming coolant fluid is compressible, are highly acceptable and they show its ability to effectively and accurately calculate the fast transient. On account of that, the results are not desirable if the sonic effect in coolant is not regarded.

Abbreviations

FT: | Fast transient
FEM: | Finite element method
PWR: | Pressurized water reactor
CI: | Channel integral
SV: | Single mass velocity
MI: | Momentum integral
SC: | Sectionalized, compressible fluid
NPP: | Nuclear power plant
FSAR: | Final safety analysis report
LOCA: | Loss of coolant accident
SHC: | Single heated channel
HFA: | Hot fuel assembly
ms: | Millisecond
μs: | Microsecond
IAEA: | International Atomic Energy Agency.

Conflicts of Interest

All authors have no conflicts of interest to declare.

Acknowledgments

This study was performed as part of long-term research into nuclear safety supported by the Nuclear Engineering Department of Science and Research Branch, Islamic Azad University. The valuable researches into various accidents of VVER-1000 reactor led us into performing this study.

References

P. K. Chan, “Comments on “Asymptotic Waveform Evaluation for Timing Analysis”¹,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 10, no. 8, pp. 1078-1079, 1991.

C. Ooi, K. Seetharamu, Z. Alauddin, G. Quadir, K. Sim, and T. Goh, “Fast transient solutions for heat transfer [FEM],” in Proceedings of the IEEE TENCON 2003. Conference on Convergent Technologies for the Asia-Pacific Region, pp. 469–473, Bangalore, India, 2003.

P. Liu, H. Li, L. Jin, W. Wu, S. X.-D. Tan, and J. Yang, “Fast thermal simulation for runtime temperature tracking and management,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 25, no. 12, pp. 2882–2892, 2006.

S. Ham and K. Bathe, “A finite element method enriched for wave propagation problems,” Computers & Structures, vol. 94-95, pp. 1–12, 2012.

M. S. Ranaa, K. Jeevanb, R. Harikrishnanc, and A. W. Rezad, “A well-condition asymptotic waveform evaluation method for heat conduction problems,” Advanced Materials Research, vol. 845, pp. 209–215, 2014.

K. Ishii, T. Nagayoshi, S. Takahashi, Y. Wada, and N. Tanikawa, “Advanced simulation technologies for nuclear power plants,” Hitachi Review, vol. 58, no. 2, pp. 94–101, 2009.

K. N. Proskuryakov, “Scientific Basis for Modeling and Calculation of Acoustic Vibrations in the Nuclear Power Plant Coolant,” Recent Advances in Petrochemical Science, vol. 2, no. 1, 2017.

J. C. M. Leung, K. A. Gallivan, and R. E. Henry, Critical Heat Flux Predictions During Blow down Transient, Argonne National Laboratory, Argonne, IL60439, U.S.A, 1981.

N. E. Todreas and M. S. Kazimi, “Thermal Hydraulic Fundamentals,” Massachusetts Institute of Technology, 1990.

G. Forti and E. Vincenti, “The codes costanza for the dynamics of liquid cooled nuclear reactor, joint nuclear research center Ispra stablishment - Italy,” Reactor Physics Department Reactor Teory and Analysis, 1967.

A. V. Levchenko, I. N. Leonov, and S. V. Kovalchuk, “Verification of the three dimensional dynamic DYNCO code on the basis of international test cases for pressurized water reactors,” Simulation Systems Ltd, 2011.

S. M. Khaleda and F. Al Mutairic, “A Numerical Technique For Solving Coupled Thermal hydraulic and Multi-Energy Group Neutron Diffusion Equations,” A Department of Studies and Basic Sciences, Community College, University of Tabuk, Saudi Arabia, 2011.

M. Hosseini, H. Khalafi, and S. Khakshournia, “Two group, three dimensional diffusion theory coupled with single heated channel model: A study based on xenon transient simulation of Tehran research reactor,” Progress in Nuclear Energy, vol. 85, pp. 108–120, 2015.

A. H. Shapiro, The Dynamics and Thermodynamics of Compressible Fluid Flow, vol. 1, Ronald Press Co., New York, NY, USA, 1953.

J. E. Meyer, “Hydrodynamic models for the treatment of reactor thermal transients,” Nuclear Science and Engineering, vol. 10, pp. 269–277, 1961.

G. Bar-Meir, “Fundamentals of Compressible Fluid Mechanics, 2007,” Version 0.4.4.2 aka 0.4.4.1j, May 21. Minneapolis, MN 55414-2411.

N. Khola and M. Pandey, “Numerical computation of one-dimensional unsteady two-phase flow using hem model and IAPWS IF-97 equations of state,” in Proceedings of the 2013 21st International Conference on Nuclear Engineering, ICONE 2013, China, August 2013.

AEOI., Reactor Final Safety Analysis Report VVER-1000 Bushehr, Chapter 4, Atomic Energy Organization of Iran, 2005.

A. H. Fadaei and S. Setayeshi, “Control rod worth calculation for VVER-1000 nuclear reactor using WIMS and CITATION codes,” Progress in Nuclear Energy, vol. 51, no. 1, pp. 184–191, 2009.

J. Plesek, R. Kolman, and D. Gabriel, “Estimation of The Critical Time Step for Explicit Integration,” in Proceedings of the 18th International Conference, Engineering Mechanics, Svratka, Czech Republic, 2012.

IAEA-TECDOC-1361, Assessment and management of ageing of major nuclear power plant components important to safety Primary piping in PWRs, Scientific and Technical Publications, Vienna, Austria, 2003.

J. M. Gonzalez-Santalo and R. T. Lahey Jr., An Exact Solution of Flow Decay Transients in Two-Phase Systems by The Method of Characteristics, The American Society of Mechanical Engineers, New York, NY, USA, 1972.

G. Marleau, A. Hébert, and R. A. Roy, “A user guide for DRAGON Version 4,” Institute of Genius Nuclear, Department of Genius Mechanical, School Polytechnic of Montreal, 2011.

D. Sekki, A. Hébert, and R. Chambon, “A User Guide for Donjon 4 Version 4,” Institute of Genius Nuclear, Department of Genius Mechanical, School Polytechnic of Montreal, 2011.

W. Wagner and H.-J. Kretzschmar, International Steam Tables, Faculty of Mechanical Engineering Chair of Thermodynamics, 2nd edition, 2007.

G. Lavialle, “The CATHARE 2 V2.5 code: Main features, CATHARE-NEPTUNE International Seminar,” Grenoble, May 2004.

TRACE-V5.0, Theory Manual, 2007, TRACE V5.0 User Manual, 2007, TRACE V5.0 Assessment Manual, 2007.

RELAP5/SCDAP//MOD3.2 Code Manuals, A Computer Code for Best-Estimate Transient Simulation of Light Water Reactor Coolant Systems During Severe Accidents, Prepared for the U.S. Nuclear Regulatory Commission, Idaho National Engineering and Environmental Laboratory, NUREG/CR-6150, 1997.

M. Rahgoshay, Encyclopedia of Advanced Subjects on Nuclear Heat Transfer, vol. 1, Chapter 3, Islamic Azad University Central Tehran Branch, Iran, 2016, ISBN: 978-964-10-4389-8.