Research Article
Print
Research Article
Evaluation of DNBR with neutronics calculation in LWR systems
expand article infoSuhail Ahmad Khan, Umasankari Kannan
‡ Bhabha Atomic Research Centre, Mumbai, India
Open Access

Abstract

The heat flux in a Light Water Reactor (LWR) system is used to estimate the Departure from Nucleate Boiling Ratio (DNBR) of the system which is an important thermal hydraulic parameter for nuclear reactors from heat removal point of view. The DNBR signifies an operational safety limit i.e. the nuclear power plant has to be operated with sufficient margin from the specified DNBR limit for assuring its safety. The DNBR is evaluated using a thermal hydraulic analysis code using inputs from neutronics calculation. The present paper presents the evaluation approach of minimum DNBR (MDNBR) during standard neutronics calculation. The DNBR calculation is performed using a core physics analysis code and burnup variation of MDNBR is studied for the full cycle length. The results of calculation are presented using the equilibrium core of 2700 MWth/900 MWe Indian Pressurized Water Reactor (IPWR). The calculations are performed using VISWAM-TRIHEXFA code system. The few group lattice parametric library for IPWR is generated by lattice analysis code VISWAM. The core follow up calculation for the equilibrium core configuration has been performed using core analysis code TRIHEXFA. A first order thermal hydraulic feedback model has been introduced into the 3D finite difference core simulation tool TRIHEXFA. The critical heat flux calculation, required for estimation of DNBR, has been performed using W-3 Tong and OKB-Gidropress correlations implemented in TRIHEXFA.

Keywords

DNBR, Heat Flux, TRIHEXFA, W-3 Tong, Indian PWR

Introduction

The power output of the reactor is limited by three dependent thermal and hydrodynamic variables: the DNBR, the maximum fuel temperature, and the core pressure drop (Δp) (Huda and Rahman 2004). The first design limitation is related to the DNB heat flux, which can be computed by means of a suitable correlation. When the heat flux become sufficiently large, the small bubble formed results in nucleate boiling and coalesces into a vapor film that covers the surface. Since the heat-transfer efficiency drops dramatically, the clad surface temperature will rise by several hundred degrees (Huda and Rahman 2004). This situation is termed as the DNBR and is defined as the ratio of the heat flux needed to cause departure from nucleate boiling to the actual local heat flux on the surface of a fuel rod. The core power distribution should be designed so as to prevent the DNBR from dropping below a chosen value under a high heat flux transient condition for the most adverse set of mechanical and coolant conditions. The DNBR evaluation is traditionally performed by thermal hydraulics analysis code using the neutronics input of core power distribution. The principal objective of this study is to evaluate the DNBR during the neutronics calculation. The calculations are performed using the VISWAM (Jagannathan et al. 2013a; Khan et al. 2016)-TRIHEXFA (Jagannathan et al. 2013b) code system.

The general practice is to estimate the critical heat flux using empirical correlations that have been validated by extensive experimental testing. Most of these correlations are proprietary and only limited information is available in literature. A methodology to calculate the critical heat flux using the published W-3 (Todreas and Kazimi 1990; Jagannathan et al. 2013b; Mathur et al. 2023) and OKB-Gidropress (Mozafari and Faghihi 2013) correlations has been incorporated into the 3D core analysis code TRIHEXFA. These correlations are traditionally popular with LWR designs. The results of the analysis are presented using the equilibrium core calculation of the Indian Pressurized Water Reactor (IPWR) (Raj et al. 2015).

The Indian Pressurized Water Reactor (IPWR) is being designed to develop indigenous capability for LWR technology. The design operating power for the IPWR is 2700 MWth with the targeted discharge burn up of 45 GWD/tU. Enriched uranium oxide is used as fuel and batch refueling scheme is adopted for fuel management and optimized to get a cycle length of 410 Full Power Days (FPDs). The fuel assembly of the IPWR uses gadolinium as burnable absorber for excess reactivity suppression. In the present analysis, W-3 and OKB-Gidropress correlations are deemed to be applicable for the hexagonal lattice geometry of the IPWR. The MDNBR for the equilibrium core of IPWR, as a function of fuel cycle burnup, has been evaluated using VISWAM–TRIHEXFA code system.

