RELAP5-3D^^

# INTRODUCTION

This paper presents and illustrates some of the key features of the
RELAP5-3D^^ code being developed at the Idaho National Engineering and
Environmental Laboratory (INEEL) for the US Department of Energy (DOE).
The purpose of the paper is to inform potential users of the code of its
unique capabilities that extend its range of applicability beyond that
presently available with any other thermal-hydraulics systems code. The
paper also discusses areas of ongoing development, future plans, and the
availability of the code to the international community.

# BACKGROUND

The RELAP5-3D^^ code is an outgrowth of the RELAP5/MOD3 code developed
at the Idaho National Engineering & Environmental Laboratory (INEEL) for
the U.S. Nuclear Regulatory Commission (NRC). Development of the RELAP5
code series began at the INEEL under NRC sponsorship in 1975 and
continued through several released versions, ending in October 1997 with
the soon to be released RELAP5/MOD3.3. The U.S. Department of Energy
(DOE) began sponsoring additional RELAP5 development in the early 1980s
to meet its own reactor safety assessment needs. Following the accident
at Chernobyl, DOE undertook a re-assessment of the safety of all of its
test and production reactors throughout the United States. The RELAP5
code was chosen as the thermal-hydraulic analysis tool because of its
widespread acceptance. Systematic safety analyses were carried out for
the DOE that included the N reactor at Hanford, the K and L reactors at
Savannah River, the Advanced Test Reactor (ATR) at INEEL, the High Flux
Isotope Reactor (HFIR) and Advanced Neutron Source (ANS) at Oak Ridge,
and the High Flux Beam Reactor (HFBR) at Brookhaven. DOE also chose
RELAP5 for the independent safety analysis of the New Production Reactor
(NPR) proposed for Savannah River before that program was canceled in
the wake of the end of the cold war.

The application of RELAP5 to these various reactor designs created the
need for new modeling capabilities. For example, the analysis of the
Savannah River reactors necessitated a three-dimensional flow model and
heavy water properties be added to the code. ATR required a new critical
heat flux correlation applicable to its unique fuel design. All
together, DOE sponsored improvements and enhancements have amounted to a
multimillion-dollar investment in the code.

Toward the end of 1995, it became clear that the efficiencies realized
by the maintenance of a single source code for use by both NRC and DOE
were being overcome by the extra effort required to accommodate
sometimes conflicting requirements. The code was therefore "split" into
two versions, one for NRC and the other for DOE. The DOE version
maintained all of the capabilities and validation history of the
predecessor code, plus the added capabilities that had been sponsored by
the DOE before and after the split.

**RELATIONSHIP OF RELAP5-3D^^ TO PRIOR VERSIONS**

At the outset of the decision to split the code into NRC and DOE
versions, the INEEL recognized the importance of retaining the pedigree
stemming from the extensive validation history of RELAP5/MOD3.
Consequently, the developmental activities with respect to RELAP5-3D^^
since the split have been carefully integrated so as not to compromise
this legacy validation. In fact, virtually all of the enhancements in
RELAP5-3D^^ are *optional* and *supplemental* to the proven performance
of RELAP5/MOD3.2. Consequently, users of RELAP5-3D^^ can confidently
apply the code using existing, one --dimensional RELAP5/MOD3.2 input
decks and expect their results to be the same or improved.

To ensure that the code remains true to its validation history (or
provides improved results), developmental versions of the code are
periodically tested using a subset of cases from the RELAP5/MOD3
validation library. Two such cases are the G.E Level Swell and THTF
experiments.

# The G.E. Level Swell experiments^1^, conducted in the late 1970's provided excellent data for assessing flashing and interphase drag models. Two tanks, one 1 ft. in diameter and the other 4 ft. in diameter were used. In the experiments, the tanks were partially filled with water and heated to near the saturation temperature. They were then depressurized through a blowdown valve. Differential pressure measurements made along the vertical length of the tank allowed for computing average void fractions as a function of elevation inside the tank.

Figure 1 shows the calculated and measured void profiles in the 4 ft.
diameter tank of Test 5801-13 at four times during the blowdown (5, 10,
15, and 20 sec.). RELAP5/MOD3.2 and RELAP5-3D^^ produced identical
results, both being in excellent agreement with the data.

The Thermal Hydraulic Test Facility (THTF) at Oak Ridge National
Laboratory (ORNL) was designed to simulate conditions in a PWR core. The
test section contained 64 electrically heated rods with internal
dimensions typical of a 17 by 17 rod bundle. Differential pressure
measurements positioned along the axial length of the bundle were used
to compute void fractions, and thermocouples were employed to measure
steam temperature and heater rod surface temperature.

In Test 3.09.10i^2^, the bundle was maintained in a partially uncovered
condition at a pressure of 4.5 MPa, an inlet mass flux of 29.8
kg/m^2^-s, and an inlet subcooling of 57.6 K. The heater rod heat flux
(uniform axial and radial) was 7.44x10^4^ W/m^2^. Figure 2 compares
RELAP5/MOD3.2 and RELAP5-3D^^ results with measured data for bundle
void profile, heater rod temperature profile, and vapor temperature
profile. Again, the two code versions are virtually identical and agree
well with the data.

The consistency of results between RELAP5/MOD3.2 and RELAP5-3D^^ is not
surprising, in light of the fact that the constitutive models for
one-dimensional flow have not been significantly altered.

# ![](media/image1.png)Figure 1. Calculated and Measured Void Profiles During GE Level Swell Test 5801-13 

# ![](media/image2.wmf)Figure 2. Measured and Calculated Results from THTF Test 3.09.10i 

# RELAP5-3D^^ CODE FEATURES

The most prominent attribute that distinguishes the DOE code from the
NRC code is the fully integrated, multi-dimensional thermal-hydraulic
and kinetic modeling capability in the DOE code. This removes any
restrictions on the applicability of the code to the full range of
postulated reactor accidents. Other enhancements include a new matrix
solver, new water properties, and improved time advancement for greater
robustness. Together with the existing modeling capabilities of
RELAP5/MOD3.2, these enhancements make the code the most powerful tool
of its kind available.

