3D space-time analysis of ATWS in compact high temperature reactor with integrated thermal-hydraulic model in code ARCH-TH

Annals of Nuclear Energy 110 (2017) 306–316

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

3D space-time analysis of ATWS in compact high temperature reactor with integrated thermal-hydraulic model in code ARCH-TH D.K. Dwivedi ⇑, Anurag Gupta, Umasankari Kannan Reactor Physics Design Division, Bhabha Atomic Research Centre, Mumbai 400085, India

a r t i c l e

i n f o

Article history: Received 2 February 2017 Received in revised form 23 May 2017 Accepted 26 June 2017

Keywords: High temperature reactors Coupled neutronics with thermalhydraulics Temperature feedbacks AER benchmark BeO moderator LBE coolant

a b s t r a c t The Compact High Temperature Reactor (CHTR) is being designed as a technology demonstrator for comprehensive Indian high temperature reactor programme for hydrogen production and similar process heat applications. The 100 kWth CHTR core consists of 233U-Th based TRISO coated fuel particles, BeO moderator and LBE coolant with natural circulation in vertical prismatic fuel assemblies. The CHTR being a new conceptual design, necessitates comprehensive integrated Neutronics/Thermal-Hydraulics (N-TH) analyses and study. For design and safety studies, 3D space-time analyses of anticipated transient without scram (ATWS) in CHTR have been carried out in detail with indigenous code system with various cases of temperature feedbacks. One of the striking features observed in the analysis is that in case of weak Doppler feedback due to high fissile content and high temperature core conditions in CHTR, the BeO moderator plays a crucial role to limit the rise in nuclear power as well as peak fuel and coolant temperatures during such transient. For these studies, 1D-radial heat conduction in multi-channel based TH module has been developed in 3D space-time code ARCH. The viability of neutronics and adiabatic Doppler feedback capability of the code has also been examined with AER benchmark problems (AERDYN01 & 02) and results are discussed. The analysis shows that the transient peak fuel and coolant temperatures are limiting at values much below the fuel safety criteria of TRISO particle (
1600 °C) and boiling point of LBE coolant (1670 °C) even with scram failure. The significance of temperature feedback effects of BeO moderator in CHTR is seems to be a first of its kind observation and is not reported in the literature. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction India is pursuing a comprehensive high temperature reactor programme (Dulera and Sinha, 2008; Sinha, 2011). The target of the progarmme is to fulfill national needs for several crucial nuclear energy applications other than those for grid based electricity generation, such as hydrogen production by means of thermo-chemical splitting water, industrial process heat applications, unattended nuclear power packs etc. The design and development of 100 kWth CHTR as a technology demonstrator and critical facility is a key technological step to Indian high temperature programme. The reactor is being designed on the basis of guidelines such as use of proliferation resistant thorium fuel cycle, passive core heat removal by natural circulation of liquid metal coolant, compact design to minimize weight of the reactor and high discharge burn-up (Sinha and Dulera, 2010; Dwivedi et al., 2010). ⇑ Corresponding author. E-mail address: [email protected] (D.K. Dwivedi). http://dx.doi.org/10.1016/j.anucene.2017.06.051 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.

As CHTR core is a new conceptual design, it requires extensive multi-physics simulations of coupled Neutronics-Thermal hydraulics phenomena including postulated LORA, LOCA, and power setback transients. In view of that, cases of unprotected anticipated transient (i.e. or Anticipated Transient Without Scram (ATWS)) in operating core configuration at beginning of the cycle (BOC) have been analyzed. It was observed that the beryllium oxide (BeO) moderator temperature feedback is very significant especially in longer duration transients. The delayed effect of moderator temperature feedback is found to be advantageous particularly when moderator is thermally coupled to heat producing fuel with less negative Doppler feedback as in CHTR. The beryllium oxide has good neutronics properties and thermal conductivity (Manly, 1964; Konings, 2012), hence it would be a better choice for moderator over graphite in a high temperature thermal reactor. It has been earlier considered as neutron reflector in various nuclear systems such as in KAMINI reactor (Mohapatra et al., 2000), in ADS (Mirvakili et al., 2016), in space reactor (Zhang et al., 2016) etc. The significance of BeO moderator temperature feedback in high temperature core configured with thorium based TRISO fuel and

307

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

LBE coolant is a first of its kind observation as per the best knowledge of authors and is discussed in the present article. The present study has been carried out with indigenous 3D space-time multi-group diffusion theory code ARCH coupled with adiabatic heating and 1D heat conduction based TH modules. The detailed analytical studies for CHTR is carried out following the current International trend for innovative reactor configurations to perform integrated/coupled Neutronics (N)-Thermalhydraulics (T-H) studies under larger multi-physics multi-scale framework (Mylonakis et al., 2014). In literature, different code systems are reported on the coupled neutronics-thermal hydraulics simulation for various rector systems. The code package IRTRAN (Aghaie et al., 2012) is developed from neutronics code WIMS + CITATION coupled with RELAP5 to analyze reactivity transients of PWRs with thermal–hydraulic feedbacks. The coupled TRACE/PARCS code (Downar, 2006) has been developed for the stability analysis of BWRs (Xu et al., 2009). Coupled neutronics / thermal hydraulics calculations are reported for High Temperature Reactors with the DALTON-THERMIX code system (Brian Boer et al., 2008). Code DYN3D-MSR, discusses MSRs analysis with coupled N-TH simulation (Krepel et al., 2007). A multi-physics simulation tool FRENETIC is developed for the quasi-3D analysis of a leadcooled fast reactor core with the hexagonal fuel element configuration (Bonifetto et al., 2013). The open-source software code DRAGON and OpenFOAM have been coupled with integrated tight coupling approach for TH feedbacks in neutronics (Wu and Rizwan-uddin, 2016). The present analysis has been carried out with an indigenous multi-physics code system currently being developed for the design and safety analyses of CHTR with temperature feedbacks. A simple adiabatic heating as well as 1D-radial heat conduction in multi-channel based TH module have been incorporated in 3D space-time analysis code ARCH. Code ARCH has earlier been upgraded with IQS (Improved Quasi-Static, Ott and Meneley, 1969) transient scheme for purpose of efficient coupled simulations (Gupta, 2012; Dwivedi et al., 2013; 2015a). At recent the conduction based TH module has been incorporated in the code for temperature feedbacks. The incorporated TH module is for single phase lumped coolant and limited to simplified model of flow in every channel with fixed inlet temperature and mass flow rate. The further development in this module is underway for simulation of natural circulation of liquid metal coolant. The natural circulation of coolant is being considered in close-loop of multichannel connected only through lower and upper plenum and will be reported in near future. The plan of the paper is as follows; schematics and conceptual core design of CHTR has been described in Section 2. Section 3 contains the details of the computer codes used for the transient analysis of CHTR. The methodology and model adopted for temperature feedbacks in 1D heat conduction based TH module of ARCH-TH has been described in Section 4. The benchmark qualifications of developed code and the analysis of ATWS in CHTR have been presented in Section 5. Section 6 summarizes the study and discusses future work plan.

