Development and Verification of an Efficient Spatial Neutron Kinetics Method for Reactivity-Initiated Event Analyses

A space/time nodal diffusion code based on the nodal expansion method (NEM), EPISODE, was developed in order to evaluate transient neutron behavior in light water reactor cores. The present code employs the improved quasi-static (IQS) method for spatial neutron kinetics, and neutron flux distribution is numerically obtained by solving the neutron diffusion equation with the nonlinear iteration scheme to achieve fast computation. A predictor-corrector (PC) method developed in the present study enabled to apply a coarse time mesh to the transient spatial neutron calculation than that applicable in the conventional IQS model, which improved computational efficiency further. Its computational advantage was demonstrated by applying to the numerical benchmark problems that simulate reactivity-initiated events, showing reduction of computational times up to a factor of three than the conventional IQS. The thermohydraulics model was also incorporated in EPISODE, and the capability of realistic reactivity event analyses was verified using the SPERT-III/E-Core experimental data.     

核动力装置自然循环及其过渡过程计算模型的建立

为准确分析含反应性反馈的核动力装置自然循环及其过渡过程中重要参数的响应特性,以核动力装置瞬态最佳估算程序RELAP5/MOD3为基础,采用两群三维时空中子动力学模型替代RELAP5/MOD3的点堆模型,并建立三维空间内中子物理与热工水力的耦合模型,编制相应的计算程序。利用所研制的程序对实际核动力装置的自然循环及其过渡过程进行分析计算,并与试验结果进行比较。结果表明:本文建立的时空中子动力学计算模型克服了点堆方程不能准确计算反应性反馈的缺点,计算精度高,研制的程序可作为核动力装置强迫循环与自然循环及其过渡过程的计算分析工具。     

Fractional neutron point kinetics equations for nuclear reactor dynamics

The fractional point-neutron kinetics model for the dynamic behavior in a nuclear reactor is derived and analyzed in this paper. The fractional model retains the main dynamic characteristics of the neutron motion in which the relaxation time associated with a rapid variation in the neutron flux contains a fractional order, acting as exponent of the relaxation time, to obtain the best representation of a nuclear reactor dynamics. The physical interpretation of the fractional order is related with non-Fickian effects from the neutron diffusion equation point of view. The numerical approximation to the solution of the fractional neutron point kinetics model, which can be represented as a multi-term high-order linear fractional differential equation, is calculated by reducing the problem to a system of ordinary and fractional differential equations. The numerical stability of the fractional scheme is investigated in this work. Results for neutron dynamic behavior for both positive and negative reactivity and for different values of fractional order are shown and compared with the classic neutron point kinetic equations. Additionally, a related review with the neutron point kinetics equations is presented, which encompasses papers written in English about this research topic (as well as some books and technical reports) published since 1940 up to 2010.     

A two-fluid code for the thermohydraulic transient analysis of CICC superconducting magnets

We present here a finite element computer model (Mithrandir) for the transient thermohydraulics of compressible helium in a Cable-In-Conduit Conductor (CICC) with central cooling hole, as presently envisaged for superconducting magnets of the International Thermonuclear Experimental Reactor (ITER). In the model the He in the hole and that in the cable bundle are treated as separate fluids, each characterized by its own flow and thermodynamic properties, coupled by exchanges of mass, momentum and energy. Results for the simulation of a quench both with and without a wall delimiting the central cooling hole are discussed. Time and space convergence of the code are demonstrated numerically.     

Numerical Solution to Space-Angle Energy-Dependent Neutron Integral Transport Equation

A numerical approach to the steady-state, space-, angle- and energy-dependent neutron transport equation is presented for neutron shielding calculations. The scattering integral, with anisotropic treatment of elastic scattering and isotropic treatment of inelastic scattering, is evaluated by the use of Gaussian and straightforward quadratures. A system of coupled one-group integral equations for all the energy meshes of interest, converted from the energy-dependent integral transport equation, is calculated by performing a line integration along the neutron path in the direction of motion. For this purpose the direction of neutron motion is represented by discrete-ordinate directions Ω_pq on the unit sphere.

The final presentation of the integral transport equation is derived in a difference form convenient for machine computation. A computation program PALLAS has been written in Fortran IV for IBM 360-75 computer to perform neutron transport calculations based on this approach.

