NUREG/IA-0146

International

Agreement Report

Implementation and Assessment of

Improved Models and Options in

TRAC-BF1

Prepared by

G. Th. Analytis

Laboratory for Thermal Hydraulics

Paul Scherrer Institute

CH-5232 Villigen PSI

Switzerland

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

October 1998

Prepared as part of

The Agreement on Research Participation and Technical Exchange

under the International Code Application and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission
AVAILABILITY NOTICE

Availability of Reference Materials Cited in NRC Publications

NRC publications in the NUREG series, NRC regulations, and Title 10, Energy, of the Code of Federal

Regulations, may be purchased from one of the following sources:

1. The Superintendent of Documents

U.S. Government Printing Office

RO. Box 37082

Washington, DC 20402-9328

<http://www.access.gpo.gov/su-docs>

202-512-1800

2. The National Technical Information Service

Springfield, VA 22161-0002

<http://www.ntis.gov/ordemow>

703-487-4650

The NUREG series comprises (1) technical and administrative reports, including those prepared for

international agreements, (2) brochures, (3) proceedings of conferences and workshops, (4) adjudications and other issuances of the Commission

and Atomic Safety and Licensing Boards, and

(5) books.

A single copy of each NRC draft report is available

free, to the extent of supply, upon written request

as follows:

Address: Office of the Chief Information Officer

Reproduction and Distribution

Services Section

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

E-mail: <GRW1 @NRC.GOV>

Facsimile: 301-415-2289

A portion of NRC regulatory and technical information is available at NRC's World Wide Web site:

<http://www.nrc.gov>

All NRC documents released to the public are available for inspection or copying for a fee, in paper,

microfiche, or, in some cases, diskette, from the

Public Document Room (PDR):

NRC Public Document Room

2121 L Street, N.W., Lower Level

Washington, DC 20555-0001

<http://www.nrc.gov/NRC/PDR/pdrl .htm>

1-800-397-4209 or locally 202-634-3273

Microfiche of most NRC documents made publicly

available since January 1981 may be found in the

Local Public Document Rooms (LPDRs) located in

the vicinity of nuclear power plants. The locations

of the LPDRs may be obtained from the PDR (see

previous paragraph) or through:

<http://www.nrc.gov/NRC/NUREGS/

S R1350/V9/lpdr/html>

Publicly released documents include, to name a

few, NUREG-series reports; Federal Register notices; applicant, licensee, and vendor documents

and correspondence; NRC correspondence and

internal memoranda; bulletins and information notices; inspection and investigation reports; licensee event reports; and Commission papers and

their attachments.

Documents available from public and special technical libraries include all open literature items, such

as books, journal articles, and transactions, Federal Register notices, Federal and State legislation,

and congressional reports. Such documents as

theses, dissertations, foreign reports and translations, and non-NRC conference proceedings may

be purchased from their sponsoring organization.

Copies of industry codes and standards used in a

substantive manner in the NRC regulatory process

are maintained at the NRC Library, Two White Flint

North, 11545 Rockville Pike, Rockville, MD

20852-2738. These standards are available in the

library for reference use by the public. Codes and

standards are usually copyrighted and may be

purchased from the originating organization or, if

they are American National Standards, fromAmerican National Standards Institute

11 West 42nd Street

New York, NY 10036-8002

<http://www.ansi.org>

212-642-4900

DISCLAIMER

This report was prepared under an international cooperative

agreemert for the exchange of technical informatio. Neither

the United States Government nor any agency thereof, nor any

of their employes, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for any third

party's use, or the results of such use, of any intormation, apparatus, product, or process disclosed in this report, or represents

that its use by such third party would not infringe privately

owned rghts.
NUREG/IA-0146

International

Agreement Report

Implementation and Assessment of

Improved Models and Options in

TRAC-BF1

Prepared by

0. Th. Analytis

Laboratory for Thermal Hydraulics

Paul Scherrer Institute

CH-5232 Villigen PSI

Switzerland

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

October 1998

Prepared as part of

The Agreement on Research Participation and Technical Exchange

under the International Code Application and Maintenance Program (CAMP)

Published by

U.S. Nuclear Regulatory Commission

ABSTRACT

A summary of modifications and options introduced in TRAC-BFl is presented and is

shown that the predicting capabilities of the modified version of the code are greatly

improved. These changes include the introduction of a different heat transfer package

during reflooding, the implementation of a simple single-phase limit procedure for forcing

the two phases to acquire the same velocity if one phase disappears, a close assessment of

the annular flow interfacial shear correlation, implementation of a simple radiation model

which seems to alleviate some numerical-oscillations problems induced by the existing

highly complex model. Furthermore, different options were introduced and tested like

upwinding some terms of the momentum equations (which seems to solve a number of

problems reported in the past), the second upwind scheme for the convective terms of

the phasic momentum equations and the implementation and assessment of a completely

different annular flow interfacial shear correlation.

The modified TRAC-BF1 is assessed against some bottom-flooding separate-effect

experiments, a "benchmark" top flooding simulation as well as against the TLTA test Nr.

6423. In the process of this task, the different options are assessed and discussed and

is shown that the predictions of the modified code are physically sound and close to the

measurements, while almost all the predicted variables are free of unphysical spurious

oscillations. The modifications introduced solve a number of problems associated with

the frozen version of the code and result in a version which can be confidently used for

LB-LOCA analyses.

iii

CONTENTS

Page

A bstract ........................................................................... iii

I Introduction .................................................................. 1

2 Summary of Modifications and Options in TRAC-BF1 ................................ 2

3 Assessment of the Modified Code and Comparison with the Frozen Version .............. 19

3.1 The NEPTUN Bottom Flooding Tests ...................................... 19

3.2 A Top Flooding Benchmark Test ......................................... 20

3.3 The TLTA Test 6423 ................................................... 21

4 Concluding Discussions ....................................................... 23

References .................................................................. 26

Figures

1 NEPTUN bottom flooding exp. Nr. 5036. Collapsed Liquid Level (m) ................... 28

2 NEPTUN bottom flooding exp. Nr. 5036. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B) ................................................... 28

3 NEPTUN bottom flooding exp. Nr. 5036. HTC to the liquid (W/m2/K) at

axial elevations 0.714 m (A) and 0.946 m (B) and void fractions near the

same elevations (C,D) ......................................................... 29

4 NEPTUN bottom flooding exp. Nr. 5050. Collapsed Liquid Level (m) .................... 30

5 NEPTUN bottom flooding exp. Nr. 5050. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B) ................................................... 30

6 NEPTUN bottom flooding exp. Nr. 5052. Collapsed Liquid Level (m) .................... 31

7 NEPTUN bottom flooding exp. Nr. 5052. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B) ................................................... 31

8 NNEPTUN bottom flooding exp. Nr. 5052. HTC to the liquid (W/m2/K) at

axial elevations 0.714 m (A) and 0.946 m (B) and void fractions near the

same elevations (C,D) ......................................................... 32

9 NEPTUN bottom flooding exp. Nr. 5036. Collapsed Liquid Level (m) .................... 33

10 NEPTUN bottom flooding exp. Nr. 5036. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B) ................................................... 33

11 NEPTUN bottom flooding exp. Nr. 5052. Collapsed Liquid Level (m) .................... 34

12 NEPTUN bottom flooding exp. Nr. 5052. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B) ................................................... 34

13 Benchmark top flooding test. Collapsed Liquid Level (m) ............................ 35

14 Benchmark top flooding test. RSTs (K) at 6 axial elevations .......................... 36

v
Figures (Cont)

Page

15 Benchmark top flooding test. Liquid velocities at 6 axial elevations .................... 37

16 TLTA Test Nr. 6423. Comparison of measured (- - -) and predicted (-) "maximum" RST histories with the frozen version of the code and the

fine-mist activated during reflooding ............................................ 38

17 TLTA Test Nr. 6423. Comparison of measured (--- -) and predicted

maximum RST histories. Andersen's annular flow shear model ....................... 39

18 TLTA Test Nr. 6423. Comparison of predicted RST histories.

Andersen's annular flow shear model ............................................ 39

19 TLTA Test Nr. 6423. Comparison of measured (- - -) and predicted

maximum RST histories with the Wallis annular flow interfacial shear correlation ......... 40

vi
1 INTRODUCTION

TRAC-BF1 is a best-estimate transient analysis thermohydraulics code for BWRs and its

original development at Idaho National Engineering Laboratories started at the early 80s in

collaboration with General Electric (GE). The original versions of the code were TRACBD1 and TRAC-BDI/MODl and were both using the semi-implicit solution scheme

for the hydraulic equations. TRAC-BF1 is an improved version of the aforementioned

previous versions and among a number of improvements which make the code both

physically more sound and numerically more robust, is the Courant - violating numerics

for the 1-Dimensional (1D) components. For a more detailed reading of the different

models in the code as well as the numerical solution techniques, the interested reader is

referred to Taylor et al., (1984) and Borowski and Wade (1992ab).

A systematic assessment of some models in TRAC-BF1 started at PSI in 1990 and

one of the areas to which attention was paid was the modelling of heat transfer during

reflooding. The code version used for this work (and also used for the work to be reported

here) was obtained from INEL in the summer of 1990 and is the UNICOS 5.1 version

named lU~test G2W2. We shall refer to this as the frozen version of the code. Since the

models which we have changed and on which we shall elaborate in this work are still the

same in the present frozen version of the code, we shall still be referring to the unmodified

version as the frozen version, although we shall be actually referring to an older version

than the one on which the term frozen would be appropriate.

In TRAC-BF 1, the code developers had already adopted the Bestion interfacial shear

correlation (Bestion, 1985) for bubbly/slug flow in rod bundles as recommended in the

past (Analytis, 1995) who showed that the corresponding model in the code is not suitable

for rod bundles. In trying to assess the code's capabilities to model reflooding, a number

of separate-effect bottom flooding tests at the heater rod bundle NEPTUN test facility

(GrUtter et al., 1980) at PSI were analysed and a number of interesting observations

were made which resulted in a set of model changes (Analytis, 1992). One of the

major model changes implemented in the code was a special wall-to-liquid heat transfer

package activated during reflooding which expresses the heat transfer coefficient (HTC)