A brief description of the IPWR is given in Section 2. Section 3 gives the description of the code systems and correlations used. Section 4 gives the result of the analysis followed by the conclusions in Section 5.

List of acronyms

BAR Burnable Absorber Rod

BOC Beginning Of Cycle

DNBR Departure from Nucleate Boiling Ratio

EFPD Effective Full Power Days

FA Fuel Assembly

FPDs Full Power Days

IPWR Indian Pressurized Water Reactor

LWR Light Water Reactor

MDNBR Minimum DNBR

MWe Mega Watt Electrical

MWth Mega Watt Thermal

TH Thermal Hydraulic

VVER Vodo Vodyanoi Enyergeticheskiy Reaktor

Description of IPWR

The IPWR is a pressurized light water reactor with installed power of 2700 MWTh/900 MWe. The IPWR design aims to demonstrate the indigenous capability in the design, safety and operational aspects of the Pressurized Water Reactor technology (Raj et al. 2015). The IPWR fuel assembly has 331 lattice locations arranged in hexagonal geometry. The 331 locations are distributed into 311 fuel pin locations, one central water tube, 18 guide tubes for control rods and one instrumentation tube. The fuel pins in the IPWR fuel assembly have profiled enrichment distribution. The burnable absorber rods are used in the fuel assembly to manage the initial excess reactivity in the core.

The IPWR core consists of 151 fuel assemblies arranged in a hexagonal lattice array. The significant parameters of the equilibrium IPWR core are given in Table 1. Soluble boron in the form of boric acid in moderator is used for burn up and Xenon, Samarium reactivity management. The cycle length achieved is 410 effective full power days (EFPD). A three-batch mode of refueling scheme is used for fuel cycle management. The equilibrium core configuration is shown in Fig. 1. The 151 fuel assemblies are distributed as 51/50/50 for three cycles as shown in Fig. 1.

Table 1.

Design Parameters of the IPWR Equilibrium Core

Rated Power MWth 2700
No. of fuel assemblies (FA) 151
FA configuration Hex
Average Linear Heat Generation rate in a pin kw/m 15.96
Power density MW/m3 87.4
Core periphery dia m 3.3
Core baffle dia m 3.46
Material of the core baffle SS + water
System pressure MPa 15.7
Active core height m 3.6
Fuel Temperature °C 625
Coolant inlet temperature °C 292
Coolant outlet temperature °C 323
Coolant Flow (m3/hr) 76700
Reactivity control Soluble boron(H3BO3 in water)
Shutdown and Control Rod clusters in fuel assembly
Control rod material B4C and Dy2O3.TiO2
Figure 1. 

Equilibrium Core of IPWR.

Description of calculation tools

Lattice analysis code VISWAM

A lattice burnup code VISWAM has been developed to cater to the current and future requirements of reactor core design computations (Jagannathan et al. 2013a; Khan et al. 2016). Currently, there are two models available in VISWAM for fuel assembly burnup calculations.

One model is based on a combination of 1-D multi group transport and 2-D few group diffusion theory. A typical fuel assembly cell consists of fuel pins of different enrichments and various heterogeneous cells like control rod, water rod or burnable absorber rod (BAR) cells. In this model, the fuel pins are classified into various pin cell types based on the enrichment and Dancoff factors. The pin cells are treated using a series of 1-D transport calculations in multi-groups using the first flight collision probability method (Pij). The square or hexagonal cell boundary is cylindricalised to allow 1-D treatment of the Wigner-Seitz cell. Heterogeneities present in the fuel assembly are treated using appropriate 1-D supercell simulations. The pincell homogenized cross sections are collapsed to few groups using appropriate supercell spectra. For non-fuel cells, the few group cross sections are obtained from the respective supercell calculation of a given heterogeneity. The fuel assembly cell is treated by 2-D few group diffusion theory using centre-mesh finite difference method.