The balance of this paper focuses on the capabilities of the
three-dimensional hydrodynamic model, the multi-dimensional kinetics
model, and the new BPLU matrix solver. Other features unique to
RELAP5-3D^^ are briefly described.

## Three-Dimensional Hydrodynamic Model

The development of the three-dimensional hydrodynamic model in RELAP5
began in 1990 under funding from the Savannah River National Laboratory
(SRNL). Since then development and testing has continued from both
planned improvements and responses to user requests. This has resulted
in a reliable model. Progress has been documented in the literature and
at various meetings (References 3-12).

The multi-dimensional component in RELAP5-3D^^ was developed to allow
the user to more accurately model the multi-dimensional flow behavior
that can be exhibited in any component or region of a LWR system.
Typically, this will be the lower plenum, core, upper plenum and
downcomer regions of an LWR. However, the model is general, and is not
restricted to use in the reactor vessel. The component defines a one,
two, or three-dimensional array of volumes and the internal junctions
connecting them. The geometry can be either Cartesian (x, y, z) or
cylindrical (r, θ, z). An orthogonal, three-dimensional grid is defined
by mesh interval input data in each of the three coordinate directions.

### Model Verification

Verification of the model has been performed by using conceptual
problems that have exact solutions. These types of problems are used to
demonstrate that the equations have been correctly coded, and are a
precursor to model validation using experimental data. Three such
problems are the "Rigid Body Rotation", "Pure Radial Symmetric Flow",
and R-θ Symmetric Flow" problems. Each of the problems is based on a
cylindrical, multidimensional component with eight rings, six sectors,
and one axial level. All six sectors are symmetrical, with non‑uniform
radial spacing. Six time dependent volumes are attached to the six outer
sectors by time dependent junctions for inlet flow specification. In
addition, six time dependent volumes are attached to the six inner
sectors by a multiple junction component. Figure 3 shows the
nodalization for 1 sector.

![](media/image3.emf)

### Figure 3. Nodalization of One Sector of Test Problem

In each of the simulations, losses due to friction and body force terms
are deactivated to create problems with exact solutions. Also, the flow
is assumed to be steady and incompressible. The azimuthal flow pattern,
where necessary, was imposed by setting the outer ring azimuthal
velocities to the desired value. All problems have analytic solutions
for velocity from the continuity and θ‑momentum equations, and analytic
solutions for pressure from the r‑momentum equations.

The Rigid Body Rotation problem represents a hollow cylinder with a
symmetric flow pattern in the azimuthal direction. Flow in the radial
direction does not exist. This problem tests only the azimuthal momentum
flux terms. The test conditions and boundary conditions are: The
azimuthal velocity at the 6.5 m radius position is 1 m/s, pressure at
the 1 m radius position is 5x10^5^ Pa, and all radial velocities are 0.0
m/s.

Comparisons between the RELAP5-3D^^ calculated results and the analytic
solution for the Rigid Body Rotation problem are shown in Figures 4 and
5. The calculated velocity and pressure profiles are seen to exactly
match the analytical solutions.

![](media/image4.png){width="3.140277777777778in" height="3.19375in"}

##### Figure 4. Azimuthal Velocity Profile for Rigid Body Rotation Problem

![](media/image5.png){width="3.1805555555555554in"
height="3.173611111111111in"}

#### Figure 5. Radial Pressure Profile for Rigid Body Rotation Problem

The Pure Radial Symmetric Flow problem represents a hollow cylinder with
a symmetric flow pattern in the radial direction. This problem tests
only the radial momentum terms. The test conditions and boundary
conditions are: All azimuthal velocities are 0 m/s, pressure at the 1 m
radius position is 5x10^5^ Pa, and radial velocity at the 7.5 m radius
position is 0.8667 m/s (from outside to inside of the ring). Comparisons
between the calculated results and the analytic solution for the Pure
Radial Symmetric Flow problem are shown in Figures 6 and 7. Again, the
RELAP5-3D^^ results are in agreement with the exact solution.

![](media/image6.png){width="2.9902777777777776in"
height="3.213888888888889in"}

##### Figure 6. Radial Velocity Profile for Pure Radial Symmetric Flow Problem

![](media/image7.png){width="3.203472222222222in"
height="3.213888888888889in"}

###### Figure 7. Radial Pressure Profile for Pure Radial Symmetric Flow Problem

The R‑θ Symmetric Flow problem represents a hollow cylinder with a
symmetric flow pattern in both the radial and azimuthal directions. The
azimuthal velocity at the 6.5 m radius position is 1 m/s, pressure at
the 1 m radius position is 5x10^5^ Pa, and radial velocity at the 7.5 m
radius position is 0.8667 m/s (from outside to inside of the ring). The
calculated and analytic solution results are seen to be in agreement
(Figures 8, 9, and 10).

![](media/image8.wmf){width="3.2291666666666665in"
height="2.8833333333333333in"}

##### Figure 8. Radial Velocity Profile for R-θ Symmetric Flow Problem

![](media/image9.png){width="2.9069444444444446in"
height="3.216666666666667in"}

##### Figure 9. Azimuthal Velocity Profile for R-θ Symmetric Flow Problem

![](media/image10.png){width="3.2270833333333333in"
height="3.216666666666667in"}

#### Figure 10. Radial Pressure Profile for R-θ Symmetric Flow Problem

Both the analytic and calculated velocity profiles come from the
continuity and θ‑momentum equations; they are each a function of radius
and outer ring velocity. The velocity profile from the Rigid Body
Rotation problem is linear (Figure 4), while the velocity profiles from
the other problems are inversely proportional to radius (Figures 6, 8
and 9). The pressures are calculated from the r‑momentum equations
(Figures 5, 7, and 10). The pressure profile for the Rigid Body Rotation
problem is concave up (Figure 5), while the others are concave down.
Additionally, the pressure profile from the R-θ Symmetric Flow problem
is influenced by the 1‑D outlet connections. There, the azimuthal flow
contributes to the pressure solution up to the point of connection.

These simple test problems verify the correct implementation of the
three dimensional continuity and momentum equations. This is a
necessary, but certainly not sufficient, test of the 3D model. Work has
begun on performing a number of validation cases using experiments
exhibiting multidimensional flow behavior.