by an empirical correlation exhibiting its experimentally observed strong dependence as a

function of the distance from the quench front(s) (QF(s)) as in the French code CATHARE

(Juhel, 1984; Bestion, 1991). Furthermore, an upwinding scheme of some terms of the

finite-differenced phasic momentum equations was implemented in the code (Analytis

and Coddington, 1992), the aim of which was to solve some problems arising from the

cell-centered way that these terms were finite-differenced in the frozen version of the

code and in particular, the way they were used in the wall shear modelling. This was

successfully accomplished. The modified version of the code was also used to analyse a

hypothetical LB-LOCA in a BWR-4 operating at 110% power and under the assumption

of limited ECCS (Analytis, 1995).

In this work, we shall summarise all the model changes and options we have implemented in TRAC-BF1. Some of these changes are a little different to the ones reported in

the past (Analytis, 1992), while some others we have found not necessary (eg the reduction

of the interfacial heat transfer from the vapour to the droplets). Additionally, we have

1
now implemented some new changes and options not implemented earlier in the code

(Analytis, 1992; Analytis and Coddington, 1992), which nevertheless we found helpful

and necessary. After discussing in some details all these modifications and options in

section 2, in section 3, we shall analyse a number of bottom flooding experiments in the

NEPTUN test facility (which is a heater rod bundle) both with the frozen and modified

version of the code and discuss the differences between them as well as the effect of

the annular flow interfacial shear correlation on the predictions. Additionally, we shall

analyse a top flooding benchmark test with different versions of the code and we shall

look in more detail on the effect of different interfacial and wall shear models on the code

predictions. Furthermore, we shall report on the analysis of the TLTA test Nr. 6423 with

different versions of the code. Finally, in section 4, we shall conclude this work by making

some recommendations.

2 SUMMARY OF MODIFICATIONS AND OPTIONS

IN TRAC-BF1

We shall now outline the code modifications and model changes we have implemented in

TRAC-BFL. Some of these changes are a little different to the ones reported in the past

(Analytis, 1992); additionally, a number of additional options have been added in the code

which we shall summarise.

(a) The main modification implemented in TRAC-BF1 was that during reflooding, a

special wall-to-liquid heat transfer package is used, similar to the one in an older

version of the French code CATHARE (Juhel, 1984; Bestion, 1991). In TRAC-BF1,

as in other existing thermohydraulic system codes like RELAP5IMOD3 (Carlson

et al., 1990), there is only one post-CHF wall-to-liquid heat transfer package independently of the physical process that is supposed to be modelled. This package is

based on modelling the wall-to-liquid film boiling HTC in subroutine HTCOR by

(Taylor et al., 1984; Borowski and Wade, 1992a)

hwL(FB) = (1 - ct)h(mod,BR) (2.1)

where hmod,BR is the modified Bromley correlation and a is the void fraction. This

wall-to-liquid HTC is used if the wall temperature T. is greater than a minimum film

boiling temperature TMIN which in the code is given either by the homogeneous

nucleation temperature, or by the Shumway correlation. If now the T,,, is lower than

TMIN but greater than the wall temperature TCHF at the critical heat flux, transition

boiling is assumed to occur, and the code uses the Bjornard quadratic interpolation

between the CHF point and the film-boiling wall-to-liquid HTC hwl(FB) ) in the

usual way

h.1(TB) = (1 - 1r) h.I(FB) + F Tc - (2.2a)

Twe- Tr

where

2
= (TF _ TM IN (2.2b)

and all the symbols have their usual meaning.

Clearly, as is well-known, the aforementioned approach does not take into account

the experimentally established fact that during reflooding, there is a strong dependence (in fact, enhancement) of the HTC in the vicinity of the QF, far above the

values one would get in non-reflooding situations. Hence, while we only introduced

the simple radiation model described below in HTCOR (notice that this subroutine is also called by the subroutine CHNQ which activates the fine-mesh for the

channel-wall), we added the subroutine HTCOR1 which is called only if the reflooding is activated and the component is the core. In this subroutine, in addition

to the aforementioned radiation model, we implemented the simple empirical wallto-liquid heat transfer correlation which is similar to the one in an older version of

the French code CATHARE and reads (Juhel, 1984)

hwl(FB) = max {(1400 - 1880 AZQF)

min (1 - 0.5), 0} + hBR V' - a (2.3)

where AZQF is the distance from the top or bottom quench-fronts (QFs) and hBR

is the original Bromley correlation with hBR - (AZQF) 0 25 . The same wallto-liquid HTC is assumed independently of whether we are in the vicinity of the

bottom or top QFs. Here, we should point out that since during reflooding, we hope

that the present approach will actually predict the "knee temperature" as being the

temperature of the rod at this elevation when the QF arrives there (as is physically

the case), if the aforementioned reflooding heat transfer logic is activated, we

have assumed the actual homogeneous nucleation temperature as the minimum film

boiling temperature rather than trying to obtain this temperature from a correlation

as is done in the code if one sets in the input ITMIN=l, by which the Shumway

correlation is chosen as the "knee temperature". We shall see in the next section that

if this option is chosen, due to the fact that the Shumway correlation predicts a "knee

temperature" higher than the homogeneous nucleation temperature, one is forcing

unphysically the different rod nodes to quench at a relatively high temperature

independently of whether the QF is in the vicinity of the node or not. As a result

of this, large oscillations may appear in the calculated quantities like velocities,

void fractions etc. In this work, in the analysis of separate-effect bottom flooding

experiments, we are not allowing any top quenching since we are preventing the

code from going to transition boiling (and subsequently, to nucleate boiling) if the

distance from the bottom QF is greater than 0.1 m. For the general case, this is

clearly not so since if there is no CCFL, during spray cooling, there will certainly

be a film created at the top of the rods (and hence, a top QF) which will start

propagating downwards.

We should point out here that during reflooding, we have implemented a similar wallto-liquid heat transfer model in RELAP5/MOD3 (Analytis, 1995, 1996). There, the

3
model is a little different in as far as in addition to defining a HTC by eq. (2.3),

the final HTC used in the film-boiling regime during reflooding is given by the

maximum of eq. (2.3) and the one of Forslund-Rohsenow for contact heat transfer

between the hot walls and the droplets. We think that a number of differences

between the constitutive equations of the two codes may be responsible for this

difference (eg the interfacial shear package is quite different in RELAP5/MOD3

when compared to the one of TRAC-BF1).

(b) There are different CHF options in the code, one of them being the Biasi correlation.

In relation to reflooding, we have already shown in relation to RELAP5/MOD3

(Analytis, 1995, 1996) that although in reality the actual value of the CHF is not

of great importance, the way it is used in the wall heat transfer package can excite

numerical oscillations if such oscillations exist on the CHF. A CHF correlation which

is a very simple mathematical expression and has been shown not to exhibit any

oscillating behaviour is the simple modified Zuber correlation and during reflooding,

this was the one used in RELAP5/MOD2 and MOD2.5 (Ransom et al., 1985).

Hence, jf the reflooding is on, we use the following CHF correlation:

qcHF _- q(CHFZuber) max 1, (sG)0) (2.4)

where G is the total mass-flux and G,,i = 67.8 kg/m2s. As we have already

mentioned before, while the actual value of the CHF is of the utmost importance in

cases that the prediction of the dry-out point is required, in the case of reflooding, it

is not so important.

(c) The drift-flux interfacial shear package in the code (subroutine FRCIF) was derived

(Andersen et al., 1983), based on the work by Ishii (1977) who computed the driftflux parameters Vgj and Co for different flow regimes. For annular flow, the Vij and

Co obtained by Ishii (1977) were based on the Wallis interfacial shear correlation

and the expressions derived by Andersen and used in TRAC-BF1 (Taylor et al.,

1984, Borowski and Wade, 1992a, 1992b) are

Co 1 + (1- a)(1 - e) (2.5a)

(a + a)

where

a 11+ 75(1 - ce) (2.5b)

VGP

e is the fraction of the entrained liquid and is calculated as in TRAC-BFl in the

following way:

e = max( (xe - 0.03),0.0)/ V1 + (0.1 + xe)2 (2.5c)

where

xe = 10-6 g( ; o D.s Rei). (2.5d)

Re, = IG1I DH / 1 (2.5e)

4
G1 = (1 - a), pt, (2.5f)

D" = DH 9.81 (p, - p9) / (2.5g)

J= = IJg/I V/a 9.81 (pP - p9) (pg / (p, - pg)).6 666 / pg, (2.5h)

and

J= Vg. (2.5i)

Then, we shall have

= 0.015 p, a (a + a) 2 (2.6a) DH

and the interfacial shear per unit volume will be

fi = 0.015 pi a (a + a) 2 (C1V9 - COVI) 2. (2.6b)

DH

where

1- -a O (2.6c)

Notice that the aforementioned annular flow interfacial shear correlation it was assumed that the continuous phase is the liquid and hence, the liquid density p, appears

as multiplier. This makes the interfacial shear up to three orders of magnitude higher

than the one of Wallis which is based on the assumption that the continuous phase

is the vapour and reads

0 p(1 + 75(1e)) - V) 2. (2.7a)

(In the Wallis correlation, Co = C, = 1). As a matter of fact, the drift-flux parameters Vgj and Co computed by Ishii for annular flow were based on this correlation;

should the assumption had been made that in annular flow, the continuous phase is

the vapour, equation (2.6b) would read

0.015 pg a (a + a)2 (CIV9 _ CoV1)2 " (2.7b) fi = DH

and this is "equivalent" to the Wallis correlation Hence, by changing the pi to

p, in the annular flow interfacial shear correlation in TRAC-BF1, one should get

results equivalent to the ones obtained by using explicitly the aforementioned Wallis

correlation. In the course of this work, we have realized that although for bottom

flooding, if the code assumes annular flow, the interfacial shear correlation as it

stands in the code is too high and always results in excessive liquid carry-over, if the

water is entering from the top, it gives quite reasonable results and can even "predict"

a kind of falling film. The opposite is true if the continuous phase is assumed to

be the vapour (as in the Wallis correlation); then, the interfacial shear is sufficiently

small to prevent excessive carry-over during bottom flooding but in contrast, if the