Comparisons are given of the numerical solutions with analytical solutions for unscattered fluxes in various geometries such as plane, spherical and two-dimensional cylindrical, for volume sources with self absorption, and with experimental spectra for angular neutron fluxes in graphite-, polyethylene- and water-shield. Excellent agreement is obtained between the present calculations and analytical or experimental results.     

Some dynamic A-n,S-2n analytic calculations

A-n, S-2n dynamic equations in neutron transport theory are given an analytic solution, suitable for numerical safety code validation. Some numerical results are presented for plane geometry with periodic boundary conditions, within the monokinetic isotropic scattering case. They turn out to be very interesting and physically well interpretable and seem to give a deep insight into the transport effects connected with the finite velocity space migration of localized neutron signals. Balance equations are solved by means of a space Helmholtz series expansion connected to a time-variable Laplace transformation technique. The equations that are to be solved for the considered slab geometry, although perfectly equivalent to S-2n equations, are actually the ones corresponding to the typical second-order A-n model. The given fully analytic solution can be used to yield standard reference results for peculiar space transients, against which any numerical safety code exploiting the same discrete ordinate model can be validated.Dipartimento di Energetica, Politecnico di Torino, Corso Duca degli Abruzzi, 24-10129 Torino, Italy. Published in Atomnaya Énergiya, Vol. 62, No. 5, pp. 320–326, May, 1987.     

Adaptive matrix formation (AMF) method of space–time multigroup reactor kinetics equations in multidimensional model

Adaptive matrix formation (AMF) method has been developed for the numerical solution of the transient multigroup neutron diffusion and delayed precursor equations in two- and three-dimensional geometry. The method is applied to a general class of two- and three- dimensional problems. The results of numerical experiments, as well as comparison with space–time experimental results indicate that the method is accurate and that the two- and three-dimensional calculations can be performed at “reasonable” computer costs. Moreover, the AMF method offers the flexibility of using smaller time steps between flux shape calculations to achieve a specified accuracy and capability, without encountering numerical problems that occur in the other conventional methods. There is a large considerable saving in computer time and costs due to the partitioning of the matrix adopted in the presented AMF method. The two- and three-dimensional problems were analyzed with the present calculations model to illustrate the accuracy and stability of the method. Furthermore, the stability of the investigated method has been tested for sinusoidal, ramp, and step-change reactivity insertions. The results are in a good agreement with those of the other less approximate methods, including the problems in which the reflector zone is perturbed.     

临界反应堆阶跃正反应性输入时中子密度响应的近似修正解

通过修正单组缓发中子先驱核衰变常量λ值，使点堆中子动力学方程单组缓发中子模型在正反应性阶跃输入时的数值计算结果趋近于六组缓发中子模型数值计算结果。在此基础上，用修正后的单组模型解析方法进行计算。结果表明：采用修正后的单组解析方法计算阶跃正反应性输入的中子密度响应，计算结果与六组的接近，满足工程计算精度要求，同时计算简便，避免了刚性问题，可以实现快速计算。     

Numerical methods for the generation of the spectrum of the multigroup slab-geometry discrete ordinates operator in neutron transport theory

We describe two numerical methods applied to the first-order form of the multigroup slab-geometry discrete ordinates equations modelling fixed-source neutron transport problems with anisotropic scattering. The numerical methods described in this article generate the spectrum and a vector basis for the null space of the multigroup slab-geometry discrete ordinates operator defined in a homogeneous domain. The first method is a more general approach of a numerical method described in a recent work by others. We then come to consider numerical and computational aspects of the first method and we propose a second method. The second method is a multigroup extension of a numerical method described in a more recent work by the present author. In order to provide those interested in implementing either method with a reference set, we present numerical results for some multigroup slab-geometry model problems with anisotropic scattering. We conclude this article with a discussion and directions for future work.     

Reactor core thermohydraulic transients with thermohydrodynamic coupling

A study of the reactor core thermohydraulics in an LMFBR has been performed for the strongly coupled thermo-hydrodynamic transients. A numerical method to solve the coupled energy-momentum equations among multichannels in a core is presented and the computer code ORIFS-TRANSIENT has been developed.The results of sample calculations for a flow coastdown transient to natural circulation following a reactor scram in a typical loop-type LMFBR are as follows: (1) the inter-subassembly coolant flow redistribution due to buoyancy forces is significant under the low flow condition, such as natural circulation; (2) the maximum coolant temperature was decreased by about 80°C (corresponding to about 22% in terms of hot channel factor) due to the flow redistribution; (3) due to thermohydrodynamic coupling between upper plenum and other regions, the maximum coolant temperature was decreased by about 9°C; (4) due to inter-subassembly heat redistribution, the maximum coolant temperature was increased by about 7°C.     