####### Model Validation

One such case recently completed is a simulation of LOFT large break
loss-of-coolant-experiment (LOCE) L2-5^13^. The LOFT facility was a
50 MW pressurized water reactor (PWR) that was designed to simulate the
response of a commercial PWR during a loss-of-coolant accident (LOCA).
Test L2-5 simulated a double-ended offset shear of a cold leg. The
experiment was selected because it was judged to provide the most
challenging test of the multi-dimensional model of all the experiments
in the existing developmental assessment test matrix used for
RELAP5/MOD3.

The analysis was accomplished in several steps. First, the original
one-dimensional model was upgraded to be consistent with current user
guidelines and to better represent Test L2-5. The upgraded
one-dimensional model was then run with the current version of
RELAP5-3D^^ to obtain baseline results.

The RELAP5/MOD3 model that was used in the developmental assessment of
LOFT Test L2-5 is illustrated in Figure 11. The model, hereafter
referred to as the one-dimensional model, contains 131 volumes, 142
junctions, and 77 heat structures.

The three-dimensional model of the LOFT reactor vessel was developed
using two multi-dimensional components. Multidimensional component 700
represented the downcomer region as shown in Figure 12. Multidimensional
component 200 represented the lower plenum, core inlet, core, and upper
plenum regions as shown in Figure 13. The vessel was divided into four
90^o^ azimuthal sectors and four radial rings. The four azimuthal
sectors corresponded to the four nozzles connecting the loops and the
vessel. One radial ring (multidimensional component 700) represented the
downcomer while the other three rings (multidimensional component 200)
corresponded to high-, average-, and low-powered regions of the core.
The axial nodalization of each multi-dimensional component was based on
that of the one-dimensional model, resulting in 6 levels for the
downcomer and 21 levels for the lower plenum, core inlet, core, and
upper plenum regions. A multiple junction (Component 709) connected the
bottom of the downcomer to the top of the third ring in the lower
plenum. The three-dimensional model of the reactor vessel was then
inserted into the one-dimensional model, with the resulting model
hereafter referred to as the three-dimensional model.

The core fuel rods were modeled with twelve heat structure geometries,
each representing the fuel rods located in a given ring and sector. Each
fuel rod heat structure geometry contained twelve heat structures,
corresponding to the number of axial levels in the hydraulic
nodalization of the core. Each fuel rod heat structure geometry was
identical except for its hydraulic connections and the radial power
peaking factor. The radial power peaking factors were 1.31, 1.13, and
0.81 for the inner, middle, and outer rings respectively. For
comparison, the radial power peaking factors for the high-powered and
low-powered rods in the one-dimensional model were 1.31 and 0.94,
respectively. The maximum axial power peaking factor was 1.58 and was
located in the bottom third of the core.

![](media/image11.wmf)

![](media/image12.wmf){width="5.989583333333333in"
height="7.534722222222222in"}

Figure 12. Nodalization of the 3D Model of the LOFT Downcomer

![](media/image13.wmf){width="6.114583333333333in"
height="8.203472222222222in"}Figure 13. Nodalization of the 3D Model of
the LOFT Lower Plenum, Core Inlet, Core, and Upper Plenum Regions

Steady-state calculations were performed for LOFT Test L2-5 with both
the one-dimensional and three-dimensional models. Table 1 shows that the
results of the steady-state calculations were in excellent agreement
with the measurements.

Table 1. A comparison of calculated and measured initial conditions in
LOFT Test L2-5.

+-------------------------------+------------+-----------+-----------+
| Parameter                     | Measured   | One-di    | Three-di  |
|                               | Value      | mensional | mensional |
|                               |            | Model     | Model     |
+-------------------------------+------------+-----------+-----------+
| Intact loop                   |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Mass flow (kg/s)              | 192.4 ±    | 192.4     | 192.4     |
|                               | 7.8        |           |           |
+-------------------------------+------------+-----------+-----------+
| Hot leg pressure (MPa)        | 14.94 ±    | 14.92     | 14.92     |
|                               | 0.06       |           |           |
+-------------------------------+------------+-----------+-----------+
| Cold leg temperature (K)      | 556.6 ±    | 556.6     | 556.7     |
|                               | 4.0        |           |           |
+-------------------------------+------------+-----------+-----------+
| Hot leg temperature (K)       | 589.7 ±    | 590.4     | 590.5     |
|                               | 1.6        |           |           |
+-------------------------------+------------+-----------+-----------+
| Pressurizer liquid level (m)  | 1.14 ±     | 1.14      | 1.14      |
|                               | 0.03       |           |           |
+-------------------------------+------------+-----------+-----------+
| Average pump speed (rad/s)    | 131.5 ±    | 130.7     | 130.7     |
|                               | 1.2        |           |           |
+-------------------------------+------------+-----------+-----------+
| Pump differential pressure    | 73.3 ± 9.2 | 63.4      | 63.5      |
| (kPa)                         |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Reactor vessel                |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Power (MW)                    | 36.0 ± 1.2 | 36.0      | 36.0      |
+-------------------------------+------------+-----------+-----------+
| Maximum linear heat           | 36.0 ± 2.7 | 34.3      | 34.3      |
| generation rate               |            |           |           |
|                               |            |           |           |
| (kW/m)                        |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Maximum fuel centerline       | 1660 ± 57  | 1710      | 1712      |
| temperature (K)               |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Differential pressure (kPa)   | 28.0 ± 1.4 | 28.1      | 28.1      |
+-------------------------------+------------+-----------+-----------+
| Steam generator secondary     |            |           |           |
| side                          |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Pressure (MPa)                | 5.85 ±     | 5.85      | 5.86      |
|                               | 0.06       |           |           |
+-------------------------------+------------+-----------+-----------+
| Mass flow (kg/s)              | 19.1 ± 0.4 | 19.1      | 19.1      |
+-------------------------------+------------+-----------+-----------+
| Feedwater temperature (K)     | 482.0 ±    | 482.0     | 482.0     |
|                               | 1.2        |           |           |
+-------------------------------+------------+-----------+-----------+
| Accumulator                   |            |           |           |
+-------------------------------+------------+-----------+-----------+
| Pressure (MPa)                | 4.29 ±     | 4.29      | 4.29      |
|                               | 0.06       |           |           |
+-------------------------------+------------+-----------+-----------+
| Temperature (K)               | 303.2 ±    | 303.2     | 303.2     |
|                               | 6.1        |           |           |
+-------------------------------+------------+-----------+-----------+
| Liquid level above standpipe  | 1.17 ±     | 1.16      | 1.16      |
| (m)                           | 0.01       |           |           |
+-------------------------------+------------+-----------+-----------+

