The multi-dimensional neutron kinetics model in RELAP5-3D^©is based on the NESTLE¹⁴ code developed by Paul Turinsky and co-workers atNorth Carolina State University under an INEEL initiative. The NESTLE code solves the twoor four group neutron diffusion equations in either Cartesian or hexagonal geometry usingthe nodal expansion method. Three, two, or one dimensional models may be used. Severaldifferent core symmetry options are available including quarter, half, and full coreoptions for Cartesian geometry and 1/6, 1/3, and full core options for hexagonal geometry.Zero flux, non-reentrant current, reflective, and cyclic boundary conditions areavailable. The steady-state eigenvalue and time dependent neutron flux problems can besolved by the NESTLE code as implemented in RELAP5-3D^© .

The few group neutron equations are spatially discretized using the Nodal ExpansionMethod (NEM). Quartic or quadratic polynomial expansions are used for the transverselyintegrated fluxes in the Cartesian and hexagonal geometries respectively. Discontinuityfactors are used to correct for homogenization errors.

Transient problems employ a user-specified number of delayed neutron precursor groups.The time discretization is done with a fully implicit method. The delayed neutronprecursor equations are integrated analytically assuming a linear variation of the neutronflux.

An outer-inner iteration strategy is employed to solve the matrix of equationsresulting from the temporal and spatial discretization. Outer iterations use Chebyshevacceleration. Inner iterations use a point or line SOR iteration scheme. The non-lineariterative strategy associated with the NEM is employed. Using the non-linear iterativestrategy means that the system of equations is equivalent to the finite difference methodusing the box scheme. This means that the code can solve either the nodal or finitedifference representation of the few-group neutron diffusion equations.

Thermal-hydraulic (TH) feedback is accomplished through the neutron cross sections. Thethermal-hydraulic conditions in the neutron kinetics nodes are computed by RELAP5-3D^© using the hydraulic model of the system specified by the userusing normal RELAP5 input. The conditions in the RELAP5 volumes are mapped onto thekinetics nodes by map functions using user-supplied data. The neutron cross sections inthe kinetics nodes are computed from RELAP5 calculated TH conditions and user suppliedcomposition data. The cross section model is quite general and the user can choose betweendifferent sets of TH conditions for use in the neutron cross section model, e.g., eithervoid fraction or coolant density as one of the independent variables in the cross sectionmodel. Neutron cross sections are computed for each kinetics node from the TH conditionsmapped onto the individual nodes and the composition data mapped onto the same nodes.

The compositions are mapped onto the kinetic nodes using user supplied compositionmaps. In this way a small number of compositions and a small number of RELAP5 volumes andheat structures may be mapped to a large number of kinetics solution nodes. Discontinuityfactors are modeled using the neutron cross section model so that they may change as theTH conditions in the kinetics nodes change during the course of the transient.

A control rod model has been implemented and each neutron kinetics node may contain anynumber of control rods. Control rods may be inserted from either the top or bottom of thereactor. The control rods may be grouped into banks, with the banks having the sameinsertion depth and velocity. The position of a control rod bank can also be determinedfrom a control variable or a general table. Scram signals may be initiated by RELAP5trips.

The reactor power and its distribution computed by the kinetics module are supplied tothe RELAP5 TH solution through the same maps that are used to compute the neutron crosssections. As the neutron flux and its distributions change during the course of thetransient, so does the power supplied to the individual RELAP5 volumes and heat structureschange accordingly.

The decay heat model as implemented in previous code versions for the point kinetics isused to compute the decay heat in the multi-dimensional neutron kinetics model. The decayheat is calculated in each kinetics node based on the fission power in that node.

The implementation of the NESTLE neutron kinetics has been verified by the simulationof the NEACRP¹⁵ three dimensional benchmark problems¹⁶. A series ofthree PWR rod ejection accidents from Hot Zero Power and Hot Full Power were proposed as abenchmark by the NEACRP. Series A, the ejection of a central rod, and series B, theejection of peripheral rods, were chosen for simulation using the spatial kinetics option.The location of the ejected control rods and the initial core configuration are shown inFigure 24. For the series A and B transients one-quarter geometry was adequate.