解析基函数展开方法求解二维三角形几何中子扩散方程

提出了一种在二维三角形几何内数值求解中子扩散方程的节块方法.节块内的各群通量分布用解析基函数近似展开,节块之间采用面偏流零次矩和一次矩进行耦合;给出了三角形几何下的节块扫描方案;采用响应矩阵技术进行迭代求解,开发了二维三角形组件中子扩散计算程序ABFEM-T.通过基准问题的校验计算,表明该方法能准确地给出有效增值系数及节块功率分布,可用于复杂的非结构几何区域的中子扩散问题的求解.     

非结构网格中子输运方程的球谐函数解法研究

从新的二阶自共扼角通量密度(SAAF：Self-Adjoint Angular Flux)中子输运方程出发．利用球谐函数对角度变量进行展开,导出了一组关于空间变量的偏微分方程组,中子通量密度的各个分量相互耦合,应用一定的迭代策略进行迭代求解。针对每一个方程,应用有限元方法对非结构网格进行离散求解。据此编写了二维球谐函数方法输运计算程序,对一系列基准题进行校算的数值结果表明,该方法具有较高的计算精度,克服了射线效应,并能用于非结构网格。     

Spherical harmonics method for neutron transport equation based on unstructured-meshes

Based on a new second-order neutron transport equation, self-adjoint angular flux (SAAF) equation, the spherical harmonics (PN) method for neutron transport equation on unstructured-meshes is derived. The spherical harmonics function is used to expand the angular flux. A set of differential equations about the spatial variable, which are coupled with each other, can be obtained. They are solved iteratively by using the finite element method on unstructured-meshes. A two-dimension transport calculation program is coded according to the model. The numerical results of some benchmark problems demonstrate that this method can give high precision results and avoid the ray effect very well.     

Safety analysis of Ulchin Nuclear Power Plant against Nihonkai-Chubu Earthquake Tsunami

A coupled numerical model based on the shallow-water equations is employed to simulate the Nihonkai-Chubu Earthquake Tsunami event which occurred in 1983. The model consists of a transoceanic propagation and an inundation model both described by finite difference methods. A special moving boundary treatment is implemented in the inundation model to track accurately a transient motion of shoreline. Maximum run-up heights along the Eastern Coastline of the Korean Peninsula are predicted and compared to field measurements. In particular, the maximum run-up height at the Ulchin Nuclear Power Plant is estimated and compared to the observed height to analyze the safety of the Plant. Although the numerical model slightly overestimates the maximum run-up heights, the model still provides reasonable predictions.     

A variational treatment for the time dependent boltzmann equation as a basis for numerical solutions conserving neutrons

A maximum principle for the time-dependent first-order Boltzmann equation is established in two independent ways:- by a generalized least squares method and by a method based on the properties of an appropriate bi-linear form. The second derivation suggests a metric for a Hilbert space which provides a geometrical interpretation of the variational principle. This interpretation leads to a Petrov-Galerkin weighted residual method in contrast to the Galerkin method of Martin for time dependent transport.

The maximum principle is used to define a figure of merit for the global error of any numerical solution for time dependent transport. The principle is used also to demonstrate the neutron conservation property of optimized numerical solutions, and the convergence of finite element methods based on the variational principle.

A direct use of the maximum principle give a conservative three level scheme for transients in the angular flux. A less massive calculation for a conservative solution proceeds in two stages. The first stage obtains a preliminary solution based on a sequence of tailored steady state calculations for the even-parity angular flux. The relevant equations are derived by making residual terms vanish in the variational principle. They are equivalent to the finite element - finite difference equations used by de Oliveira and Wood for the analysis of oil-well logging by means of a neutron pulse. For the second stage the preliminary solution is weighted at each time step. The weights are determined by the maximum principle to yield a conservative solution which is continuous in time. The solution can be arranged to have also a continuous time derivative. Thus the maximum principle can be used both as a way of predicting and correcting solutions.     