Transient simulations of LOFT Test L2-5 were performed using the
one-dimensional and three-dimensional models. Comparisons between
calculated and measured results, including the sequence of events, the
overall system response, and the core thermal response, are next
discussed.

The measured sequence of events for LOFT Test L2-5 is presented in
Table 2. The test was initiated at 0.0 s when the quick-opening blowdown
valves began to open. A reactor trip signal was generated at 0.02 s on
low hot leg pressure, and the reactor scrammed shortly thereafter. The
operators tripped the primary coolant pumps at 0.94 s. Flow from the
accumulator, HPIS, and LPIS began at 16.8 s, 23.90 s, and 37.32 s,
respectively. The accumulator emptied at 49.6 s. The fuel rod cladding
temperature first departed from near the saturation temperature at 0.91
s. The peak cladding temperature occurred at 28.47 s. All of the core
cladding had been quenched by 65 s, concluding the interesting portion
of the test.

Table 2 also presents the calculated event times with the
one-dimensional and three-dimensional RELAP5 models. The calculated
event times with both RELAP5 models were generally in reasonable
quantitative agreement with the measured values. One exception was that
the calculated peak cladding temperatures with both models occurred near
6 s, compared to about 28 s in the test. As will be shown later, the
measured cladding temperatures increased slowly between 5 s and 28 s
while the calculated temperatures decreased slowly, so the effect of the
difference in timing is not as significant as might be inferred from
Table 2. The core cladding also quenched 10 s to 15 s earlier in the
calculations than in the test.

The results of the assessment calculations are presented in the form of
comparison plots, which show measured values and the corresponding
calculations with the one-dimensional and three-dimensional RELAP5
models.

A comparison of calculated and measured primary system pressures is
presented in Figure 14. The calculated primary system pressures were
similar with both models and in reasonable agreement with the data. The
measured curve contains an inflection point at about 16 s, roughly
corresponding to the initiation of accumulator injection that was more
pronounced than in the calculations. The calculated pressures were also
slightly less than the measured values after 40 s.

Table 2. Calculated and measured sequence of events for LOFT Test L2-5.

  ---------------------------------- ------------ ----------------- -------------------
                                     Time after                     
                                     rupture (s)                    

  Event                              Test         One-dimensional   Three-dimensional
                                                  model             model

  Test initiated                     0.0          0.0               0.0

  Reactor trip signal                0.02 ± 0.01  0.02              0.02

  Quick opening blowdown valves      0.04 ± 0.01  0.04              0.04
  fully opened                                                      

  Cladding temperatures initially    0.91 ± 0.2   0.84              0.28
  deviated from saturation                                          

  Primary coolant pumps tripped      0.94 ± 0.01  1.60^1^           1.60^1^

  Subcooled break flow ended (cold   3.4 ± 0.5    3.7               3.9
  leg)                                                              

  Steam control valve fully closed   9.38 ± 0.05  9.38              9.38

  Pressurizer emptied                15.4 ± 1.0   15.5              15.5

  Accumulator injection initiated    16.8 ± 0.1   15.0              14.3

  HPIS injection initiated           23.90 ± 0.02 23.90             23.90

  Maximum cladding temperature       28.47 ± 0.02 6.0               6.3
  reached                                                           

  LPIS injection initiated           37.32 ± 0.02 37.32             37.32

  Accumulator emptied                49.6 ± 0.1   50.0              50.5

  Core cladding quenched             65 ± 2       49                55
  ---------------------------------- ------------ ----------------- -------------------

1\. The measured voltage and current to the pumps did not drop
instantaneously to zero following the trip of the pumps in the test.
Since the code assumes an instantaneous trip, the pump trip in the
calculations was delayed to coincide with the measured decrease in pump
speed and differential pressure.

![](media/image14.wmf){width="6.019444444444445in"
height="5.017361111111111in"}

#### Figure 14. Calculated and measured primary system pressure for LOFT Test L2-5

Figures 15 and 16 show calculated and measured mass flow rates in the
broken loop cold leg and broken loop hot leg. The measured flow in the
broken loop cold leg was substantially larger than the flow in the
broken loop hot leg, particularly during the first 5 s of the transient.
The fluid upstream of the cold leg break was subcooled for several
seconds while the fluid upstream of the hot leg break was at the
saturation temperature almost immediately after the break, leading to
higher critical flow rates on the cold leg side. Both RELAP5 models
predicted this trend.

![](media/image15.wmf){width="4.961111111111111in"
height="4.108333333333333in"}

Figure 15. Calculated and measured mass flow rates in the broken loop
cold leg for LOFT Test L2-5

![](media/image16.wmf){width="4.81875in" height="3.8152777777777778in"}

Figure 16. Calculated and measured mass flow rates in the broken loop
hot leg for LOFT Test L2-5

The calculated results are generally within the uncertainty of the
measurements and are thus considered to be in excellent agreement with
the data.

A comparison of calculated and measured mass flow rates in the intact
loop hot leg is shown in Figure 17.

![](media/image17.wmf){width="4.634027777777778in"
height="3.3881944444444443in"}

Figure 17. Absolute value of the calculated and measured mass flow rates
in the intact hot leg for LOFT Test L2-5