2. Description of CHTR core In the conceptual core design of 100 kWth CHTR, it is a U-Thorium fuelled, BeO moderated and liquid metal (Lead Bismuth Eutectic, LBE) cooled prismatic block type vertical reactor (Dulera and Sinha, 2008; Sinha and Dulera, 2010; Sinha, 2011). The core (Fig. 1) consists of 19 hexagonal fuel assemblies (FAs). These FAs are surrounded by 150 mm thick top and bottom axial reflectors blocks and eighteen cylindrical BeO reflector blocks followed by graphite radial reflector. For power regulation and control,

233

Downcomer Fuel Tube Movable BeO block

Graphite radial reflector LBE coolant in channel

Fixed BeO reflector

Fig. 1. Schematic cross-sectional view of CHTR core.

twelve control rods (CRs) made of tantalum alloy are provided in centrally located coolant channels of twelve outer FAs. The reactor has two independent shutdown systems for scram. The primary shutdown system contains six shutoff rods similar to the CRs falls under gravity in six inner coolant channels. Twelve axially movable BeO reflectors blocks act as secondary shutdown system, which is to be withdrawn out to introduce high neutron leakage. It has also six burnup compensation rods inserted axially in fixed BeO blocks for coarse control of initial excess reactivity during long refueling interval (Dwivedi et al., 2010). The major design parameters of CHTR are given Table 1. The fuel assembly of CHTR (Fig. 2a) contains BeO moderator with centrally located 700 mm long graphite tube with fuel. In the central portion of graphite tubes, LBE coolant flows from bottom to top plenum through natural circulation with inlet temperature of 900 °C. The core average coolant mass flow rate in hot operating condition at full power is 6.7 kg/s (Dulera et al., 2008). Each graphite tube carries within it the fuel inside twelve equispaced longitudinal bores of 10 mm diameter. These bores are filled with fuel compacts made from TRISO coated particles embedded in graphite matrix. The TRISO particles (Fig. 2b) are in the form of micro-spheres of (233U + Th) carbide kernel coated with three types of layers of soft pyrolitic carbon, SiC and hard carbon layer (Price, 2012; Konings, 2012). In the present design of CHTR core, Gadolinium, added in the fuel kernel of the TRISO particles

Table 1 Major design parameters of CHTR. Reactor power Core configuration Fuel Fissile inventory Volume packing of TRISO particles in fuel compacts Fuel channels/pitch Core life cycle Core burnup Moderator Reflector Coolant Coolant inlet/outlet Fuel length Power regulation Primary SDS Secondary SDS

100 kWth prismatic block, vertical U-Th in TRISO particles fuel 33.75 wt% 42.3% 233

19/135 mm 15 effective FPYs 68,000 MWd/t BeO BeO and Graphite Lead Bismuth Eutectic 900 °C/1000 °C 700 mm 12 control rods 6 shut-off rods 12 movable BeO blocks

308

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316 Table 2 Energy group boundaries of homogenize macroscopic cross-section.

Fig. 2a. Schematic cross-sectional view of single channel/fuel assembly in CHTR.

Number of Energy Group

Upper Energy Boundary (eV)

1 2 3 4 5

1.96  107 6.06  106 9.12  103 2.60 0.40

eration of homogenized condensed 5-group parameters with ITRAN. These homogenized 5-group cross-sections have been used for whole core calculation with code ARCH. The energy group boundaries of homogenized macroscopic cross-section in 5-group are given in Table 2. 3.2. 3D space-time analysis code ARCH

Fig. 2b. Schematic of TRISO coated particle.

of central fuel assembly, is used as burnable poison for passive control of reactivity. 3. Neutronic simulation methodology and tools The CHTR simulations are being carried out in standard two stage analyses. In the first stage, the lattice cells are being simulated in minute details with integral transport theory code ITRAN (Krishnani, 1982) with multi-group data library. In this step, fewenergy-groups cell-homogenized cross section data are being generated for each cell type. The temperature and density dependent few-group homogenized cell parameters are also tabulated for feedbacks. The delayed neutron kinetics parameters are also computed. In the second stage, the whole core simulations are being carried out with 3D space-time diffusion theory code ARCH (Gupta, 2012). The code ARCH uses few-group condensed and cell homogenized data generated by ITRAN. The brief descriptions of the in-house codes ITRAN and ARCH have been given in the following sub-sections for information. 3.1. Lattice analysis code ITRAN The code ITRAN is based on integral neutron transport theory for lattice simulation of nuclear reactors. It uses first flight collision probability and interface current methods. The code generates homogenized condensed few-group cross-sections and lattice parameters of various fuel types and reflector materials in the core as well as delayed neutron kinetics parameters at different burnup. It requires multi-group nuclear data libraries in WIMS format as input. In the present analysis, ENDFB-VI, 172 group nuclear data library has been utilized. The data library sensitivity analysis for the study is being planned separately. The effect of double heterogeneity in CHTR core has been studied with RPT method (Kim et al., 2005) and found to be insignificant as high packing fraction of TRISO particles in compacts (Table 1). Therefore, for present study, the TRISO particles along with graphite matrix in the fuel compacts has been smeared in the fuel bores of CHTR fuel assembly for gen-