Cellular neural networks (CNN) simulation for the TN approximation of the time dependent neutron transport equation in slab geometry

This paper describes the application of a multilayer cellular neural network (CNN) to model and solve the time dependent one-speed neutron transport equation in slab geometry. We use a neutron angular flux in terms of the Chebyshev polynomials (T_N) of the first kind and then we attempt to implement the equations in an equivalent electrical circuit. We apply this equivalent circuit to analyze the T_N moments equation in a uniform finite slab using Marshak type vacuum boundary condition. The validity of the CNN results is evaluated with numerical solution of the steady state T_N moments equations by MATLAB. Steady state, as well as transient simulations, shows a very good comparison between the two methods. We used our CNN model to simulate space–time response of total flux and its moments for various c (where c is the mean number of secondary neutrons per collision).     

Power series solution method for solving point kinetics equations with lumped model temperature and feedback

Point kinetics equations are stiff differential equations, and their solution by the conventional explicit methods will give a stable consistent result only for very small time steps. Since the neutron lifetime in a LMFBR is very short, the point kinetics equations for LMFBRs become even stiffer. In this study the power series solution (PWS) method is applied for solving the point kinetics equations for a typical LMFBR. A Fortran program is developed for accident analysis of LMFBRs with the PWS method for solving the point kinetics and a lumped model for solving the heat transfer equations. A new technique is developed with fixing factor to find out the average temperature at the peak power node (PPN) without performing temperature calculations at all axial nodes in a reactor fuel pin. The temperature at PPN also decides whether the reactor is within the design safety limit (DSL) or it has entered a serious transient that may lead to an accident. The coupled heat transfer and point kinetics models for a peak power node give the average fuel, clad and coolant temperatures. For the transient over power accidents (TOPA), this is the best way for calculating the temperature, with minimum amount of computations. TOPA analyses are carried out with PWS method. It is found that the PWS methodology uses a small number of numerical operations, while the computational time and the accuracy are comparable with the available fast computational tools. This methodology can be used in nuclear reactor simulation studies and accident analysis.     

Reference model parameter identification of space–time dependent reactivity in a CANDU-PHWR

The problem of estimating reactivity transients from an observed neutron flux transient is considered. This is relevant, for example, to analyzing a power rundown test or to estimating reactivity variations associated with some computer codes that do not specifically compute individual reactivity components. A method is presented which utilizes inverse space–time kinetics and optimal state estimators to extract the components of the reactivity transient from observed neutron flux measurements. The approach takes into account geometric characteristics and composition of the reactor core, as well as reactor operating conditions. Measurements from a limited number of in-core neutron flux detectors are the inputs used to extract reactivity components that fit a modal model of the reactor, referred to as the “reference model”. An improved solution for the reactivity components is then generated using the modal approximation solution for the neutron transport equation in conjunction with optimal estimation techniques. The method has been applied to a reactivity initiated accident in which a transient is initiated by a non-uniform loss-of-coolant. This results in a spatially varying neutron overpower transient that is terminated by the asymmetric insertion of shutoff rods. In this paper the Joint Extended Kalman Filter and Rauch–Tung–Striebel smoother is employed to estimate the neutron flux distribution in the core and identify the reactivity components of the reference model. The reference model in the state space and the Kalman filter algorithm are shown. Results of numerical simulations of the reactor transient and the optimal estimation of the reactivity components are presented to demonstrate the capabilities of the method.     

Coupled dynamics in the physics of molten salt reactors

The present paper is devoted to the analysis of the coupled thermo-fluid and neutronic dynamics of fast fluid-fuel multiplying nuclear systems. A completely coupled model is needed since in some fast reactors designs, the velocity pattern could be very complicated and strongly affected by the neutron dynamics via the heat source from fission reactions. Furthermore, the neutron dynamics is strongly affected by the thermohydrodynamics via the motion of precursors and by feedback effects. The methods typical of solid fuel reactors of previous generations are not sufficient to handle these more highly coupled concepts. In the preset paper, we consider the coupling between neutronics and thermohydrodynamics with simple but realistic hypotheses assumed to model the evolution of all the variables involved in the calculation. The numerical scheme used represents the current state of the art in the solution of non-linear systems: the Newton–Krylov algorithm. Several calculations are presented to demonstrate the ability of the methods described here to study the behavior of molten salt reactors in both steady state and transient situations.