The instrument measured only the magnitude of the flow rate and not its
direction. Consequently, the absolute values of the flow rates, which
are presented in Figure 17, provide a more direct indication of the
agreement between the calculated and measured results. In the
calculations, the flow in the hot leg was generally towards the steam
generator until 5 s. The flow then reversed, going towards the reactor
vessel due to the pump trip and the corresponding flow coastdown. The
maximum negative flow occurred at about 10 s and was caused by vapor
generation in the steam generator u-tubes, which forced flow towards the
reactor vessel. The draining of the pressurizer also contributed to the
flow from the hot leg to the reactor vessel. Based on the comparisons
shown in Figure 17, a similar flow reversal probably occurred in the
experiment. The trends in the calculations were similar to those
observed in the test except that the measured results were more
oscillatory, particularly between 35 s and 60 s, which roughly
corresponded to the reflooding of the core.

In general, there was little difference between the results from the
one-dimensional and three dimensional models insofar as loop behavior
was concerned. That was not unexpected. Differences were obviously
expected in the vessel.

Figure 18 presents calculated and measured fuel centerline temperatures
in the central, high-powered bundle. The measured temperature decreased
rapidly during the first 5 s of the transient, then remained nearly
constant until 63 s, when the temperature decreased rapidly to near the
saturation temperature of the fluid. The rapid temperature decrease at
63 s is referred to as a quench and indicates that the fuel rod surface
had been wetted by the advancing mixture level during the reflooding of
the core. The calculated temperatures were similar to the measurement
except that the quench occurred earlier than in the test, particularly
in the calculation with the one-dimensional model, and both calculated
temperatures decreased at a faster rate than in the test between 5 s and
the quench time.

![](media/image18.wmf){width="4.4534722222222225in"
height="3.2534722222222223in"}

Figure 18. Calculated and measured fuel centerline temperatures in ring
1, sector 3, level 8 for LOFT Test L2-5

Comparisons between calculated and measured cladding temperatures as a
function of elevation are presented in Figures 19, 20, and 21. The
results correspond to axial levels 6, 7, and 8 of ring 1, sector 2 of
the three-dimensional model. Ring 1 represents most of the central,
high-powered fuel rod bundle, and sector 2 represents the quadrant
connected to the broken loop hot leg. The one- and three-dimensional
models produced similar results, underpredicting the peak cladding
temperature and quenching earlier than in the experiment. However,
quenching behavior in the three-dimensional model more closely matched
the data. This is mainly attributed to the fact that the high-powered
fuel rods were attached to a relatively hot fluid channel in the
three-dimensional model whereas the high-powered fuel rods were attached
to a single, average channel in the one-dimensional model.

Significant radial variations in cladding temperatures were observed in
both the test and the calculation with the three-dimensional model as
shown in Figures 22 and 23. Figure 22 shows measured temperatures in
rings 1 through 3 at level 9. The thermocouples referred to in the
figures were all located at the same elevation and within the same
sector so that the measured differences in results were due to radial
effects. Figure 23 shows the corresponding calculated results with the
three-dimensional model. The calculated and measured temperatures both
show the influence of the radial power profile, with the highest
temperatures in ring 1, the high-powered ring, and the lowest
temperatures in ring 3, the low-powered ring. The radial variation in
cladding temperatures was more pronounced in the test than in the
calculation, primarily because the time of CHF varied more in the test.

The hydraulic responses of the primary and secondary coolant systems
calculated with RELAP5-3D^^ and the three-dimensional input model were
generally in reasonable agreement with LOFT Test L2-5. The results were
generally as good as or better than the results obtained using the
RELAP5/MOD3 or RELAP5-3D^^ one-dimensional models. The calculated
thermal response of the core fuel rods with the three-dimensional model
was also generally similar to that observed in the test. The calculated
peak cladding temperature was 990 K while the measured peak cladding
temperature was 1078 K. The most significant deviations between the
calculated and measured thermal responses were that the calculated peak
cladding temperature occurred earlier than in the test and that the
top-down rewet that was observed near 15 s in the test was not
predicted.

![](media/image19.wmf){width="3.3534722222222224in"
height="2.566666666666667in"}

Figure 19. Calculated and measured cladding temperatures in ring 1,
sector 2, level 6 for LOFT Test L2-5

![](media/image20.wmf){width="3.3472222222222223in"
height="2.386111111111111in"}

Figure 20. Calculated and measured cladding temperatures in ring 1,
sector 2, level 7 for LOFT Test L2-5

![](media/image21.wmf){width="3.1840277777777777in"
height="2.3986111111111112in"}

Figure 21. Calculated and measured cladding temperatures in ring 1,
sector 2, level 8 for LOFT Test L2-5

![](media/image22.wmf){width="4.590972222222222in"
height="3.7270833333333333in"}

Figure 22. Measured cladding temperatures showing radial effects in
sector 2, level 9 for LOFT Test L2-5

![](media/image23.wmf){width="4.68125in" height="3.3333333333333335in"}

Figure 23. Calculated cladding temperatures showing radial effects in
sector 2, level 9 for LOFT Test L2-5

## Multi-Dimensional Neutron Kinetics Model

The multi-dimensional neutron kinetics model in RELAP5-3D^^ is based on
the NESTLE^14^ code developed by Paul Turinsky and co-workers at North
Carolina State University under an INEEL initiative. The NESTLE code
solves the two or four group neutron diffusion equations in either
Cartesian or hexagonal geometry using the nodal expansion method. Three,
two, or one dimensional models may be used. Several different core
symmetry options are available including quarter, half, and full core
options 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 are available. The steady-state eigenvalue
and time dependent neutron flux problems can be solved by the NESTLE
code as implemented in RELAP5-3D^^.

The few group neutron equations are spatially discretized using the
Nodal Expansion Method (NEM). Quartic or quadratic polynomial expansions
are used for the transversely integrated fluxes in the Cartesian and
hexagonal geometries respectively. Discontinuity factors 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 neutron precursor equations are integrated
analytically assuming a linear variation of the neutron flux.

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

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

The compositions are mapped onto the kinetic nodes using user supplied
composition maps. In this way a small number of compositions and a small
number of RELAP5 volumes and heat structures may be mapped to a large
number of kinetics solution nodes. Discontinuity factors are modeled
using the neutron cross section model so that they may change as the TH
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 any number of control rods. Control rods may be inserted
from either the top or bottom of the reactor. The control rods may be
grouped into banks, with the banks having the same insertion depth and
velocity. The position of a control rod bank can also be determined from
a control variable or a general table. Scram signals may be initiated by
RELAP5 trips.

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

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