The second model uses interface current method based on 2D collision probability. In this model, a lattice cell may be a fuel pincell, water rod cell or an absorber rod cell. The geometry of the cell is not changed, i.e., the outermost region of the cell is retained as the square or hexagonal shape without cylindricalisation. The collision probabilities are calculated for single lattice cell in 2D geometry. For fuel assembly calculation, the lattice cells are linked using interface currents by using the double P2 (DP2) expansion of angular flux at the pincell boundary. In this method, the neutron transport equation for group ‘g’ (the group index ‘g’ is omitted for simplicity), when discretized over a region consisting of NV zones and NS surfaces reduces to linear flux and current equations (assuming flat flux approximation)

ΣjVjϕj=α=1NSν=0NνPjανSαJ-,αν+i=1NVPjiqi. (1)

SαJ+,αν=β=1NSμ=0NμPαβνμJ-,βμSβ+i=1NVpαiνqi. (2)

J-,αν=β=1NSμ=0NμAαβνμJ+,βμ. (3)

The summation over ν in above equations represents the order of expansion of angular flux at pincell boundary. Here qi = ΣsiViϕi + SiVi is the total source in region i, Σsi is the self scattering cross section within the group and Si is the fission and scattering source in a group given by

Si=g=1ggGΣsiggϕig+χgkg=1GνΣfigϕig. (4)

Here Pji gives the probability of a neutron emitted uniformly and isotropically in region i and having its first collision in region j and Pν gives the probability of neutron entering through surface α uniformly in mode ν and having first collision in region j. Pναi is the probability that neutrons emitted uniformly and isotropically in region i will escape through surface α in mode ν and Pνμαβ is the probability that neutrons entering through surface β uniformly in mode μ will be transmitted through the cell and out through surface α in mode ν without making a collision.

The fuel assembly cell calculation in VISWAM generates several parameters like the infinite multiplication factor ‘k, power distribution, homogenized and collapsed five or two group cross sections of the entire assembly cell, and kinetics parameters like delayed neutron fraction ‘β’ and prompt neutron mean life time ‘l’. Average energy yield per fission is also calculated as the weighted mean value from the fissions in all actinides. The fuel pins are categorized into several burnup zones. Depletion equations are solved in these burnup zones using spectrum averaged reaction cross sections to determine the evolution of isotopic densities as a function of burn up. The burnup equations are solved using 4th order Runge-Kutta method.

The results reported here have been obtained using the few group cross section lattice data generated using first model of 1D transport and 2D diffusion. This is because this model is computationally efficient and the few group lattice library with parametric variation of fuel temperature, coolant temperature, coolant density, Xenon and Samarium concentration is generated in reasonably small time.

Core analysis code TRIHEXFA

TRIHEXFA is a 3D diffusion theory code developed indigenously for performing the full 3D core calculation with burnup. It has been extensively validated against the many benchmarks and experimental data (Jagannathan et al. 2013b). TRIHEXFA has been made compatible to the few group parametric cross section library generated by lattice analysis code VISWAM. While performing the 3D core calculation the flux, power and corresponding temperature distributions should be iteratively determined for internal self-consistency. The few group parametric library is generated using lattice code VISWAM for homogenised fuel assembly in controlled/uncontrolled condition for a set of reactor state parameters of boron concentration, fuel temperature, coolant temperature, coolant density, saturated xenon concentration and samarium concentration. The value of these parameters pertaining to nominal operating condition is called reference state. In TRIHEX-FA the perturbation of cross sections due to difference in the local values of temperature, Xenon and Samarium from the nominal values is updated periodically (Jagannathan et al. 2013b). These perturbations are applied to all the cross sections in five groups. The ratio of the five group cross section values for the perturbed state of a parameter with respect to the five group cross sections of its reference value is evaluated at the input grid values. The local coolant temperature for kth fuel mesh in axial direction is obtained as

Tc (k) = Tc (k − 1) + ΔT (k). (5)

where ΔT (k) for kth mesh is calculated using heat balance in the following expression.

ΔT(k)=PMESH*PTOCOOL*POW(k)MeshFlow*Cp*ρin*1000.0. (6)

where:

PMESH – Average Power produced in single mesh (w)

PTOCOOL – power to coolant (0.97)

POW (k) – relative power of the mesh

MeshFlow – Flow in single mesh area (m3/hr)

