Parallelised Krylov subspace method for reactor kinetics by IQS approach

Nuclear reactor kinetics involves numerical solution of space–time-dependent multi-group neutron diffusion equation. Two distinct approaches exist for this purpose: the direct (implicit time differencing) approach and the improved quasi-static (IQS) approach. Both the approaches need solution of static space-energy-dependent diffusion equations at successive time-steps; the step being relatively smaller for the direct approach. These solutions are usually obtained by Gauss–Seidel type iterative methods. For a faster solution, the Krylov sub-space methods have been tried and also parallelised by many investigators. However, these studies seem to have been done only for the direct approach. In the present paper, parallelised Krylov methods are applied to the IQS approach in addition to the direct approach. It is shown that the speed-up obtained for IQS is higher than that for the direct approach. The reasons for this are also discussed. Thus, the use of IQS approach along with parallelised Krylov solvers seems to be a promising scheme.     

A preconditioned Krylov method for solution of the multi-dimensional,two fluid hydrodynamics equations

A preconditioned Krylov method is introduced for reducing the computational burden in the solution of the multi-dimensional, two fluid hydrodynamics equations. The BILU3D preconditioned BICGSTAB method was applied to the linearized continuity equation in the inner iteration of the fully implicit solution of the mass, energy, and momentum equations in the EPRI code VIPRE-02. A nuclear reactor thermal-hydraulics test problem was performed with 104 multi-dimensional flow channels in the vessel. For the simulation of a typical pressurized water Reactor steam line break transient, the overall execution time was reduced by more than 50% compared to the existing solution techniques which utilize stationary alternate direction implicit (ADI) iterative methods for the inner iteration.     

Generalized power series method with step size control for neutron kinetics equations

Based on the power series method (PWS), a generalized power series method (GPWS) has been introduced for solving the point reactor kinetics equations. The stiffness of the kinetics equations restricts the time interval to a small increment, which in turn restricts the PWS method within a very small constant step size. The traditional PWS method has been developed using a new formula that can control the time step at each step while transient proceeds. Two solvers of the PWS method using two successive orders have been used to estimate the local truncation errors. The GPWS method has employed these errors and some other constraints to produce the largest step size allowable at each step while keeping the error within a specific tolerance. The proposed method has resolved the stiffness point kinetics equations in a very simple way with step, ramp and zigzag ramp reactivities. The generalized method has turned out to represent a fast and accurate computational technique for most applications. The method is seemed to be valid for a time interval that is much longer than the time interval used in the conventional numerical integration, and is thus useful in reducing computing time. The method constitutes an easy-to-implement algorithm that provides results with high accuracy for most applications where, the reactor kinetics equations are reduced to a differential equation in a matrix form convenient for explicit power series solution. Results obtained by GPWS method: attest the power of the theoretical analysis, they demonstrate that the convergence of the iteration scheme can be accelerated, and the resulting computing time can be greatly reduced while maintaining computational accuracy. The point kinetics equations have been solved as a preliminary simple case aimed at testing the applicability of the GPWS method to solve point kinetics equations with feedback or, space kinetics problems.     

Analysis and improvements of the DPN acceleration technique for the method of characteristics in unstructured meshes

Algebraic preconditioners, renumbering techniques and a two-level algebraic multigrid method have been implemented to speed up the Krylov iterations of the DP_N equations used for the acceleration of the method of characteristics in unstructured meshes. These algorithms were customized to take advantage of the cell-based structure of the DP_N equations. Moreover, two techniques to speed up the solution of the multigroup eigenvalue MOC equations have been implemented. A solution of the multigroup eigenvalue DP_N equation has been developed to provide a first guess for the external transport iterations. Next, a multigroup DP_N acceleration method has been developed to accelerate the thermal iterations. This latter development has been particularly useful because our standard multigroup rebalancing acceleration was counterproductive in the presence of heavy absorbents. All these acceleration techniques have been incorporated in the spectral code APOLLO2. Numerical examples and comparisons are given for the 6-group eigenvalue Atrium benchmark problem. Our best calculation, an initialized ILU0-preconditioned DP₁ scheme with thermal acceleration, was 7.7 times faster that the free iteration calculation, while the total number of transport iterations was divided by 17.     

Code validation for the magnetohydrodynamic flow at high Hartmann Number based on unstructured grid