### Verification

The implementation of the NESTLE neutron kinetics has been verified by
the simulation of the NEACRP^15^ three dimensional benchmark
problems^16^. A series of three PWR rod ejection accidents from Hot Zero
Power and Hot Full Power were proposed as a benchmark by the NEACRP.
Series A, the ejection of a central rod, and series B, the ejection 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 in Figure 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 parallel pipes as shown in Figure 25. Each pipe was described using a
series of heat structures and control volumes to model the fuel and
coolant from a single assembly. A separate inlet reservoir was provided
for each pipe and was set at a constant pressure, temperature, and boron
concentration. Each pipe was connected to its reservoir using a
time-dependent junction that provided identical, constant flow to each
assembly. The outlet of each pipe was connected to an outlet reservoir
using a series of branches. Since a maximum of nine junctions per branch
is allowed, additional branches were used to accommodate the total of 47
pipes of the quarter core model.

During the course of the analysis, various thermal‑hydraulic mesh
structures were examined. The finest axial mesh used for the pipes
corresponds to that prescribed in the NEACRP problem with the exception
that the three smallest nodes at the bottom and top of the core were
combined into a single thermal‑ hydraulic node. The original NEACRP mesh
structure for the neutronic solution was retained. This provided for a
total of 14 axial thermal‑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 temperature rise
anticipated in the reflector.

![](media/image24.wmf){width="4.603472222222222in"
height="3.490972222222222in"}

![](media/image25.wmf){width="4.7340277777777775in" height="3.49375in"}

###### Figure 24. Initial Core Configuration for Series A and B Rod Ejection Transients

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 RELAP5 consistent with that prescribed in the
NEACRP problem,

$$T = (1 - \alpha)T_{F,C} + {\alpha T}_{F,S}$$

a very small central and peripheral fuel pellet mesh was used such that
the area of these regions corresponds to the weighting of the surface
temperature, *T~F,S~*, and centerline fuel temperature, *T~F,C~*
specified in the benchmark problem. All other fuel pellet meshes were
equidistant.

![](media/image26.wmf){width="4.429861111111111in"
height="4.095833333333333in"}

### Figure 25. RELAP5-3D Core Nodalization for Test Problem

Because the benchmark problem specifies gap conductance and RELAP5 only
accepts gap

conductivity, the following relation was used:

$$k = hDr$$

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

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

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

#### Table 3. Time Step Sizes Used in the Analyses

  ----------------------------- -------------------- ---------------------
  Case                          0 to 1 second        1 to 5 seconds

  Zero Power (A1,B1)            0.001                0.01

  Full Power (A2,B2)            0.01                 0.05
  ----------------------------- -------------------- ---------------------

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

Table 4. Summary of the RELAP5 NEACRP Benchmark Calculation Results

  ------------- ------- ------- ------- ------- ------- ------- ------- -------- ------- ------- --------
                A1                      B1                      A2                       B2      

  Steady State                                                                                   

  Critical      563.4   -4.6    1253            -1.4    1154            -6.8     1183            -4.8
  Boron Conc.                                                                                    
  (ppm)                                                                                          

  Assembly      2.865   0.3%    1.926           -0.2%   2.203           -0.8%    2.101           -0.4%
  Peaking                                                                                        
  Factor                                                                                         

  Max. Fuel     286.0   0.0%    286.0           0.0%    1612            -3.6%    1528            -3.1%
  Temp. (deg.                                                                                    
  C)                                                                                             

  Rod Worth     820     -0.2%   836             0.6%    91              1.7%     99              -0.1%
  (pcm)                                                                                          

  @ Time of                                                                                      
  Maximum Power                                                                                  

  Time (sec)    0.611   0.046   0.519           0.004   0.100           -0.020   0.110           -0.010

  Relative Core 0.911   -22%    2.252           -7.8%   1.085           0.5%     1.063           0.0%
  Power                                                                                          

  Maximum Fuel  337.2   -2.3%   325.5           -0.9%   1613            -3.5%    1528            -3.1%
  Temp. (deg.                                                                                    
  C)                                                                                             

  @ Time = 5                                                                                     
  sec                                                                                            

  Relative Core 0.191   -2.2%   0.317           -1.1%   1.036           0.0%     1.038           0.0%
  Power                                                                                          

  Max. Fuel     650.0   -3.2%   549.9           -1.5%   1632            -3.5%    1539            -3.1%
  Temp. (deg C)                                                                                  
  ------------- ------- ------- ------- ------- ------- ------- ------- -------- ------- ------- --------

As indicated in Table 4, the RELAP5 steady state results for the hot
zero power cases are in reasonably good agreement with the reference
results. RELAP5 does show a slight negative bias in the prediction of
the critical boron concentration for both cases A1 and B1. Also note
that RELAP5 is consistently lower than the reference in its prediction
of the power in the rodded locations for the steady‑state. This is
consistent with the negative bias in the critical boron prediction since
a greater estimation of the total rod worth would lead to a lower
critical boron concentration. The negative bias is larger in case 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 the reference result. As shown in Figure 26, agreement of
the transient results is excellent for 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 magnitude
22.1% lower than the reference calculation.

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

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

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

![](media/image27.wmf){width="6.219444444444444in"
height="3.782638888888889in"}

![](media/image28.wmf){width="6.110416666666667in"
height="4.189583333333333in"}

###### Figure 26. Transient Core Power for Rod Ejection from Hot Zero Power

##### 

![](media/image29.wmf){width="5.997916666666667in"
height="3.5729166666666665in"}

## ![](media/image30.wmf){width="5.997222222222222in" height="4.118055555555555in"}

#### Figure 27. Transient Core Power for Rod Ejection from Hot Full Power

#######  {#section-2 .unnumbered}

####### Steady-State

