Investigation of a LOCA in a typical MTR by a novel best-estimate code

Progress in Nuclear Energy 86 (2016) 141e161

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

Investigation of a LOCA in a typical MTR by a novel best-estimate code Roozbeh Vadi a, b, *, Kamran Sepanloo b a b

Department of Energy Engineering and Physics, Tehran Polytechnic University, 424 Hafez Ave., P.O. Box: 15875-4413, Tehran, Iran Nuclear Science and Technology Research Institute, End of North Karegar Ave., P.O. Box: 14359-836, Tehran, Iran

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 April 2015 Received in revised form 28 July 2015 Accepted 24 October 2015

In this study, a new CFD-based thermo-hydraulic code is introduced and tested on some standard verification problems. This is a time-dependent and three-dimensional code that has no restriction with regard to the details of geometry of solution domain. The CFD-code is also capable of modeling both single-phase and two-phase flow. This code is specifically devised for combining with neutronic codes by a novel scheme. For this study, the Point-Kinetics equations with six groups of precursors are chosen as neutronic part. Methods of solving this set of equations are surveyed with the criterion of best compatibility with the CFD code. The resulting combined code is applied to investigate a loss of coolant accident in a typical MTR. This accident involves reactivity insertion, flow reversal and natural convection simulation. All models required for these phenomena are also developed and implemented into the combined code. The results show considerable improvements in comparison with the conventional conservative approach. © 2015 Elsevier Ltd. All rights reserved.

Keywords: LOCA CFD Point-kinetics equations MTR Best-estimate approach

1. Introduction Nuclear energy office of the US department of energy established an Advanced Modeling and Simulation Office (AMSO) in 2008 (AMSO, 2014). AMSO has initiated two parallel programs (NEAMS, 2014; CASL, 2014) to develop a new generation of integrated modeling and safety software with the capability of “predictive simulation” of key nuclear energy problems. The concept of predictive simulation is to remove or minimize the use of empirical correlations (mathematical based methods) to define simulation tools by incorporating the first-principles effects of the governing physics and transcending them to the traditional “engineering” scale (Banfield, 2013). After this program, similar efforts have been started in the European Union (F-BRIDGE, 2014), Russia (3D flow, 2014), South Korea (SMART project, 2014) and Japan (IVNET, 2014; ASRC, 2014). This study is a part of an academic program that has been started with the same motive. Considering the purpose, developing a high-fidelity CFD code is of crucial importance. Therefore, all of the mentioned programs (NEAMS, 2014; CASL, 2014; Banfield,

* Corresponding author. Nuclear Science and Technology Research Institute, End of North Karegar Ave., P.O. Box: 14359-836, Tehran, Iran. E-mail addresses: [email protected] (R. Vadi), [email protected] (K. Sepanloo). http://dx.doi.org/10.1016/j.pnucene.2015.10.013 0149-1970/© 2015 Elsevier Ltd. All rights reserved.

2013; F-BRIDGE, 2014; 3D flow, 2014; SMART project, 2014; IVNET, 2014; ASRC, 2014), use CFD approach for thermo-hydraulic simulation, although the applied numerical methods are different. Section 2 gives a brief review of the methods and treatments applied to the CFD code introduced in this study. Algorithm of this code is handpicked in a way that it can be used for a novel scheme which solves thermo-hydraulic and neutronic equations simultaneously. For this study, the Point-Kinetics equations with six groups of precursors are chosen for neutronic analysis. Criterion for selecting method of solving these equations is the capability of combining it with the CFD code in an efficient manner. Four premier methods are surveyed in Section 3 to assess which one fulfills the criterion better. Then algorithm of combining the selected method and the CFD code is presented in Section 4. The resulting combined code is used to investigate a loss of coolant accident in Tehran Research Reactor (TRR) that is a typical MTR. Section 5, gives a brief description of this reactor and Section 6, explains the accident scenario. This accident involves reactivity insertion, flow reversal and natural circulation. Additional models required for simulating these phenomena are presented in Section 7.2.2 to Section 7.2.4. Two models are developed to simulate the accident. A simple model that uses conservative approach (Section 7.1) and a comprehensive model that uses best-estimate approach (Section 7.2). Since the code has no restriction as to details of the geometry, it is compatible with both approaches. The simple model is

142

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Nomenclature

beff ¼

P6

i¼1 bi

f

Effective delayed neutron fraction Reactivity coefficient Time at which flow reversal occurs after scram Power level at which flow reversal occurs Time it takes for flow velocity to reach stagnation (zero value) from initial velocity Time it takes for flow velocity to reach equilibrium (VE) from stagnation (zero value) Equilibrium velocity of pure (steady-state) natural convection mode An arbitrary dependent variable

Subscripts F M D o n i c

Fuel Moderator Density Old value of variable New value of variable Number of delayed neutron group (from 1 to 6) Cell

a Abbreviations LOCA Loss of Coolant Accident MTR Material Test Reactor CFD Computational Fluid Dynamics AMSO Advanced Modeling and Simulation Office FA Fuel Assembly PPF Power Peaking Factor Variables n(t) Power level (neutron concentration) r Reactivity bi Effective delayed neutron fraction of delayed group i Ci Precursor concentration of delayed group i li Decay constant of delayed group i L Prompt neutron life-time T Temperature D Density H Active height (length) of fuel plates

presented to compare with the results of former studies on this accident (Aghaie et al., 2012; AEOI, 2001) that used the same approach because direct comparison of the results of a conservative study and a best-estimate one is either wrong or meaningless due to the fact that the basic assumptions are completely different. On the basis of the step-by-step verifications of Section 2.2, Section 3, Section 7.1 and Section 7.2.1, the results of the comprehensive model are presented in Section 7.2.5. The results show that the comprehensive model provides more accurate prediction of the accident, thereby refining the countermeasures that should be taken. 2. Introducing the new thermo-hydraulic code

tR PR tz tE VE

Bonhus, 1994), the accuracy of this method is comparable to the node-based ones while it retains considerable portion of the computational cost reduction. In this method the solution is assumed to vary linearly. In Fig. 1, the change in cell values between cell co and ci along the vector Dri from the centroid of cell co to cell ci, can be expressed as (Anderson and Bonhus, 1994):


 ðVfÞco $Dri ¼ fci  fco

(1)

If similar equations for each cell surrounding the cell co are written, the following system (written in compact form for convenience) will be obtained (Anderson and Bonhus, 1994):

½JðVfÞco ¼ Df

(2)