The code ARCH has been developed as a comprehensive computer code for various time-dependent neutronics analyses in both rectangular as well as triangular-Z lattice based core designs (Gupta, 2012). The code solves the time-dependent multi-group diffusion equation along with delayed neutrons precursor’s equations with Fully Implicit as well as with efficient Improved Quasi Static (IQS) method. For presented feedback studies, the IQS option has been used. The spatial discretisation is carried out by Finite Difference Methods (FDM). Apart from the conventional schemes, the code uses advanced Krylov subspace method based algorithms. For the presented analysis, the 1D-radial heat conduction in multi channel based TH module has been developed in code ARCH. This version of the code is named as ARCH-TH (Analysis of Reactor kinetics in Cartesian and Hexagon geometry with ThermalHydraulic feedbacks). 3.2.1. The IQS module of ARCH-TH Reactor transient simulations become computationally extensive in case of coupled N-TH analysis with direct methods. Thus an efficient Improved Quasi-Static (IQS) technique has also been included in code ARCH (Dwivedi et al., 2013). The IQS method is based on the factorization technique of the neutron flux into the product of the amplitude and shape functions (Ott and Meneley, 1969). These functions can be obtained by solving the amplitude and shape equations, respectively, which are decomposed from the original neutron diffusion equation. In most of the rector transients, the shape of the neutron flux distribution in the core changes slowly. The numerical calculation for the weakly time dependent shape equation is time-consuming, which is calculated at relatively large time steps in IQS method than that is used in the fully implicit scheme of direct method. The rapidly varying amplitude function can be easily solved using smaller time steps in point kinetics equations, due to its inexpensive computation costs. Computational efficiency of the IQS method becomes more significant for reactivity-initiated transients, where the temporal variation of the shape function is relatively slower than that of the amplitude function. In IQS method used here, the time derivative of shape function is replaced by a backward difference of first order (Ott and Meneley, 1969). The two–tier scheme has been adopted for numerical solution of time varying functions. The shape function is supposed to be weakly time dependent; therefore, it is being computed at larger time intervals (i.e. macro time steps  102 to 101 sec.). In this approach delayed source term of precursors is calculated directly by implicit time integration from the flux history (Stacey, 2001) and treated as inhomogeneous source in shape equation. The shape equation is being solved in IQS module with Bi-conjugate Gradient Stabilized (BiCGStab) technique (Van der

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

309

Vorst, 2002; Saad, 1996). The rapidly varying amplitude function is being computed at smaller time intervals (micro steps  104 to 103 sec.) in point kinetic like equations in every macro time interval. The point kinetic parameters (ie. dynamic reactivity, integrated precursors’ concentration) are recomputed periodically to account for changes in the material properties (absorber rod movements and feedbacks from changes in temperature profile) between the successive macro steps. The fourth order Runge-Kutta method has been adopted to compute Point kinetics equations for better accuracy. Linear interpolation of dynamic reactivity has been considered within each macro time interval. 4. Methodology and incorporation of conduction based TH module For temperature feedbacks in steady-state as well as during transients; a thermal hydraulic module based on 1D-radial heat conduction in multi-channel model has been developed in code ARCH and is referred as ARCH-TH. The nuclear heat generated in the fuel present in twelve bores of each fuel channel/assembly in CHTR is being conducted in the adjacent region of graphite fuel tube and BeO moderator and ultimately transferred to LBE coolant flowing through the center of each fuel tube (Fig. 3a). Solving 1Dradial heat equation in hexagonal axial meshes of the channels, outer surface of each hexagonal axial mesh is assumed to be in circular form and twelve fuel bores are considered in annular shape conserving cross-sectional area (Fig. 3a). In the present model, heat transfer in the axial direction is ignored due to large ratio of height to diameter of the fuel assembly. Additionally, the single phase coolant is assumed to be flowing in each cylindrical coolant channel with given inlet temperature and mass flow rate (Fig. 3b). The net heat flux at outer surface of the BeO moderator is also neglected by imposing adiabatic boundary condition, after taking advantage of relatively small conductivity as well as temperature gradient at these boundaries. In the model, time dependent radial heat conduction equations are being solved for the nuclear heat transfer from fuel to adjacent region and finally to coolant in each axial mesh of every fuel assembly/channel independently. These axial meshes in each FA are simulated through the lumped coolant flow model as outlet coolant temperature of every mesh is taken as the inlet for the next elevated axial mesh in the channel. The whole core calculation model for 3D space time analysis of CHTR has been illustrated in Fig. 4. In the model, the time dependent radial heat equation has been numerically solved by the corner-mesh finite difference scheme. The corner mesh scheme has been preferred as it gives better information about nodal temperatures (Jain, 1989). The time dependent heat conduction equation in one dimensional cylindrical geometry (Todreas and Kazimi, 2012) is given as follows:

Fig. 3b. Schematic of coolant flow in multi-channel model with known inlet temperature (Tci) and mass flow rate (Cw).

qm Cp


 @Tðr; tÞ 1 @ @Tðr; tÞ ¼ q000 ðr; tÞ þ rk @t r @r @r

ð4:1Þ

where, qm , Cp and k are material density, specific heat at constant pressure and thermal conductivity at position ‘r’ respectively. q000 ðr; tÞ is the power density at position ‘r’ and time ‘t’. Tðr; tÞ represents the temperature at point ‘r’ and time ‘t’. At the interface of two conducting material A and B, the following boundary condition is applied to temperature Tðr; tÞ:




kA

@T @r






¼ kB A

@T @r



ð4:2Þ B

The derivatives are taken at the two sides of the interface. At the interface of conducting material and convective (coolant) material (i.e. at inner wall of fuel tubes in CHTR core), boundary condition is as follow:




kA

@T @r



¼ hðTw  Tb Þ

ð4:3Þ

rw

where, derivative is taken at the boundary of the conductive and convective materials. The variable ‘h’ is the heat transfer coefficient; ‘Tw ’ and ‘Tb ’ are the temperature of the surface of the coolant channel wall and average bulk temperature of the coolant in the mesh respectively. The heat transfer coefficients ‘h’ is determined as follow:

h¼

kc Nu Dh

ð4:4Þ

where, ‘kc ’ and ‘Dh ’ are the thermal conductivity of coolant and hydraulic diameter of the coolant channel respectively. ‘Nu ’ is the Nusselt number. Nusselt number has been computed by Cheng & Tak correlation for low Prandtl number fluid flow of LBE coolant in CHTR. This correlation is proposed in Cheng and Tak (2006) and given as follow,

Cheng & Tak correlation : Nu ¼ A þ 0:018  P0:8 e where, ‘A’ is defined as,

8 if Pe 6 1000 > < 4:5 4 A ¼ 5:4  9  10  Pe if 1000 < Pe 6 2000 > : 3:6 if Pe > 2000

Fig. 3a. Schematic cross-sectional view of channel/ fuel assembly in CHTR and its equivalent 1D radial representation.

where, ‘Pe’ is the Peclet number, defined as:

ð4:5aÞ

310

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

Core map as at axial mid-plane Fuel assembly (FA) FA with Gadolinium Movable BeO reflector Fixed BeO reflector block Graphite reflector Reactor vessel

Finer details such as downcomers and other under lying materials in hexagons have been considered by filling the individual triangles. Each hexagon axial mesh is divided in 24 triangles (2nd order division) in the simulation with code ARCH-TH as below:

Fig. 4. The whole core simulation model for 3D space time analysis of CHTR.

Pe ¼ Re  Pr

ð4:5bÞ

In the present study of CHTR transient, Cheng & Tak correlation has been considered over the other correlations option available in the TH module. The other option of correlations for computation of heat transfer in case of liquid metal coolant such as Lubarski & Kaufman, Ibragimov, Notter & Sleicher and Kirrilov & Ushako correlations (Borgohain et al., 2011) are also made available in the module. In Eq. (4.5b), Reynolds number (Re) and Prandtl number (Pr) are defined as:

Re ¼

Pr ¼

Gc Dh

l

;

Cp l kc

ð4:6aÞ

ð4:6bÞ

where, ‘Gc ’ is the mass flux, ‘m’ is the dynamic viscosity of the coolant. In case of CHTR, ‘Gc ’ and ‘Dh ’ are taken as follow:

Cw Gc ¼ pðr2c  r2r Þ

ð4:7aÞ

Dh ¼ 2ðrc rr Þ

ð4:7bÞ