+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 5     | PAN   |       |       |       | 0.293 | 0.354 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | REL   |       |       |       | 0.286 | 0.359 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | \% D  |       |       |       | -2.4  | 1.4   |       |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 6     |       |       |       | 0.752 | 0.533 | 0.497 | 0.285 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.753 | 0.534 | 0.501 | 0.290 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.1   | 0.2   | 0.8   | 1.8   |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 7     |       |       | 0.545 | 0.757 | 0.393 | 0.380 | 0.206 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | 0.529 | 0.757 | 0.382 | 0.382 | 0.202 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | -2.9  | 0.0   | -2.8  | 0.5   | -1.9  |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 8     |       | 0.964 | 0.867 | 1.000 | 0.745 | 0.301 | 0.294 | 0.226 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.964 | 0.867 | 1.000 | 0.745 | 0.292 | 0.296 | 0.230 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.0   | 0.0   | 0.0   | 0.0   | -3.0  | 0.7   | 1.8   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 9     | 0.533 | 0.793 | 0.575 | 0.945 | 0.951 | 0.527 | 0.214 | 0.285 |
|       |       |       |       |       |       |       |       |       |
|       | 0.516 | 0.792 | 0.557 | 0.943 | 0.951 | 0.528 | 0.209 | 0.289 |
|       |       |       |       |       |       |       |       |       |
|       | -3.2  | -0.1  | -3.1  | -0.2  | 0.0   | 0.2   | -2.3  | 1.4   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|       | I     | J     | K     | L     | M     | N     | O     | P     |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+

####### At Time of Maximum Power

+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 5     | PAN   |       |       |       | 0.128 | 0.150 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | REL   |       |       |       | 0.125 | 0.151 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | \% D  |       |       |       | -2.3  | 0.7   |       |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 6     |       |       |       | 0.362 | 0.242 | 0.214 | 0.120 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.362 | 0.242 | 0.214 | 0.121 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.0   | 0.0   | 0.0   | 0.8   |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 7     |       |       | 0.316 | 0.390 | 0.188 | 0.169 | 0.088 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | 0.307 | 0.390 | 0.183 | 0.170 | 0.086 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | -2.8  | 0.0   | -2.7  | 0.6   | -2.3  |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 8     |       | 0.790 | 0.562 | 0.540 | 0.371 | 0.140 | 0.126 | 0.093 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.790 | 0.561 | 0.540 | 0.370 | 0.136 | 0.127 | 0.094 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.0   | -0.2  | 0.0   | -0.3  | -2.9  | 0.8   | 1.1   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 9     | 1.000 | 0.778 | 0.390 | 0.513 | 0.474 | 0.248 | 0.093 | 0.117 |
|       |       |       |       |       |       |       |       |       |
|       | 1.000 | 0.777 | 0.378 | 0.513 | 0.474 | 0.248 | 0.090 | 0.118 |
|       |       |       |       |       |       |       |       |       |
|       | 0.0   | -0.1  | -3.1  | 0.0   | 0.0   | 0.0   | -3.2  | 0.9   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|       | I     | J     | K     | L     | M     | N     | O     | P     |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+

####### At Final Time t = 5 seconds

+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 5     | PAN   |       |       |       | 0.143 | 0.168 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | REL   |       |       |       | 0.140 | 0.172 |       |       |
|       |       |       |       |       |       |       |       |       |
|       | \% D  |       |       |       | -2.1  | 2.4   |       |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 6     |       |       |       | 0.392 | 0.266 | 0.239 | 0.135 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.393 | 0.268 | 0.241 | 0.138 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       |       | 0.3   | 0.8   | 0.8   | 2.2   |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 7     |       |       | 0.333 | 0.417 | 0.205 | 0.188 | 0.099 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | 0.323 | 0.418 | 0.199 | 0.189 | 0.097 |       |
|       |       |       |       |       |       |       |       |       |
|       |       |       | -3.0  | 0.2   | -2.9  | 0.5   | -2.0  |       |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 8     |       | 0.802 | 0.581 | 0.569 | 0.397 | 0.153 | 0.142 | 0.106 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.802 | 0.581 | 0.570 | 0.398 | 0.149 | 0.143 | 0.108 |
|       |       |       |       |       |       |       |       |       |
|       |       | 0.0   | 0.0   | 0.2   | 0.3   | -2.6  | 0.7   | 1.9   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
| 9     | 1.000 | 0.785 | 0.403 | 0.540 | 0.505 | 0.269 | 0.104 | 0.134 |
|       |       |       |       |       |       |       |       |       |
|       | 1.000 | 0.785 | 0.391 | 0.541 | 0.506 | 0.271 | 0.102 | 0.136 |
|       |       |       |       |       |       |       |       |       |
|       | 0.0   | 0.0   | -3.0  | 0.2   | 0.2   | 0.7   | -1.9  | 1.5   |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|       | I     | J     | K     | L     | M     | N     | O     | P     |
+-------+-------+-------+-------+-------+-------+-------+-------+-------+

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

## BPLU Matrix Solver

The Border Profiled Lower Upper (BPLU) matrix solver is used to
efficiently solve sparse 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 most efficient
for solving systems that correspond to networks, such as pipes, but is
efficient for any system that it can permute into border-banded form.

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

![](media/image1.wmf)Figure 29. 3dflow problem

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

Table 5. Comparison of Run Times for 3dflown.i Using Default and BPLU
Solvers

+-----------------+-----------------+-----------------+-----------------+
| Case            | Default Solver  | BPLU Solver     | Ratio           |
|                 |                 |                 |                 |
|                 | (CPU secs.)     | (CPU secs.)     |                 |
+-----------------+-----------------+-----------------+-----------------+
| 1               | 7.179600        | 2.437199        | 2.94584         |
+-----------------+-----------------+-----------------+-----------------+
| 2               | 7.142498        | 2.110301        | 3.38459         |
+-----------------+-----------------+-----------------+-----------------+
| 3               | 6.903399        | 2.718000        | 2.53988         |
+-----------------+-----------------+-----------------+-----------------+
| 4               | 6.141501        | 2.421701        | 2.53603         |
+-----------------+-----------------+-----------------+-----------------+
| 5               | 5.512910        | 2.117097        | 2.60399         |
+-----------------+-----------------+-----------------+-----------------+
| 6               | 5.817507        | 2.698403        | 2.15591         |
+-----------------+-----------------+-----------------+-----------------+
| 7               | 6.166502        | 2.432401        | 2.53515         |
+-----------------+-----------------+-----------------+-----------------+
| 8               | 7.405095        | 2.116199        | 3.49924         |
+-----------------+-----------------+-----------------+-----------------+
| 9               | 6.396306        | 2.696500        | 2.37208         |
+-----------------+-----------------+-----------------+-----------------+

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

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