In order to analyze the magnetohydrodynamic (MHD) effect in liquid metal fusion blanket, a parallel and high performance numerical code was developed to study MHD flows at high Hartmann Number based on the unstructured grid. In this code, the induced current and the Lorentz force were calculated with a current density conservative scheme, while the incompressible Navier–Stokes equations with the Lorentz force included as a source term was solved by projection method, a set of method were used to improve the computing performance such as Krylov subspace method and AMG method. To validate this code, three benchmarks of MHD flow at high Hartmann Number were conducted. The first benchmark was the case of Shercliff fully development flow, the second benchmark was the MHD flow in a circular pipe with changing external magnetic field, and the third benchmark was the MHD flow in a pipe with sudden expansion. In these cases the Hartmann Numbers were from 1000 to 6000. The code good computing performance, and numerical results show matched well with the analytical and experimental results.     

Explicit Time Integration Scheme Using Krylov Subspace Method for Reactor Kinetics Equation

The spatial discretization form of the space-dependent reactor kinetics equation is a first-order simultaneous ordinary differential equation in time. Conventional numerical methods of the space-dependent kinetics equation, i.e., the generalized Runge-Kutta method, the implicit method (backward Euler method), and the Theta method, are based on the time difference approximation. However, the present study adopts the analytical solution of the space-dependent kinetics equation expressed by the matrix exponential and no time difference approximation is used. In this context, our present approach is classified as an explicit method in which no iteration calculation on space and energy is necessary. The Krylov subspace method is used to evaluate the matrix exponential observed in the solution of the spatially discretized space-dependent kinetics equation. The Krylov subspace method is implemented into a space-dependent kinetics solver. In order to examine the effectiveness of the Krylov subspace method, the TWIGL benchmark problem is analyzed as a verification calculation. The calculation results show the effectiveness of the present method especially in the step reactivity perturbation.     

三维六角形节块多群中子扩散程序NDHEX

本文介绍用DIF3D (NOD)求解二、三维六角形几何系统下中子扩散方程的理论模型及数值计算方法。六角形节块内的中子通量密度分布采用高次多项式近似表示,最后导出通量矩方程及偏流的响应矩阵方程。应用粗网再平衡和渐近源外推方法加速收敛。参考此方法编制了计算程序NDHEX,并对一些六角形基准问题进行了计算。结果表明:NDHEX的计算结果与DIF3D(NOD)的计算结果符合很好;与差分程序相比,具有更高的精度与计算效率。它可用于快堆计算。     

Charge evolution of swift-heavy-ion beams explored by matrix method

The approach to charge equilibrium of a beam of swift ions can be expressed by a matrix with rows and columns expressing the initial and instantaneous charge state, respectively. We have explored the capability of this matrix formalism in comparison with the standard method, i.e., numerical solution of linear rate equations as implemented in the ETACHA code. The matrix method is computationally more efficient, and the acceptable numerical error can be pre-defined. The agreement with predictions based on ETACHA is generally good. Significant discrepancies are primarily due to differences in the applied cross sections for electron capture and loss. Comparisons have been made also with experimental results for sulphur in carbon.     

注量展开节块法求解三维六角形几何中子扩散方程

提出了一种在三维六角形几何节块内数值求解中子扩散方程的节块法该方法把节块内各群中子注量分布用解析基函数近似展开为了改善节块耦合关系．提出了,一种新的节块边界条件：面平均偏流零次矩和一次矩同时保持连续。此外．将响应矩阵技术应用于迭代求解过程,使得该方法具有较高的计算效率基于本文提出的模型,发展了三维六角形组件中子扩散计算程序FEMHEX。通过对二维、三维VVER基准问题校验计算表明,该方法能高效．准确的给出有效增殖系数以及节块功率分布。     

PWS: an efficient code system for solving space-independent nuclear reactor dynamics

The reactor kinetics equations are reduced to a differential equation in matrix form convenient for explicit power series solution involving no approximations beyond the usual space-independent assumption. The coefficients of the series have been obtained from a straightforward recurrence relation. Numerical evaluation is performed by PWS (power series solution) code, written in Visual FORTRAN for a personal computer. The results are applied to the step reactivity insertion, ramp input, zigzag input, and oscillatory reactivity changes. When the reactivity is given, including the case in which the feedback reactivity is a function of neutron density, the developed method can provide a straightforward procedure for computing reactor dynamics problems. The solution of this method was compared to some other analytical and numerical solutions of the point reactor kinetics equations; the results proved that the approach is both efficient and accurate to several significant figures.     

