Research Article 
Corresponding author: Suhail Ahmad Khan ( suhailak@barc.gov.in ) Academic editor: Georgy Tikhomirov
© 2023 Suhail Ahmad Khan, Umasankari Kannan.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Khan SA, Kannan U (2023) Evaluation of DNBR with neutronics calculation in LWR systems. Nuclear Energy and Technology 9(2): 7783. https://doi.org/10.3897/nucet.9.98452

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 VISWAMTRIHEXFA 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 W3 Tong and OKBGidropress correlations implemented in TRIHEXFA.
DNBR, Heat Flux, TRIHEXFA, W3 Tong, Indian PWR
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) (
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 W3 (
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 MW_{th} 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, W3 and OKBGidropress 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.
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
The IPWR is a pressurized light water reactor with installed power of 2700 MW_{Th}/900 MW_{e}. The IPWR design aims to demonstrate the indigenous capability in the design, safety and operational aspects of the Pressurized Water Reactor technology (
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
Rated Power MW_{th}  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/m^{3}  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 (m^{3}/hr)  76700 
Reactivity control  Soluble boron(H_{3}BO_{3} in water) 
Shutdown and Control  Rod clusters in fuel assembly 
Control rod material  B_{4}C and Dy_{2}O_{3}.TiO_{2} 
A lattice burnup code VISWAM has been developed to cater to the current and future requirements of reactor core design computations (
One model is based on a combination of 1D multi group transport and 2D 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 1D transport calculations in multigroups using the first flight collision probability method (P_{ij}). The square or hexagonal cell boundary is cylindricalised to allow 1D treatment of the WignerSeitz cell. Heterogeneities present in the fuel assembly are treated using appropriate 1D supercell simulations. The pincell homogenized cross sections are collapsed to few groups using appropriate supercell spectra. For nonfuel cells, the few group cross sections are obtained from the respective supercell calculation of a given heterogeneity. The fuel assembly cell is treated by 2D few group diffusion theory using centremesh 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 N_{V} zones and N_{S} surfaces reduces to linear flux and current equations (assuming flat flux approximation)
$\Sigma j{V}_{j}{\varphi}_{j}=\sum _{\alpha =1}^{{N}_{S}}\sum _{\nu =0}^{{N}_{\nu}}{P}_{j\alpha}^{\nu}{S}_{\alpha}{J}_{,\alpha}^{\nu}+\sum _{i=1}^{{N}_{V}}{P}_{ji}{q}_{i}$. (1)
${S}_{\alpha}{J}_{+,\alpha}^{\nu}=\sum _{\beta =1}^{{N}_{S}}\sum _{\mu =0}^{{N}_{\mu}}{P}_{\alpha \beta}^{\nu \mu}{J}_{,\beta}^{\mu}{S}_{\beta}+\sum _{i=1}^{{N}_{V}}{p}_{\alpha i}^{\nu}{q}_{i}$. (2)
${J}_{,\alpha}^{\nu}=\sum _{\beta =1}^{{N}_{S}}\sum _{\mu =0}^{{N}_{\mu}}{A}_{\alpha \beta}^{\nu \mu}{J}_{+,\beta}^{\mu}$. (3)
The summation over ν in above equations represents the order of expansion of angular flux at pincell boundary. Here q_{i} = Σ_{si}V_{i}ϕ_{i} + S_{i}V_{i} is the total source in region i, Σ_{si} is the self scattering cross section within the group and S_{i} is the fission and scattering source in a group given by
${S}_{i}=\left(\sum _{\underset{{g}^{\u2019}\ne g}{{g}^{\u2019}=1}}^{G}{\Sigma}_{si}^{{g}^{\u2019}\to g}{\varphi}_{i}^{{g}^{\u2019}}+\frac{{\chi}_{g}}{k}\sum _{{g}^{\u2019}=1}^{G}\nu {\Sigma}_{fi}^{{g}^{\u2019}}{\varphi}_{i}^{{g}^{\u2019}}\right)$. (4)
Here P_{ji} gives the probability of a neutron emitted uniformly and isotropically in region i and having its first collision in region j and P^{ν}_{jα} 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 4^{th} order RungeKutta 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.
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 (
T_{c} (k) = T_{c} (k − 1) + ΔT (k). (5)
where ΔT (k) for k^{th} mesh is calculated using heat balance in the following expression.
$\Delta T\left(k\right)=\frac{PMESH*PTOCOOL*POW\left(k\right)}{MeshFlow*{C}_{p}*{\rho}_{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 (m^{3}/hr)
C_{p} – Specific heat (joules/g/°C)
ρ_{in} – density of inlet water in kg/m^{3}.
The fuel temperature for k^{th} fuel mesh is obtained as
T_{f} (k) = T_{c} (k) + (T_{rated} − T_{c} (k)) * POW (k). (7)
where T_{rated} 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 (
${\rho}_{\text{pert.}}=\frac{\Delta k}{k}=\frac{\int dr{\sum}_{g=1}^{G}{\varphi}_{g}^{\u2020}\left[\nabla \xb7\Delta {D}^{g}\nabla {\varphi}_{0}^{g}\Delta {\Sigma}_{t}^{g}{\varphi}_{0}^{g}+{\sum}_{{g}^{\u2019}=1}^{\mathrm{G}}\Delta {\Sigma}^{{g}^{\u2019}\to {g}^{\u2019}}{\varphi}_{0}^{{g}^{\u2019}}+\frac{{\chi}^{\mathrm{g}}}{{k}_{0}}{\sum}_{{g}^{\u2019}=1}^{\mathrm{G}}\Delta \left(\nu {\Sigma}_{f0}^{{g}^{\u2019}}\right){\varphi}_{0}^{{g}^{\u2019}}\right]}{\int dr{\sum}_{g=1}^{G}{\varphi}_{g}^{\u2020}{\chi}^{\mathrm{g}}{\sum}_{{g}^{\u2019}=1}^{\mathrm{G}}\nu {\Sigma}_{f0}^{{g}^{\u2019}}{\varphi}_{0}^{{g}^{\u2019}}}$. (8)
where the variables have their usual meaning. This method has been implemented in TRIHEXFA 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.
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: W3 correlation (
The W3 correlation, developed at Westinghouse by Tong (
For a channel with axially uniform heat flux, the correlation, giving the uniform critical heat flux. (q″_{cr}), is given by the following equation (
q″_{cr} = [(2.022 − 0.06238p) + (0.1722 − 0.01427p)e^{(18.177 − 0.5697p)xe}]
× [(0.1484 − 1.596x_{e} + 0.1729x_{e}x_{e})2.326G + 3271]
× [1.157 − 0.869x_{e}] × [0.2664+0.8357e^{−124.1Dℎ}]
× [0.8258 + 0.0003413(ℎ_{f} − ℎ_{in})]. (9)
where the symbols used and their valid ranges are as follows
q″_{cr} is in kW/m^{2}
p = Pressure in MPa (5.5 to 16)
x_{e} = Steam quality (–0.15 to 0.15)
G = Coolant mass flux in kg/m^{2}s (1356 to 6800)
D_{h} = 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 (
$h\left(z\right)={h}_{in}+\frac{1}{G}{\int}_{0}^{{l}_{c}}{q}^{\u2019\u2019}\left({z}^{\u2019}\right)\pi Dd{z}^{\u2019}$. (10)
where
l_{c} – Axial length of DNB location in nonuniformly heated channel
h (z) – Enthalpy at axial location z
The coolant quality in axial mesh k is calculated using following expression
$x\left(k\right)=\frac{{h}_{k}{h}_{f}}{{h}_{fg}}$. (11)
where h_{fg} is the enthalpy of evaporation in kJ/kg.
The axially nonuniform heat flux (q″_{cr}_{,n}) is obtained by applying a corrective F factor to the uniform critical heat flux
${q}_{cr,n}^{\u2019\u2019}=\frac{{q}_{cr}^{\u2019\u2019}}{F}$ (12)
where the factor F is given by
$F=\frac{C{\int}_{o}^{l}{q}^{\u2019\u2019}\left({z}^{\u2019}\right){e}^{c(l{z}^{\u2019})}d{z}^{\u2019}}{{q}^{\u2019\u2019}\left(l\right)[1{e}^{Cl}]}$ (13)
Here l is the distance extant to DNB as predicted by uniform q″_{cr} model and C is an experimental coefficient given by
$C=\frac{4.23\times {10}^{6}[1{x}_{e}(l){]}^{7.9}}{{G}^{1.72}}{m}^{1}$ (14)
The OKBGidropress correlation was obtained on the basis of results of experimental studies of DNB in the rod bundles for the conditions of VVER1000 type reactors and reported in the final safety analysis report of Bushehr Nuclear Power Plant (
q″_{cr}_{,n} = q″_{cr}_{,unif} F (15)
where q″_{cr}_{,unif} is the uniform critical heat flux distribution and is given by
q″_{cr}_{,unif} = 0.795(1 − x_{e})^{0.105p−0.5} × G^{0.184−0.311xe} × (1 − 0.0185p) (16)
and the factor F is given by
$F={\left[\frac{{\int}_{{L}_{CHF}}^{{L}_{CHF}}\Delta {L}^{q\left({z}^{\u2019}\right)d{z}^{\u2019}}}{q\left(z\right).\Delta L}\right]}^{n}$ (17)
$n=3.7919.61\frac{p}{{p}^{*}}+17.88{\left(\frac{T}{{T}^{*}}\right)}^{2}$ (18)
where ΔL = 55.D_{h} 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=\frac{{q}_{cr,n}^{\u2019\u2019}}{{q}^{\u2019\u2019}}$ (19)
Where
q″_{cr}_{,n} = Critical heat flux calculated as described above in each axial mesh (in kW/m^{2})
q″ = Prevailing heat flux in each axial mesh (in kW/m^{2}).
The minimum DNBR and its location in the axial mesh are searched from calculated 3D distribution of DNBR.
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 (
Fig.
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.
As seen from Figs
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 W3 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.
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.