Two considerations lead to further efficiencies. First, the block
structure discussed above leads to a profile rather than a solid band
structure for which zero operations can be avoided by using a variable
row length during row reduction. Profiled-structures also arise from
connectivity relationships among the variates. Therefore, the BPLU
algorithm uses a profile, rather than a banded, solver. Second, a
nearly-banded matrix can frequently be decomposed into the sum of a
banded and an outrigger matrix in more than one way. Selection of the
bandwidth for the decomposition is tied to reducing the amount of work
associated with solving the reorganized system.

## OTHER IMPROVEMENTS IN RELAP5-3D^^

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

-   Graphical User Interface (beta-test)

-   More implicit momentum equation permitting larger time step sizes

-   Multidimensional heat conduction model for modeling graphite
    structures in RBMK reactors

-   Windows '95/98/NT version

Ongoing and planned DOE funding will provide additional improvements:

-   Ability to couple to other codes (e.g. CFD)

-   Random-access plot/restart file

-   Fuel swelling model

## RELAP5-3D^^ AVAILABILITY

## The INEEL will make the RELAP5-3D^^ code available to domestic and international organizations that join the International RELAP5 Users Group (IRUG). This is a subscriber-funded group that takes part in guiding the future development and maintenance of the code. Organizations interested in obtaining the code will be given the opportunity 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 levels of participation. More information on joining IRUG can
be found at http://www.inel.gov/relap5/irug.htm.

######## **REFERENCES**

1.  J. A. Findley and G. L. Sozzi, "BWR Refill-Reflood Program -- Model
    Qualification Task Plan," EPRI NP-1527, NUREG/CR-1899, GEAP-24898,
    October 1981.

2.  T. M. Anklam, R. J. Miller, M. D. White, "Experimental Investigation
    of Uncovered-Bundle Heat Transfer and Two-Phase Mixture Level Swell
    Under High-Pressure and Low Heat Flux Conditions," NUREG/CR-2456,
    ORNL-5848, March 1982.

3.  K. Carlson, R. Riemke, R. Wagner, J. Trapp, "Addition of
    Three-Dimensional Modeling," RELAP5/TRAC-B International Users
    Seminar, Baton Rouge, LA, November 4-8, 1991.

4.  R. Riemke, "RELAP5 Multi-Dimensional Constitutive Models,"
    RELAP5/TRAC-B International Users Seminar, Baton Rouge, LA, November
    4-8, 1991.

5.  K. Carlson, R. Riemke, R. Wagner, "Theory and Input Requirements for
    the Multi-Dimensional Component in RELAP5 for Savannah River Site
    Thermal-Hydraulic Analysis," EGG-EAST-9878, Idaho National
    Engineering Laboratory, July, 1992.

6.  K. Carlson, C. Chou, C. Davis, R. Martin, R. Riemke, R. Wagner,
    "Developmental Assessment of the Multi-Dimensional Component in
    RELAP5 for Savannah River Site Thermal-Hydraulic Analysis,"
    EGG-EAST-9803, Idaho National Engineering Laboratory, July, 1992.

7.  K. Carlson, C. Chou, C. Davis, R. Martin, R. Riemke, R. Wagner, R.
    Dimenna, G. Taylor, V. Ransom, J. Trapp, "Assessment of the
    Multi-Dimensional Component in RELAP5/MOD2.5, Proceedings of the
    5^th^ International Topical Meeting on Nuclear Reactor
    Thermal-Hydraulics, Salt Lake City, Utah, USA, September 21-24,
    1992.

8.  P. Murray, R. Dimenna, C. Davis, "A Numerical Study of the Three
    Dimensional Hydrodynamic Component in RELAP5/MOD3, " RELAP5
    International Users Seminar, Boston, MA, USA, July, 1993.

9.  G. Johnsen, "Status and Details of the 3-D Fluid Modeling of
    RELAP5," Code Application and Maintenance Program Meeting, Santa Fe,
    NM, October, 1993.

10. E. Tomlinson, T. Rens, R. Coffield, "Evaluation of the RELAP5/MOD3
    Multidimensional Component Model, RELAP5 International Users
    Seminar, Baltimore, MD, August 29 -- September 1, 1994.

11. A. Shieh, "1D to 3D Connection for the Nearly-Implicit Scheme," INEL
    Report, June, 1997.

12. K. Carlson, "1D to 3D Component Connection for the Semi-Implicit
    Scheme," INEL Report, June, 1997.

13. C. Davis, "Assessment of the RELAP5 Multi-Dimensional Component
    Model Using Data from LOFT Test L2-5," INEEL Report LDRD 3101,
    September, 1996.

14. R. M. Al-Chalabi, et al, "NESTLE: A Nodal Kinetics Code,"
    Transactions of the American Nuclear Society, Volume 68, June, 1993.

15. H. Finnemann and A. Galati, "NEACRP 3-D LWR Core Transient Benchmark
    -- Final Specifications," NEACRP-L-335 (Revision 1), January, 1992.

16. J. L. Judd, W. L. Weaver, T. Downar, J. G. Joo, "A Three Dimensional
    Nodal Neutron Kinetics Capability for RELAP5," Proceedings of the
    1994 Topical Meeting on Advances in Reactor Physics, Knoxville, TN,
    April 11-15, 1994, Vol. II, pp 269-280.

17. H. Finnemann, et al, "Results of LWR Core Transient Benchmarks,"
    Proceedings of the Joint International Conference on Mathematical
    Methods and Supercomputing in Nuclear Applications, Vol. 2, pg. 243,
    Kernforschungszentrum, Karlsruhe, Germany, April, 1993.

##