Subchannel analysis on advanced computer architectures

Math1 Comput. Modelling, Vol. 11, pp. 47-52, Printed in Great Britain

1988

MODELLING

ON

SUBCHANNEL ANALYSIS ARCHITECTURES

John

A. Turner,

0895-7177/88 $3.00 + 0.00 Pergamon Press plc

J. Michael

PARALLEL

COMPUTERS

ON ADVANCED

COMPUTER

Doster

North Carolina State University Department of Nuclear Engineering Box 7909, Raleigh, N.C. 27695-7909

Abstract. Heat removal in light water nuclear reactor cores involves turbulent two-phase flow parallel to small diameter fuel rods. This situation results in flow which is highly three-dimensional, and accurate modelling is essential to the design and safety analysis of nuclear reactor systems. The present work is part of an ongoing effort in the Nuclear Engineering Department at North Carolina State University to evaluate the use of advanced computer architectures for nuclear engineering problems including neutronics and thermal-hydraulics. The work presented here concentrates on the solution of the seven-banded block-tridiagonal linear systems arising from the modelling of multi-phase flow in fuel rod bundles. Keywords.

Array

ear systems; two-phase

processors;

parallel

banded

processing;

flow; vector

linear

subchannel

The BWR system sisting of a nuclear

lin-

thermal-hydraulics;

Subchannel analysis codes such as COBRA and VIPRE are often used to help evaluate nuclear reactor core safety limits in both PWRs and BWRs during normal operating conditions, anticipated transients, narios. These codes

to BWRs.

is a closed steam cycle core located in a reactor

conves-

and assumed accident are used in conjunction

scewith

reactor system thermal-hydraulic codes such as TRAC and RELAP, and both types of codes require significant run times even on large computers. Vectorized and multi-tasked versions of the systems codes have reduced run times on both single and multi-CPU supercomputers, and some have achieved much-better-than-realtime simulation (Moore, 1982; Makowitz, 1986; Ishiguro and Harada, 1986; Sills and Doster, 1986a). The objective of our work is to develop and evaluate new algorithms to speed up subchannel analysis codes. The motivation is twofold: 1) solution on large machines will come cheaper, and 2) problems which could previously only be solved on large machines might become solvable on advanced engineering work-

sel, a turbine-generator set, a condenser, and a feedwater supply system. System pressure is approximately 1000 psi (6900 kPa), and output can be as much as 1400 MW electric (4000 MW thermal). Water enters the reactor vessel subcooled, flows upward parallel to 12 ft (3.66 m) long, small diameter fuel rods, and exits the core as saturated steam with a quality of about 14%, and a void fraction of about 75%. The steam then goes through moisture separators and dryers, exits the reactor vessel, proceeds to the turbine-generator set, the condenser, and is returned to the reactor vessel by the feedwater supply system. MCM

block-tridiagonal

It should be noted that flow in the core is timedependant, highly turbulent, three-dimensional, and two-phase (steam and water).

Two types of light water nuclear reactor systems, Pressurized Water Reactors (PWRs) and Boiling Water Reactors (BWRs), are commonly in operation. Though the models and algorithms developed here are applicable to either system, will be limited

analysis;

pipelines.

INTRODUCTION

discussion

systems;

cmt-c 47

48

Proc. 6th Int. Conf. on Mathematical

stations (microcomputers and/or vector pipelines cessors). These machines price/performance ratio

A drift-flux subchannel analysis code called SWIRL has been developed to allow the evaluation of various approaches to implementation of this problem on advanced computer architectures.

HYDRODYNAMIC

The microscopic fluids core can be written:

MODEL

equations

for flow in the

Continuity