for annulus channel.where, Cw is the average coolant mass flow rate in each channel. The change in given inlet temperature and mass flow rate of coolant in each of the channels is assumed to be negligible during reactor transient in CHTR. ‘rc ’ is the radius of the coolant channel centrally located in the fuel tube and ‘rr ’ is radius of the absorber rod partially inserted in the channel for reactor control. The change in properties of Pb-Bi eutectic coolant such as isobaric specific heat (Cp ), mass density (qm ), viscosity (m) and thermal conductivity (kc Þ with temperature (Tc ) are being considered as per given equations in OECD NEA-6195 (2007). The 1D-radial heat conduction equations is being solved by discretizing each conducting region in several fine radial divisions/ meshes in each axial mesh of the channels. For feedbacks, volume weighted average temperatures of each type of material (fuel/moderator) are being computed in every axial mesh of each fuel assemblies in the core. The heat conduction equations have been solved with fully implicit method. The Bi-conjugate Gradient Stabilized technique has been utilized for solution. The 5-group homogenized cross sections data-set generated at various temperatures of the

materials at small intervals with lattice code ITRAN is utilized for temperature feedbacks in the study. The linear interpolation scheme is adopted for evaluations of multi-group cross-section data at core material temperature during transient for reactivity feedback. The interpolation scheme with generated data set at small temperature intervals is admissible and computationally efficient. 5. Results and analyses The ATWS transient simulation analysis in CHTR has been carried out with indigenously developed code system (ITRAN + ARCH-TH). The neutronics of developed code has been qualified with benchmark analysis of transient without any temperature feedback as well as with adiabatic Doppler feedback. The results of benchmark problem with code and subsequent analysis of ATWS in CHTR have been presented in the following sub-sections. 5.1. AER benchmark qualification of ARCH-TH The viability of the code ARCH-TH has been examined with well known AER benchmark problems (AER-DYN-001 & 002). The predictions of code ARCH-TH and comparison with reference results have been presented as the following. 5.1.1. AER benchmark without feedback (AER-DYN-001) The benchmark problem deals with the case of asymmetric control rod ejection transient without any feedback in VVER type geometry (Keresztúri et al., 1992). Initial power is close to zero and power rise during transient is not very large. In the transient, asymmetric control rod ejection is followed by shut-down system activation. A control rod is ejected from 0.0 to 0.08 s.; a constant velocity is assumed (Fig. 5a). Steady state simulation with ARCH code shows that the reactivity worth of ejected control rod is 0.00484 as compared to 0.00482 with KIKO3D. The power variations during transient has been computed and compared with code KIKO3D and DYN3D (Keresztúri et al., 2000). Results obtained with IQS module of code ARCH-TH are found to be matching well (Fig. 5b). The variations of dynamic reactivity during transient have been predicted and are found to be in good agreement with KIKO3D (Fig. 5c). The axial power distribution in given fuel assembly at various moments of transient with ARCH-TH has also been predicted and compared (Fig. 5d). The normalized radial power

311

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

10

Axial Power at position 203

Fig. 5a. Horizontal map of half-reactor in AER-Dyn001 benchmark, the ejected rod is marked with No. 26.

ARCH-TH at 0 s ARCH-TH at 6 s KIKO3D at 0 s KIKO3D at 6 s DYN3D at 0 s DYN3D at 6 s 5

0 2

3

4

5

6

7

8

9

10

11

No. of Axial Layer Fig. 5d. Axial power distribution in the 2nd neighbor FA of the ejected rod at t = 0 and t = 6 s.

2.0

ARCH-TH KIKO3D DYN3D

Relative Power

1.5

1.0

0.5

Fig. 5b. Variation of power with time in AER-Dyn001.

0.0 422

427

432

437

442

FA

0.006

ARCH-TH KIKO3D

0.003

Fig. 5e. Normalized radial power distribution at t = 0 s in 3rd axial plane (ejected rod in No. 426).

Reactivity

0.000 -0.003

8 -0.006 -0.009

ARCH-TH KIKO3D DYN3D

6

-0.015 -0.018 0

1

2

3

4

5

6

Time(sec)

Relative Power

-0.012

4

2

Fig. 5c. Variation of core dynamic reactivity during transient in AER-Dyn001.

distributions have been estimated with ARCH-TH and compared for a particular axial plane during transient shown in Figs. 5e and 5f.

0 422

427

432

437

442

FA Fig. 5f. Normalized radial power distribution at t = 6 s in 3rd axial plane.

5.1.2. AER benchmark with adiabatic Doppler feedback (AER-DYN002) The validation of fuel temperature feedback capability of ARCHTH has been carried out with VVER-440 based kinetics benchmark AER-DYN-002 for adiabatic heating of the fuel (Grundmann and

Rossendorf, 2000). In the benchmark problem, the transient is initiated by the ejection of the eccentric rod (with reactivity worth of 2.0$) in 0.16 s at hot zero power (Fig. 6a). The control rod removal

312

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

2000

ARCH-TH

o

Max Fuel Temperature ( C)

2500

Fig. 6a. Core map with the fuel types in AER-Dyn002 benchmark. Absorbers of bank K6 at positions of types 4/2.

KIKO3D

1500

1000

500

0 0.00

0.25

0.50

0.75

1. 00

1.25

1.50

1.75

2.00

Time (Sec) Fig. 6c. Variation of peak fuel temperature in AER-Dyn002.

6000

5000

Integral Power (MW-Sec)

speed of 12.5 m/s is considered. The initial reactor power is 1.375 kW. The feedback mechanism is based on the adiabatic increase of fuel temperature from the initial value of 260 °C. The excess heat produced is not considered to be removed from the fuel during the transient. The rising power is arrested with fuel temperature feedback only. The transient is simulated up to 2 s. The predictions of code ARCH-TH have been examined for this benchmark and it is found to be in good agreement with other results (Grundmann and Rossendorf, 2000). The variation of nuclear power (Fig. 6b), peak fuel temperature (Fig. 6c), time integral power (Fig. 6d) and core reactivity (Fig. 6e) have been computed by ARCH-TH and are found to be matching very well with code KIKO3D (Kereszturi et al., 2009). The variation of power peaking factor predicted with ARCH-TH has been computed and was found to be matching well with DYN3D (Fig. 6f). The normalized radial averaged axial power distribution at steady state (t = 0 s, Fig. 6g) as well as after 2 s of initializing the transient (t = 2.0 s, Fig. 6h) have been compared and found to be in good agreement with other benchmarked codes.

4000

3000

ARCH-TH

1000

0 0.00

5.2. Benchmark qualification of 1D-radial heat conduction based TH module

KIKO3D

2000

0.25

0.50

0.75

1.00

1.25

1.50

1.75

2.00

Time (Sec) Fig. 6d. Variation of integral power with time.

The 1D-radial conduction based TH module in ARCH-TH is primarily developed for the transient analyses of single phase molten metal Lead Bismuth Eutectic (LBE) cooled reactor core e.g. CHTR. The TH module integrated with a point kinetics solver has been examined earlier (Dwivedi et al., 2017) with benchmark analysis

100000

ARCH-TH Nuclear Power (MW)

80000

KIKO3D

60000

40000

20000

Fig. 6e. Variation of core dynamic reactivity.

0 0.0

0.1

0.2

0.3

0.4

Time (Sec) Fig. 6b. Nuclear Power variation with time in AER-Dyn002 benchmark.

0.5

based on beam interruptions in LBE-cooled subcritical system (D’Angelo and Gabrielli, 2003). The same point kinetics equation solver is being utilized for solution of point kinetics in IQS module of ARCH-TH. The parameters and input data for simulation of steady state as well as beam interruptions transients in LBE cooled,

313

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

Fig. 6f. Variation of power peaking factor. Fig. 6h. Normalized axial (radially averaged) power distribution at t = 2.0 s.

Results of conduction based TH module of code, coupled with Point Kinetics Fuel centerline temperature Clad surface temperature Fuel surface temperature Coolant temperature

1073 1023 973

o

Temperature ( K)

923 873 823 773 723 673 623 573 0

15

30

45

60

75

90

Axial position(cm)

Fig. 7a. Steady state axial temperature distribution results of benchmark with TH module of code, coupled with point kinetics.

Fuel centerline temperature predicted with TH module of code, coupled with point kinetics beam trip transient interruption for 12s interruption for 6s beam interruption for 3s beam interruption for 1s

Fig. 6g. Normalized axial (radially averaged) power distribution at t = 0.0 s.

1023

5.3. Analyses of anticipated transient without scram (ATWS) in CHTR The ATWS analysis in hot operating core configuration of CHTR has been carried out with the present 1D-radial heat conduction based TH module of ARCH-TH. The transient case has been studied with several feedback cases viz. with only fuel temperature feedback, with fuel and moderator temperature feedbacks as well as with fuel, moderator and coolant temperature feedbacks taken together. The transient is assumed to be initiated due to inadvertent withdrawal of a control rod, consequently about 1.0 mk positive reactivity in 2.3 s has been added in the hot critical core of CHTR at full power during BOC at zero burnup. The effective delayed neutron fraction (beff) in CHTR has been found to be 4.6 mk with lattice code ITRAN. The fuel temperature (Doppler) coeffi-

973 923

Temperature (K)

MOX fuelled assembly of XADS-type Accelerator-driven System have been taken as provided in the benchmark (D’Angelo and Gabrielli, 2003). The predictions of the TH module in steady state (Fig. 7a) and in cases of beam interruptions transients (Figs. 7b and 7c) have been found in good agreement with reference results (D’Angelo and Gabrielli, 2003).

873 823 773 723 673 623 573 0

5

10

15

20

25

Time (s)

Fig. 7b. Fuel centerline temperature at the fuel zone mid-plane transients induced by beam interruptions of different duration: 1 s, 3 s, 6 s, 12 s and by a definitive beam trip with TH module in code, coupled with point kinetics.

cient of reactivity in hot operating core of CHTR is found to be 5.6  106 k/°C. Due to higher mass, density and also the macroscopic thermal scattering cross-section of BeO as compare to gra-

314

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316 Coolant outlet temperature predicted with TH module in code, coupled with point kinetics beam trip transient, intrruption for 12s , interruption for 6s beam interruption for 3s , and beam interruption for 1s

1110 1100

683 673

1090

Temperature ( C)

663

o

Temperature ( K)

1080

o

653 643 633 623 613 603 593 583

1070 1060

Peak Fuel Temperature Peak Moderator Temperature Peak Coolant temperature

1050 1040 1030

573 0

5

10

15

20

1020

25

Time (s)

0

phite in the lattice, major contribution in neutron moderation is by BeO moderator. The CHTR core is found to be under moderated with average neutron energy 1 eV. At operating temperature, when temperatures of both BeO and graphite are being changed, the moderator temperature coefficient (MTC) is found to be 1.8  105 k/°C as compare to 1.49  105 k/°C with only BeO temperature changes. The MTC is found to be more negative than the fuel temperature coefficient. The coolant temperature and density feedback coefficient of reactivity in hot operating core of CHTR is negative but very small and estimated about to be 3.8  107 k/°C in hot operating core. The aforesaid temperature coefficients are computed with ITRAN + ARCH code system. For ATWS scenario, it is assumed that no shutdown system is activated during transient and power rise is being arrested by feedbacks from core temperature rise (Fig. 8a). The study shows that the coolant temperature and density feedbacks are almost negligible where as moderator temperature feedback effect is very significant (Dwivedi et al., 2015b). The ATWS with fuel, moderator and coolant temperature feedbacks has been studied in details. The coolant inlet temperature in hot operating core of CHTR is considered as 900 °C. In the steady-state operations at full power, core average coolant mass flow rate is taken as 6.7 kg/s. Therefore, the average mass flow rate in each of nineteen coolant channel is

with only Fuel temperature feedback with fuel and moderator temperature feedback with fuel, moderator and coolant temperature feedback

900 800

Power (kW)

700 600 500 400 300 200 100 0 0

100

200

300

400

500

600

Time(sec) Fig. 8a. Power variation with time during ATWS in CHTR in various temperature feedbacks.

100

200

300

400

500

600

Time(sec) Fig. 8b. Temperature variation with time during ATWS in CHTR with fuel, moderator and coolant temperature feedbacks in simulation.

Net core dynamic reactivity (k) with fuel, moderator and coolant temperature feedback

0.0012 0.0010 0.0008 0.0006

Reactivity

Fig. 7c. Variations of outlet coolant temperature during transients induced by beam interruptions and by trip with TH module in code, coupled with point kinetics.

0.0004 0.0002 0.0000 -0.0002 -0.0004 -0.0006 0

100

200

300

400

500

600

Time(sec) Fig. 8c. Net core dynamic reactivity variation during ATWS with fuel, moderator and coolant temperature feedbacks.

considered 0.353 kg/s/channel as input to TH module. The Chen and Tak correlation is considered for heat transfer coefficients estimation in the present study. In the analysis, an ideal heat sink condition has been assumed i.e. temperature of coolant coming after passing through the heat sink is constant at core inlet during transient. The steady state simulation has been carried out before initializing the transient with temperature feedbacks by code ARCH-TH. The axial fuel and coolant temperature profile in the hottest channel in steady state have been shown as predictions at time t = 0 in Figs. 8d and 8e respectively. In the steady state, radial temperature profile in the maximum powered axial mesh of hottest channel has also been predicted and shown as graph for t = 0 in Fig. 8f. Radially averaged normalized axial power distribution under steady state condition in the core has been shown in Fig. 8g. The variations of nuclear power, temperature and reactivity in the core have been predicted up to 600 s of the transient. The variation of nuclear power during ATWS has been predicted considering various cases of temperature feedbacks in the simulation (Fig. 8a). It was found that with feedbacks from fuel and moderator temperature, power is rising up to 4.7 times of the initial power. It was also observed that the coolant temperature and density feedback is not playing any significant role (Fig. 8a). The variation of

315

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

Axial fuel temperature pr ofile during transient in hottest channel

0.06 0 sec 20 sec 35 sec 50 sec 100 sec 150 sec 200 sec 400 sec 600 sec

o

Fuel Temperature ( C)

1075 1050 1025 1000 975 950 925 0

15

30

45

60

75

Radially averaged power

1100

0.05 0.04

0.03

0.02 0.01

90

Axial position (cm)

0.00

Fig. 8d. Axial fuel temperature profile in hottest fuel assembly during transient.

20

30

40

50

60

70

80

90

Axial position (cm) Axial coolant temperature profile during transient

1100 1080

o

Coolant Temperature ( C)

in hottest channel

0 sec 20 sec 35 sec 50 sec 100 sec 150 sec 200 sec 400 sec 600 sec

1060 1040 1020 1000 980 960 940 920 900 880 0

15

30

45

60

75

90

Axial Position (cm) Fig. 8e. Axial coolant temperature profile in hottest channel during transient.

Radial temperature profile during transient at max powered mesh in hottest channel 1100 0 sec 10 sec 20 sec 35 sec 50 sec 100 sec 150 sec 200 sec 300 sec 400 sec 500 sec 600 sec

1060

o

Temperature ( C)

1080

1040 1020 1000

Fig. 8g. Radially averaged normalized axial power distribution in the core in steady state condition.

The change in core material temperature profiles have also been predicted with ARCH-TH in the study. The variation of axial temperature profile of fuel and coolant of the hottest channel in steady-state as well as during the transient are shown in Figs. 8d and 8e. The plots in Fig. 8d show the fuel temperature profile during heating of the core at various moments from beginning to 150 s of the transient. It is observed that peak fuel temperature point is shifting towards the center of the fuel assembly and it is due to the center peaked power production profile in the channel. Whereas, after about 150 s, the peak fuel temperature point is again shifting towards the top of the core and it is stabilizing after 600 s. It is due to higher heat removal rate compared to the production in this duration. The coolant temperature profile is also predicted (Fig. 8e) at various moments during transient in the hottest channel of the core. The radial temperature profile during the transient in peak powered axial hexagon mesh (i.e. near the middle of the hottest channel) has also been shown in Fig. 8f. The steady state radial temperature profile (at beginning of the transient), the maximum temperature point is found to be at outer surface of the annular fuel region and the adjacent BeO moderator region is at the same temperature. It is due to adiabatic boundary conditions assumed at outermost radial surface. After the transient has been initiated and power starts to rising, the peak temperature point is in the annular fuel region with asymmetric parabolic shape of temperature profile in fuel (i.e. heat source) region from beginning to nearly 150 s of transient. It is found that the core heat is eventually transferred to the coolant of each channel, and peak temperature point is in the BeO moderator region after around 150 s of the transient.

980 0

1

2

3

4

5

6

7

8

Radial position (cm) Fig. 8f. Radial temperature profile in peak powered axial mesh of the hottest channel during transient.

core reactivity and peak temperature of the core component during transient are presented and analyzed in details with the case of fuel, moderator and coolant temperature feedbacks in the simulation (Figs. 8b and 8c).

6. Conclusions and discussions The 100 kWth Compact High Temperature Reactor (CHTR) is being designed and developed as technology demonstrator for comprehensive Indian high temperature reactor programme. The physics design of the core necessitates a comprehensive integrated neutronics and thermal-hydraulics analysis for which a capability has been developed and 3D space-time analysis of ATWS in CHTR has been carried out with various temperature feedbacks using the developed code ARCH-TH. The analysis is done with more efficient IQS time integration scheme. The neutronics of the code

316

D.K. Dwivedi et al. / Annals of Nuclear Energy 110 (2017) 306–316

without feedbacks has been qualified with AER-DYN001 benchmark. The Doppler feedback capability of the code has been validated with AER-DYN002 benchmark. The predictions of code ARCH-TH are found to be in very good agreement with other benchmarked codes as in case of no feedback and adiabatic Doppler feedback in simulations. The 1D-radial heat conduction based TH module has also been introduced in the code and examined with a benchmark problem simulating beam interruptions in an accelerator-driven sub-critical system. A detailed study has been carried out to simulate the ATWS in hot operating CHTR during BOC with code ARCH-TH with various feedbacks conditions. The transient power, peak fuel temperature, peak coolant temperature, moderator temperature and radial as well as axial temperature profile in each channel/fuel assembly of CHTR have been predicted and analyzed. The analysis shows that the BeO moderator temperature feedback along with Doppler feedback is very significant and plays a vital role in arresting the power rise at much lower values as compared with case of only fuel temperature feedback. The LBE coolant temperature and density feedback in CHTR was found to be not much of significance. The development of conduction based TH module in code ARCH-TH is in progress for the simulation of natural circulation of liquid metal coolant in close-loop of multichannel connected only through lower and upper plenum. It is also planned to augment the module with essential modifications for considering two phases of coolant to enhance applicability of the code for simulation of water cooled reactor system. These results will be reported in the near future. Acknowledgements We express our acknowledgement to Dr. P.D. Krishnani, Raja Ramanna Fellow, BARC. We are very grateful to Director, Reactor Design and Development Group, BARC for encouragement and support. We are very grateful to Dr. H.P. Gupta, Raja Ramanna Fellow, ThPD, BARC, for fruitful discussion and suggestions. The authors are sincerely thankful to Shri A Borgohain, RED, BARC, for discussion regarding various correlations used in low Prandtl number fluid flows of liquid metal coolant. References Aghaie, M., Zolfaghari, A., Minuchehr, A., 2012. Coupled neutronic thermal– hydraulic transient analysis of accidents in PWRs. Ann. Nucl. Energy 50, 158– 166. Boer, B., Lathouwers, D., Ding, M., Kloosterman, J.L., 2008. Coupled neutronics/ thermal hydraulics calculations for High Temperature Reactors with the DALTON - THERMIX code system. International Conference on the Physics of Reactors, ‘‘Nuclear Power: A Sustainable Resource”. Casino-Kursaal Conference Center, Interlaken, Switzerland, September 14–19. Bonifetto, R. et al., 2013. A full-core coupled neutronic/thermal-hydraulic code for the modeling of lead-cooled nuclear fast reactors. Nucl. Eng. Des. 261, 85–94. Borgohain, A., Maheshwari, N.K., Vijayan, P.K., Sinha, R.K., 2011. CFD analysis on heat transfer in low Prandtl number fluid flows. NURETH 14, 457. Cheng, X., Tak, N.-I., 2006. Investigation on turbulent heat transfer to lead–bismuth eutectic flows in circular tubes for nuclear applications. Nucl. Eng. Des. 236, 385–393. D’Angelo, A., Gabrielli, F., 2003. Benchmark on Beam Interruptions in an Accelerator driven System- Final Report on Phase I Calculations. NEA/NSC/DOC(2003) 17, Italy. Downar, T.J., 2006. PARCS n2.7 US NRC Core Neutronics Simulator. School of Nuclear Engineering, Purdue University, W. Lafayette, Indiana.

Dulera, I.V., Sinha, R.K., 2008. High temperature reactors. J. Nucl. Mater. 383, 183– 188. Dwivedi, D.K., Gupta, A., Krishnani, P.D., 2010. Challenges in core reactivity management and control optimization in physics design of Compact High Temperature Reactor. BARC Newsletter 317, 22–28. Dwivedi, D.K., Gupta, A., Srivastava, A.K., Krishnani, P.D., 2013. Incorporation of Improved Quasi Static (IQS) scheme module in 3-D space-time analysis code ARCH. ARP 2013 – Advance in Reactor Physics: Simulation Techniques and Analysis Methodologies, Anushakti Nagar, Mumbai, India, October 23–25. Dwivedi, D.K., Gupta, A., Krishnani, P.D., 2015a. 3D space time analysis of anticipated transient without scram in CHTR with fuel temperature feedback. Thorium Energy Conference 2015 (ThEC15), Mumbai, India, October 12–15, 2015. Dwivedi, D.K., Gupta, A., Krishnani, P.D., 2015b. Significance of moderator temperature feedback in 3-D space time safety analysis of Compact High Temperature Reactor (CHTR). Proceedings of CANSAS – 2015 & NRTHS – 2015, Anushaktinagar, Mumbai, India, 8–11 Dec. 2015. Dwivedi, D.K., Gupta, A., Kannan, U., 2017. Benchmark qualification of 1D-radial heat conduction based TH module of point kinetics code PATH., report n.-RPDD/ HTR/76/2017. Grundmann, U., Rossendorf, FZ, 2000. AER Benchmark Specification Sheet- Test ID: AER-DYN-002. http://aerbench.kfki.hu/aerbench/Dyn002.pdf. Gupta A., 2012. ARCH: A 3D space time analysis code in Cartesian and Hexagon geometries. NSRP-19, Mamallapuram, TN, India, 12–15, Dec. Jain, V.K., 1989. Point-kinetics model with one-dimensional (radial) heat conduction formalism. report no. BARC-1452. Keresztúri, A., et al. 1992. A Three-Dimensional Hexagonal Kinetic Benchmark Problem. 2nd AER Symposium, Paks Hungary. Keresztúri, A., et al. 2000. AER benchmark Specification Sheet. aerbench.kfki.hu/ aerbench/Dyn001.doc. Kereszturi, A., et al., 2009. http://aerbench.kfki.hu/aerbench/dyn002_solaeki.txt. Kim, Y., Kim, K., Noh, J.M., 2005. Reactivity-Equivalent Physical Transformation for Homogenization of Double-Heterogeneous Fuels. Transactions of the Korean Nuclear Society Autumn Meeting, Busan, Korea, October 27–28, 2005. Konings, R., 2012. Comprehensive Nuclear Materials, vol. 2. Elsevier. Krepel, J., Rohde, U., Grundmann, U., Weiss, F., 2007. DYN3D-MSR spatial dynamics code for molten salt reactors. Ann. Nucl. Energy 34, 449–462. Krishnani, P.D., 1982. CLUB-A multi group integral transport theory code for analysis of cluster lattices. Ann. Nuclear Energy 9, 255–260. Manly, W.D., 1964. Utilization of BeO in reactors. J. Nucl. Mater. 14 (1964), 3–18. Mirvakili, S.M. et al., 2016. Neutronic and thermo hydraulic analysis of a modeled subcritical uranyl nitrate aqueous reactor driven by 30-MeV protons. Ann. Nucl. Energy 97, 171–178. Mohapatra, D.K. et al., 2000. Moderator temperature effect on reactivity in light water moderated experimental reactors. Ann. Nucl. Energy 27, 969–983. Mylonakis, A.G., Varvayanni, M., Catsaros, N., Savva, P., Grigoriadis, D.G.E., 2014. Multi-physics and multi-scale methods used in nuclear reactor analysis. Ann. Nucl. Energy 72, 104–119. OECD/NEA, Nuclear Science Committee, Working Group on Lead-Bismuth Eutectic, 2007. Handbook on Lead-Bismuth Eutectic Alloy and Lead Properties, Materials Compatibility, Thermal hydraulics and Technologies, ISBN 978-92-64-99002-9. Ott, K.O., Meneley, D.A., 1969. Accuracy of the quasistatic treatment of spatial reactor kinetics. Nucl. Sci. Eng. 36, 402–411. Price, M.S.T., 2012. The dragon project origins, achievements and legacies. Nucl. Eng. Des. 251, 60–68. Saad, Y., 1996. Iterative Methods for Sparse Linear Systems. PWS Publishing, New York. Sinha, Ratan Kumar, 2011. Advanced nuclear reactor systems – An Indian perspective. Energy Procedia 7, 34–50. Sinha, R.K., Dulera, I.V., 2010. Carbon based materials- applications in high temperature reactors. Indian J. Eng. Mater. Sci. 17, 321–326. Stacey, W.M., 2001. Nuclear Reactor Physics. Wiley, New York. Todreas, N.E., Kazimi, M.S., 2012. Nuclear Systems, Vol-1: Thermal Hydraulic Fundamental. CRC Press, Taylor & Francis Group. Van der Vorst, H., 2002. Iterative Krylov Methods for Large Linear Systems. Cambridge Press. Wu, H., Rizwan-uddin, 2016. A tightly coupled scheme for neutronics and thermal– hydraulics using open-source software. Ann. Nucl. Energy 87, 16–22. Zhang, Wenwen et al., 2016. Transient thermal–hydraulic analysis of a space thermionic reactor. Ann. Nucl. Energy 89 (2016), 38–49. Xu, Y., Downar, T., Walls, R., Ivanov, K., Staudenmeier, J., March-Lueba, J., 2009. Application of TRACE/PARCS to BWR stability analysis. Ann. Nucl. Energy 36, 317–323.