The RELAP5 core model for the benchmark problem consisted of a sequence of parallelpipes as shown in Figure 25. Each pipe was described using a series of heat structures andcontrol volumes to model the fuel and coolant from a single assembly. A separate inletreservoir was provided for each pipe and was set at a constant pressure, temperature, andboron concentration. Each pipe was connected to its reservoir using a time-dependentjunction that provided identical, constant flow to each assembly. The outlet of each pipewas connected to an outlet reservoir using a series of branches. Since a maximum of ninejunctions per branch is allowed, additional branches were used to accommodate the total of47 pipes of the quarter core model.

During the course of the analysis, various thermal-hydraulic mesh structures wereexamined. The finest axial mesh used for the pipes corresponds to that prescribed in theNEACRP problem with the exception that the three smallest nodes at the bottom and top ofthe core were combined into a single thermal- hydraulic node. The original NEACRP meshstructure for the neutronic solution was retained. This provided for a total of 14 axialthermal-hydraulic meshes in each pipe. For purposes of the thermal-hydraulic calculation,the reflector was modeled as a separate, single pipe because of the very small temperaturerise anticipated in the reflector.

For the fuel pin model 8, 1, and 2 meshes in the fuel pellet, gap, and cladding,respectively, were used. In order to generate an effective Doppler temperature in RELAP5consistent with that prescribed in the NEACRP problem,

a very small central and peripheral fuel pellet mesh was used such that the area ofthese regions corresponds to the weighting of the surface temperature, T_F,S,and centerline fuel temperature, T_F,C specified in the benchmarkproblem. All other fuel pellet meshes were equidistant.

Because the benchmark problem specifies gap conductance and RELAP5 only accepts gapconductivity, the following relation was used:

where k and h are the gap conductivity and conductance, respectively, andD r is the gap width. This expression is valid for asmall gap width as used in the benchmark problem.

The axial mesh structure used for the neutronic solution was identical to thatspecified in the benchmark problem and four nodes per fuel assembly were used in theradial plane. The partial cross sections prescribed in the NEACRP benchmark were processedinto an equivalent set of cross section multipliers to coincide with the modeling used inRELAP5. Some minor discrepancy persisted since it is not possible to construct acompletely consistent set of partial cross section data for the diffusion coefficient usedin RELAP5 based on the partial transport cross section data specified in the benchmarkproblem.

The steady state solution was obtained by running a null transient with RELAP5 for 20seconds using time steps of 0.05 secs. The current algorithm provides for athermal-hydraulics update only at the completion of each k-effective solution. Theconvergence criteria for the neutronic solution were set at 10^-4 and 10^-5for the global fission source and k-effective, respectively. The coupling coefficients inthe NEM nonlinear iteration method were updated every five outer iterations. For thetransient solution, the global fission source convergence criterion was increased to 10^-5.The time step sizes for the transients are given in Table 3.

The results of the four transient cases are summarized in Table 4. For each of the fourcases the predicted values and their deviations from the reference are shown. Thereference results are those reported at the Karlsruhe conference¹⁷. Thetransient power for each of the four cases is compared with the reference in Figure 26 forA1 and B1 and in Figure 27 for A2 and B2. The RELAP5 radial power distribution at axialplane 6 (out of 16) for problem A1 is compared to the reference in Figure 28 for thesteady state (t=0), for the maximum power condition (t=0.611 secs.), and for theasymptotic core state (t=5 secs.).

As indicated in Table 4, the RELAP5 steady state results for the hot zero power casesare in reasonably good agreement with the reference results. RELAP5 does show a slightnegative bias in the prediction of the critical boron concentration for both cases A1 andB1. Also note that RELAP5 is consistently lower than the reference in its prediction ofthe power in the rodded locations for the steady-state. This is consistent with thenegative bias in the critical boron prediction since a greater estimation of the total rodworth would lead to a lower critical boron concentration. The negative bias is larger incase A1 because there are several more rods inserted than in case B1.

The RELAP5 ejected rod worth prediction for both cases is in good agreement with thereference result. As shown in Figure 26, agreement of the transient results is excellentfor case B1, however, discrepancies exist in the RELAP5 and reference results for case A1.RELAP5 predicts the maximum power will occur 0.046 seconds later and have a magnitude22.1% lower than the reference calculation.

The asymptotic core state predicted by RELAP5 is in reasonably good agreement with thereference result as shown in Table 4, although RELAP5 shows a slight negative bias in theasymptotic core power.