Cp – Specific heat (joules/g/°C)

ρin – density of inlet water in kg/m3.

The fuel temperature for kth fuel mesh is obtained as

Tf (k) = Tc (k) + (TratedTc (k)) * POW (k). (7)

where Trated is the average fuel temperature at nominal power.

After calculating local parameters of the mesh, the cross section ratios for the mesh are obtained using interpolation for that parameter’s local value in the mesh. The perturbed cross section for a mesh is then obtained by multiplying the reference few group cross sections with the ratios for each type of perturbation, viz., fuel and coolant temperature, xenon and samarium loads.

A model to estimate reactivity coefficients using first order perturbation theory has been implemented in TRIHEXFA. This model is based on the method described in Stacey (Stacey 2007). In this method, it is assumed that for small perturbations the direct and adjoint flux of the reference and perturbed state are nearly same. With this approximation, the reactivity change due to perturbation of any of the parmeters listed aboved can be determined from the expression:

ρpert. =Δkk=drg=1Gϕg·ΔDgϕ0g-ΔΣtgϕ0g+g=1GΔΣggϕ0g+χgk0g=1GΔνΣf0gϕ0gdrg=1Gϕgχgg=1GνΣf0gϕ0g. (8)

where the variables have their usual meaning. This method has been implemented in TRIHEX-FA code. The direct and adjoint fluxes were always caculated during core follow up simulations for estimating core kinetics parameters viz., effective delayed neutron fraction ‘βeff’ and prompt neutron life time ‘l’. These flux distributions were used to estimate each type of reactivity coefficients. There is no necessity to perform separate 3D simulations for each type of perturbed states, as existed in the old ‘2K’ approach. The flow chart for neutronics and thermal hydraulic (TH) calculation in TRIHEXFA is shown in Fig. 2.

Figure 2. 

Flow Chart of Neutronics and TH Calculations in TRIHEXFA.

DNBR calculation methodology

For calculating DNBR, critical heat flux needs to be evaluated. The heat fluxes are determined from the neutronic power at each and every defined volume of the region in the core. The ratio between the predicted heat flux and the actual operating heat flux is called the departure from nucleate boiling ratio (DNBR). This ratio changes across the length of fuel bundle and reaches a minimum value called MDNBR and is one of the important thermal hydraulic safety parameter which needs to be evaluated for safe operation of PWR reactors. The critical heat flux is calculated using the look up tables or empirical correlations developed for various reactor types having different operating parameters and flow channel geometry. We have used two correlations for calculating the MDNBR margin for equilibrium core of IPWR: W-3 correlation (Todreas and Kazimi 1990) and the OKB correlation (Mozafari and Faghihi 2013).

The W-3 correlation, developed at Westinghouse by Tong (Todreas and Kazimi 1990; Mathur et al. 2023), is developed for axially uniform heat flux, with a correcting factor for non-uniform flux distribution (Todreas and Kazimi 1990). For predicting the DNB condition in a non-uniform heat flux channel, the following two steps are to be followed:

  1. The uniform critical heat flux. (qcr) is computed with the W-3 correlation, using the local reactor conditions.
  2. The non-uniform DNB heat flux qcr , n) distribution is then obtained (assuming a flux shape similar to that of the reactor) by dividing qcr by the F factor.

For a channel with axially uniform heat flux, the correlation, giving the uniform critical heat flux. (qcr), is given by the following equation (Todreas and Kazimi1990; Mathur et al. 2023)

qcr = [(2.022 − 0.06238p) + (0.1722 − 0.01427p)e(18.177 − 0.5697p)xe]

× [(0.1484 − 1.596xe + 0.1729xe|xe|)2.326G + 3271]

× [1.157 − 0.869xe] × [0.2664+0.8357e−124.1Dℎ]

× [0.8258 + 0.0003413(fin)]. (9)

where the symbols used and their valid ranges are as follows

qcr is in kW/m2

p = Pressure in MPa (5.5 to 16)

xe = Steam quality (–0.15 to 0.15)

G = Coolant mass flux in kg/m2s (1356 to 6800)