Krylov sub-space methods for K-eigenvalue problem in 3-D neutron transport

The K-eigenvalue problem in nuclear reactor physics is often formulated in the framework of Neutron Transport Theory. The fundamental mode solution of this problem is usually obtained by the power iteration method. Here, we are concerned with the use of a Krylov sub-space method, called ORTHOMIN(1), to obtain a more efficient solution of the K-eigenvalue problem. A Matrix-free approach is proposed which can be easily implemented by using a transport code which can perform fixed source calculations. The power iteration and ORTHOMIN(1) schemes are compared for two realistic 3-D multi-group cases with isotropic scattering: an LWR benchmark and a heavy water reactor problem. In both the schemes, within-group iterations over self-scattering source are required as intermediate procedures. These iterations are also accelerated using another Krylov method called conjugate gradient method. The overall work is based on the use of Sn-method and finite-differencing for discretisation of transport equation.     

Simulation of a small loss of coolant accident by using RETINA V1.0D code

RETINA has been developed for modeling of two-phase flow situations in full-scope simulators of nuclear power plants. A special feature of RETINA is that both RETINA V1.0D (drift-flux — 5 equations) and RETINA V1.0-2V (two-fluid — 6 equations) approach are available in the code and the same constitutive relations are used in both cases. The governing equations are discretized implicitly, and an automatic derivation algorithm determines the Jacobian matrix, which is partitioned taking into account the special structure of nuclear power plants. Partitioned inverse formula is used to solve the global equation system providing the possibility of multi-level parallelization. Heat structures are modeled in two dimensions and are coupled to the flow equations explicitly. Since the code will be used in real-time simulators, we paid special attention to time-effective solution. In this paper, we demonstrate the ability of our code by simulating a small loss of coolant accident Paks Model Circuit (PMK). The simulation results are compared to real measurements obtained by Paks Model Circuit.     

Modifications to COBRA-TF to model dispersed flow film boiling with two flow,four field Eulerian–Eulerian approach – Part 2

In the first part of the paper, the modifications performed to improve the dispersed flow film boiling model in COBRA-TF have been described. The improvements were achieved by adding a small droplet field to the code’s solution scheme. The conservation equations, the source terms for the equations and the models developed were summarized. In this paper, the effects of spacer grids on the dispersed flow heat transfer and COBRA-TF modifications for the spacer grid models are presented. The results of the code predictions are presented by comparing the experimental data from Rod Bundle Heat Transfer experiments with the results of code simulations performed with original and modified code. Measurements and calculations for the spacer grid temperature have been compared. The results of the analysis performed with the modified code indicate the improvement in code predictions for the spacer grid temperature.     

Nodal kinetics model upgrade in the Penn State coupled TRAC/NEM codes

The Pennsylvania State University currently maintains and does development and verification work for its own versions of the coupled three-dimensional kinetics/thermal-hydraulics codes TRAC-PF1/NEM and TRAC-BF1/NEM. The subject of this paper is nodal model enhancements in the above mentioned codes. Because of the numerous validation studies that have been performed on almost every aspect of these codes, this upgrade is done without a major code rewrite. The upgrade consists of four steps. The first two steps are designed to improve the accuracy of the kinetics model, based on the nodal expansion method. The polynomial expansion solution of 1D transverse integrated diffusion equation is replaced with a solution, which uses a semi-analytic expansion. Further the standard parabolic polynomial representation of the transverse leakage in the above 1D equations is replaced with an improved approximation. The last two steps of the upgrade address the code efficiency by improving the solution of the time-dependent NEM equations and implementing a multi-grid solver. These four improvements are implemented into the standalone NEM kinetics code. Verification of this code was accomplished based on the original verification studies. The results show that the new methods improve the accuracy and efficiency of the code. The verification of the upgraded NEM model in the TRAC-PF1/NEM and TRAC-BF1/NEM coupled codes is underway.     

Development of efficient Krylov preconditioning techniques for multi-dimensional method of characteristics