The RELAP5 steady state results for the hot full power cases show a slight negativebias in the prediction of the critical boron concentration.

As shown in Figure 27, the RELAP5 transient results for case A2 shows a positive biasin the prediction of the maximum power, whereas case B2 is in close agreement with thereference results. As shown in Table 4, there is good agreement between RELAP5 and thereference in the prediction of the time of the maximum power, as well as in the predictionof the asymptotic core power.

Figure 28. Core Radial Power Distribution at Axial Plane 6 for Case A1

The Border Profiled Lower Upper (BPLU) matrix solver is used to efficiently solvesparse linear systems of the form AX = B. BPLU is designed to take advantage of pipelines,vector hardware, and shared-memory parallel architecture to run fast. BPLU is mostefficient for solving systems that correspond to networks, such as pipes, but is efficientfor any system that it can permute into border-banded form.

Speed-ups are achieved for RELAP5-3D^© running with BPLUover the default solver. For almost one-dimensional problems, there is no speed-up;however, for problems with wider bandwidths, especially those with three-dimensionalregions, significant speed-ups may be achieved. One of the standard installation problems,"3dflown.i" illustrates the reduction in run time that can be achieved. Theproblem is a simple cube subdivided into a 3x3 region in each of the Cartesian coordinatedirections (Figure 29). There are nine cases examined with this model, comprised of flowin each coordinate direction (x,y,z) of vapor only, liquid only, and a two-phase mixture.

Figure 29. 3dflow problem

Table 5 compares the CPU times required for the nine cases for the default solver andthe BPLU solver. Results show an average reduction in CPU time of 63%.

The BPLU algorithm is designed to solve "nearly-banded" coefficient matricesdirectly. It is a variation of banded Gaussian Elimination with partial pivoting thatincorporates variable reordering. A nearly-banded matrix is a sparse matrix that can bewritten as the sum of a banded matrix and a matrix of outriggers or non-zero entries thatlie outside the band. The algorithm restructures the linear system into a form that can besolved efficiently on a vector-parallel computer, does not waste operations on zeroentries, and controls the creation of non-zero elements.

The algorithm is suitable for and has been tested with the linear systems of thefollowing variety: the RELAP5 semi-implicit time step matrix, the RELAP5 nearly-implicittime step matrix, and the NEWEDGE field equations matrix. All these matrices are nearlybanded with outriggers. Further, the structure of the nearly-implicit matrix is similar tothat of the semi-implicit matrix, but with the elements replaced by 2x2 blocks ofelements. The NEWEDGE matrix replaces the 2x2 blocks with 5x5 blocks.

Two considerations lead to further efficiencies. First, the block structure discussedabove leads to a profile rather than a solid band structure for which zero operations canbe avoided by using a variable row length during row reduction. Profiled-structures alsoarise from connectivity relationships among the variates. Therefore, the BPLU algorithmuses a profile, rather than a banded, solver. Second, a nearly-banded matrix canfrequently be decomposed into the sum of a banded and an outrigger matrix in more than oneway. Selection of the bandwidth for the decomposition is tied to reducing the amount ofwork associated with solving the reorganized system.

Other improvements in the RELAP5-3D^© code include:

More accurate water properties derived from the National Bureau of Standards and National Research Council

Enhanced robustness from intelligent recovery from advancement failure

Ongoing and planned DOE funding will provide additional improvements:

Graphical User Interface

More implicit momentum equation permitting larger time step sizes

Multidimensional heat conduction model for modeling graphite structures in RBMK reactors

FORTRAN 90

Windows ‘95/NT version

 

The INEEL will make the RELAP5-3D^© code available todomestic and international organizations that join the International Users Group. This isa subscriber-funded group that takes part in guiding the future development andmaintenance of the code. Organizations interested in obtaining the code will be given theopportunity to participate at one of two levels:

Participant - Organizations that wish to receive the code, future updates, and some on-call assistance.

Member - Organizations that wish to receive the code, future updates, a guaranteed level of on-call assistance, and participate in voting on improvements.

Both domestic (U.S.) and non-domestic organizations may choose between these two levelsof participation. Subscription costs for one year (two year minimum) are shown below. Moreinformation on the structure and function of the users group will soon be available on thehome page: relap5.inel.gov.

Questions or problems regarding this web site should be directed to danidl@inel.gov.

Last modified: Thursday November 12, 1998.