Dh = equivalent heated diameter in m (0.005 to 0.018)

f = Saturated liquid enthalpy in kJ/kg

in = Inlet enthalpy in kJ/kg (>930)

The enthalpy at DNB location is calculated using following expression (Mathur et al. 2023)

h(z)=hin+1G0lcq(z)πDdz. (10)

where

lc – Axial length of DNB location in non-uniformly heated channel

h (z) – Enthalpy at axial location z

The coolant quality in axial mesh k is calculated using following expression

x(k)=hk-hfhfg. (11)

where hfg is the enthalpy of evaporation in kJ/kg.

The axially non-uniform heat flux (qcr,n) is obtained by applying a corrective F factor to the uniform critical heat flux

qcr,n=qcrF (12)

where the factor F is given by

F=Colq(z)e-c(l-z)dzq(l)[1-e-Cl] (13)

Here l is the distance extant to DNB as predicted by uniform qcr model and C is an experimental coefficient given by

C=4.23×106[1-xe(l)]7.9G1.72m-1 (14)

The OKB-Gidropress correlation was obtained on the basis of results of experimental studies of DNB in the rod bundles for the conditions of VVER-1000 type reactors and reported in the final safety analysis report of Bushehr Nuclear Power Plant (Mozafari and Faghihi 2013). Experimental data on critical heat fluxes were obtained both under uniform and non-uniform distribution of heat flux along the bundle length on the basis of experimental data on studying DNB on 7-rod bundles and 19-rod bundles. The non-uniform critical heat flux distribution (qcr,n) is determined by the following correlation:

qcr,n = qcr,unif F (15)

where qcr,unif is the uniform critical heat flux distribution and is given by

qcr,unif = 0.795(1 − xe)0.105p−0.5 × G0.184−0.311xe × (1 − 0.0185p) (16)

and the factor F is given by

F=LCHFLCHF-ΔLq(z)dzq(z).ΔLn (17)

n=3.79-19.61pp*+17.88TT*2 (18)

where ΔL = 55.Dh is the relaxation length in m, p* = 22.1 Mpa is the critical pressure and T* is the critical temperature at p*.

The 3D DNBR distribution in the axial and radial meshes of equilibrium core of IPWR is calculated as

DNBR=qcr,nq (19)

Where

qcr,n = Critical heat flux calculated as described above in each axial mesh (in kW/m2)

q″ = Prevailing heat flux in each axial mesh (in kW/m2).

The minimum DNBR and its location in the axial mesh are searched from calculated 3D distribution of DNBR.

Results and discussion

In order to evaluate the DNBR, the equilibrium core configuration is obtained by performing the successive core follow up calculations using a predetermined loading pattern. This calculation is performed by considering the space dependent feedbacks due to fuel temperature, coolant temperature, saturated Xenon and Samarium (Khan and Kannan 2017). Figs 3, 4 show a typical volume averaged axial variation of the fuel and coolant temperature in the fuel assembly. Due to the axial variation of these temperatures, it is necessary to perform the equilibrium core calculation with the space dependent feedbacks. The DNBR is evaluated during the equilibrium core calculation.

Figure 3. 

Typical Axial Fuel Temperature Profile.

Figure 4. 

Typical Axial Coolant Temperature Profile.

Fig. 5 gives variation of minimum DNBR obtained using W-3 and OKB Gidropress correlations with cycle burn up for equilibrium core of IPWR. It is seen that the minimum DNBR value obtained using W-3 correlation stays close to 2.0 throughout the cycle length under nominal operating conditions. This value of DNBR is higher than the internationally accepted safety limit of 1.3 for the pressure water reactors (Hortal 2003) implying a safe operation during the cycle operation. The minimum DNBR value predicted using OKB correlations is greater than 2.1 throughout the cycle.

Figure 5. 

Variation of MDNBR with Cycle Burnup.