Two techniques are proposed in the preconditioning for the Krylov sub-space method called the Generalized Minimal RESidual (GMRES) method to accelerate inner iterations based on the method of characteristics (MOC). The GMRES method is an iterative method to solve a linear algebraic system byminimizing the norm of the residual vector. The proposed preconditioning technique is based on the first flight collision probability which is efficiently made by the multi-dimensional MOC code. To simplify the preconditioner, slight couplings among regions are ignored by considering the mean free path. And another proposed technique makes simplified preconditioner by the scaling matrix which can homogenize and de-homogenize the fuel region and the cladding region. The scaling technique reduces the size of the matrix and also reduces the calculation time of inverse matrix. Numerical results show that the preconditioner simplified by the mean free path efficiently reduces the number of iterations for the GMRES algorithm. And the scaling technique keeps the efficiency of preconditioner even in the multi-dimensional geometry. The total calculation time is found to be reduced when these techniques are employed.     

An analytical solution of the point kinetics equations with time-variable reactivity by the decomposition method

In this work, we report an analytical solution for the point kinetics equations by the decomposition method, assuming that the reactivity is an arbitrary function of time. The main idea initially consists in the determination of the point kinetics equations solution with constant reactivity by just using the well-known solution results of the first-order system of linear differential equations in matrix form with constant matrix entries. Applying the decomposition method, we are able to transform the point kinetics equations with time-variable reactivity into a set of recursive problems similar to the point kinetics equations with constant reactivity, which can be straightly solved by the mentioned technique. For illustration, we also report simulations for constant, linear and sinusoidal reactivity time functions as well comparisons with results in literature.     

An analytical solution of the point kinetics equations with time-variable reactivity by the decomposition method

In this work, we report an analytical solution for the point kinetics equations by the decomposition method, assuming that the reactivity is an arbitrary function of time. The main idea initially consists in the determination of the point kinetics equations solution with constant reactivity by just using the well-known solution results of the first-order system of linear differential equations in matrix form with constant matrix entries. Applying the decomposition method, we are able to transform the point kinetics equations with time-variable reactivity into a set of recursive problems similar to the point kinetics equations with constant reactivity, which can be straightly solved by the mentioned technique. For illustration, we also report simulations for constant, linear and sinusoidal reactivity time functions as well comparisons with results in literature.     

An improved numerical scheme with the fully-implicit two-fluid model for a fast-running system code

A new computational method is implemented in the FISA-2 (Fully-Implicit Safety Analysis-2) code to simulate the thermal-hydraulic response to hypothetical accidents in nuclear power plants. The basic field equations of FISA-2 consist of the mixture continuity equation, void propagation equation, two phasic momentum equations, and two phasic energy equations. The fully-implicit scheme is used to eliminate a time step limitation and the computation time per time step is minimized as much as possible by reducing the matrix size to be solved. The phasic energy equations written in the nonconservation form are solved after they are set up to be decoupled from other field equations. The void propagation equation is solved to obtain the void fraction. Spatial acceleration terms in the phasic momentum equations are manipulated with the phasic continuity equations so that pseudo-phasic mass flux may be expressed in terms of pressure only. Putting the pseudo-phasic mass flux into the mixture continuity equation, we obtain linear equations with pressure variables only as unknowns. By solving the linear equations, pressures at all the nodes are obtained and in turn other variables are obtained by back-substitution. The above procedure is performed until the convergence criterion is satisfied. Reasonable accuracy and no stability limitation with fast-running are confirmed by comparing results from FISA-2 with experimental data and results from other codes.     

Exact Error Estimation for Solutions of Nuclide Chain Equations

The exact solution of nuclide chain equations within arbitrary figures is obtained for a linear chain by employing the Bateman method in the multiple-precision arithmetic. The exact error estimation of major calculation methods for a nuclide chain equation is done by using this exact solution as a standard. The Bateman, finite difference, Runge-Kutta and matrix exponential methods are investigated.

The present study confirms the following. The original Bateman method has very low accuracy in some cases, because of large-scale cancellations. The revised Bateman method by Siewers reduces the occurrence of cancellations and thereby shows high accuracy. In the time difference method as the finite difference and Runge-Kutta methods, the solutions are mainly affected by the truncation errors in the early decay time, and afterward by the round-off errors. Even though the variable time mesh is employed to suppress the accumulation of round-off errors, it appears to be nonpractical. Judging from these estimations, the matrix exponential method is the best among all the methods except the Bateman method whose calculation process for a linear chain is not identical with that for a general one.