water enters from the top, one obtains rather unphysical results, the entering liquid

is pushed to the bottom and no film is created at the top. Recently *

* P. Coddington, Private communication (1996)

5
it was suggested that this pathological behaviour may be due to the particular

form of the wall shear model used in the code (which is also responsible for the

pathological behaviour reported in the past (Analytis and Coddington, 1992), which

led us to the upwinding of some terms of the momentum equations) and attempts

were made to investigate this possibility. Under item (j) below, we shall further

report on this problem and show that if one uses a wall shear model like the one of

TRAC-PF1 (Liles et al., 1986; Gao and Leslie, 1989) which is not based on the total

pressure-drop (hence requiring subsequent partitioning of the wall shear between

the two phases) but rather defines the wall friction individually for each phase (or if

one tries somehow not to "weight" the wall shear between the liquid and the vapour

with the liquid fraction), most of the aforementioned problems are solved.

Finally, in the code, the following values of CIO are used in the phasic momentum

equations:

If1 - a < 0.1,

CIO = 10 CIO (1 - a). (2.8a)

Ifl - e < 0.1,

CIO = 10 CIo. (2.8b)

Ifl - e > 0.1,

Cio - 1 - e (2.8c)

As an option, we introduced in the code the original Wallis correlation (which, in a

way, is equivalent to changing p, to p2 in the correlation in TRAC-BF1).

(d) As an additional option, we implemented in the code the Bharathan annular flow

interfacial shear correlation used in RELAP5/MOD2, the interfacial shear coefficient

of which reads:

CIO = San,,n Pg (0.0025 + 0.1375 10.0(9.07 R) F(1" + 4.74 R')) (2.9)

where the symbols are defined as follows:

- 6.1192775 ý/-b ci(1

SD-D (2.10a)

F = max(lO-8 , 0.5 (1 - ago) D*), (2.10b)

R* = min 30, 4 (2.10c)

Rý = (2.10d)

D* is defined by (2.5g) and

ab = max(O, (1 - a) Vgbb), (2.10e)

6
igs =V11 - C•b (2.10f)

where

Vgbb = e-S max(O, (1 - I0- 4 Vi(i - a) Pi Did)) (2.1Og)

and

2.5 I -g) -(p (2.10h)

Pg

where

S= 4 10-5 V') 6 (2.10i)

The distribution coefficients C, and C, in the interfacial shear per unit volume fi are

equal to 1 and ad = 10-. The employment of this correlation during the analysis

of different tests alleviated many problems encountered with the aforementioned

other two correlations. This we implemented in a new subroutine called FRCIFR

which is called if this option is required.

(e) The drift-flux coefficients Co and C1 are modified in the interfacial shear subroutine

FRCIF if V, < 0 and V9 > 0; if this is the case, CCFL is assumed in FRCIF, even if

the actual CCFL flag in the code (ICCFL) is switched-off. Though, one can easily

have a situation in which at one node, V, < 0 and V. > 0 (for example, during

bottom flooding), without having CCFL, but simply because the droplet is too heavy

to be lifted by the vapour and falls back. Hence, we removed the modifications of

Co and C1 which were responsible for some numerical oscillations, since clearly,

one can have at one time step V, < 0 and at the next, V, > 0. Though, these

modifications are still activated but now only if for that node ICCFL 3& 0. To

achieve this, the flag ICCFL is made an argument of the subroutines FRCIF and

FRCIFR.

(f) In order to let the liquid velocity V, approach the vapour velocity Vg as a- > 1, in

a previous work (Analytis, 1992), we have utilised an approach similar to the one

used in the old version of the code, TRAC-BD1. In this work, we use a different

approach and assume that the interfacial shear calculated by the code is valid up

to a = 0.9999999999; for a > 0.99999999999, a value given by a(1 - a)1020

was assumed while between the two limits, we used an interpolation. We did this

modification both for the ID and 3D vessel component. For doing this, some limits

defined in the solution subroutines TFIE and TF3E of the 1 D and 3D momentum

equations had to be modified. This approach is different to the one reported by

Analytis (1992).

(g) Similarly to other system thermohydraulic transient analysis codes like RELAP5, in

TRAC-BF1, in the semi-implicit numerical solution scheme, (subroutine TF1E),

for each edge of a computational cell and for each time-step, after linearising the

interfacial (and wall) shear terms (for more details, see (j) below), the 2 phasic discretized momentum equations are solved as a system of 2 linear algebraic equations

7
with the new phasic velocities at the edges (k + 1) of the computational volumes

as unknowns, e.g.

( M 11 M 12 k JjR 2 1 a

M21 M22 ,+1R

where a superscript indicates the time-step and

1______ cto (a, k+1⁄2)

(fw)k+, n (2.1 lb)

is a function of a number of old-time variables. The subscript "p" refers to the phase

p (g or 1) while f,, is the wall friction, and C1o (-, k+ i) the flow-regime dependent

interfacial shear coefficient which depends on the junction void fraction a, k+1

while •n a twiddle over indicates that is upwinded. An upwinded phasic quantity

Yp, k+ L is defined by

_n Y,.k+1⁄2 = (WP) (YP)n + (WPI) (Yp)•+ (2.12a)

where 1 VZ ,k+1⁄21 > 0

(WP) = (2.12b)

0~ 0n

(WP1) = 1 - (WP). (2.12c)

(Yp)[ and (YP)I are the cell-centred phasic quantities up-stream and down-stream

the junction (k + 1), respectively. Finally, in eq.(2.1 la), N are

Dn'+1 Dn

f .tk+l -- PLk

01P =,PP}k+1⁄2 -(PP)k+2

n°+2 n V V")k+1⁄2, ... (2.13)

-(Cty PP)k L..

where P is the pressure and (Vp V , the convective phasic velocity terms. k~•y

In the aforementioned formulation, quantities like 1/(a, Pp)k+" and

1(pp,)'k+1⁄2 are cell-length averaged (indicated by a "bar" over them) as follows:

Axk (ap p,)k + Axk+1 (a, p)•+, (2.14)

AXk + AXk+1

8
where AXk, AXk+l are the cell lengths of the adjacent cells.

As has already been shown elsewhere (Analytis and Coddington, 1992), the numerical modelling of the phasic wall friction (f,,)' in eq.(2.1 lb) and in particular, the

incorporation of the discretized form of this term into the general finite-difference

formulation of the phasic momentum equations is of the utmost importance. For

our purpose, it is sufficient to note that this term is a function of the square of the

total mass flux Gn+I (evaluated at the junction k + 1), i.e. } 2 ~

= (akgkt ,k+f I, +'

k+1 21

where all the symbols have their usual meaning, and the terms (ap Pp)k+I are

cell-length averages defined by eq.(2.14).

If we now inspect (fn)k+1⁄2 as given by eq.(2.15) and consider the case in which

pure steam (ag,k = 1) from cell k (left) is injected into the cell (k + 1) (to the

right) containing pure water (ag,k+l = 0). The void fraction ag, k+1⁄2 at the junction

(k + 1) between the two volumes is the upwinded void fraction and for our case

with V., k+ > 0, will be equal to 1. Additionally, at this junction one will usually

have I,+ k+1 > 0 and hence, the wall friction given by eq.(2.15) will attain a high

value due to the cell-length averaging of the term (at pt)' + (cf. eq.(2.14)) which

will dominate. This leads to an unphysical and high pressure drop between the two

adjacent cells. Clearly, for this very simple case, this problem will be most severe

when a., k = 1 and ag, k+1 = 0 but will always show up when a., k >> ag, k+1

while for ag, k - ag, k+l, it will disappear. For a large multi-component system

calculation, one may come across a similar situation on a number of occasions in

the course of a transient, but it would be very difficult to quantify.

The most obvious way of eliminating this problem is by modifying the cell-length

averaging procedure of the terms (ap pp)n+1⁄2 in eq.(2.1 la) and employing a donorcell differencing (upwinding) approach instead. Hence, we now define

(a,, PP)n L = (WP)(a,, pp)'• + (WP1)(a,, Pp)'+i (2.16)

where (WP) and (WP1) are given by eq.(2.12b) and (2.12c). By upwinding these

terms (and hence, the total mass flux G'+ L), one can see that for our simple steam

injection problem, Gk+ L will be given by (for Vn +,1⁄2 > 0)

2- 7

Gk+I.= (a p) k + 1 . (2.17)

There is now no liquid contribution to this term, since (a, p,)n = 0 and hence, the wall

friction will not attain a very high value. Finally, we should mention that although

the origin of the problem was traced to the cell-length averaging procedure used

to define the total mass flux in the wall friction terms, we found it "numerically"

9
advantageous to extend these modifications and upwind all (ap pp)' +, terms in

eqs.(2.1 lb) and (2.13) as shown in eq.(2.12a). More discussion of this problem is

given elsewhere (Analytis and Coddington, 1992). In the present work, we have also

included this upwinding option in the discretization of the 3D phasic momentum

equations for the vessel in subroutine TF3E. We shall further discuss this point in

the following section, where we shall present some TLTA calculations; in particular,

we shall show that the two different discretization schemes result in different rod

surface temperature (RST) histories. Finally, we should mention here that since as

we shall discuss later in this section, we found that in general the wall shear model

of the code was highly problematic and we decided to change it to the one of TRACPFl, it is possible that the reasons for which we had to introduce the upwinding

option (Analytis and Coddington, 1992) for some terms of the momentum equations

are no longer there. In this work, we did not try to answer this question.

(h) As an option, we introduced the second upwind scheme for the convective terms of the

1D phasic momentum equations (subroutine TF1E). When this option is activated,

it is only used for the PIPE and CHAN components. Employment of this scheme

for all components resulted in some problems during the analysis of the TLTA test

Nr. 6423 whose origin we had no time to investigate. In any case, as we shall

discuss in the following section, the employment of this scheme resulted in a little

different RST histories for the aforementioned TLTA test when compared to the

ones predicted when the donor-cell differencing scheme of TRAC-BFl is used. The

second upwind scheme which is also used in RELAP5, is numerically less diffusive

than the one in TRAC-BF1. For any phasic velocity V, we define the quantities

= 0.51Vk_1⁄2. + Vk+l (2.18a)

and

d,= 0.5jVk+, + Vk+11 (2.18b)

and the cell-centered velocities

Vk = 0.5(Vk_.. + Vk+L) (2.18c)

and Vk+1 = 0.5(Vk+I + Vk+!). (2.18d)

Subsequently, we define

VkVk = VkVk + 0.5C(Vk-1 - Vk+L) (2.19a)

and

Vk+lVk+l = Vk+lVk+l + 0.5f+(Vk+1⁄2 - Vk+2). (2.19b)

An alternative, higher-order accurate scheme which we also implemented (and

should be even less diffusive) is one in which c and c+ are defined by:

1

f- = 0.51Vk_ + Vk+1⁄2 - 1(Vk+i + Vk_.1⁄2 - 2Vk+L)I (2.20a)

10
and1 and =0.5lVk+L + Vk+!- I(Vk+8 + V_.L - 2Vk+1⁄2)I (2.20b)

and subsequently, for the cell-centered velocities

1

Vk = 0.5(Vk_1⁄2 + Vk+i 4 1 (V +i + Vk_1⁄2 - 2Vk+1⁄2)) (2.20c)

and

Vk+= 0.5(Vk+1⁄2 + Vk+1 -(Vk+8 + VkL - 2Vk+L)), (2.20d)

while finally,

VkVk = VkVk +

0.56-(Vk_. - Vk+. - (Vk+i + Vk_1⁄21 - 2Vk+i.)) (2.21a)

and

Vk+lV÷+l = Vk+aVk+l +

0.5e+(Vk+L Vk+! - (Vk+i + Vk_.. - 2Vk+1⁄2)). (2.21b)

Finally, in either case, the convective term can be written as

1 a(V 2 ) k+I =_ON (O.5(Vk+l Vk+1 - Vkvk))(22a k= CONV - (2.22a)

where

Axk - Axk + AXk+l (2.22b) 2

and Ax,, and Ax,,+l are the lengths of the meshes upstream and downstream the

junction k+ 1. There have been some cases for which a smoother and "less diffusive"

behaviour of the code predictions was indeed observed when this scheme was used

for the convective terms of the phasic momentum equations. We shall not discuss

this issue any further, but we think it is worth investigating.

Before concluding this section, we should say that there are different ways of

finite-differencing the spatial part of the virtual mass term in the phasic momentum

equations; these we also introduced as an option. The spatial part of the virtual

mass term is defined by:

= k. pca).VV - V) (2.23a) ( apg (2.23a)

for the vapour momentum equation and

fvmI = kv. ( 1 Pc.(a ) .VDV(Vt - V,) (2.23b) k(1 -e otpI

11
for the liquid, where VD = (1 - a)V9 + aV1, k,,,, is the virtual mass coefficient

and the subscripts C and D are referring to the continuous and dispersed phase,

respectively (Taylor et al., 1984, Borowski and Wade, 1992b). In the code, the a

in Pc at the junction j + 1 is defined by 0.5(a 2 _i. + a+i). Instead, we used the

upwinded j+I . Furthermore, the finite-differenced term VD.V(Vg - VI) is defined

in the code as

VD.V(Vg - Vl + =

VDoJ+1⁄2 ((1 + Ski,)(V,,j+ - Vij+1⁄2. - Vgj_1⁄2 + Vtd_1⁄2)

+l- Si)gj!- Vji-g'+- + V1j0 Aj1(2.23c)

where

Ski. = sign(1, VDj+L) (2.23d)

Instead, we finite-differenced this term as follows:

VD.V(V-= (v- +4 [(wv(V, V 2

(1 - (WV))(Vg,,+1 - Vgj+1⁄2) - (WL)(V ,.+1⁄2 - V1,j-_) +

(1 - (WL))(V~j+1 - Vij+)])/Axi+1⁄2 (2.23e)

where WL and WV are 1 or 0, depending one whether V, +. and V- i• are greater

or less than 0, respectively. Notice that with the definition of VD, during finitedifferencing, some of the terms in VD.V(Vg - VI) can be expressed by the second

upwind formulation discussed in this section.

(i) A simple radiation model similar to the one in RELAP5 was implemented in TRACBF1 in the subroutine HTCOR holding the heat transfer correlations. The radiation

model in the code was found to induce some oscillations and was de-activated from

the input in all our runs. Comparison of this very simple model with the sophisticated

model in the code revealed very similar temperature histories during the analysis

of bottom flooding low flooding rate experiments without the oscillations observed

when the original model of the code was used.

(j) As we have already mentioned before, although the employment of the Wallis annular

flow interfacial shear correlation ( or of the one used in the code but with p,

in the place of p.) is necessary if one wants to reproduce results from bottom

flooding separate-effect test experiments, this seems to create some problems for

some transients. In particular, for the TLTA test Nr. 6423 it leads to a very small

time-step approximately 65 seconds after the initiation of the transient. In contrast,

the correlation in the code (derived by assuming that in annular flow, the continuous

phase is the liquid), always over-predicts the liquid carry-over in bottom flooding

tests, but it performs quite reasonably for system calculations. During a number of

discussions, it was suggested * that the origin of this problem may

lie on the modelling of the wall shear which was also responsible for the pathological * P. Coddington, Private communication (1996)

12
behaviour reported elsewhere (Analytis and Coddington, 1992). Hence, in order to

investigate this possibility, we implemented in the code the wall shear models of

TRAC-PF1 (Liles et al., 1986; Gao and Leslie, 1989), a "heuristic" model which

does not "weight" the shear between the wall and the liquid with the liquid fraction,

as well as an other simple model similar to the one of TRAC-P in which the

wall shear is also separately defined for the two phases. These modifications we

implemented in subroutine FRCW; though, the ID phasic momentum equations

solution subroutine TF1E was also appropriately modified. Before explaining the

models, we shall first outline the discretized phasic momentum equations for the

ID components of TRAC-BF1 and the way we had to modify them in order to

implement the TRAC-PFI-type wall shear models. The corresponding equations

for the vapour and the liquid are written as:

(V;n+'- vn)k+1⁄2L

At

+ (Pi+¶ - PO)

Azk+1⁄2. (Pg)k+ 4

+ (Vnb6V)k+l +

conv, term

Co(Cv "g,k+ "+l ( -2)+ - C" 0 v,k+1⁄2, ,,+1 Y)q

F-+-/-p)•

Cwg(vn+l\)2

-+ g cosB =0

(aPg )k+1⁄2

(2.24a)

and

V")k. IAt

+ (PT+1 - Pk)

Axk+4 (Tk+1⁄2

I1 c'o (CcV1

£7• + (v1:#SvI)k+! + - -

conv. term

(f7in+

(Vn+l )2 cW,, I k+1⁄2)

+ - +2osO = 0

(( - a)pl)k+L

f.+ I,, j1(( 10 )oPI+1⁄2 - -n

4-+I. - cov-') - Ct)PI)kL1⁄2

I((-,a)P,) '+i1⁄2

(2.24b)

where overlined quantities indicate that they are cell-length averaged. For the sake

of notational convenience, we have ignored the virtual mass terms. As is wellknown, in the above equations, the terms with the flow-regime dependent factor

IYo are due to the interfacial shear (with q = 2 or 4 in the different interfacial shear

correlations (Taylor et al., 1984), while the ones with the C_,g and c.at factors are due

to the wall shear. Notice that in the momentum equations (2.24a) and (2.24b), the

exponent q is defined by

q = (1 - X2)qbubb + X2 ((1 - (ent))qann + (ent)qdrp) (2.24c)

where qbubb = 4 ( = 2 for bundles), q•nn = 2 and qdrp = 4. X2 = 0 at the void

fraction corresponding to the end of the bubbly/slug flow regime and is equal to

13
1 at the annular-dispersed flow, while (ent) is the entrainment fraction. Hence,

in general, the exponent q is not an integer quantity and is defined as above for

smoothing purposes.

For each time-step and each junction k + 1, the above system of algebraic equations

can be solved as outlined under item (g) above only after linearising the interfacial

and wall shear terms so that the new-time phasic velocities appear linearly. We do

this by approximating (Taylor et al., 1984)

C -(q - 1) (V;)' + q(V;)'-1 (V;+')). (2.25a)

for the interfacial shear term (where VR = (C 1Vg - Co VI)), and

fl _ c(-- (Vn) 2 --+ 2(Vn) (V/n+1)), (2.25b)

ffl •_• (fl.) 2 + 2(V g)(14+)). (2.25c)

for the wall shear. Though, the aforementioned linearization of the wall shear terms

is not used in TRAC-BF1. The general modelling of the wall shear in the code is

done as follows: We first define

dpn

- p,-n(Vn,)2 dz •f (2.26a)

and

Idpn

Cfg- dz ((2.26b)

where the subscript wf means "wall friction". At the end-of-time step, the total

pressure gradient due to the wall friction is divided between the phases using the

void fraction as a measure of the surface area in the cell "melted" by each phase;

we obtain

(&fl) 2 p'•cfgVgflVg"+ + a"(l - anpfctla/&TJl"+I (2.26c)

and dp ~ I~dpln+l

- (1-a)-@ dz dz

(a•)(1 - vvg•+ g")pc 1 + I (1 - ") 2

pc 1 V7vf+1 (2.26d)

14
In principle, one could easily linearise the above wall shear model along the lines

of eqs (2.25b) and (2.25c) by approximating the implicit phasic velocity products

Vn+1 V"+ 1 by (2VnVn+1 - (Vn) 2 ) rather than by V'V"+' as done in eqs. (2.26c)

and (2.26d). As we said before, these terms are not linearised in the code, although

we have included such a linearization as an option if the original wall shear model

of the code is used. Due to the special form of this model, both the liquid and vapour

wall shear contributions are appearing in each of the phasic momentum equations

(Borowski and Wade, 1992b); in particular, the former is the wall shear entering the

vapour and the latter, the liquid momentum equations, respectively. In assessing

the different wall shear models, we also corrected a logic error in the subroutine

TFHE solving the phasic momentum equations, in which the calculated wall shear

coefficients were not re-set to 0 if there were back-ups. *

In order to investigate the possible relation of the aforementioned problems with

the wall shear model in the code, we implemented in TRAC-BF1 the wall shear

models of TRAC-PF1 (Liles et al., 1986; Gao and Leslie, 1989). In these models,

the wall shear of the liquid appears in the liquid momentum equation and the one

of the vapour in the vapour momentum equation. The TRAC-PFI models are the

homogeneous and annular flow models and can be summarised as follows:

I. Homogeneous Model

The pressure drop due to the wall friction is generally calculated from

dzPw = c 9V2 + cV 2 (2.27)

where c,• and c,,t are defined by

= &PgCfg (2.28a)

DH

and

cwt = (1 - DHa)plcf (2.28b)

Then, the frictional pressure drop can be written as

dp _ ap-cfg (1 - t)pic1 IV1 2

dp D Fit Wg 2 + + (2.29)

dz DH DH

In the homogeneous model, the friction coefficients cfg and cf I are calculated as

follows:

cf L,hom = 2f (2.30)

Cfg,hom = 2f, if a < 0.9

= (10a - 9)2(21 - 20a)Cfl,hom, if a > 0.9 (2.31)

P. Coddington, Private communication (1996)

15
Finally, f is defined by

f = 0.032, if Re < 500

= 0.032 - 5.25 10-6(Re - 500), if 500 < Re < 5000

= 0.046(Re)- 2 , if Re > 5000. (2.32a)

where

Re = GDH G = PmVm. (2.32b)

The mixture viscosity Pm is calculated by

S+ (2.32c)

Pm Pg PM

and the flow quality x by

1

= = (1 - a)pIVI (2.32d)

cpgVg

Finally, the mixture density and velocity are defined by

Pm = apg + (1-a)pI (2.32e)

and

V,, apV 9 + (1 - a)pI VI (2.32f)

H. Annular Flow Model

By this, it is not actually implied that flow is annular, but it is just a "label" given to

this particular model. In this model, the wall friction coefficients are given by

cf, = 2f, if a•<0.9

= (1 - wf)2f + WjCfl,h., if 0.9 < a < 0.9995

= cllthomo, if a > 0.9995 (2.33a)

and

cfg = 0, if a<0.9

= (lO -_ 9)2(21 - 20a)cf1 , if 0.9 < a < 0.9995

= Cfg,homo, if a > 0.9995 (2.33b)

where

= (a - 0.9) (2.33c) 0.0995

16
Then, the friction factor f is defined as

f = foP (2.34a)

where b

f-P = a + (Re)c (2.34b)

with

a = 0.026 k + 0.133 k (2.34c)

b = 22 ( k). 4 4 (2.34d)

c =1.62 (2.34e)

where k = 5 10', while the Reynolds number is the same like the homogeneous

model. Finally, qS2 is given by

0= PL (2.34f)

PM

In implementing the aforementioned TRAC-PF1 wall shear models in the code, we

also linearised the wall shear term (cf. eq. (2.25b), (2.25c))

M. No Weighting of Wall-Friction to the Liquid

Although even in the TRAC-PFI model the wall shear is "weighted" with the

corresponding phasic fraction, in reality, this is not the case and eg in annular flow,

the wall shear force acts only on the liquid film and not on the vapour or the droplets.

For this, in addition to the TRAC-PF1 models, we implemented in the code a logic

according to which eq. (2.31) reads

cf, = 2f, if a<0.96

= (10a - 9.6)2(21 - 20a)cji, if a > 0.96 (2.35)

while eq. (2.29) is modified as

CgV 2 (WAWET) pcfIV1 2 dp _+ (WAWET)+ (2.36)

dz DH DH

where WA WET is internally defined in the code and represents the fraction of the

particular volume which is in contact with the liquid: WAWET = 1 indicates that

the wall of the volume is wet, while WAWET = 0 indicates that it is dry. The

aforementioned logic does actually somehow captures the fact that if the wall is wet,

the wall shear is only between the wall and the liquid and is not actually "weighted"

with the liquid fraction. In implementing all the three aforementioned wall shear

models in the code, we also linearised the wall shear term (cf. eq. (2.25b), (2.25c)).

17
Concluding, we should say that in addition to the models above, we have recently

tried and implemented in the code two more wall friction models with equally successful

results. In the first model, cf and cfg in eqs. (2.28a) -(2.28b) are given by

cit = Af(e)

cfg = Agf (Reg),

(2.37a)

(2.37b)

where Rep denotes the Reynold's number of the phase p; eg for the vapour,

Re. = 'pgVgDH

11g

(2.37c)

The ramp function A is defined by

\9 = 0, ifa<0.999

a - 0.999

- 0.0009 ' if 0.999 < a < 0.9999

= 1, ifa > 0.999 (2.37d)

At = 1 - Ag. (2.37e)

The friction factor f(Rep) is for smooth pipes and is defined by

f(Re = 1) 1 max(64 0.0055+ 0.55 S Rep ( 5 (Rep)O-333) (2.37f)

while for turbulent flow (Rep > 5000), we use

f (R p) = 1 (0.046,

f(Re~2 - ReO.) (2.37g)

In the second model which is the Churchill's model, * one first defines

ff = ((Re) .0

+ (aa + bb)1.5' (2.38a)

where aa, bb and cc are defined by

* D.R. Page, Penn State University, Private communication

(1997)

18
'166

aa -- (2.457 log o(7)'9 .276 (2.38b)

Ree +.2a

bb = (37530 , (2.38c) Re)(23c

cc = 1 (2.38d)

where Ed is the surface roughness divided by the hydraulic radius (half the hydraulic

diameter) and the Reynold's number Re is defined by eq. (2.32b). Then, cfI and cfg in

eqs. (2.28a) -(2.28b) are given by

ef1 = 4 ff, (2.39a)

cfg = C I, (2.39b)

A number of runs of different problems was also made with the aforementioned wall

shear models; in all cases, the results were rather similar to the ones obtained when the

TRAC-PF1 homogeneous wall shear model was activated; hence, in this work, we shall

not report on these results.

3 ASSESSMENT OF THE MODIFIED CODE AND COMPARISON WITH THE FROZEN VERSION

3.1 The NEPTUN bottom flooding tests

In this section, we shall analyse 3 bottom flooding tests in the NEPTUN test facility at

PSI. Some of the important boundary conditions of the tests are summarised below:

Nr. 5036 (P = 4.1 bar, AT, = 10 K, VIN = 0.015 m/s),

Nr. 5050 (P = 4.1 bar, AT, = 78 K, VIN = 0.15 mis),

Nr. 5052 (P = 4.1 bar, AT. = 78 K, VIN = 0.025 m/s),

In the analysis of the bottom flooding experiments presented here with the modified

version of TRAC-BF1, we have used for annular flow the interfacial shear correlation in

the code but with p. in place of pi (i.e. equivalent to the original Wallis correlation; see

19
Section 2 under (g)). Additionally, we have made all the runs by using the upwinding

option described in section 2 under (e). In a previous work (Analytis and Coddington,

1992), we have shown that for these experiments, the differences between the predictions

obtained by using the upwinding and the standard code cell-length average approach

are almost non-observable. Furthermore, we shall present some results obtained with

the modified version of the code but with the Bharathan annular flow interfacial shear

correlation of RELAP5/MOD2 and compare them with the results obtained by using the

Wallis correlation.

In Fig. 1 - 8, we compare for these 3 tests the collapsed liquid level (CLL), RST

histories and HTCs to the liquid and void fractions obtained by the frozen and modified

versions of the code with the corresponding measurements. Clearly, for all 3 experiments,

the predictions of our modified version are very close to the measurements. In particular,

for test Nr. 5036, the large differences between the measured and predicted RSTs are not

only due to modelling of the wall-to-liquid HTC during reflooding, but also due to the

excessively high liquid carry-over predicted by the frozen version of the code, as a result

of the high annular flow interfacial shear, a problem we already discussed in the previous

section (the runs with the modified code were made with the annular flow interfacial shear

in which the continuous phase is assumed to be the vapour eg, we have used p. instead

of pj). All the NEPTUN runs were made with the original wall shear model of the code;

analysis of one case with the TRAC-PF1 wall shear models discussed in this work (and

the upwinding option switched-off) did not show any significant differences and the very

good agreement between measurements and predictions prevailed. Hence, we do not

show these runs here.

In Fig.9 - 12, we show the same comparisons for the NEPTUN tests Nr. 5036 and

5052, but with the Bharathan annular flow interfacial shear correlation activated. If we

compare these predictions with the corresponding ones shown before, one can clearly see

that for the tests in which the code does assume annular flow for some time during the

transient, the CLLs are now lower than before. This indicates that the Bharathan annular

flow interfacial shear is higher than the one of Wallis, but still lower than the one used in

the frozen version of the code. Consequently, the different elevations quench later than

before, since now there is less liquid in the bundle.

3.2 A top flooding benchmark test

In this section, by using the input deck of the NEPTUN exp. Nr. 5036, we shall simulate

a top flooding benchmark test with the same initial and boundary conditions as the ones

of the real, bottom flooding test. The aim of this is to assess the way the code performs

under "clean" top flooding conditions. In doing this, we shall point out all the problems

arising from the modelling of the annular flow interfacial shear as well as the effect of the

wall shear model on the code predictions. The CCFL flag was activated at the top node

of the test-section.

In our analysis, we shall compare different predicted variables of interest by analysing

the case with 3 different code versions: The modified one with the Wallis annular flow

20
interfacial shear correlation and the TRAC-PF1 "annular flow" wall shear option, the

modified one with Wallis correlation and the wall shear model of the frozen version, and

finally the frozen version of the code. In Fig. 13, we compare the predicted CLLs, in Fig. 14

the RSTs at 6 different axial elevations and in Fig. 15, the predicted liquid velocities, also

at 6 different axial elevations. By inspecting these figures, one can readily draw the

following important conclusions:

(a) Similarly to the bottom flooding tests, the frozen version of the code excites oscillations

which result in spontaneous CLL depletions. The result of this is that the RSTs are

reaching very high values due to liquid deficiency in the test section and the fact

that before the liquid reaches the bottom of the rod bundle and starts accumulating,

is being sporadically expelled. Similarly to the corresponding real bottom flooding

test, the elevations in the middle of the rods do not quench. There are various

reasons for this behaviour, one of them being the wrong modelling of the HTC

during reflooding as discussed in this work.

(b) With the frozen version, the liquid velocities are exhibiting large-amplitude spikes

resulting in excessive liquid carry-over from the test-section. In the upper parts

of the core, there exists a small negative liquid velocity (indicating a kind of

falling film) which is suddenly interrupted spontaneously (and reverses sign) by

the aforementioned large liquid velocity spikes. The creation of the "film" at the

higher elevations is due to the very high annular flow interfacial shear used in the

frozen version of the code as we discussed in this work, which results in creating

and holding-up the film at the top of the rods.

(c) When the modified code is used with the Wallis annular flow interfacial shear correlation but with the standard wall friction package, the liquid velocities are now

everywhere large and negative and there is hardly any film formed at the top of the

rods due to the low annular flow interfacial shear which pushes the liquid to the

bottom of the test-section. Clearly, the wall shear is too low to hold the liquid film.

We have already discussed this problem in this work.

(d) When the modified code with the Wallis annular flow interfacial shear correlation but

also with the wall shear correlation of TRAC-PFl is used, the predictions look much

more reasonable and physically sound. In particular, there is a slowly downwardsprogressing film from the top of the rods and all predicted quantities are free of

unphysical oscillations. In particular, the liquid velocities are well-behaved and

free of spurious and unphysical oscillations.

3.3 The TLTA Test 6423

The description of the TLTA facility and tests are described elsewhere (Taylor et al., 1984),

while the input deck used for the analysis reported in this work was the one which was

used by the code developers and was sent to us together with the code. In this work, we

shall report on the analysis of the TLTA test Nr. 6423. The same test was analysed in

the past (Analytis, 1992; Analytis and Coddington, 1992); though, since our final code

21
version is now a little different to the one we used in the previous work, we shall repeat

here the calculations by using different options, and compare the predictions with the ones

of the frozen version of the code as well as with the measurements. The results of our

comparisons in the past showed that for this particular test, the predicted RST histories

are sensitive to the discretization scheme used in the solution of the 1D and 3D phasic

momentum equations (i.e., to whether some of the terms are upwinded or not).

In Fig. 16, we show the maximum RST history predicted by the frozen version of the

code with the fine-mesh switched-on, and we compare it to the measured one. Clearly,

as also shown and discussed in the past (Analytis, 1992), the predicted RST history is not

agreeing with the measurements: The peak RST is approximately 320 K higher than the

measured one and quenching occurs much later (approximately up to over 350 s later than

the measurements show). Furthermore, the axial elevation exhibiting the peak RST is

different (lower) than the corresponding ones when the test is analysed with our modified

versions. It is worth reminding the reader that as has already been pointed out (Analytis,

1992), in contrast to these predictions, if the upwinding scheme described under (g) above

is used, even with the frozen version of the code the differences between measurements

and predictions are no longer so large. We shall not discuss this point any further since

as we have already extensively elaborated in this work, independently of whether the

predictions of the frozen version for this test are much better if this scheme is used, there

are a number of problems with some of the physical models in the code.

In Fig.17, we compare the measured and predicted maximum RST history for the

TLTA test Nr. 6423 by using the modified version of the code, both with the standard and

upwinding options in the ID phasic momentum equations. The CCFL modifications of

the distribution coefficients C0 and C1 in the interfacial shear subroutines are activated (as

explained in the previous section) only if the CCFL model is activated for that junction.

One can see that although the modified code over-predicts the maximum RST by 50 K

when the upwinding option is used, the predicted RST(s) are close to the measured one.

Notice that the modified code does not predict the first dry-out.

The upwinding used in the case shown in Fig.17 was only for the ID components.

In Fig.18, we show a comparison of the predicted maximum RST history for the same

TLTA test by using the modified version in which the upwinding option is activated both

for the ID and 3D components, with the one in which the upwinding option is used only

in the ID components (as in Fig.17). Clearly, although the peak RST is for both cases

almost the same, when the upwinding option is used in all components, quenching occurs

later, following a plateau-type region of the RST. In other axial elevations (not shown

here), the peak RST of the nodes as well as the dry-out times are also a little different.

Finally, in Fig.19, we compare the predictions of the modified code with the Wallis

annular flow interfacial shear correlation (i.e. by setting p. in the place of p, in the

Andersen annular flow interfacial shear correlation in subroutine FRCIF), but with three

different wall shear models: The homogeneous and annular flow model ones of TRACPF1 and a rather "heuristic" one in which the shear between the wall and the liquid is not

"weighted" with the liquid fraction, but rather with the fraction of the volume which is

wet. The upwinding option was not activated in these runs. One can clearly see that all

the runs are completed without any problems. In contrast, with the original wall shear

22
model and the Wallis annular flow interfacial shear, it is very difficult to run this case due

to the fact that as we said before, after 65 s, the code was selecting very small time-steps.

Concluding this section we should mention that we have also analysed this test by

activating the second upwind scheme in the phasic momentum equations (for the CHAN

and PIPE components; see previous section under (h)). The predictions by analysing the

test with this option were very close to the ones obtained by using the standard scheme

in the code and in actual fact, exhibited some characteristics indicating less numerical

diffusion as compared to the results obtained by using the standard model in the code.

In any case, the main message conveyed from the analysis of this test by the frozen and

modified versions of the code is that as it is patently clear from Fig. 16, the frozen version,

for a number of reasons, cannot even remotely predict the peak RST history.

4 CONCLUDING DISCUSSIONS

In this work, we reported on a series of model changes and options implemented in TRACBF1. Some of these changes are a little different to the ones reported in a previous work

(Analytis, 1992) while some of the implemented options are new and serve to assess not

only physical models, but are also related to numerics. The main conclusions of our work

can be summarised as follows:

(a) The implemented empirical wall-to-liquid HTC activated in the code during reflooding

seems to perform very well. In particular, analysis of a number of NEPTUN bottom

flooding experiments with the new version of the code showed very good agreement

between measurements and predictions. Though, to achieve this good agreement,

one had to use the Wallis annular flow interfacial shear correlation which assumes

that the continuous phase is the vapour, instead of the one already in the code which

although based on the same correlation, is derived by assuming that the continuous

phase is the liquid, which results in a factor p, as a multiplier. This makes the shear

excessively high and in the case of bottom flooding, can greatly deplete the liquid

inventory.

(b) Two different annular flow interfacial shear correlations were implemented and tested

in the code and in addition to &he Wallis correlation (being the same one like the one

already in the code but with p. instead of pl) which seems to give good predictions

for bottom flooding experiments (in which this flow regime is activated in the

code), we also implemented as an option the Bharathan annular flow interfacial

shear correlation as used in RELAP5/MOD2. Clearly, due to the particular simple

flow-regime map of TRAC-BF1 (bubbly/slug flow on the one end, annular/mist on

the other and interpolation between), if during a particular transient the annular/mist

flow interfacial shear is activated (which is the case, if the wall is wet and the liquid

fraction is rather small), the code predictions are bound to be very sensitive to the

actual functional form of that particular interfacial shear correlation. We have shown

that within the framework of the TRAC-BF1 flow-regime map, for bottom flooding,

while the Wallis annular flow interfacial shear correlation tends to be a little too low

23
resulting in a small over-prediction of the liquid inventory, the one of Bharathan

exhibits the opposite effect, both of them being in general quite acceptable.

(c) The new approach for letting the liquid velocity approach the one of the vapour as

the liquid disappears seems to work well and is superior to the one we used in a

previous work (Analytis, 1992).

(d) We have shown that the RST histories obtained for the TLTA test Nr. 6423 are very

close to the measured ones whether the standard cell-length averaging or upwinding

schemes are used for some terms of the phasic momentum equations. Also a number

of calculations for this test were performed with the modified code including the

new wall shear package of TRAC-PF1, and the same good agreement between

measurements and predictions resulted. The second upwind scheme used in the

finite-differencing of the convective terms of the ID phasic momentum equations

was also tested (only if the component was a CHAN or a PIPE) and the code

predictions were also very close to the ones obtained when the first upwind scheme

is used. The second upwind scheme should be numerically less diffusive than the one

used in the code. Furthermore, we have shown that for this test, the RST histories

predicted by the unmodified code are not even remotely near the measurements,

indicating the inevitable general inability of the frozen version of the code to model

reflooding, but also a number of other problems from which the frozen version is

plagued.

(e) A number of outstanding problems related to the modelling of the annular flow interfacial shear were investigated by changing the wall shear model of the code and

implementing the one of TRAC-PF I. By doing this, the Wallis annular flow interfacial shear correlation could be used independently of whether we were modelling

separate-effect bottom flooding tests, top flooding tests or "system" calculations

like the TLTA. In particular, a number of problems which under certain conditions

manifested themselves when the Wallis annular flow interfacial shear correlation

was activated in the code were actually due to the wall shear model; most of these

problems were eliminated when the wall shear model was changed. A top-flooding

benchmark test was also analysed with different versions of the code (including the

frozen version) and the problems with the wall shear model of the code were clearly

identified. In view of all these we would strongly recommend changing the wall

shear model in the code which seems to create a number of problems in favour of

an alternative one.

Concluding, we should say that although at least as far as physical models are concerned, this work is restricted to a rather limited area in the code (wall heat transfer during

reflooding and interfacial and wall shear), there are certainly a number of other problems

which require attention as well as models whose improvement would be highly desirable.

Our general recommendations in relation to the problems analysed and discussed in this

work can be summarised as follows:

(a) Implementation of a wall-to-liquid HTC during reflooding which takes into account the distance from the upper or lower QFs. Such a correlation was

24
successfully used in this work and its physical and realistic behaviour clearly

demonstrated. This model is already implemented in the code by the Penn

State University team and delivered to the US Nuclear Regulatory Commission (NRC).

(b) During reflooding, since the CHF is only used indirectly and we are not actually

interested in the exact value of this quantity, it would be beneficial to use a

simple CHF correlation similar to the one used in this work, which exhibits a

smooth behaviour. This will greatly damp oscillations excited during reflooding

whose origin can sometimes be traced back to the complex functional form of

the CHF used in the transition boiling.

(c) Implementation of the high void fraction single-phase limit in the momentum

equations as explained in this work.

(d) Modification of the Co and C1 coefficients in the interfacial shear subroutine

FRCIF ONLY IF there is really CCFL.

(e) Implementation of the right annular flow interfacial shear in the code by changing

the liquid density in the correlation to the vapour density (Wallis correlation

as in TRAC-PF1 and RELAP5/MOD3).

(f) Changing the wall friction model of the code which seems to induce a number of

problems in favour of a model similar to the one in TRAC-PF which does not express the wall friction in terms of total pressure drop (and hence, subsequently

requiring partitioning of the shear between the two phases), but separately as

wall friction on the vapour and the liquid. Ideally, one would need a model

which eg in annular flow, would account for the fact that the wall friction Is

ONLY between the wall and the liquid. This modification may also render the

upwinding scheme of some terms of the phasic momentum equations reported

in this work unnecessary. Both items (e) and (f) will be parts of the new version

of the code.

(g) A further investigation of the advantages of using a second-upwind type scheme

for the convective terms of the phasic momentum equations is worth pursuing. Furthermore, the possible extension of the code for accommodating an

additional liquid field would be desirable.

Acknowledgements

We would like to thank Paul Coddington with whom we had a number of discussions

in the past. In particular, the investigation of the effect of upwinding some terms of the

phasic momentum equations was performed jointly with him, while a number of recent

discussions led to the study of the shortcomings of the wall shear model in the code.

25
References

G. TH. ANALYTIS, Trans. Amer. Nucl. Society, 52, 481 (1986).

G. TH. ANALYTIS, Developmental Assessment of TRAC-BF1 with Separate-Effect and Integral

Reflooding Experiments. NURETH-5 Conference, 21-24/9/1992, Salt Lake City, Utah, U.S.A., Vol. I,

pp 287 - 293 (1992).

G. TH. ANALYTIS, Trans. Amer. NucL. Soc., 73, 512 (1995).

G. TH. ANALYTIS, Development and assessment of a modified version of RELAP5/MOD3, NURETH7 Conference, 10- 15/9/1995, Troy, New York (1995).

G. TH. ANALYTIS, Developmental Assessment of RELAP5/MOD3.1 with Separate-Effect and Integral

Test Experiments: Model Changes and Options, PSI Report Nr. 96-09 (April 1996); Nucl. Engng and

Dsgn, Vol. 163, 125 (1996).

G. TH. ANALYTIS and P. CODDINGTON, The effect of donor cell differencing of some terms of the

momentum equations of TRAC-BF1 on code predictions, NURETH-5 Conference, 21- 24/9/1992, Salt

Lake City, Utah, U.S.A., Vol. H, pp 2 87 - 293 (1992).

J.G.M ANDERSEN, K.H. CHU AND J.C. SHAUG, BWR Refill-Reflood Program Task 4.7 - Model

Development: Basic Models for the BWR Version of TRAC, NUREG/CR-2573 (1983).

D. BESTION, Interfacial Friction Determination for the 1D-6-Equations Two-Fluid Model used in the

CATHARE Code. Paper presented at the European Two-Phase Flow Group Meeting, Southampton,

U.K., June 3-7, 1985.

D. BESTION, The Physical Closure Laws in the CATHARE code, nucL Engng and Dsgn, 124, pp 229 -

245 (1990).

J. A. BOROWSKI AND N. L. WADE (Editors), TRAC-BFi/MODI: An Advanced Best-Estimate

computer Program for BWR Accident Analysis: Model Description, NUREG/CR-4356, EGG-2626, Vol.

1 (1992a).

J. A. BOROWSKI AND N. L. WADE (Editors), TRAC-BF1/MOD1 Models and Correlations,

NUREG/CR-4391, EGG-2680, (1992b).

K. E. CARLSON et. al., RELAP5/MOD3 Code Manual Volume IV: Models and correlations (Draft),

June 1990.

S. GAO and D. C. LESLIE, Two-fluid Modeling of Vertical Annular Two-Phase Flow, Report AERER13633, AERE Harwell (1989).

H. GROTER, F. STIERLI, S. N. AKSAN and G. VARADI, NEPTUN Bundle reflooding experiments:

Test Facility description, EIR Report Nr. 386 (1980).

26
M. ISHII, One-Dimensional Drift-Flux Model and Constitutive equations for Relative Motion Between

Phases in Various Two-Phase Flow Regimes, ANL-77-47 (1977).

D. JUHEL, A Study on Interfacial and Wall Heat Transfer downstream a Quench Front, in "The First

International Workshop on Fundamental Aspects of Post-Dryout Heat Transfer," NUREGICP-0060, pp

179-203 (1984).

D. R. LILES et. al., TRAC-PF1I/MOD1: An Advanced Best-Estimate computer Program for Pressurized

Water Reactor Thermal-Hydraulic Analysis NUREG/CR-3858, LA-10157-MS, 1986.

V. H. RANSOM et. al., RELAP5/MOD2 Code Manual Volume 1: Code Structure, System Models and

Solution Methods, NUREG/CR-4312, (1985).

D. D. TAYLOR et. al., TRAC-BDI: An Advanced Best-Estimate Computer Program for Boiling Water

Reactor Transient Analysis, NUREG/CR-3633 (1984).

27
2.00- I I I.

1.•0-

1.25-

120 1.o0- -J

0.. 2 .. 00. .. .... .. .5. 

92001 A900.

100A0

IO..j 700. 2

Fig00. Pt .N 50.-

m(A)ad096m() -): Modified TRAC-BF1; ( ...... ): frozen -)mesren. version; 

o aue et

O0D . "1 I.. 0• !o .•... ..

00.

001:1

.. I I I,-.IOI

0. =50. 100. "150. 200. 2S 300. 0. 0. 00. 150 20. 258. 300.

TIMr (s) TIME (s>

Fig.2 NEPTUN bottom flooding exp. Nr. 5036. RSTs (K) at axial elevations 0.714

m (A) and 0.946 in (B). (--): Modified TRAC,-BF1; ( .... ): frozen version; (- --

measurements.

28
a0

0

0

to

I

z

0

0

too.-

7S0.-

600.-

4S0.-

300.-

150.-

. . .I - -I - - i. . .

750.-

A0

0

0

-j

0

0:

........ H ......... a ... ......... ...................

00.-

300.-

ISO.-

0.

0.

B -

PI AN*PF 504 10oIO 20 0 . . . .I. . . . II I I 0. 50. 1n0. 0. 1 200. 250. -00.

TIME ($) so. 1oT. ISO. 2C0. TIME: (s) 250. 300.

OD

0

0

0

0

•0

ID

0.

..J

4

0

0

0

CL

4O

I

0.

Fig.3 NEPTUN bottom flooding exp. Nr. 5036. HTC to the liquid (W/m 2/K)

at axial elevations 0.714 m (A) and 0.946 m (B) and void fractions neax the same

elevations (C,D). (-): Modified TRAC-BF1; ( ...... ): frozen version.

29
t. .!

;,; 5•,• .4

0.25-

0.00 00 r ... .. .

0. m0. 20. z0. - 30. 40. 50. 

TIME (a)

Fig.4 NEPTUN bottom flooding exp. Nr. 5050. Collapsed Liquid Level (m).

-): Modified TRAC-BF1; ( ...... ): frozen version, ITMIN=0; (-.-.-): frozen version,

ITMIN=1, (- -- ): measurements.

200122. a 1 00. .. .. I .

1100. 11 500-

. A 1 B

goo.- .. 0 0. 0.. "....... . 5000.

0 I

000. ' ~ -%7.0.

o *\.?o-

- 80. %" ""..-. •. 1 1 "

o I 800

50w. 500.

I ;I

1 I l

N . __ I-. 

0. 50 20.,0 40. 50. 60. I

TIM "( ). I0. 20. "30. 40. SO. 40. TIE TIME )

Fig.5 NEPTUN bottom flooding exp. Nr. 5050. RSTs (K) at axial elevations

0.714 m (A) and 0.946 m (B). (---): Modified TRAC-BFi; ( ...... ): frozen version,

ITMIN=0; (-.-.-): frozen version, ITMIN=1, (--- ): measurements.

30
I., ]" •ir i pii

0.25-

0.00 I

0. 20. 40. • o. :. 0. 20. a0. Io.

TIME (s)

Fig.6 NEPTUN bottom flooding exp. Nr. 5052. Collapsed Liquid Level (m)

(-): Modified TRAC-BF1; (......): frozen version; -- -): measurements.

1200. .

oA '000- B

30900.1

A B

too. 3O0.--

S 800.

000. 1.. 700

800 •. :0 I/ 1i. a. 1o Io. lo 3. to

Fi.7NPTN oto fodig x. r 5002.]P ~ K taileeain .1

.00.

300. I. 0. .0I 

0. 20. 0.0 . . Sol. 3010. 120. 340 18.20 . 80 TM 80 *0. 30 4. 30

TIME (s) TM

Fig.7 NEPTUN bottom flooding exp. Nr. 5052. RSTs (K) at axial elevations 0.714

m (A) and 0.946 m (B). (-): Modified TRAC-BF1; (...... ): frozen version; (--

measurements.

co

0

0

030

In

31
" " 1 l - - - i - • . [ . . , I v'J' l . . 1 ` . . . . . . .

• I I I I I

00.-

?so.-

450.-

0

0

Co

0I

0

Cr

300-

ISO.-

A

I i * A.I

0

0

0

Io

-I

0

h1 1 0. 20. 40. 40. 80. 100. 120. 140. 160.

TIME (s)

100

0

0

0

I

0(

-J

0

0

0

-J

Fig.8 NEPTUN bottom flooding exp. Nr. 5052. HTC to the liquid (W/mn2 /K)

at axial elevations 0.714 m (A) and 0.946 m (B) and void fractions near the same

elevations (C,D). (-): Modified TRAC-BF1; ( ...... ): frozen version.

32
I.J

-j 0.73-

0.00

0. so. 100. 'ISO. ZOO. 25o. 300.

TIMIt C()

Fig.9 NEPTUN bottom flooding exp. Nr. 5036. Collapsed Liquid Level (m)

(-): Modified TRAC-BF1; ( ...... ): Modified TRAC-BF1 with the annular flow

interfacial shear of RELAP5/MOD2; (- - -): measurements.

1100.-

1000.-

.00.-

SOC.-

700.-

go

0

0

0

1zD

S.-

0

0

0

ID

B -

IX

I -

I I -

1~

600-

SOO.-

.00.

300. SI I

10. 0lc. "11. 100. 250. 300.

TIME (s)

Fig.10 NEPTUN bottom flooding exp. Nr. 5036. RSTs (K) at axial elevations 0.714

m (A) and 0.946 m (B). (-): Modified TRAC-BF1; ( ...... ): Modified TRAC-BF1

with the annular flow interfacial shear of RELAP5/MOD2; --- -): measurements.

33
w td /' if /- -.

0.75-

0.00

TIMI (S)

Fig.11 NEPTUN bottom flooding exp. Nr. 5052. Collapsed Liquid Level (m)

(-): Modified TRAC-BF1; ( ...... ): Modified TRAC-BF1 with the annular flow

interfacial shear of RELAP5/MOD2; (- - -): measurements.

I I I I I U

0

0

0

(C

I-

,..

1100.-

1000.-

I0O0.-

700.-

G00.-

A

%%\

0

0

to•0

U.-

SO".-

"00. 0. _.11_..i..

I;. -A.

TIME (S) TIME (s)

Fig.12 NEPTUN bottom flooding exp. Nr. 5052. RSTs (K) at axial elevations 0.714

m (A) and 0.946 m (B). (-): Modified TRAC-BF1; ( ...... ): Modified TRAC-BF1

with the annular flow interfacial shear of RELAP5/MOD2; (-- -): measurements.

34
.U ...... ... 1 - - U.....1 1

-

I I I I I I

1.75-

1.50-

1.25-

LLJ

-j

1.00-

0.75-

0.50-

0.25-

0.00

a.

P

I .1- 11- ... . . .. V

50. 100. 150.

TIME

I

200.

(S)

I I I I I I I ] I I .................

250. 300. 350.

Fig.13 Benchmark top flooding test. Collapsed Liquid Level (m). (-): Modified

TRAC-BF1 with Wallis interfacial shear correlation and the annular flow wall shear

of TRAC-PF1; ( ..... ): Modified TRAC-BF1 with Wallis interfacial shear correlation;

(-.-.-): frozen version.

35
t2 ...

0

0

0

%0

0

0

0

2200o.-

o00.-

SM0.-

300.

C ,00. 2. 000. 250 300 300.

TIME (s)

12, ,.

001

0

0

UIC

000.-

--- ----------------------------

0

a

300.

0. 00. o00. '30. no. 250. 300. 35

TIME (s)

O.

1ao-. .

1100-

i I -" |

N

0

0

*0

0

0o

000..

500.

700.-

600.-

500.-

.. 1, 30&.1- i .

0. 0 200. 200. 200. 050. 300. U Sa30. e . 2. 300. 206C 300• 350.

TIME (s) TIME (s)

Fig.14 Benchmark top flooding test. RSTs (K) at 6 axial elevations. (-): Modified

TRAC-BF1 with Wallis interfacial shear correlation and the annular flow wall shear

of TRAC-PF1; ( ..... ): Modified TRAC-BF1 with Wallis interfacial shear correlation;

(-.-.-): frozen version.

36
0' 0.0 MAJ - 11U l

S !.i p] }: ' t. 1 ";5l .F 1 , -.. "IJ S ~ • iv 

•Eif.j. I '

F i - ;

T0MEoo) .. 0 0 2. ) j

! --

2.00

ii ' i *il ,

>D >

0. ZSQ!0. $0 .

TIM ()TIM

. I 1

ZI *Z l•Eji•:

ii

J 

J

H -0* ~I L' 

|IlO lll-'SO

0. SO. 7-0. ". 0. 0. 0 .

TIME (z) TIME (s)

Fig.15 Benchmark top flooding test. Liquid velocities at 6 axial elevations. (--):

Modified TRAC-BF1 with Wallis interracial shear correlation and the annular flow

wall shear of TRAO-PF1; (....): Modified TRAC-BF1 with Wallis interfacial shear

correlation; (-.-.-): frozen version.

37
0

"" - 800.

0

Ia. 700.-"

(//

600.-

400.-

300.0 0 . .. . .. I. . . . .I. . . . . . . . .. I......... ......... . . . .

0. 100. 200. 300. 400. S00. 600. 700.

TIME (s)

Fig.16 TLTA Test Nr. 6423. Comparison of measured (- - -) and predicted (-)

"maxdmum" RST histories with the frozen version of the code and the fine-mesh

activated during reflooding.

38
.11 C

0qN

U'

900.-

Boo.i

700..

4100.-

SOO.

500.-

*00.-

1 4 I

I

I -, '

I I

3 / I

a ~ / ii /

/ I

'S

I I I . ..

I I 4 4 300. 0. 10. 00o. Is6. 200. 210.

TIME (a)

I I . -. 300. iSO. AM.

Fig.17 TLTA Test Nr. 6423. Comparison of measured (- - -) and predicted maximum RST histories. Andersen's annular flow shear model. (-): Modified TRACBF1; ( ...... ): Modified TRAC-BF1 with upwinding option.

SOO.-

700.

40.-

In

0•

I .. ..... i . . . . .i TI ..... . • ... . .. 1 . ... .. i . .. .. i ... ..

0. $0. t00. SSG. .200. 250. 300. 3110.

TIME (a)

4OO.

Fig.18 TLTA Test Nr. 6423. Comparison of predicted RST histories. Andersen's

annular flow shear model. (-): Modified TRAC-BF1, upwinded both in 1D and

3D components; ( ...... ): Modified TRAC-BF1 upwinded only in 1D components.

39
9 0 0 .. ......... I. .. .. .. ... ...... I......... I. ........ . . ......... . . . .........

900.- T-rrr~ 7-~100,r

800.-

o •I.-. !.

o 600o.7ii•;/ /

-- 500.-

II

400.-

300........ r

0. 50. 100. 150. 200. 250. 300. 350. 400.

TIME (s)

Fig.19 TLTA Test Nr. 6423. Comparison of measured (- - -) and predicted maximum RST histories with the Wallis annular flow interfacial shear correlation. (-):

Modified TRAC-BF1 with homogeneous TRAC-PF1 wall shear model; (..... ): Modified TRAC-BF1 with annular flow TRAC-PF1 wall shear model; (-.-.-): Modified

TRAC-BF1 with no liquid fraction weighting of the wall shear to the liquid.

40
NRC FORM 335 U.S. NUCLEAR REGULATORY COMMISSION 1. REPORT NUMBER

(2419) (Assigned by NRC, Add V.a, Supp., Rev.

NRCM 110Z and Addendum Numbers. If any.) 32oI.320 BIBILIOGRAPHIC DATA SHEET

2. TIJE AND SUBTrm.E NUREGAIA-0146

Implementaton and Assessment of Improved Models and Options in TRAC-BF1 3. DATE REPORT PUBLISHED

MONTH YE

October 1998

4. FIN OR GRANT NUMBER

W6667

5. AUTHOR(S) 6. TYPE OF REPORT

G. Th. Analytis Technical

7. PERIOD COVERED P-custve Dates)

8. PERFORMING ORGANIZATION - NAME AND ADORESS gVA ,C pmvide DLvivko OftR oreRoiom U.& NucMIowRegubbxy Ca on and mileing address,; fonacor.

provide name and maiing dus.)

Laboratory for Thermal Hydraulics

Paul Scherrer Institute

CH-5232 Villigen PSI

Switzerland

Q. SPONSORING ORGANIZATION - NAME AND ADORESS #,RC, tye 'Sam. as abow, ifonkackr provide AMC Dvaix OfRce Regior n U.S. Nutmar ftgulaty Commisaiw

and mulang -d'&-.)

Office of Nuclear Regulatory Research

Division of Systems Technology

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

10. SUPPLEMENTARY NOTES

D.D. Ebert NRC Project Manager

11. ABSTRACT (20o wods or bu)

A summary of modifications and options introduced in TRAC-BF1 Is presented a nd is shown that the predicting capabilities

of the modified version of the code are greatly Improved. These changes includ e the Introduction of a different heat transfer

package during reflooding, the implementation of a simple single-phase limit procedure for forcing the two phases to acquire

the same velocity if one phase disappears, a close assessment of the annular fI ow Interfacial shear correlation, implementation

of a simple radiation model which seems to alleviate some numerical-oscillation problems induced by the existing highly

complex model. Furthermore, different options were Introduced and tested like upwinding some terms of the momentum

equations, the second upwind scheme for the convective terms of the phasic mome ntum equations and the implementation

and assessment of a completely different annular flow Interfacial shear correlation.

Modified TRAC-BF1 Is assessed against some bottom-flooding separate-effect expe riments, a "benchmark" top flooding

simulation as well as against the TLTA test No. 6423. In the process of this task, the different options are assessed and

discussed and is shown that the predictions of the modified code are physicall y sound and close to the measurements, while

almost all the predicted variable are free of unphysical spurious oscillations. The modifications introduced solve a number of

problems associated with the frozen version of the code and.result in a version which can be confidently used for LBLOCA

analyses.

12. KEY WORDS/DESCRIPTORS (L#Mwords aprsses w assWist ws Inw ocaftr iMe a'l) IS. AVALAB•JTY STATEMENT

TRAC-BF1 unlimited

LOCA 14. SECURrnY CLASSaFAnON

LBLOCA (This PNG,)

reflood unclassified

TLTF MThs Rep

assessment unclassified

interfacial shear 15. NUMBER OF PAGES

radiation model

16. PRICE

NRC FORM 335 (2-")
Federal Recycling Program

NUREG/IA-0146 IMPLEMENTATION AND ASSESSMENT OF

IMPROVED MODELS AND OPTIONS IN TRAC-BF1

OCTOBER 1998

UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, DC 20555-0001

OFFICIAL BUSINESS

PENALTY FOR PRIVATE USE. $300

FIRST CLASS MAIL

POSTAGE AND FEES PAID

USNRC

PERMIT NO. G-67

I SYSTEM'S SR

Dc 05