The DNBR, for a given set of operating parameters, varies along the length of the flow channel. At the coolant entrance, it is high indicating low heat flux, and correspondingly high critical heat flux, thus a system is well within safety limits. As fuel temperature increases, the DNBR decreases until it reaches a minimum. Fig. 6 gives the axial variation of DNBR in the hottest channel at beginning of cycle (BOC). As seen from Fig. 7, which gives the radial power peaking distribution at BOC along with fuel assembly numbers, the hottest channel is 137. The minimum DNBR is seen at axial mesh 15 in channel 137 using both the correlations. Fig. 6 shows a smooth axial variation of DNBR and is as per the expected physical phenomenon.

Figure 6. 

Axial Variation of DNBR in hottest Channel.

Figure 7. 

Radial Power Distribution at BOC.

As seen from Figs 5, 6, there is a difference between the DNBR calculated from the two correlations. Since the inputs such as relative power of channels and dimensional inputs are same for the two correlations, this difference can be attributed to different critical heat fluxes and axial heat flux distributions calculated using two correlations. It is seen that DNBR values in bottom and top axial meshes are comparatively large in Fig. 6 using both the correlations. This is due to the non validity of correlations in these meshes.

Conclusion

The variation of minimum DNBR during core cycle has been evaluated for the equilibrium core of IPWR. The calculation is performed using the VISWAM–TRIHEXFA code system. The thermal hydraulic criterion is being evaluated along with neutronic simulation in a internally coupled formalism. The equilibrium core configuration is obtained with multiple core follow up calculations with the space dependent feedbacks due to fuel temperature, coolant temperature, saturated Xenon and Samarium. To evaluate the minimum DNBR, the critical heat flux is calculated using the W-3 and OKB Gidropress correlations available in the TRIHEXFA code system. The minimum DNBR is seen to be of the order 2.0 throughout the cycle duration under nominal operating conditions which is within the accepted limit for safe operation.

Acknowledgements

The first author is thankful to Dr. V. Jagannathan and Arvind Mathur (Both Ex RPDD, BARC) for codeveloping and valuable discussions. The authors are thankful to Shri Gopal Mapdar of RPDD, BARC for plotting some figures appearing in the present paper.

References

  • Hortal J (2003) Safety Margins: Deterministic And Probabilistic Views. In: Safety margins of operating reactors–Analysis of uncertainties and implications for decision making, IAEA–TECDOC–1332 (January 2003).
  • Jagannathan V, Khan SA, Mathur A (2013a) VISWAM–A Modular Code System for State of Art Analysis of Nuclear Reactors–Developmental Ideas and Current Status”. Proc. of ARP 2013–Advance in Reactor Physics: Simulation Techniques and Analysis Methodologies, BARC Training School Hostel, Anushaktinagar, Mumbai, India, October 23–25, 2013.
  • Jagannathan V, Khan SA, Mathur A (2013b) Recent Advances in 3D Core Diffusion Code TRIHEX–FA. Proc. of ARP 2013–Advance in Reactor Physics: Simulation Techniques and Analysis Methodologies, BARC Training School Hostel, Anushaktinagar, Mumbai, India, October 23–25, 2013.
  • Khan SA, Kannan U (2017) Equilibrium Core calculation of IPWR with space dependent feedbacks. RPDD Note No. RPDD/IPWR/41/2017, March 30, 2017.
  • Khan SA, Mathur A, Jagannathan V (2016) Incorporation of Interface Current Method Based on 2D CP Approach in VISWAM Code System for Hexagonal Geometry. Annals of Nuclear Energy 92: 161. https://doi.org/10.1016/j.anucene.2016.01.041
  • Mathur AK, Jagannathan V, Khan SA, Singh S (2023 [under publication]) Calculation of DNBR Margin in TRIHEXFA code using W-3 Correlation. Report under publication.
  • Mozafari MA, Faghihi F (2013) Design of annular fuels for a typical VVER-1000 core: Neutronic investigation, pitch optimization and MDNBR calculation. Annals of Nuclear Energy 60: 226. https://doi.org/10.1016/j.anucene.2013.04.035
  • Raj D, Sarkar A, Mapdar G, Kannan U (2015) Technical Design report of IPWR : Reactor Physics chapter. RPDD Internal Report, November 23, 2015.
  • Todreas NE, Kazimi MS (1990) Nuclear Systems I – Thermal Hydraulic Fundamentals. Taylor & Francis, New York.
login to comment