2.1. Methodology Table 1 gives a brief review of numerical method and special treatments implemented into the new CFD code along with the references in which they are explained. Many modifications are applied to the original methods mentioned in the table. These modifications are mainly focused on two points; reducing computational resources required and applying these methods to an unstructured and hybrid mesh grid. The most noteworthy modifications are:
Improving method of applying the least-square scheme for gradients evaluation. The finite-volume method employed in this study is essentially a cell-based technique. It means that discrete values of variables are preserved at the cell centers in each numerical step and only they are available for the calculations in the next one. It stands in contrast to the node-based methods such as finite-element method where the values are preserved at every nodes of the cell. As a result, computational cost of finite-volume method is intrinsically lower than node-based ones. However the accuracy of this method is highly dependent on the technique used for evaluation of the variables gradients. In this regards, the “least squares cell-based” scheme has been touted as the one of the best available methods because according to the mathematical proof of Ref. (Anderson and

where [J] is a geometric-based coefficient matrix. The objective here is to determine the cell gradient by solving the minimization problem for the system of the non-square coefficient matrix in a least-squares sense. This linear-system of equation is to be solved during the solution whereas in this study, the decomposition process proposed in Ref. (Anderson and Bonhus, 1994) has been used. This decomposition yields a matrix of weights for each cell. Thus for the cell-centered scheme this means that the three components of the weights (Wxio, Wyio, Wzio) are produced for each of the faces of cell co. The key feature of current study modification is that these weights are pure functions of geometry (numerical grid). Therefore instead of solving a linear-system of equations in every numerical step, it is possible to calculate these weights prior to solution at the stage of reading the numerical grid into the code and then during the solution, the components of the gradient at the cell center are computed by following simple multiplication. n X
  j
Wio $ fci  fco Vfj c ¼ o

j ¼ x; y; z

(3)

i¼1

In this study, to save even more computational resources, these weights are not preserved in dynamic memory but they are written to a binary file prior to the solution and during the solution, they are extracted from the file at the appropriate time.

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

143

Table 1 Brief review of numerical method and special treatments implemented into the new CFD code. Title

Applied method

Numerical method

Finite volume (Leveque, 2004)

Solver approach

Single-phase flow Two-phase flow

Turbulence model

Simple Advanced

Linearization shape function Method of evaluation of dependant variables on faces of a numerical cell Pressure interpolation scheme (including the momentum equations)

General Wall modeling General Wall modeling

Standard k-ε (Launder and Spalding, 1972) Standard wall-function (Launder and Spalding, 1972) Realizable k-ε (Shih et al., 1995) Enhanced wall-treatment (Chen and Patel, 1988)

Second-order linear function (Leveque, 2004) Convection terms (except for pressure term in the momentum equations) Diffusion terms

Central-difference (Leonard, 1991)

Natural-convection dominated problems Other problems

Body force weighted model (Patankar, 1980) Central-difference (Leonard, 1991)

Gradients evaluation method

Second-order upwind (Barth and Jespersen, 1989)

Least squares cell-based (Anderson and Bonhus, 1994) Differential type (Wang, 2000)

Gradient (slope) limiter Pressure-evelocity coupling scheme (pressure-ecorrection equation)

Pressure-based (Segregated) (Chorin, 1968) Pressure-based (coupled) (Chorin, 1968)

Steady-state problems Time-dependant problems

SIMPLE (Rhie and Chow, 1983) PISO (Issa, 1986)

Under-relaxation method

Variables (explicit) relaxation scheme (Leveque, 2004)

Solution accelerator

Algebraic Multi-Grid (AMG) with V-cycle recursive procedure (Hutchinson and Raithby, 1986)

solution method of Linear Algebraic Equations (LAE) set

LU decomposition with partial implicit pivoting (Press et al., 1992)

General multiphase model

Drift flux (Mixture) model (Ferzieger and Peric, 1996) Two-flow (Eulerian) model (Ferzieger and Peric, 1996)

Wall boiling models in the context of the Eulerian multiphase model

Sub-cooled boiling regime Departed nucleate boiling regime


Modifying the PISO scheme for pressureevelocity coupling. The pressure-based solver used in this study, employs an algorithm which belongs to a general class of methods called the projection method (Chorin, 1968). In the projection method, the “pressure correction equation” is used instead of the continuity equation which is obtained by manipulating continuity and momentum equations (Chorin, 1968). As noted in Table 1, in this study, two algorithms are employed as pressureevelocity coupling schemes for the pressure correction equation; the “SIMPLE” algorithm (Rhie and Chow, 1983) for the steady-state problems and the “PISO” algorithm (Issa, 1986) for the transient ones. Both algorithms use a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field. This

Fig. 1. Control volumes used to illustrate the method of evaluating variables gradients.

RPI nucleate boiling model (Kurul and Podowski, 1991) Extended RPI model (Lavieville et al., 2005)

relation is of the first-order for SIMPLE algorithm and it is of the second-order in the case of PISO algorithm. Although the SIMPLE algorithm works properly for most of the problems but one of its limitations is that new velocities and corresponding fluxes do not satisfy the momentum balance after the pressure correction equation is solved due to the fact that in the pressure-based solver, the pressure correction equation is solved after the momentum equation (see Fig. 2). As a result, the calculation must be repeated until the balance is satisfied. It can dramatically increase the number of iterations required for convergence, especially for transient problems. To improve the efficiency of this calculation, the PISO algorithm performs two additional corrections: momentum correction (Ferzieger and Peric, 1996) and skewness correction (Issa, 1986). The geometry of solution domain pertaining to a nuclear reactor is often complicated due to existence of high aspect ratio (see Section 7.2.1). This matter results in existence of some degree of skewness in the generated numerical grid. In this case, simultaneous coupling of the momentum and skewness corrections at the same pressure correction equation source may cause divergence or a lack of robustness. Therefore in this study, method of handling the momentum and skewness corrections inside the PISO algorithm is modified to apply two or more iterations of skewness correction for each separate iteration of neighbor correction. This technique allows a more accurate adjustment of the face mass flux correction according to the normal pressure correction gradient and reduces the total number of iterations required for achieving convergence.
Using geometric-based method for AMG accelerator instead of mathematic-based one (Hutchinson and Raithby, 1986).

144

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

between adjacent cells (Hutchinson and Raithby, 1986). In Algebraic Multigrid (AMG) scheme selected for this study, all of the required steps including creation of the coarse levels meshes are conducted mathematically (Hutchinson and Raithby, 1986). This feature makes the AMG scheme particularly attractive for use on unstructured meshes. In this study, all instructions of the original scheme have been followed except for the method of creating the coarse levels meshes. In this regard, due to the following reasons, coarse levels meshes are created geometrically. 1. The original scheme uses a mathematical algorithm called “grouping” in order to merge the fine level cells into the coarse level ones (Hutchinson and Raithby, 1986). If the numerical grid, besides being unstructured, is hybrid as well, it will increase number of the required calculations for each iteration in order to implement the algorithm and as a result the computational cost is high in this case. 2. The grouping algorithm attempts to collect cells into groups of fixed size, typically two or four and all efforts in this algorithm have been devoted to ensure that these groups are unique. In other words, no regards have been paid to resolution variation of numerical grid or type of the problem at hand.

Fig. 2. Algorithm of the CFD-based thermo-hydraulic code.

As noted in Table 1, LU decomposition with partial implicit pivoting is used for solving linear equations set. Although this solver like other solvers of the “point-implicit” category rapidly removes local (high-frequency) errors in the solution, global (lowfrequency) errors are reduced at a rate inversely related to the mesh size. Thus, for a large number of nodes, the solver “stalls” and the residual reduction rate becomes prohibitively low (Press et al., 1992). Multigrid techniques allow the global errors to be addressed by using a sequence of successively coarser meshes. This method is based upon the principle that global errors existing on a fine mesh can be represented on a coarse mesh where it again becomes accessible as the local error; because there are fewer coarse cells overall, the global corrections can be communicated more quickly

Therefore, in this study, in addition to the original numerical grid, two other set of grid are prepared beforehand which are obtained from successive coarsening of the original one. It is obvious that any kind of manipulation can be applied to these coarse level meshes according to the resolution of numerical grid in any specific part of it or type of the problem at hand. Moreover, the data pertaining to these grids are saved in a form of binary file prior to the solution and during the solution they are extracted from the file at the appropriate time and no calculation is required. Therefore, by this method, the computational cost also reduces considerably. This code is based on the time-dependent conservation equations for single-phase and two-phase flow. The form of the conservation equations used in the current study can be found in Ref. (Ferzieger and Peric, 1996). As noted above, the code is developed in order to work on unstructured and hybrid cells of two or three dimensional numerical grid, for any free geometry. The algorithm of this code (shown in Fig. 2) is handpicked in a way that it can be used for a novel scheme which solves thermo-hydraulic and neutronic equations simultaneously, in contrast to the conventional coupling scheme. Although this scheme is practically used in this study (see Section 4) but we will explain it in another paper because the real applications and advantages of this scheme can be

Fig. 3. Comparison of Nusselts number variation along the heated wall (current study and Ref. (Baughn et al., 1984)).

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

demonstrated more clearly when neutronic equations have spatial dependency which is not the case in the current study.

2.2. Verification Full process of verification and validation of the new CFD code is beyond the scopes of this study. Instead the code was tested on 56 standard verification problems from Ref. (ANSYS Inc, 2013). The results of a few of these attempts are presented in Table 2. Figs. 3e6 are referred in the table to compare the results with the references. These examples should cover physical concepts of the thermohydraulic phenomena which are relevant to the current study, such as turbulence simulation (Baughn et al., 1984), natural convection (Kuehn and Goldstein, 1978) and conjugate heat transfer (Incropera and Dewitt, 2006). Note that the maximum error is less than three percent, in all 56 cases. This should be reliable enough for most of the cases.

2. An open-integral method of order 5 from branch of explicit integral method (Burden and Faires, 2010). 3. A closed-integral method of order 4 from implicit branch which is well known as the Gear method (Kinard, 2003; Gear, 1971). 4. A semi-analytical method which is well known as Spectral or Stairway method (Kinard, 2003). These methods are tested on all standard verification problems from Ref. (Kinard, 2003) but specifications of the one that decided the best method is presented in Table 3. The results are tabulated in Table 4. Note that all methods work fairy well at time-step of 0.001 [sec] but the only method that reaches acceptable results at 0.01 [sec] is the last one. It should be mentioned that, in contrast to Runge-Kutta method, both integral methods did not diverge at this time-step but amount of the errors are off the chart. As a result the Spectral method is selected for this study. 4. Algorithm of combing thermo-hydraulic and neutronic analyses

3. Selecting method of solving point-kinetics equations Point-Kinetics equations are declared by following set of ordinary differential equations (Ott and Neuhold, 1985).

8 6 X > > > dnðtÞ ¼ rðtÞ  beff nðtÞ þ > li Ci ðtÞ < dt L i¼1 > > > > : dCi ðtÞ ¼ bi nðtÞ  li Ci ðtÞ; i ¼ 1; …; 6 dt L

145

(4)

There are many methods to solve this set of coupled equations (Kinard, 2003). However to choose the best method for the current study, following remarks should be considered.
It is an established fact that neutronic responses are generally faster than thermo-hydraulic ones (Ott and Neuhold, 1985). Therefore neutronic equations such as Eq. (4) generally require shorter time-step to grasp correct variation of the variables.
In this study, thermo-hydraulic equations and the method of solving them (Section 2) are more complicated and require much more computational resources than the neutronic ones. For these reasons, the criterion of selecting the best method can be stated in one sentence. “It should work with longest possible time-step while it keeps accuracy and stability”. Following prominent methods are surveyed to assess which one is capable of fulfilling this criterion better. 1. Runge-Kutta method of order 4 from branch of differential methods (Burden and Faires, 2010).

Algorithm of combining Spectral method (for solving PointKinetics equations) and the CFD code is represented in Fig. 7. The mathematical stability of the iterative part of this algorithm was proven in Ref. (Chorin, 1968). The thermo-hydraulic feedbacks are included by following equation whereby the reactivity required for solving Eq. (4), is evaluated at every time-step.

rðtÞ ¼ r0 ðtÞ þ aTF ðTF ðtÞ  TFo Þ þ aTM ðTM ðtÞ  TMo Þ   DM ðtÞ  DMo  aMD DMo

(5)

where ro(t), denotes the reactivity induced by the movement of control rod(s) and D and T refer to average density and temperature respectively that are weighted by the volume of the relevant region (fuel or moderator). Subscripts F and M refer to fuel and moderator respectively and subscript o refer to initial (steady-state) value of the relevant variable and aTF, aTM, aMD are reactivity coefficients of fuel temperature, moderator temperature and moderator density respectively. 5. Brief description of Tehran Research Reactor (TRR) TRR is a typical pool-type MTR. The relevant specifications of this reactor are presented in Table 5. A schematic diagram of TRR is illustrated in Fig. 8. In the normal operational condition, coolant flows downward through the core by the primary cooling circuit pump. However when shutdown or scram commences, the direction of flow reverses after a specific amount of time (hence it flows

Table 2 Results of testing the new CFD-code on three standard verification problems. Test case no.

Description

Reference

1

Turbulent heat transfer in a pipe expansion

(Baughn et al., 1984)

2

Natural convection in a concentric annulus

(Kuehn and Goldstein, 1978)

3

Conjugate heat transfer in a composite solid block

(Incropera and Dewitt, 2006)

4

Boiling in a pipe - reaching critical heat flux and simulation of post-dryout conditions

(Hoyer, 1998)

Ref. result(s)

Current study result(s)

Error (maximum value) [%]

Fig. 3(a)

Fig. 3(b)

1.79

Fig. 4(a)

Fig. 5(a)

Fig. 4(b)

Fig. 5(b)

1.07

1.91

The cooled wall 378 [ K]

The adiabatic wall 413 [ K]

The cooled wall 378.15 [ K]

The adiabatic wall 413.12 [ K]

The cooled wall 0.04

The adiabatic wall 0.03

Fig. 6(a)

Fig. 6(b)

2.09

146

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 4. Comparison of temperature distribution on the bottom wall of symmetry (current study and Ref. (Kuehn and Goldstein, 1978)).

Fig. 5. Comparison of temperature distribution on the top wall of symmetry (current study and Ref. (Kuehn and Goldstein, 1978)).

Fig. 6. Comparison of temperature distribution on the tube wall (current study and Ref. (Hoyer, 1998)).

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

147

Table 3 Conditions of a standard test case with ramp reactivity insertion from Ref. (Kinard, 2003). Quantity

Value

b1 b2 b3 b4 b5 b6 l1 l2 l3 l4 l5 l6

0.000266 0.001491 0.001316 0.002849 0.000896 0.000182 0.0127 s1 0.0317 s1 0.115 s1 0.311 s1 1.4 s1 3.87 s1 0.00002 s Ramp, 0.1 $/s

L

r

upward) and the residual heat would be removed by pure natural convection. 6. The accident scenario Safety analysis report of TRR (AEOI, 2001), categorized the loss of coolant accidents as follows: 1. Breaks below the grid level that covers all the piping and instruments of the primary cooling system outside the reactor pool (see Fig. 8). 2. Breaks above the grid level, involving the irradiation beam-tube facilities (AEOI, 2001). The first category is of the interest to this study. If such an accident should happen, obviously, it would be intended to isolate the pool as soon as possible by closing the control valves, otherwise the pool would be losing water and water level would be dropping consistently. Furthermore, continuous leakage of water to the damaged area could escalate the cause of initiating event or even results in secondary accident(s) (AEOI, 2001). However if the isolation of the pool happens immediately after the scram, it might put the core in danger because in that case, the direction of flow reverses immediately and the consequent jump in temperature due to instantaneous stagnation of flow might violate the safety margins of temperature. For these reasons, a compromise should be established. From one hand, one intends to isolate the pool as soon as possible, from other hand, fuel plates should have enough time to cool down. Therefore the main target here is to find the efficient time interval after the accident for isolating the pool. 7. Simulation, results and discussion 7.1. Simple model of the accident (conservative approach) This accident had been investigated by a conservative approach in Refs. (Aghaie et al., 2012; AEOI, 2001). In this section, all

Fig. 7. Algorithm of combing thermo-hydraulic and neutronic analyses.

simplifying assumptions of these references are applied to the combined-code of this study as much as possible. The results are compared with their counterparts from Ref. (AEOI, 2001) so as to verify current study and also as a basis to assess effects of revoking these assumptions in the following section. In order to comply with these assumptions, domain shown in Fig. 9 is chosen. The upper and lower parts are added in order to include cooling effects of flow on the plate from bottom and top and also to apply inlet velocity in a manner that is consistent with following sections of this study. Boundary conditions and thermo-physical properties that are applied to this model are presented in Table 6. List of these assumptions and conditions are as follows (Aghaie et al., 2012; AEOI, 2001):

Table 4 Comparison of different methods for solving point-kinetic equations with the conditions of Table 3. Runge-Kutta method of order 4 Time [sec] 2 4 6 8 9

Dt ¼ 0.001 s 1.337871 2.228519 5.59001 42.77891 478.4113

Dt ¼ 0.01 s Diverged Diverged Diverged Diverged Diverged

Open integral method of order 5

Dt ¼ 0.001 s 1.3382 2.227591 5.57899 42.786 478.524

Dt ¼ 0.01 s 1.8741 4.3716 12.9001 101.054 1351.231

Gear method of order 4

Spectral method

Dt ¼ 0.001 s

Dt ¼ 0.01 s

Dt ¼ 0.001 s

Dt ¼ 0.01 s

1.33819 2.227593 5.57912 42.786 478.521

1.8699 4.3602 12.8813 100.742 1347.092

1.337932 2.228331 5.581545 42.78140 487.4531

1.338046 2.228029 5.581859 42.78544 487.5113

Target results (Kinard, 2003) 1.3379 2.2283 5.5815 42.781 487.45

148

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161 Table 5 Specifications and main operating conditions of TRR (AEOI, 2001). General Thermal power

5 MW

Core dimensions (first operating core)

40.5  38.54  89.7 cm

Grid plate locations for fuel assemblies

96

Each location cross-sectional dimensions

7.71  8.1 cm

Reflector material

Graphite/Light water

Control elements

4 Ag-In-Cd shim safety rods 1 stainless steel regulating rod

Thermo-hydraulics Coolant material

Light water

Cooling method

Operational condition Full shutdown condition

Operating pressure

0.171 Mpa

Primary coolant flow rate

500 m3/h (2200 gpm)

Coolant inlet temperature (full power)

37.8  C (100  F)

Coolant outlet temperature (full power-first core)

47.1  C (115.7  F)

Forced flow primary loop Natural circulation

Fuel assembly Fuel

Plate type (MTR), LEU 20%, Al Clad

Fuel assembly dimensions (Max.)

7.6  8.01  87.8 cm

No. of fuel plates in standard fuel assembly

19

No. of fuel plates in control fuel assembly

14

Coolant passage thickness

0.27 cm

Fuel meat material

U3O8-Al

Fuel plate dimensions

0.15  7.05  65.5 cm

Fuel meat dimensions

0.07  6.00  61.5 cm

1. Only a single channel that consists of a plate and a passage between two plates are chosen as simulation domain. 2. The model of simulation and the solver are basically onedimensional in axial direction and in this direction only 10 numerical grids are considered. The numbers after word “T” in Fig. 10 refer to the location of these grids from top to bottom in an increasing manner. 3. Shape of power distribution is assumed to be symmetric cosine function in axial direction and constant in the other directions. 4. Shutdown process is simulated by instantaneous drop of power level from 5 MW to 200 kW and this level, remains constant afterwards. 5. Velocity reversal process is also assumed to occur instantaneously at 10 s after shutdown by exchanging location of flow inlet and outlet and enforcing velocity of 0.02 [m/s] to the new inlet as a constant value. In order to compensate any possible consequences of these assumptions, a large conservative power peaking factor is selected (Aghaie et al., 2012; AEOI, 2001). But precise value of this factor is not directly stated in neither of these references. For the results presented in Fig. 10, PPF is equal to 3.0. 7.2. Comprehensive model of the accident (best-estimate approach) All parts of this section are prerequisites that are required for removing five simplifying assumptions and conditions, noted in Section 7.1 and therefore for attaining an accurate simulation of the accident, based on best-estimate approach. Four cases are considered for study in this section. These cases

are defined on the basis of time intervals between commencing scram and isolating the reactor pool (see Section 6). The intervals are 10, 1, 0.25 and 0.0 [sec] for case 1 to 4 respectively. Other specifications of these cases are discussed in Section 7.2.3. 7.2.1. Geometry of comprehensive model and verification of this model Geometry of the simple model (Fig. 9) only includes single fuel plate and its cooling passage. Domain of solution is expanded here to include whole fuel assembly. Scattered view of essential consisting parts and a view of the final assembled solution domain are shown in Fig. 11. The central part consists of fuel meat and clad plates as solid regions and passages between these plates and upper and lower plenum as fluid regions. The side part is a flow passage for cooling outer side of the side plates. As Fig. 11 shows by four thin legs, this region continues to the end of narrow passages that have been implemented in the grid plate for passing this cooling flow (AEOI, 2001). The top and bottom part are actually parts of the pool and plenum medium respectively. These parts are included to simulate effects of flow mixing and separation between central and side part and also for answering this question that how flow-rate divides between these two passages. To verify this model and also to assess importance of including all of these geometrical details, normal operational conditions of the first core configuration (Fig. 12) are selected to test on this model. Boundary conditions and material properties that are used for this simulation are same as the initial values mentioned in Table 6. Results of temperature rise and pressure drop presented in Table 7 are sorted in a descending manner, from left to right, based on the numerousness of the error imposed from different

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

149

Fig. 8. Schematic diagram of Tehran Research Reactor (AEOI, 2001).

geometrical and physical assumptions. The complete model is the result of discarding all of these assumptions. The results are compared with the counterpart data from the SAR of TRR (AEOI, 2001) which are presented in the last column of Table 7. These results show that the largest amount of error is obtained when effects of flow mixing and separation are ignored by omitting side, top and bottom parts (Fig. 11). From total flow rate of 8.521 [kg/s] that enters domain, share of side part is 0.0107 [kg/s] which means 0.126%. This amount is 27% lower than 0.0172 [kg/s] that is obtained from fraction of cross sectional area of these two passages. Therefore pressure loss and temperature rise in the side part are always higher than the central part. For this reason, passage of the side part is the hot channel. This matter can be observed more clearly in Fig. 13 in which counters of temperature distribution on some virtual cross-sections, vertical to the axis of the fuel assembly, are shown. Second recognized source of error is applying the simple turbulence model instead of the advanced one (see Table 1). Third source of error roots in disregarding 38 very narrow passages that exist at the both end of each plate. Cross sectional dimension of these passages is only 0.75  1.5 [mm]. Results of fifth column of Table 7 prove importance of including temperature-dependent material properties.

quantity (condition 4). Whereas it is an established fact that after a full shutdown process, due to the large negative reactivity worth of control elements, power level should drop continuously and also any changes in temperature of fuel and in temperature and density of moderator would cause consequent reverse reactivity feedbacks (Ott and Neuhold, 1985). These effects can be included by solving Point-Kinetic equations (Eq. (4)). This set of equations has been extensively validated and used for evaluating power level of MTRs such as TRR during transients (Kazeminejad, 2006, 2008; Salama and El-Morshedy, 2011). By using these equations, one assumes that the general shape of power distribution is fixed over a transient and only magnitude of this shape that is the main output of point-kinetics equations could be variant in every time-step (Ott and Neuhold, 1985). Therefore, precise estimation of three following prerequisites is crucial to use point-kinetics equations in an accurate manner.

7.2.2. Requirements of using point-kinetics equations for power excursion In the simple model (Section 7.1), power variations are simulated by instantaneous drop of power level at the shutdown moment and other than that, power level is an ever-constant

Configuration of the first operating core of TRR (Fig. 12) is selected for this study because it is the most compact one (AEOI, 2001), hence it provides the highest average power per fuel element. The CITVAP code from MTR_PC package (Villarino and

1. Estimating spacial shape of power distribution that depends on the core configuration. 2. Evaluating kinetic parameters that are denoted by b and L in Eq. (4). 3. Calculating reactivity at every time-step, using Eq. (5) that requires evaluation of reactivity coefficients.

150

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

procedure instructed in Ref. (IAEAIAEA-TECDOC-233, 1978). The results are presented in Fig. 15 and Table 8. It should be mentioned that by solving the Point-Kinetics equations in the form Eq. (4) (with six groups of precursors) all the instantaneous and delayed energy productions are included due to the fact that variations of the concentrations of all isotopes are accounted for, in the calculations (Ott and Neuhold, 1985). This fact is not only true about normal operational condition but also includes the shutdown mode. Only the none-fission reactions due to the excess neutrons plus b- and g-decay energy of (n,g) products are left out from these calculations (Ott and Neuhold, 1985). According to general analysis of Ref. (IAEAIAEA-TECDOC-233, 1978) and the specific study of Ref. (AEOI, 2001) for TRR, these reactions comprise a fraction between 2 and 3.5 percent of the total energy at every time-step. The upper limit of this range is added to the power that is calculated by the Point-Kinetics equations at every time-step after commencing scram. According to Ref. (Ott and Neuhold, 1985), when the accident starts from normal operational (critical) condition, the ratios of initial concentration of delayed neutron precursors to the initial power are calculated as follows:

C0i b ¼ i n0 Lli

i ¼ 1; …; 6

(6)

After first time-step, initial concentration of precursors is known from the previous one. 7.2.3. Flow reversal model The equilibrium velocity of natural circulation in steady-state condition can be evaluated by following equation (Todreas and Kazimi, 2001).

DPloss ¼ DPBuoyancy

(7)

In the case of a pool-type reactor, each side can be estimated as follows.

DPloss ðVE Þ ¼ DPfriction loss þ DPform loss þ DPflow acceleration Fig. 9. Geometry of solution domain for the simple model of the accident (symbolic).

Carlos, 1993) is applied to estimate spatial shape of power for this configuration. The code WIMSD-5B (NEA, 2003) is also employed to generate macroscopic cross-sections required for these calculations. The simulation is performed for two positions of control rods; when they are completely out of the core and when they are halfway in the core. Results are shown in Fig. 14. Note that in both cases, hot fuel assembly locates at central position but the power level is higher for the first condition. Therefore power shape of this case, Fig. 14 (a2), is selected. It means that, for the comprehensive simulation, the total power produced in the hot fuel assembly is divided equally among the 19 consisting fuel plates using this shape of power. The CITVAP code is also used for evaluating the kinetic parameters and reactivity coefficients according to the standard

(8)

I DPBuoyancy ðVE Þ ¼

Dgvl loop

Z

¼


b a

Z Dgdz 

c b

Z Dgdz 

 ¼ DPool  Da ðVE Þ gH

d c

Z Dgdz 

a

Dgdz d

(9) The convention of direction for the loop integration is presented in Fig. 16. If it is possible to state average density (Da) and pressure drop in Eqs. (8) and (9), as a function of equilibrium velocity (VE), this velocity can be easily calculated by solving Eq. (7). For this purpose, in this study, results of detailed simulation of a FA in the natural convection and steady state mode are used. These data are tabulated in Table 9 for the four test cases of this study. Then a free-

Table 6 Boundary conditions and material properties applied to the simple model. Initial inlet temperature Initial (steady-state) inlet velocity Operating Pressure Initial heat source Coolant material properties (r, m, k, Cp) Fuel & clad material properties (r, k, C)

37.8  C 1.1705 [m/s] 0.171 MPa Cosine function of axial direction with the average value of configuration of the first operating core (see Fig. 12) Temperature-dependant polynomials that are fitted to the data extracted from steam table at the operating pressure. From reference (IAEA IAEA-TECDOC-949, 1997). (fuel's thermal conductivity should be defined as a function of Temperature)

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

151

Fig. 10. Comparison of average cross-sectional temperature during the accident using the simple model (current study (CFD) and SAR (AEOI, 2001) - see Section7.1 for definition of T#).

152

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 11. Geometry of solution domain for the comprehensive model of the accident.

Fig. 12. Configuration of the first operating core of Tehran Research Reactor.

order polynomial curve is fitted to each one of them, as Fig. 17 represents an example for the case 1. If these polynomial functions are used in Eq. (7), the unknown velocity will be obtained. This treatment is extended here to achieve a method for evaluating time-dependant profile of velocity reversal. This is done by performing the loop integral on time-dependant and simplified

form of coolant momentum equation as follows.

dðDVÞ vP ¼ Dg  losses dt vl

(10)

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

153

Table 7 Pressure drop and temperature rise in the first operating core of TRR in comparison with numerical results by enforcing different geometrical and physical assumptions (the complete model is the result of discarding all of the assumptions - each assumption is enforced separately). Applied assumptions and simplifications

Discarding the side, bottom and top parts from Fig. 11

Using the simple turbulence model (see Table 1)

Discarding passages at ends of fuel plates from Fig. 11

Using constant thermophysical properties

Complete model

Target results (AEOI, 2001)

Dp [kPa] DT [ C]

12.033 10.524

18.815 9.4457

13.778 10.501

14.572 9.1929

15.870 9.2796

z16 9.3

Fig. 13. Counters of temperature [ K] distribution on some virtual cross-sections, vertical to the axis of fuel assembly.

Fig. 14. Comparison of power distribution in the first operating core of TRR based on different position of control rods; (a) when they are completely out of the core, (b) when they are halfway in the core.

H /H


 d
V DPool  Da ðV Þ ¼ DPool  Da ðV Þ  gH  DPlosses ðV Þ dt (11)

0H


 d ðVDa ðVÞÞ ¼ DPool  Da ðVÞ  gH  DPlosses ðVÞ dt

(12)

By putting the same polynomial relations that were obtained from CFD simulation into differential Eq. (12) and solving it, the

154

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 15. Variation of reactivity coefficients against the relevant variable (the fitted correlations are to be substituted in Eq. (5)).

Table 8 Kinetic parameters of the first operating core of TRR (Calculated using the CITVAP code (Villarino and Carlos, 1993)). Quantity

b1

b2

b3

b4

b5

b6

L [ms]

Value

0.000226

0.001772

0.001590

0.003211

0.000974

0.000345

45.1

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

155

There are two approaches to perform last step: 1. The Boussinesq approximation. In this approach, density is assumed to be a constant value (D0) in all conservation equations, except for the buoyancy term in the momentum equation. In this term, density is defined as a linear function of temperature.

DgxD0 bðT  T0 Þg

(13)

The constant density, D0 is evaluated at a presupposed average temperature, T0. Slope of the linear function is thermal expansion coefficient (b) that is also an approximate value due to unknown range of temperature variation. Fig. 16. The convention of direction for the loop integration of Eq. (9).

results presented in Fig. 18 and Table 10 have been obtained for the four considered cases of this study. Method of solving, initial condition and applied time-step are presented in Table 11. Two noteworthy facts about these results are as follows: 1. Fig. 18 shows that for all four cases of this study, velocity profiles converge to the equilibrium velocities (last column of Table 10) that are calculated by Eq. (7) with the same pressure and density correlations. 2. By increasing power and the consequent increases of established buoyancy force, time that is required to overcome initial coolant momentum and reach stagnation (zero value of velocity), decreases. But due to the same reason after stagnation (in negative mark range of velocity), it takes more time to reach equilibrium velocity for higher initial power level than a lower one. So the gap between the total times reduces.

7.2.4. Natural convection model When heat is added to a fluid, temperature gradient causes density difference in that fluid. The buoyancy force formed by the density difference can induce a flow in the fluid. This buoyancyinduced flow is termed “natural convection”. Following steps should be taken to include effects of natural convection in a CFD simulation (Launder and Spalding, 1972) or any similar approach that solves general conservation equations on a free geometry.
Solving the energy conservation equation.
Including the buoyancy term (Dg) in the momentum conservation equation.
Defining density (D) as a function of temperature in a proper manner.

2. If a thermohydraulic table is available for the fluid, it is possible to extract variation of density against temperature directly from the table at the operating pressure. Obviously, this approach is expected to have a better accuracy. The method used in this study is to fit a polynomial function through these data and then use the resulting correlation instead of density in all conservation equations.

7.2.5. Results of comprehensive model and discussion All requirement of the comprehensive model of the accident and the scenarios selected as case studies are explained in the preceding sections. In the current section, results of simulating these cases with the combined-code of this study are presented. Similar to the simple model (Section 7.1, Fig. 10) and Ref. (Aghaie et al., 2012; AEOI, 2001), it is assumed that reactor operates normally for 30 s at full power and then the accident occurs. Figs. 19 and 20 shows fuel, clad and coolant temperature during accident, averaged over surfaces at fuel middle and top position, respectively. It should be confirmed that this is a three-dimensional simulation on the geometry depicted in Fig. 11. Therefore any point on the diagrams represents an averaged value over a crosssectional surface that virtually cuts through all relevant (fuel, clad and coolant) passages at a specific location such as middle of the fuel plate. Note that only coolant temperature at top of fuel exceeded its maximum initial value during the accident (considering whole passage). Also among all temperature jumps, only in the case 4, where it is assumed that scram and flow reversal occur simultaneously, a very mild safety margin is violated. This margin (noted in Fig. 20 (c) by 71.85  C line) belongs to 10% lower range of absolute saturation temperature at the operating pressure of reactor (AEOI, 2001). Note that 4.33  C was subtracted from the evaluated temperature (90% of the absolute saturation temperature) as the upper margin of the sensitivity analysis (AEOI, 2001). In this range, a certain amount of sub-cooled boiling could be expected. However in view of the magnitude and the duration of

Table 9 Variations of average density and pressure drop against equilibrium velocity of natural convection. VE [m/s] Case 1 Case 2 Case 3 Case 4

DP [pas] Da [kg/m3] DP [pas] Da [kg/m3] DP [pas] Da [kg/m3] DP [pas] Da [kg/m3]

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

36.466 991.854 43.994 984.143 47.363 981.372 49.105 979.955

54.021 993.177 56.521 988.261 58.372 986.577 59.305 985.698

75.242 993.889 73.689 990.315 74.312 989.077 74.924 988.438

100.083 994.321 96.253 991.583 95.754 990.640 95.824 990.153

128.415 994.603 122.954 992.391 121.942 991.637 121.652 991.251

159.274 994.802 154.392 992.947 152.388 992.324 150.960 992.008

194.767 994.949 189.499 993.355 187.506 992.820 186.034 992.548

235.220 995.063 228.732 993.665 226.017 993.198 224.483 992.963

156

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 17. Variations of average density and pressure drop against equilibrium velocity of natural convection mode for the case 1 (the fitted correlations are to be substituted in Eqs. (8), (9) and (12)).

Fig. 18. Profile of velocity reversal for the four cases of the comprehensive model.

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161 Table 10 Data of flow reversal process for the four cases of the comprehensive model (for definition of the variables see nomenclature). Case no.

tR [sec]

PR [KW]

tz [sec]

tE [sec]

VE [m/s]

1 2 3 4

10 1 0.25 0

176.168 453.032 542.013 586.562

1.48 1.21 1.10 1.03

0.580 0.670 0.695 0.735

0.0161 0.0267 0.0293 0.0310

Table 11 The method and data required for solving differential Eq. (12). Method of solving Initial condition Applied time-step

Runge-Kutta method of order 4 (Incropera and Dewitt, 2006). Velocity of 1.1705 [m/s] at time zero 0.005 [Sec]

157

coolant temperature in the range, no considerable amount of vapor is expected and the detailed results of simulation with the single phase and the two-phase model (Fig. 21) show that even in the worst-case scenario (case 4) the temperature of coolant never reaches the saturation value corresponding to the operating pressure of system, therefore as expected the void fraction is equal to zero in the whole domain and for all of the cases. Following remarks should be noted with regard to these results;
According to SAR of TRR (chap. 8 and 9) (AEOI, 2001), 0.25 s of electrical delay should be included in any calculations that involves signal transfer. It means that even if it is attempted to isolate the pool immediately after scram (hence the conditions

Fig. 19. Average cross-sectional temperature during the accident at the middle of fuel plate for the four cases of the comprehensive model.

158

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 20. Average cross-sectional temperature during the accident at the top of fuel plate for the four cases of the comprehensive model.

of case 4), it will happen 0.25 s later (hence the conditions of case 3).
Considering the type of control valves that isolate the pool (see Fig. 8) and pressure at which they should operate (AEOI, 2001), 0.75 s of mechanical delay has to be added to the electrical delay mentioned (Spellman and Drinan, 2001). Therefore even if it is attempted to reach the conditions of case 4, it will end up at the conditions of case 2. As a result, in practice, no safety margin will be violated.

However, true advantage of the code introduced in this study, lies in three-dimensional and point by point simulation. Fig. 21 shows a temporal sequence of temperature distribution contours of coolant on some virtual sections parallel to the fuel plate in the hot channel for the case 4. Comparing these results with Figs. 19 and 20 show that temperature in the thermal boundary layer of the hot channel is at least 5  C higher than the counterpart average value. Therefore by replacing average values with extensive, point

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

159

Fig. 21. Temporal sequence of temperature distribution contours of coolant on three virtual sections parallel to the fuel plate in the hot channel for the case 4 (X refers to distance from the fuel plate).

by point results, an improved accuracy and prediction and also better physical understanding should be expected. Fig. 22 demonstrates variation of power level during the accident for the four cases of this study that are the main output of solving point-kinetic equations. This diagram is started from a small fraction of a second after 30 s to avoid the almost instantaneous drop of power from 5 MW to less than 450 kW. Otherwise the main part of the diagram (after 30 s) would be indiscernible.

Large negative reactivity induced by control rods fall (10beff) surpassed relatively small reactivity contributions due to fuel temperature and coolant density and temperature, which are positive or negative depending on the different phases of the transient. For this reason, differences among the power levels are negligible. However, the order of positioning of diagrams is consistent with the results of Figs. 19 and 20. For example, in the specific time period that is magnified in Fig. 22, the order of

160

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161

Fig. 22. Variation of power level during the accident for the four cases of the comprehensive model (this diagram is started from a small fraction of a second after 30 s to avoid almost instantaneous drop of power from 5 MW to less than 450 kW).

positioning of the four cases is compatible with thermo-hydraulic feedbacks that are expected from Figs. 19 and 20, in the same specific interval. In the end, it should be mentioned that, in practice, the actual time interval between the scram and isolation of pool is instructed to be 30 [sec] (AEOI, 2001). It is three times of the interval that was evaluated to be safe by the over-conservative approach of SAR (see Section 7.1). This is a typical case of conservative design. According to the results of this study, the interval can be reduced to a fraction of a second safely and thereby avoiding large amount of coolant loss and potential secondary accidents because of this leakage to the damaged area of the primary cooling circuit (see Section 6). 8. Conclusion A novel CFD-code was introduced in this study. It is a timedependent and three-dimensional code that does not have any limitation on the details of geometry. The code is also capable of modeling both single-phase and two-phase flow. These virtues make it possible to use the code for thermo-hydraulic simulation of different types of reactors. The algorithm of this code is chosen in order to combine it with neutronic codes efficiently. Point-kinetics equations with six groups of precursors were selected as the neutronic part for this study and four well-known methods of solving them were surveyed to choose the best method for combining with the CFD-code. The combined-code was used to simulate a LOCA in a typical MTR. This accident involves some sophisticated physical phenomena. All required model were also devised and implemented into the combined-code, including a new semi-analytical model for flow reversal. The results showed that by using the new code, quality of simulation and prediction were improved, not only because of enhanced accuracy but also because of providing extensive, point by point, three-dimensional results. References 3D flow, ERCOSAM-SAMRA and RSQUD projects. Nuclear Safety Institute of the Russian Academy of Science (IBRAE). [Online] Russian Federation [Cited: December 17, 2014.] http://www.ibrae.ac.ru. AEOI, 2001. “Safety Analysis Report for the Tehran Research Reactor (LEU)”. Tehran, Iran. Aghaie, M., Zolfaghari, A., Minuchehr, A., Norouzi, A., 2012. Transient analysis of break below the grid in Tehran research reactor using the newly enhanced COBRA-EN code. Ann. Nucl. Energy vol. 49, 1e11. AMSO. Advanced Modeling and Simulation Office. [Online] the USA DOE/NE [Cited:

December 17, 2014.] http://energy.gov/ne/nuclear-reactor-technologies/ advanced-modeling-simulation. Anderson, W., Bonhus, D.L., 1994. An implicit upwind algorithm for computing turbulent flows on unstructured grids. Comput. Fluids 23 (1), 1e21. ANSYS Inc, 2013. ANSYS Fluid Dynamics Verification Manual. Canonsburg, Pennsylvania, November. ASRC. Advanced Science Research Center. Sector of Nuclear Science Research (SNSR). [Online] Japan Atomic Energy Agency [Cited: December 17, 2014.] asrc.jaea.go.jp. Banfield, J.E., 2013. Semi-implicit Direct Kinetics Methodology for Deterministic, Time-dependent, and Fine-energy Neutron Transport Solutions. PhD thesis. University of Tennessee, Knoxville, Tennessee. Barth, T.J., Jespersen, D., 1989. The design and application of upwind schemes on unstructured meshes. In: Technical Report AIAA-89e0366. AIAA 27th Aerospace Sciences Meeting, Reno, Nevada. Baughn, J.W., Launder, B.E., Hoffman, M.A., Takahashi, R.K., 1984. Local heat transfer downstream of an abrupt expansion in a circular channel with constant Wall heat flux. J. Heat Transf. 106, 789e796. Burden, R.L., Faires, J.D., 2010. Numerical Analysis, ninth ed. Brooks/Cole Press, Florence, Kentucky. CASL. The Consortium for Advanced Simulation of Light Water Reactors. [Online] the USA [Cited: December 17, 2014.] http://www.casl.gov. Chen, H.C., Patel, V.C., 1988. Near-Wall turbulence models for complex flows including separation. AIAA J. 26 (6), 641e648. Chorin, A.J., 1968. Numerical solution of Navier-Stokes equations. Math. Comput. 22, 745e762. F-BRIDGE. Framework of Basic Research for Innovative Fuel Design for GEn IV systems project. [Online] the European Union. [Cited: December 17, 2014.] http://www.f-bridge.eu. Ferzieger, J.L., Peric, M., 1996. Computational Methods for Fluid Dynamics. SpringerVerlag, Heidelberg. Gear, C.W., 1971. Simultaneous numerical solution of differential-algebraic equations. IEEE Trans. Circuit Theory 18 (1), 89e95. Hoyer, N., 1998. Calculation of dryout and post-dryout heat transfer for tube geometry. Int. J. Multiph. Flow 24, 319e334. Hutchinson, B.R., Raithby, G.D., 1986. A multi-grid method based on the additive correction strategy”. Numer. Heat. Transf. 9, 511e537. IAEA, IAEA's Technical Document Vol. 949, 1997. Thermo-physical Properties of Materials for Water Cooled Reactors. International Atomic Energy Agency, Vienna, Austria. IAEA, IAEA's Technical Document Vol. 233, 1978. Research Reactor Core Conversion from the Use of Highly Enriched Uranium Fuel to the Use of Low Enriched Uranium Fuels. International Atomic Energy Agency, Vienna, Austria. Incropera, F.P., Dewitt, D.P., 2006. Fundamentals of Heat and Mass Transfer. Wiley, NY, 5th Edition, pg. 117. Issa, R.I., 1986. Solution of implicitly discretized fluid flow equations by operator splitting. J. Comput. Phys. 62, 40e65. IVNET. Innovative and Viable Nuclear Energy Technology Development project. [Online] Research Group for Thermal and Fuel Energy. Japan Atomic Energy Agency [Cited: December 17, 2014.] nsectionjaea.go.jp. Kazeminejad, H., 2006. Reactivity insertion limits in a typical pool-type research reactor cooled by natural circulation. Ann. Nucl. Energy 33, 252e261. Kazeminejad, H., 2008. Thermalehydraulic modeling of flow inversion in a research reactor. Ann. Nucl. Energy 35, 1813e1819. Kinard, M., 2003. Efficient Numerical Solution of the Point Kinetics Equation. M.Sc. Thesis, Texas Tech. Univ.

R. Vadi, K. Sepanloo / Progress in Nuclear Energy 86 (2016) 141e161 Kuehn, T.H., Goldstein, R.J., 1978. An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli. J. Heat Transf. 100, 635e640. Kurul, N., Podowski, M.Z., 1991. On the modeling of multidimensional effects in boiling channels. In: Proceedings of the 27th National Heat Transfer Conference, Minneapolis, Minnesota, USA. Launder, B.E., Spalding, D.B., 1972. Lectures in Mathematical Models of Turbulence. Academic Press, London, England. Lavieville, J., Quemerais, E., Mimouni, S., Boucker, M., Mechitoua, N., 2005. NEPTUNE CFD V1.0 Theory Manual”. EDF. Leonard, B.P., 1991. The ultimate conservative difference scheme applied to unsteady one-dimensional advection. Comp. Methods Appl. Mech. Eng. 88, 17e74. Leveque, R.J., 2004. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, UK. NEA, 2003. WIMSD-5B, a Neutronic Code for Standard Lattice Physics Analysis. NEAMS. Nuclear Energy Advanced Modeling and Simulation program. [Online] the USA DOE/NE [Cited: December 17, 2014.] http://neams.ne.anl.gov. Ott, K.O., Neuhold, R.J., 1985. Introductory Nuclear Reactor Dynamics. American Nuclear Society Press, La Grange Park, Illinois. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow”. Hemisphere, Washington, DC. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes

161

in C”., second ed. Cambridge University Press, UK. Rhie, C.M., Chow, W.L., November 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J. 21 (11), 1525e1532. Salama, A., El-Morshedy, S.E., 2011. Numerical simulation of the IAEA 10 MW generic MTR reactor under loss of flow transient. Ann. Nucl. Energy 38, 564e577. Shih, T.H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k-ε eddy-viscosity model for high reynolds number turbulent flows - model development and Validation”. Comput. Fluids 24 (3), 227e238. SMART project. Integrated Modular Advanced ReacTor project. [Online] Korea Atomic Energy Research Institute [Cited: December 17, 2014.] https://www. kaeri.re.kr. Spellman, F.R., Drinan, J., 2001. Fundamentals for the Water and Wastewater Maintenance Operator Series: Piping and Valves. Technomic Publishing Company, Inc., Lancaster, Pennsylvania (USA). Todreas, N.E., Kazimi, M.S., 2001. Nuclear Systems II, Elements of Thermal Hydraulic Design”. Hemisphere Publishing Corporation, NY. Villarino, E., Carlos, A.L., 1993. CITVAP v.3.1 Reactor Calculation Code, MTR_PC Package. Nuclear Engineering Division, INVAP, Argentina. Wang, Z.J., 2000. A fast nested multi-grid viscous flow solver for adaptive cartesian/ quad grids”. Int. J. Numer. Methods Fluids 33, 657e680.