$+O.(pl;J=0 Internal

(1)

Energy _-_

apu x+

v

.(pu

;,

t

with multiple CPUs or attached array prohave a very attractive (Zee, et al., 1986).

= -P

G.

v + a.

4

(2)

Modeliing

11 V V,

= = = =

W

=

time internal energy vector velocity axial velocity component lateral velocity component

These equations could simply be applied to the two-phase mixture, in which case a homogeneous equilibrium model would result. At the other extreme, each of the above equations could be applied to both the liquid phase and the vapor phase (or the mixture and one of the phases) to obtain a full two-fluid model (Lahey and Moody, 1977). Though the two-fluid model most closely describes the physical system, it presents prodigious difficulties due to increased computational requirements and uncertainties in interfacial interaction terms. We have chosen a model between these two extremes, the four equation drift-flux model, which eliminates the explicit representation of the interfacial terms while preserving the velocity difference between the The drift-flux model uses a phases (slip). semi-empirical phenomenological model based on bubble rise velocity to model the relative velocity of the two phases. This allows adequate modelling of slip with the least additional complexity, and is a fairly common and Reed, 1978; Bessho and

Axial

Momentum

Sami,

6

.(pV,

G)

= -g+

l

(3)

v

Lateral

-

jpw

V)

= $+

0.

a,

(4)

l

-

P

= density

= UT = St = P = Q =

axial shear stress lateral shear stress gravitational constant pressure heat transferred to fluid

flow

with

axial

This allows the use of uncoupled and lateral momentum equations is clearly good except for severe flow blockages.

axial and axial

Thermodynamic

where: 6,

assumptions

is small

compared

flow

Momentum

apw at+

simplifying

oz +pgr

t i

Lateral

(Liles 1985;

1986; and others).

Thus far, the following have been made: 8PVz -+ at

approach Uchikawa,

l

equilibrium

This allows the use of a single mixture energy equation, and is valid except for rapid transients.

The Zuber-Findlay void-quality in the axial direction -

model

holds

This is the drift-flux assumption, allows the use of a single mixture

and axial

Proc. 6th Int. ConJ on Mathematical

momentum equation drift flux correlation.

along

with

Mixture

DISCRETIZATION

the

If we further assume homogeneous flow between subchannels and that turbulent mixing yields no net mass exchange between subchannels, the following equation set results:

49

ModeIling

SCHEME

The equation set is discretized staggered mesh as shown below:

spatially

on a

Continuity

(5) FIG. Mixture

Internal

$+

+ .(pff ;,

-6.

Mixture

8PVz

v-

al+

= -PG.

;

+

(6) Donor cell averaging and upwind difference schemes are used where appropriate, and the equation of state is linearized using a first order Taylor series expansion. This linearization imposes a timestep restriction due to density truncation error.

Momentum

.(pV,

G,

= -g+

G

*&

+

(7)

-i?.{~(&)~j,}+pgz

Time discretization is semi-implicit. Time derivatives are expanded in a forward difference, with pressures and velocities (or mass fluxes) evaluated

Mixture

Lateral

apw al+

v

In addition correlation, p(pu,

the

new

timestep

and

Momentum

earity

is achieved

.(PW

restricts timestep size by the material Courant limit, i.e., the timestep must be less than the minimum transport time through any cell (axial or lateral):

t;)

=

g+

;s.

;,

(8)

At<

= drift flux = total volumetric

at the past timestep,

convected

evaluated

1 denotes liquid phase g denotes vapor phase a = void fraction i

at

terms

where:

‘g 3D

Mesh

Thermodynamic properties such as pressure, density, and internal energy are defined at cell centers (k - 1, k, k + 1, etc.) and velocities (or mass fluxes) are defined at cell boundaries (k - 3, k + i, etc.).

(~(,,-.&~}+6.Y

Axial

1 - Staggered

Energy

so that

in the new time variables.

linThis

&$,I,.,, 7l~l,,,l,,,tg)

where: Ar flux

to these equations and the drift flux an equation of state of the form p = P) is required for closure.

Al AL

= lateral spacing = timestep = axial cell length

can be reduced The resulting set of equations to a single linear system in the spatial pressure

Proc.

50

6th Int. Conf. on kfathematical

distribution: ai ’k-1 I Pf:e: ’

+ bikPftAt

+ c

$ ai,k++Pf:h+:

t+At

%kp,k

=

as a minor flow blockage would all subchannels, thus requiring full problem.

t

S,,k

(10)

where: i denotes subchannel of interest j denotes an adjacent subchannel k denotes axial cell ij denotes from subchannel i to subchannel

j

each timestep

involves

of coefficient

2.

Solution

3.

Back substitution

4.

Timestep Courant

Most is

of linear

Timestep

matrix system

setup

control

tests

and back and/or involve

Solution

LINEAR

SOLUTION

Methods for solving linear systems fall naturally into two classes, direct and iterative. Direct methods give an exact solution (in the absence of roundoff error), but storage requirements can

substitution

multi-taskable. repeated

SYSTEM

lateral error)

condi-

tionals, but are also parallelizeable. Solution the linear system requires more scrutiny.

of

Consider first the size of the system. If Nchan denotes the number of subchannels and Ncell denotes the number of axial cells, the size of the system is:

(Nchan)(Ncc~l) by (Nchan)(Ncen) A typical BWR fuel assembly has sixty-four fuel rods in an 8 by 8 array, yielding eighty-one subchannels. Using twenty-five axial cells, the linear system is 2025 by 2025 (twenty-five 8 1 by 81 blocks). Symmetry to reduce

1 - System

i

variables

control tests (axial and limit and density truncation

vectorizable

TABLE

It should be noted that as other portions of the code are vectorized, this fraction will increase, so it is on the solution of the linear system that we are currently concentrating our efforts.

for remaining

of the matrix

highly

flow in of the

the following

steps:

1. Setup

perturb solution

Timing studies on a VAXstation II show that the solution of the iinear system quickly becomes a major portion of the computational effort. Table 1 shows the fraction of total CPU time spent solving the linear system in SWIRL for a fixed number of axial cells. The Gaussian elimination scheme (GE) is standard, but recognizes zeros below the diagonal and uses no pivoting, and the block-tridiagonal solver (BT) will be discussed later.

A back substitution step for the remaining system variables (p, pu,pV, pw) is thus required The coeffifollowing the pressure solution. cient matrix is clearly seven-banded and blocktridiagonal. So note that

Modelling

considerations can sometimes be used the system size, but a situation such

be greater, tages suffer

while

iterative

methods

offer advan-

in storage and roundoff control, but from slow or uncertain convergence.

can

Though iterative methods are inherently parallel, it has been shown that fluids equation solutions of this type are sensitive to matrix solution errors. Thus, tight convergence criteria are required for iterative methods, and efficient direct methods can become superior (Sills and 1986b). For this reason, we are initially ering direct methods.

Doster, consid-

Let Ak, Bk, and Ck be the klh main diagonal, sub-diagonal, and super-diagonal blocks of the coefficient matrix, respectively. Also let Xk and Sk be the Icth blocks of the solution and source vectors, respectively. Note that Ak, Bk, and Ck are Nchan by Nchan arrays and Xk and Sk are N than by 1 arrays.

Proc. 6th Int. Co&

on Mathematical

Using this notation, consider the following algorithm (an extension of the familiar tridiagonal solver, where n is the number of blocks, which is Ncch):

Dl = A;‘Ci El = A;‘& for Ic=2ton-1: Dk = (Ak - BkDk_,)-lCk Ek = (Ak - BkDk_1))‘(Sk - BIIEL-I) - E&E,-1) E, = (A, - B,D,_r)-‘(S, X* = E, for k=n-1to1: Xk

FIG.

=

Ek

-

Ddk+l

2 - Block-Tridiagonal

The following

observations

Algorithm

about

the algorithm

are significant: 1.

The algorithm generates no fill outside the off-diagonal blocks, so storage requirements are minimal,

2.

The operations vector pipeline tures,

NO.

Modelling

TABLE 2 Scalar BT TBT

F Ch. 16 25 36

0.15 0.50 1.40 3.40 7.38 14.67

c RAY BT

0.06 0.20 0.56 1.34 2.08 5.73

X-MP

vector

i

TBT

0.04 0.09 0.20 0.42 0.76 1.35

0.02 0.04 0.08 0.16 0.29 0.51

1 t\

r; lesults Speedup BT TBT

I1

3.75 5.56 7.00 8.10 9.71 10.87

3.00 5.00 7.00 8.38 9.93 11.24

vo methods, consider first In comparing the the approximate operation counts. BT requires a matrix inverse, two matrix-matrix multiplies, and three matrix-vector multiplies, for approximately 3n3 + 3n2 operations for each block. In TBT the two matrix-matrix multiplies become n2 operations instead of n3, for an approximate operation count of n3 + 57~’ for each block. The ratio of n3 terms is thus 3, and we see an approach to this limit in both scalar and vector modes. The vector/scalar speedups indicate that some chaining confirms that the algorithm CRAY-like architectures.

for both methods is occuring. This is quite suitable for

DISCUSSION

3.

required or array

are well suited to processor architec-

The diagonal blocks are five-banded matrices and the off-diagonal blocks are diagonal matrices.

RESULTS

The algorithm CPU of a CRAY forms. The first BT in Table 1. of the algorithm

has been implemented on one X-MP/48 in FORTRAN in two we have already referred to as BT is a general implementation in Figure 2. The second, called

TBT, is the same algorithm, the fact that the off-diagonal matrices.

but makes use of blocks are diagonal

Table 2 shows results for a fixed number of blocks (axial cells) with block size (number of subchannels) ranging from 16 to 81. Shown are CPU times in seconds in scalar and vector modes, and vector/scalar speedups, where: Speedup

=

Scalar

mode

CPU

time

vector

mode

CPU

time

(11)

These preliminary results of the form described here and efficiently on vector cessor architectures. The very little timization, ter results TBT. Finally,

indicate that systems can be solved rapidly pipeline or array proalgorithms underwent

in the way of architecture-specific opand we feel we can achieve even betafter fine-tuning the vector version of

if we

can

achieve

similar

results

on

smaller machines such as the VAXstation II with attached array processor(s) or other workstations, this work could have application in advanced engineering design and analysis workstations, improved thermal-hydraulic models for training simulators, and/or advanced operator aids and on-line control systems.

REFERENCES

Bessho, Y. and S. Uchikawa (1985). Subchannel Analysis Program for Boiling Water Reactor Fuel Bundles Based on Five Conservation Equations of Two-Phase Flow. Nucl. Tech., 60, 210-217.

52

Proc. 6th Int. Conf. on Mathematical Modelling

Lahey, R. T., Jr. and F. J. Moody (1977). The Thermal-Hydraulics of a Boiling Water Nuclear Reactor. Amer. Nucl. Sot. Liles, D. R. and Wm. H. Reed (1978). A SemiImplicit Method for Two-Phase Fluid Dynamics. J. Comp. Phys., 26, 390-407.

Sills, E. D. and J. M. Doster (1986). Algorithms for the Solution of the Mixture Drift-Flux Equations on Multi-CPU Vector Pipeline Machines. Trans. Amer. Nucl. sec., 53, 258.

Experiments with the RECode. w

Sills, E. D. and J. M. Doster (1986). Evaluation of Algorithms for the Solution of the Drift-Flux Equations on Advanced Computers. Proc. 2”d Intl. Conf. on Simulation Methods in Nucl. Eng., 678.

Sami, S. (1986). A Transient Two-Phase Velocity Difference Model for Drift Calculation in CANDU Thermohydraulic Codes. Nucl. Tech 75, 283-297. -3

Zee, S. K., E. D. Sills, P. J. Turinsky, and J. M. Doster (1986). Solution of Neutronic and Thermal-Hydraulic Problems on an Engineering Workstation. Trans. Amer. Nucl. sot., 53, 270.

Makowitz, H. (1986). Numerical in Concurrent Multiprocessing LAP5 Nuclear Reactor System Sci. Eng., 92, 136.