NUREG/IA-O 151

International

Agreement Report

Verification of RELAP5/MOD 3

With Theoretical and Numerical

Stability Results on Single-Phase,

Natural Circulation in a

Simple Loop

Prepared by

J. C. Ferreri, ARN

W Amnbrosini, USP

Autoridad Regulatoria Nuclear Universita degli Studi di Pisa, Facolta di Ingegneria

Av. del Libertador 8250 Dipartimento di Costruzioni Meccaniche e Nucleari

1429 Buenos, Aires, Argentina Via Diotisalvi 2, 56126 Pisa, Italy

Offce of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

February 1999

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

P.O. 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/ordernow>

703-487-4650

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

international agreements, (2) brochures, (3) pro- ceedings of conferences and workshops, (4) adju- dications 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:

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

SR1 350N9/lpdr/htmi>

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, from--

American 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

agreement for the exchange of technical information. Neither

the United States Government nor any agency thereof, nor any

of their employees, 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 information, apparatus, product or process disclosed in this report or represents

that its use by such third party would not infringe privately

owned rights.
NUREG/IA-015 1

International

~ Agreement Report

Verification of RELAP5/MO4D 3

With Theoretical and Numerical

Stability Results on Single-Phase,

Natural Circulation in a

Simple Loop

Prepared by

J. C. Ferreri, ARN

W. Afinbrosini, USP

Autoridad Regulatoria Nuclear Universita degli Studi di Pisa, Facolta di Ingegneria

Av. del Liertador 8250 Dipartiinento di Costruzioni Meccaniche e Nucleari

1429 Buenos, Aires, Argentina Via Diotisalvi 2,56126 Pisa, Italy

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

February 1999

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

VERIFICATION OF RELAP5/MOD3

WITH THEORETICAL AND NUMERICAL STABILITY RESULTS

on single-phase, natural circulation in a simple loop

ABSTRACT

The theoretical results given by Pierre Welander [1] are used to test the capability

of the RELAP5 series of codes to predict instabilities in single-phase flow. These

results are related to the natural circulation in a loop formed by two parallel

adiabatic tubes with a point heat sink at the top and a point heat source at the

bottom. A stability curve may be defined for laminar flow and was extended to

consider turbulent flow. By a suitable selection of the ratio of the total buoyancy

force in the loop to the friction resistance, the flow may show instabilities. The

solution was useful to test two basic numerical properties of the RELAP5 code,

namely: a) convergence to steady state flow-rate using a "lumped parameter"

approximation to both the heat source and sink and, b) the effect of nodalization to

numerically damp the instabilities. It was shown that, using a single volume to

lump the heat source and sink, it was not possible to reach convergence to steady

state flow rate when the heated (cooled) length was diminished and the heat

transfer coefficient increased to keep constant the total heat transferred to (and

removed from) the fluid. An algebraic justification of these results is presented,

showing that it is a limitation inherent to the numerical scheme adopted.

Concerning the effect of nodalization on the damping of instabilities, it was shown

that a "reasonably fine" discretization led, as expected, to the damping of the

solution. However, the search for convergence of numerical and theoretical results

was successful, showing the expected nearly chaotic behavior. This search lead to

very refined nodalizations. The results obtained have also been verified by the use

of simple, ad hoc codes. A procedure to a issess the effects of nodalizations on the

prediction of instabilities threshold is outlined in this report. It is based on the

experience gained with the aforementioned simpler codes.

III

CONTENTS

Abstract i

Contents v

List of Tables vii

List of Figures vii

Executive Summary xi

Acknowledgments xiii

Nomenclature xv

1. Introduction 1

11. Theoretical Analysis 2

Ill. Numerical Results 6

IV. Conclusions 14

References 15

APPENDIX A 26

Internal Memorandum on "Limitations of the Use of Heat

Exchanger Approximation for a Point Heat Source

APPENDIX B 38

Theoretical Stability Analysis for Turbulent Natural

Circulation in a Simple Loop

APPENDIX C 43

-Results From the Application of RELAP5/MOD3.x Code Plots of Selected Results -List of Cases Considered and Their Nodalizations

APPENDIX D 64 Input Deck for Welander's Problem

V

LIST OF TABLES

Table 1 Specification of a Test Problem

LIST OF FIGURES

Figure C.la TEST #1 Base Case using RELAP35/IMOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1 m, Heated Length = 1 m,

Initial Mass Flow Rate=O. 1 Kg/s)

Figure C.lb TEST #1 Base Case

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure CAla.

Figure C.2a TEST #2 Sensitivity of results to heated length using RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1im, Heated Length = 0.5m,

Initial Mass Flow Rate=0.1 Kg/s)

Figure C.2b TEST #2 Sensitivity of results to heated length

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.2a.

Figure C.3a TEST #3 Sensitivity of results to heated length using RELAPSIMOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1 m, Heated Length = 0.1 m,

Initial Mass Flow Rate=0.1 Kg/s)

Figure C.3b TEST #3 Sensitivity of results to heated length

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.3a.

Figure C.4a TEST #4 Sensitivity of results to node number using RELAP25/IMOD3.0

Mass Flow Rates in Legs and TDV junction, 25 nodes per leg,

(Length of Components = 0.4m, Heated Length = 0.1lm,

Initial Mass Flow Rate=0.1 Kg/s)

Figure C.4b TEST #4 Sensitivity of results to node number

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.4a.

Figure C.5a TEST #5 Sensitivity of results to node number using RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 50 nodes per leg,

(Length of Components = 0.2m, Heated Length = 0.l1m,

Initial Mass Flow Rate=O.1 Kg/s)

Figure C.5b TEST #5 Sensitivity of results to node number

ViI
Figure C.6a

Figure C.6b

Figure C.6c

Figure C.6d

Figure C.7a

Figure C.7b

Figure C.7c

Figure C.7d

Figure C.7e

Figure C.7f

Figure C.8.a

Figure C.8b

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.5a.

TEST #6 The recovery of unstable behavior using RELAP5/MOD3.O

Mass Flow Rates in Legs and TDV junction, 90 nodes per leg,

(Length of Components =0. 11 m, Heated Length = 0. 1 m,

Initial Mass Flow Rate=0.1 Kg/s)

TEST #6 The recovery of unstable behavior

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.6a.

TEST #6 The recovery of unstable behavior

Source Fluid Temperature vs. Mass Flow Rate,

conditions as in Figure C.6a.

TEST #6 The recovery of unstable behavior

Sink Fluid Temperature vs. Mass Flow Rate, conditions as in Figure C.6a.

TEST #7 The recovery of unstable behavior starting from (almost) rest

using RELAP5IMOD3.0

-Mass Flow Rates in Legs and TDV junction, 90 nodes per leg,

(Length of Components =0.1 1 m, Heated Length = 0.1 m,

Initial Mass Flow Rate = 0.001 Kg/s)

TEST #7 The recovery of unstable behavior starting from (almost) rest

Mass Flow Rate in Ascending Leg, 90 nodes per leg,

(Length of Components =0. 11 m, Heated Length = 0.1 m,

Initial Mass Flow Rate = 0.001 Kg/s)

TEST #7 The recovery of unstable behavior starting from (almost) rest

Fluid Temperatures at Source, Ascending Leg Cell 5, Sink and

Descending Leg Cell 5, conditions as in Figure C.7a.

TEST #7 The recovery of unstable behavior starting from (almost) rest

Fluid Temperatures at Source and Sink, conditions as in Figure C.7a.

TEST #7 The recovery of unstable behavior starting from (almost) rest

Source Fluid Temperature vs. Mass Flow Rate,

conditions as in Figure C.7a.

TEST #7 The recovery of unstable behavior starting from (almost) rest

Sink Fluid Temperature vs. Mass Flow Rate, conditions as in Figure C.7a.

TEST #14 Search for Neutral Stability using RELAP5IMOD3.O

Mass Flow Rates in Ascending Leg and toward TDV, 37 nodes per leg,

(Length of Components =0.1 1 m, Heated Length = 0.1 m,

Initial Mass Flow Rate = 0.001 Kg/s)

TEST #14 Search for Neutral Stability

conditions as in TEST #7

Fluid Temperatures at Source and Sink conditions as in Figure C.8a.

Vill
Figure C.8c:

Figure C.8d

Figure C.9a

Figure C.9b

Figure C.9c

Figure C.1 0

TEST #14 Search for Neutral Stability

Source Fluid Temperature vs. Mass Flow Rate,

conditions as in Figure C.8a.

TEST #14 Search for Neutral Stability

Sink Fluid Temperature vs. Mass Flow Rate, conditions as in Figure C.8a.

TEST #19 Sensitivity of results of TDV isolation using RELAP5IMOD3.O

Mass Flow Rates in Ascending Leg, 37 nodes per leg,

(Length of Components = 0.263m, Heated Length = 0.1lm,

Initial Mass Flow Rate=0.001 Kg/s)

TEST #19 Sensitivity of results of TDV isolation at 3000s

Flow rates to/from TDV, Test #9 with valve left open up to 9000s

Conditions as in Figure C9.a

TEST #19 Sensitivity of results of TDV isolation at 3000s

Flow rates in Ascending Leg, Comparison of Tests #9 and #19

Conditions as in Figure C9.a

TEST #14 Search for Neutral Stability, conditions as in TEST #7

Mass Flow Rates in Ascending Leg using RELAP5/MOD3.O and MVOD2,

90 nodes per leg, Length of Components =0.1l1m, Heated Length = 0.1lm,

Initial Mass Flow Rate = 0.00 1 Kg/s

Ix

EXECUTIVE SUMMARY

Theoretical results dealing with natural circulation in single-phase in a loop

formed by two parallel adiabatic tubes with a point heat sink at the top and a point

heat source at the bottom are considered. They are used to test the capability of

the RELAP5 series of codes to predict instabilities in these flows.

The solution was useful to test two basic numerical properties of the

RELAP5 code, namely:

a) Convergence to steady state flow-rate using a "lumped parameter"'

approximation to both the heat source and sink

b) The effect of nodalization to numerically damp the instabilities.

It was shown that, using a single volume to lump the heat source and sink,

it was not possible to reach convergence to steady state flow rate when the

heated (cooled) length was diminished simultaneously keeping constant the total

heat transferred to (and removed from) the fluid. This is due to the numerical

properties of the finite-differences scheme used.

Concerning the effect of nodalization on the damping of instabilities, it was

shown that "reasonably fine" nodalizations led, as expected, to the damping of the

solution. However, the search for convergence of numerical and theoretical results

was successful, showing the nearly chaotic behavior of the system. This search

led to refined nodalizations. The results obtained have also been verified by the

use of simple, ad hoc codes.

A procedure to assess the effects of nodalizations on the prediction of

instabilities threshold using RELAP5/MOD3.x is proposed. This procedure is

based on the experience gained with simpler codes used and is outlined in this

report.

X1

ACKNOWLEDGEMENTS

Work was started while one of the authors (J.C.F.) was on leave at the Universitci

di Pisa at the end of 1994. The partial support of the International Atomic Energy

Agency for his stay is gratefully acknowledged. The kind and friendly collaboration

of the Staff of the Dipartimento Construzioni Meccaniche e Nucleari is also

gratefully recognized. Subsequent work has been done in both institutions. The

results are the product of a f ruitf ul, non-formal collaboration. The CAMP provided

the intellectual frame for this continued effort.

.xtII

NOMENCLATURE

Roman Letters

A

a, b

C

Op

D

Fnod(q)

f

A

g

H

J

L

N

Q

q

R

Re

S

T

Ts

t

z

Area

Coefficients of the friction law

Courant number

Specific heat at constant pressure

Diameter

Diffusion coefficient

Source-sink heat transfer multiplier for the

nodal solution

Friction factor

Amplification matrix

Gravity

Heat transfer coefficient

Jacobian matrix

Loop length

Number of nodes

Volumetric flow rate

Dimensionless volumetric flow rate*

Friction parameter (from Welander [1])

Reynolds number based on pipe diameter

Axial coordinate along the loop

Fluid temperature

Fluid ternperature in the source

Time

Eigenvalue

[M2]

[J/(in3 K)]

[-1

[in3 s]

[-1

I-]

Greek Letters

V

p(A)

Subscripts

ss

w

Wel

0

Superscripts

n

n+1

Buoyancy parameter in dimensionless

momentum equation

Isobaric expansion coefficient

Friction parameter in dimensionless

momentum equation

Fluid dynamic viscosity

Fluid kinematic viscosity

Spectral radius

steady-state value

wall

Welander

reference or initial value

[kgl(ins)]

[M2/s]

1-1

Value at time level n

Value at time level n+1

XV
ss Steady-state value

Abbreviations

US NRC United States Nuclear Regulatory Commission

UP University of Pisa

FTUS Forward Time Upwind Space

ODE Ordinary Differential Equation

PDE Partial Differential Equation

XV'
1. INTRODUCTION

The assessment of best estimate, nuclear safety codes is a very important

stage in. their development, including several steps. These steps are embedded in

some accepted procedures like the U.S. NRC Code Scaling, Applicability and

Uncertainty verification methodology [2] or the UP Uncertainty Methodology based

on Accuracy Extrapolation [3]. Both specify the verification of the influence of the

numerical methods implemented in such codes on the results obtained as a part

of the applicability step. Then, defining suitable analytic tests and using them to

set the conditions to satisfy in order to get meaningful results, is an important part

of the verification stage.

The interest of the present work deals with the RELAP5IMOD3.x [4] series

of codes. An important contribution to this subject has been contributed by Wuliff

et al. [5]. It is one of the rare examples of reports dealing with the effects of

discretization errors in the response of thermal-hydraulic system codes in plant

behavior analysis. In a more limited context and despite of its restrictions, the

problem of natural circulation in single-phase flows is quite common in many

situations of interest in the nuclear industry. The removal of decay heat in intact

loops after pump shutdown is only one example of the possible scenarios.

In this report the results of Welander [1] are revisited. These results deal

with the natural circulation in a very simple situation, namely: flow driven by

buoyancy in a loop constructed with parallel, adiabatic, vertical tubes. The loop

has a point heat source in its bottom and a point heat sink in its top. The mode of

heat transfer is such that the product of the heat transfer coefficient in the heated

length times the heated length is kept constant when the heated length tends to

zero. This defines the heating rate. The same concept is applied in the heat sink.

In this way, the physically finite heat source (sink) becomes a point.

Welander [1] performed his analysis only for laminar flow. However, his

results may be generalized to consider turbulent flow. Chen [6], among many

other authors, performed such analysis for a more general loop considering

heated (cooled) horizontal tubes of finite length. In this report the analysis in [1] is

generalized in order to keep strictly the same hypotheses of the original derivation.

The analysis is included as APPENDIX A to this report. The results of such

analysis are applied to verify whether or not the RELAP5 results may converge

when the lengths of the source and sink are varied accordingly to the abovementioned criterion. APPENDIX B to this report gives an algebraic justification for

non-convergence of the upwind space approximation for the heat transfer. Having

defined the theoretical stability map, an unstable flow condition is defined to check

the effect of the nodalization on the appearance of oscillations. In this way it is

possible to define the limits of applicability of a coarse nodalization. It is important

to point out from the very beginning that it is possible to get convergence of the

predicted results to the expected oscillatory analytical results. However, a

nodalization much more detailed than the usually accepted for stable flows must

be used.

These aspects may be also analyzed far more easily using ad hoc codes of

similar numerical properties. Then some results are shown using in-house
developed codes. These results serve to illustrate a procedure that could be

implemented using the RELAP5 numerical approximation. This may be of

importance to quantify the effects of the nodalization in the damping of

perturbations in unstable (or marginally stable) systems.

In what follows the geometry of the loop and the equations that specify the

theoretical problem are briefly discussed. They include their generalization to

consider turbulent flow and the analysis of the lumped parameter approximation

for the point heat source and sink. The theoretical implications of the latter, from

the numerical point of view, are discussed. Then, results obtained using the

RELAP5/MOD3 and comparisons with RELAP5/MOD2 are presented, confirming

the theoretical hypotheses. The discussion of the results and the conclusions

close the report.

II. THEORETICAL ANALYSIS

Figure la,b shows the geometry of this simple hydraulic system, adapted

from [1]. It is a schematic representation of the loop, appropriate for its

discretization following the RELAP5 code specifications. These are shown in

Figure lc,d. The vertical legs are adiabatic, smooth, circular tubes of length U/2

and diameter D. The length of both the heat source and the heat sink is S. The

cross section of the tubes is A. When S tends to zero, the heat transfer f rom/to the

heat source/sink (shaded in Figure lc,d) is increased to keep constant the total

heat transferred. The role of the volume denoted as TDV (a time dependent

volume) will be clarified later. This element, of course, is not present in the

analysis in reference [1].

To deduce the governing equations, the validity of the Boussinesq

expansion of the fluid density in terms of temperature will be assumed. The friction

will be computed from the Poiseuille law for laminar flow and from the Blasius law

for turbulent flow.

Then, the governing equations are:

i) The momentum equation integrated along the loop:

dq ~a-R( D 1e 2-0 JgA f(J)2 (1)

dt 16 (v .A) _L 0

In Eq. (1), the anti-symmetry of the system has been considered.

2
Hi) The energy equation, expressed in terms of temperature:

C [.DT+q DT] {DSH(TF-T) (2)

In equation (2), the first equality holds at the source/sink, and the second indicates

the adiabatic walls of the tubes.

iii) The heat balance at the point heat source/sink:

In the previous equations, written in dimensional form, q is the volumetric

flow-rate, T is the fluid temperature, H is the heat transfer coefficient, v is the fluid

kinematic viscosity, 0 is a coefficient equal to 1 for laminar flow and to 1/4 for

turbulent flow (the exponent in the Blasius friction law), a is the coefficient in the

friction law: it equals 16 for laminar flow and 0.079 for turbulent flow, P3 is the

volumetric expansion coefficient of the fluid and g is the acceleration of the gravity.

The sub-index SI means a condition. evaluated at the input of the source, subindex S means a condition prevailing at the source and sub-index F indicates a

temperature evaluated at the source wall. In Eq. (1) there is a derived coefficient R

and, in Eq. (3), it is y. Their definitions are:

R=32 v and: iD 4

The definition of R takes into account the total frictional force that must exist in the

loop to make the results in [1] coincident with those given -considering the

Poiseullie law. The effect of the assumptions regarding the friction law will be

discussed later.

Setting the temperature at the heat source wall at AT degrees above the

average loop temperature and considering steady state, then, using Equations (1-

3), the flow-rate may be obtained from the following implicit equation:

25-q 2- = &T+B.q 2- .[i1e{ s] (5)

where 8 is defined by:

2a(D 1 4vL

IPgAL

The solution of Eq. (5) must be iterated if the flow-rate does not coincide

with the postulated flow regime. This is not the case when the Reynolds number is

adopted as the dependent variable, instead of flow rate. In non-dimensional terms,

3
Eq. (5) specifies the relation between flow-rate and steady state temperature in

the loop, namely:

2. Q2-0

S= 1.- e-"0

OS (6)

SS

where:

Q=q E)T a= gATL K= 7tDH

(KS)A ' AT' = 2(KS) 2 ' C-A

a (DKS )1-8 RL C F, 6 -v 2I S and 5-we= a (DKS)'-

When 0=1 and a=16 the previous equations reduce to the ones in [1]. The

linearized analysis of the stability of the steady state equations may be performed

in similar terms as in [1]. Now, c shows the influence of the Reynolds number in

the system through the friction law coefficients.

Let us now consider a "lumped parameter" approach to model the heat

transfer from the heat source toward the fluid. Note that the anti-symnmetrical

situation must be considered for the sink. Details of the derivation may be found in

APPENDIX B. The steady state energy equation, evaluated at the source, gives

the value of the temperature in the fluid. For this purpose the steady state energy

equation and an upwind approximation are used.

Then, the value for the fluid steady state temperature at the source is:

E) 1 (7) SS l+2Qss

Note that Qss is only a function of the quotient Wse that is a non-dimensional

parameter depending on the physical parameters of the system and on the

product (KS) that must be kept constant, it may be concluded that the results

obtained applying the r.h.s. of Eq. (2) as it stands will not change using a smaller

node length. For example: if Ws=e-2 and 0=1 (laminar flow), then Qss=O.958 and

e55s=o.479 from the theory. The value, as computed from Eq. (7) is 9ss=0.342.

There is only one way to obtain the correct steady state value of ess, it is to get

the *analytical solution of Eq. (3) averaging over the length S. Following this

procedure and using the upwind approximation for the advective term, the

expression for ess is:

2F

4
The value of ess computed with this expression is coincident with the theoretical

value.

5
Ill. NUMERICAL RESULTS

Three versions of RELAP5 [4] have been applied in the present work,

namely: MOD2, MOD3.O, and MOD3.1. Due to the simplicity of the physical

problem considered (single-phase flow, constant wall temperature and heat

transfer coefficient), no relevant differences in the behavior of the different

versions was noted in relation to the capabilities in predicting the stability

threshold. Then, the results obtained will be reported specifying the code version

just for the sake of completeness. The tests involving the MOD2 version have

been also considered to check alternative, more detailed nodalizations to the

heat/sink volumes in the early stages of this work.

The basic nodalization adopted in the analysis is reported in Figure 1 d. It

consists of two BRANCH components (100 and 130) introduced to simulate the

lower heating and the upper cooling devices. A very thin structure is connected to

each of them, with imposed external temperature, in order to simulate the actual

boundary conditions consisting in an imposed internal wall temperature. The two

BRANCH components are joined with two PIPE components (110 and 120),

simulating respectively the ascending and the descending legs. A SNGLJUN (300)

and a TMDPVOL (400) have been included to compensate for the expansion or

contraction of the fluid volume 'inside the system. However, specific calculation

tests have shown that the dynamic behavior of the system is negligibly affected by

the presence of these components.

The application of the system code to Welander's problem was concerned

with a physical system having the following main characteristics (see other details

in Table 1):

" pipe inner diameter: 0.1 m;

* length of the legs: 10 m;

" length of the heating and cooling sections: 0.1 mn to 1 m;

* overall loop length: 20 m + 2 x heating section length;

* heating surface temperature: 30 00;

* cooling surface temperature: 20 OC;

* heating and cooling surface heat transfer coefficient: 20000 W/(m2K);

* system pressure: 105 Pa;

" pipe surface roughness: 10-7 M.

Appendix C reports the list of the calculation cases considered, together

with the plots of the related results. Most of the cases where run with a time step

of 0.5 s. Various heating lengths and number of nodes were used to test the

convergence of RELAP5/MOD3.x to the theoretical solution.

The topics of main interest were:

6
1 . the treatment of the source and the sink, which in the upwind scheme adopted

by the code cannot be correctly simulated using a single node, since the

exponential trend of the fluid temperature inside the heating and the cooling

sections cannot be adequately accounted for, and

2. the effect of truncation error related to the discretization of the legs by a finite

number of nodes: this is known to introduce spurious dissipative effects that

result in damping the oscillations.

3. the effect of Courant number on scheme numerical phase lags, which may

lead to over prediction of instability for a given nodalization.

The flow was started at very low flow-rate and began to grow as time

elapsed. In damped calculations, the temperature in the loop increased up to the

point in which it reached the average value as given by their theoretical value, the

latter depending on the adopted criteria for the source/sink lumping.

The main conclusions coming from the analysis of the obtained results

were the following:

" The steady-state conditions predicted by Welander cannot be achieved, unless

an appropriate source heat transfer multiplier is included in the code.

" For a given time-step, it is possible to predict both stable and unstable

conditions, depending on the number of nodes adopted to discretize the legs.

The analysis to follow applies to the first conclusion. For these calculations

the following spatial discretization have been considered: 10 nodes in each leg

(Ax = 1 m) and S = 1.0, 0.5 and 0.1 m (correspondingly, the values for H were:

2000, 4000 and 20000 W/ M2 X). It must be pointed out that the flow was unstable

with these parameters. The code was RELA5IMOD3.O. As expected, the value of

the steady state temperature was almost the same, with small differences that

were due to the different total friction in the loop. The effect of the TDV volume

was not important as will be discussed later. Naturally, oscillations in the transition

to steady state were present. The behavior of the system may be observed in

Figure 2 for the fluid temperature at both the heat source and sink. The results

obtained confirmed the results predicted by the theoretical analysis. Then, to

improve convergence toward theoretical values, an average, lumped heat source

must be used, as suggested previously.

Let us now consider the problem of convergence of the code results to the

expected unstable results. To do this, let as suppose that we do not know that the

flow is unstable and proceed as CSAU, UMAE (and engineering judgment)

indicate: look for convergence of computed results. For this purpose the

nodalization considering S = 0.1 m and 20000 W/ M2 X was used. The number of

nodes in each leg was varied accordingly up to the point in which the expected

chaotic behavior was obtained. The latter made its first appearance when the

number of nodes per leg was increased from 25 to 50. Nearly neutral oscillations

in flow-rate were observed for 37 nodes per leg.

7
Figure 3 shows the time variation of flow-rate using 90 nodes per leg and

0.001 Kg/s of starting flow-rate. As may be observed the flow is oscillatory,

switching from positive to negative values after some cycles of oscillation. Fluid

temperature varies correspondingly.

The system behavior may be observed more clearly when the temperature

in the sink is plotted versus flow-rate in the loop. Figure 4 is a clear indication of

the chaotic behavior of the system.

Figure 5 presents results of the calculations, performed with 40 and 60

volumes per leg, still1 using the semi-implicit numerical method of the code. They

were started~close to steady-state conditions. The values of a and F_ corresponding

to the considered physical problem are around 339 and 2.3 respectively,

identifying a physically unstable system (see Figure A.3), and it may be found that

the time-step adopted in the calculation, At = 0.05 s, roughly corresponds to a

dimensionless time-step At=1 0- It can be noted that, with the two numbers of

nodes adopted, stable and unstable conditions are predicted, owing to the

different damping introduced by truncation error.

Simply closing the connecting junction in the middle of an unstable

transient controlled the effect of the TDV. The effect of this action on the flow

variables was negligible and could not be recorded. Also, the same results have

been obtained running the same cases with RELAP5/MOD2. The results were

also verified running a small, ad hoc code using an integrated equation approach

for the momentum equation and a fully up-wind, explicit approximation for the

energy equation. The equivalent, length-averaged heat source/sink approach was

used for lumping, what allowed obtaining of the exact steady state flow-rate value

when possible. The same numerical behavior was observed.

In what follows, a more detailed numerical analysis will be outlined, mak 'ing

reference to the FTUS approximation. This numerical scheme is quite similar to

the one used in RELAP5. Then the algebraic equations expressing the energy

balance are the following:

=r (1 _ C) on + C on,(=,..N

Gn+i 1 1- C cAŽ-Fn (Qn)Qon +Con1 At n

N ~ As n )N N1A o

1= (8)

q q40

O1 =(1 + C)~ E- (i=2,. ,N-1)

Gn11+C _t A FdQn nCn+At nd() A&s ) 1 2AS~ (9)

8
0 n1 _0 n+1

where C is the Courant number:

C Qn -At (10) As

and

N-i

with N being the number of nodes. The function Fnod(q) is a source-sink heat

transfer multiplier, similar to the one adopted for the nodal expansion, introduced

in order to calculate the steady-state conditions in coincidence with the exact

solution.

The momentum equation is discretized in time as follows:

Qfll ýl+._eIo_ - F (Qn)2-e ~ (12)

The steady-state conditions calculated by the method are:

=9S1 1 = ... = E)SS,N.1 = -(ESSN = OSS,leg

eSe 'Ffl(QSS) (13)

SSIg2Q ss + Fn~(Q s)

where:

Fd (QSS) =QSS (1-e-1/Qss) /e-11ass (14)

Thus, the result is:

a 1 - e-"0 2 N-2 (15) QS2v q i/OSS 1 + e-1 /c N-i1

It can be easily shown that for N - othe above equation becomes

coincident with the exact expression given by Equation (5).

The assessment of the effects of nodalizations on the onset of instabilities

will be performed using the previous expressions. However, a's will be mentioned

later, some other methodologies have been also considered. This approach can

be considered the numerical analogue of the classical linear stability analysis of

PDEs. The main reasoning behind the methodology is shortly summarized in what

follows. A finite-difference numerical method for a time-marching problem can

been written as an algebraic fl-vector equation relating the n7 values of the

unknown function at the n-Th. and (n+1)-Th. time level (yfl and yfl+l), grid

9
parameters (in the present case, As and At) and physical parameters (a, s). This

algebraic equation represents the discretized form of the original PIDEs together

with the related boundary conditions. In our specific case for the above-described

numerical methods it is:

F(yn' yn', At, As, a, C) = 0 (16)

It will be now shown that:

" Studying stability of steady-state solutions of a mathematical problem is

capable of clearly pointing out the effect of truncation error brought about by

the nodalization.

* Care must be taken in avoiding numerical instabilities or in recognising them in

the obtained stability maps.

The vector function F is generally non-linear. Therefore, determining the

steady state conditions (i.e., the fixed points) may require the iterative solution of

the equation:

E(yn = ySS Yn1 _ ptysac (17)

Once the fixed points have been determined, their stability can be studied

through linearization by perturbation. Then, considering small deviations from the

selected fixed point:

yn =yS +(5y)fl yfl+l y S + (8y)f+l (18)

Substituting into Equation (16), second order terms can be neglected and

Equation (17) can be used to reach the following relationship between

perturbations at the n-Th. and at the (n+1)-Th. time levels:

(By) n+1 = ...(Jf+l) - .l .J *(By) n (19)

where jnand J"' denote the Jacobian matrices of F with respect to yfl and yfl+l

respectively, calculated at the selected fixed point. It is cleariy understood that the

inverse of Jn+1 must exist for any meaningful time-marching numerical scheme; in

particular, Jn+' can be made equal to the identity matrix for explicit numerical

methods and boundary conditions. It is then argued that stability can be discussed

considering the eigenvalues of the matrix:

A = _(Jf+l)-l .jf (20)

As results assuming exponential growth or decay of perturbation vectors. In

particular, given the spectral radius of the matrix, p(A), it is useful to consider the

quantity:

AP = P(A-1 (21)

quantifying the margin in excess to neutral stab ility, which takes negative values

for stable conditions and positive values for unstable ones. This quantity can be

therefore used to find neutral stability conditions and to set up stability maps.

Then, it is here preferred to calculate Ap throughout a selected a-e rectangular

10
domain, thus identifying with the aid of contour plots regions with a different

degree of stability. This method is easier to implement in computer programs and

has the advantage to provide a greater deal of information, at the price of a

reasonable increase in computing effort.

Figure 6 reports the results obtained for the FTUS method with 30, 40, 50

and 100 nodes and At=10 It can be noted that a very small unstable region is

found within the addressed domain with 30 nodes, whereas increasing the detail

of discretization unstable conditions are predicted for lower and lower values of cc.

This clearly explains the above-discussed behavior and shows the dramatic

quantitative impact of truncation error on the prediction of stability.

Figure 7 illustrates the map for 100 nodes. It is interesting to compare the results

obtained for the explicit upwind method with the results of a modal solution' with a

second order term simulating numerical diffusion. With this aim, the diffusion

coefficient is defined as:

S(Q) IQ As ( 1 J AI 1t) (22)

as resulting from the analysis of truncation error for the FTUS method. Figure 8

shows the predicted flow rate variation with time. The diffusion coefficient was the

one corresponding to 1000 nodes under the FTUS approximation. The number of

modes considered in the modal expansion was 500. As may be observed, both

approximations behave similarly. Exact coincidence is precluded by the nonlinearity of the system and the approximate values given by expression (6). The

latter fact is due to the use of the nominal, steady state value of Q used for its

calculation. It may be shown that the predicted stability conditions are very similar

for the nodal and the modal solution with equivalent dissipative effects and the

agreement is improved by increasing the number of nodes. This confirms the

overwhelming importance of the second order term alone in determining the

overall truncation error effect on stability predictions.. Figure 6 also shows the

linear stability curve obtained by the modal solution with no diffusion (B(Q)=0). It

shows, as expected, its close agreement with the stability curve obtained by the

conventional linear stability analysis shown in APPENDIX A.

Finally, Figure 9 shows the stability maps for various 1st and 2nd order

methods obtained with At=1 0'4. It is clearly visible that 2nd order methods provide

relatively very accurate predictions of the stability boundary. The changes

observed in the maps increasing the number of nodes up to 100 are minimal,

supporting the conclusion that in the present case the effect of truncation error on

stability prediction is due almost exclusively to the second order dissipative term.

The low Courant number used makes the FTUS results almost as diffusive as the

ITUS ones.

This consists in a modal expansion of the governing equations and boundary conditions (1-3) in

terms of Fourier series of sinus and cosinus. Then, a coupled system of ODEs is obtained, which is

free of spatial discretization error. Considering a reasonably high number of modes, a reference

solution may be obtained. Then, the energy conservation equation may be generalized to consider

a diffusive term, with diffusion depending on flow rate.

I V
The analysis of the stability maps obtained with the FTUS method

and 0=0.8 or with second order methods and At = 10-3, pointed out the

appearance of regions of instability larger than those predicted by the reference

modal solution. The possible reasons for this behavior were mainly two:

* An effect of truncation error at large time steps on the velocity of propagation

of perturbations along the legs, related to the dispersive effects brought about

by odd-order derivatives;

" The effect of the use of a large time step in momentum equation, giving rise to

a general tendency to instability due to delayed feedback.

In what follows it will be shown that, although the first effect cannot be

excluded, the second is the leading one, being responsible for the observed

overprediction in the extent of the unstable region.

In all the cases considered (FTUS and 2nd order methods), the

adopted discretization for the momentum equation implied that both the

temperature integral along the loop and the friction term have been evaluated

making use of old-time-level quantities. The results obtained making use of this

form are shown in the Figures 15 and 16 below for both the FTUS method. It can

be noted that in the case of FTUS with C=0.8, the overprediction of the extent of

the stability region is remarkable. In the case of the MacCormack method with At=

10-3, a slight discrepancy with the modal solution was also found.

The following approach was used to discriminate among the possible

causes of this convergence "coming from outside" to the neutral curve:

" Considering a fully implicit linearized approximation of the friction term in the

momentum equation: this approach did not lead to any improvement in the

results. On the contrary, the degree of overprediction was larger, showing a

more pronounced dependence on C.

" Evaluating the friction term and the temperature integral term on the basis of

(rn-i)-Th time level values (it did not increase too much the computational

effort):

Therefore, the following form of momentum equation results:

Qfl1Qf~~a e,~ + '(Qfl)leflt

W 2

Making use of this momentum equation in the FTUS scheme, the stability

maps reported in Figures 17 and 18 were obtained. It can be noted that now no

overestimate of the unstable region was obtained, showing the usual diffusive

behavior.

12
The analysis above clearly demonstrates that of the two above envisaged

reasons for overestimating instability, that related to the time discretization of the

momentum equation is certainly dominant.

13
IV. CONCLUSIONS

The aforementioned calculations allowed arriving to the following

conclusions:

i) The analytical results given in [1] are very useful to verify how large predictive

codes, like the RELAP5 series of codes, treat lumped heat transfer structures.

This may also be the case with other codes. Its proper understanding is very

important when the lumping of heat exchanging structures is imperative. Imposing

an appropriate criterion for the lumping avoid convergence errors in these cases.

ii) The effect of numerical diffusion may damp oscillations in a system, even under

unstable flow conditions. Appropriate nodalization permits to recover the correct

physical behavior. Looking for convergence of computed results is always

rewarding in these situations.

iii) The methodology adopted in this work for setting up stability maps shows that a

linear stability analysis based on numerical methods is effective in pointing, out the

effect of truncation error. Furthermore, if the appropriate nodalization detail and/or

higher order schemes are adopted, it can even result as reliable as the usual

frequency-domain techniques.

As shown by the obtained results, the methodology adopted in the present

work is useful to quantify the effects of nodalization and time step choice (i.e.

discretization errors) for system codes like RELAP5. This will be the subject of

further research.

14'
V. REFERENCES

[1] Welander, P., On the Oscillatory Instability of a Differentially Heated Fluid

Loop, J. of Fluid Mechanics, vol. 29, p 17-30, 1967.

[2] Boyack, B.E. et at., Quantifying Reactor Safety Margins. Part 1: An Overview of

the Code Scaling, Applicability and Uncertainty Evaluation Methodology, Nuclear

Eng. and Design, vo~l. 119, p. 1-15, 1991.

[3] D'Auria, F., Debrecin, N. and Galassi, G.M., Outline of the Uncertainty

Methodology based on Accuracy Extrapolation, Nuclear Technology, vol. 109, p.

21-38, 1995.

[4] Carlson, K.E. et al., RELAP5/MOD3 Code Manual Volume 1: Code Structure,

System Models and Solution Methods, NUREG/CR-5535, 1990.

[5] Wulff, W., H.S. Cheng, A.N. Mallen and U.S. Rothagi, BWR Stability Analysis

with the BNL Engineering Plant Analyzer, NUREG/OR 5816, 1992

[6] Chen, K., On the oscillatory Instability of Closed-Loop Thermosyphons, ASME

Trans., J. of Heat Transfer, vol. 107, p. 826-832, 1985

15
PARAMETER

L

LU2

S

A

D

TAV

AT

Tw (at source)

Tw (at sink)

p

H

power to the source

ss mass flow rate

ss vol. flow rate

Reynolds number

ss (T-TAV)

D.K.S/v

U

a

Ewei

qss

VALUE

SI Units

20.00

10.00

0.10

0.00785

0.10

25.00

5.00

30.00

20.0

0.000255

997.00

4177.

0.000864

8 .6660d-07

20000.

2813.31

0.717

.71 93.d-03

10568.

0.522

2216.7

0.09158

338.94

1.4436

2.3027

4.7674

0.1045

TABLE 1 Specification of a test problem

16r
Sink

S =1 A

a)

b)=OL

F", EW- ýA

s=2

Source

Heating

FIGURE 1la,b A schematic view of the natural circulation loop, adapted from [1 ]

17
Figure 1c. The loop geometry in a RELAP5 like definition.

rEý 4 ILI -

100, 130 = Branch

110, 120 = Pipe

300 = Single Junction

400 = Time Dependent Voluir

110 120

1Tý

Figure 1 d - RELAP5 code nodalization for the analysis of Welander's problem

18
VaG,.11 -02-24-19g7

U

S

U

a

0.

Ea

I-

*0

.5

28.00

27.00

26.00

25.00

24.00

23.00

22.00

21.00

20.00

19.00

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Fig. 2 The evolution of fluid temperature at the source

and the sink vs. time using 10 nodes per leg and S=0.1 m, a = 339, c =2.3

WfW M 1 -02-24-1997

2.00

1.50

1.00

a0

2

.50

0

-.50

-1.00

-1IS0

0 1000. 2000. 3000. 4000. 5000. 6000.

Time (s)

7000. 8000. 9000.

Fig. 3 The time variation of loop flow-rate.

Number of nodes per leg: 90, S=Q-1lm, a = 339, c =2.3

19
30

=28

o 26

424

22 _

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5

Loop Flow Rate [kg/s]

Projection of the attractor in the phase space as calculated by

RELAP5/MOD3 using 90 nodes per leg, a = 339, z =2.3

Figure 4 -

0.72

*0.71

~'0.7

0.69

0.68

0.67

0.66

0 1000 2000 3000 4000 5000 6000 7000 8000

Time [s]

Figure 5 - RELAP5/MOD3 results for the dynamic behavior of a single-phase

thermosyphon loop, a = 339, F_ =2.3

20
10

8

6

C

4

2

0

0 200 400 600 800 1000

cc

Figure 6 The effect of the number of nodes on the

Neutral stability curve using the FTUS scheme.

21
Ap

5E-004

4E-004

3E-004

2E-004

IE-004

OE+000

*1E-004

-2E-004

-3E-004

-4E-00

-5E-004

-6E-M

-7E-004

-8E-004

-9E-004

-1E-003

C

100 200 300 400 500 600 700 800 900 1000

cc

Figure 7 Stability map for the FTUS scheme, 100 nodes and At=1 0 -

C

U)

'2

0

U,

2

10

5

0

-5

-10

1000 Nodes __ _ _ _ _ _ _ _. r ii

----- O00Modes + _____ ____

D(1000) ill v

r

0 5 10 15

Dimensionless Time

Figure 8 The flow rate for the FTUS scheme using 1000 nodes and its simulation

using a modal expansion of 500 modes and expression (22)

22
NEUTRAL CURVES USING 100 NODES

10

8

6

I

FE

4

2

0

0 - 200 400 600 800 1000

Figure 9 Comparison of the neutral curves obtained with different approximations

(ITUS: implicit-time, Upwind-Space, MOCOR: Mac Cormack)

23
E

500 600 700 0

(X

Figure 10 - FTUS scheme with C=0.8 and explicit momentum equation (30 nodes)

C

100 200 300 400 500 600 700 800 9 1000

(X

Figure 11 - FTUS scheme with 0=0.8 and explicit momentum equation

(100 nodes)

24
IC

6

C

Theortif SB I

s1

2

A

100 '200" '300 '400 500' 600 700 800 1000o

aX

Figure 12 - FTUS scheme with 0=0.8 and semi-implicit momentum equation

(30 nodes)

I

II - LI I TheOretical SB

4

S

2

$ 100 200-~300 400 500 600

aX

Figure 13 - FTUS scheme with 0=0.8 and semi-implicit momentum equation (100

nodes)

25
APPENDIX A

THEORETICAL STABILITY ANALYSIS FOR TURBULENT

NATURAL CIRCULATION IN A SIMPLE LOOP

In this APPENDIX, the analysis in [1] is generalized to consider turbulent

flow, keeping strictly the same hypotheses of the original derivation. The results

emerging from this analysis have been used in the main text to setup an unstable

problem in such a way that a "reasonably" detailed discretization using the

RELAP5 MOD 3.x [4] codes completely masked the physical situation, giving false

stable solutions.

In what follows the geometry of the loop and the equations that specify the

theoretical problem are briefly discussed.. They include their generalization to

consider turbulent flow. The neutral stability curves are obtained as a function -of a

new parameter, allowing the definition of a theoretical stability map.

The solution of equation (6) (see main text) must be iterated if the flow-rate

does not coincide with the postulated flow regime. In non-dimensional terms

equation (6) specifies the relation between flow-rate and steady-state temperature

in the loop, namely:

T. uniform in the leg

(A. 1-a)

2Ts, ( s -Yl

(A. 1-b)

(A. 1-c)

where the variables have the same meaning than in the main text and, again, subindex ss means a condition evaluated at flow steady-state conditions.

The linear analysis of the stability of the steady state equations may be

performed in similar terms as in [1 ]. To this end we first introduce a perturbation to

steady-state flow conditions in the form:

26
q= qs, + q'

T= T + T'

(A.2)

Introducing these definitions into the non-dimensional form of equations

(A. 1-A.2) we get:

dt SS 0 'd

(A.3)

a~' q a T'=0

(A.4)

and:

T'0+mT'l+nq = 0

(A.5)

where:

M = 1 -TS

n= TS

(A.6)

and the sub-indexes 0 and 0, 1

respectively.

1 means temperatures evaluated at s =

In order to study the growth of perturbations, an exponential form in terms

of time is adopted. Then:

q'= q .er.t

To= T er~t

(A.7)

In the A.7, r is a complex and T is the non-disturbed value of T.

Then:

,aT

T'o + m-T1i-n q= 0

(A.8)

(A.9)

and:

27
I+ me" +

n a qss (1 e- 0rq~

r c,+s qý'O) ] e 0

82-= (2-0) c

(A. 10)

where:

The neutral stability curves may be found putting:

r = iwo

Then, inserting this definition in (A. 10), it is found:

1+ m-eiWIqqq + n a qs

1(0 ~i)+E, q('-) ]

(1 - e- ko/q$) = 0

(A.1 1)

Let us introduce the following definitions:

a = n a / qr

c = F- / qs

(A. 12)

Then, equation (17) becomes:

ei- -1I

a

ro (koc + u Eq1&)

= 0

(A. 13)

Where:

(0)a [ DKS 1i~ [6 J

(A. 14)

Imposing that the imaginary part of (A.13) must equal 0, it results:

-2 2 2(0) -

(C u qSSO + co)(1m - 2a = 0

Furthermore, from equation (1-b):

(A.15 )

1 - - Ln - S

(A. 16)

Then,

28
(A. 17)

Now, equation (A. 15) may be written as:

-2

-2 Ct)

E + u 2 (qss)2(1 -0 A c 0

(A. 18)

where:

(A. 19)

Finally, equation (A. 18) is written as the equation of an ellipse, then:

[ A J2 + -2 [A~

2 (uqss)('bo)=

(A.20)

Furthermore, setting the real part of equation (A.13) equal to zero and using

(A. 18), it is found that:

- 1 1 -CL

- i-e CO -cot g(- T S 2

(A.21)

_ E~quations (A.20) and (A.21) are solved by iteration to find the values of

a and F_ Fixing a value for Tss, what implies a definition for qss, performs it. After

these values are found, the corresponding values of a and F-, as well as the values

of all other variables of interest may also be found.

29
RESULTS AND DISCUSSION

Results have been obtained taking a base case, as specified in Table 1. At

first glance it is not essential to set a physically Sound situation, because the

results may be also interpreted in non-dimensional terms. However, as mentioned

before, the use of the present results to set a test, benchmark problem for a large

thermal-hydraulic system code, led our studies. This approach was useful, so it

was decided to keep the same philosophy. The dimensions of the loop may be

considered typical of a real experimental rig. Results have been obtained keeping

the same geometry and changing the heat flow (as measured by HTC) at both the

source and the sink. The range adopted for the variation of HTC was such that the

flow changed from fully laminar to fully turbulent.

Figure A.1 shows the neutral stability curve obtained by solving equations

(A.20) and (A.21) and considering turbulent flow. Figure A.2 is a stability map,

obtained by solving equations (1 -3) by a modal decomposition method, using 100

and perturbing the resulting ODES according to a standard linear stability analysis

technique. It may be considered a reference solution for the differential problem

and the borderline separating the stable and unstable zones must be coincident

with the one in Figure A.1. In Figure A.3 the results of Figure A.1 have been replotted and are compared with the values obtained from Figure A.2. As may be

observed the values coincide fairly well. The conspicuous point in the middle of

the unstable zone corresponds to the data of Table 1.

'Figure A.4 shows two neutral curves: a) one for laminar flow, obtained by

setting HTC=1 000 and b) one for turbulent flow, obtained by setting HTC=20000

in correspondence with the parameters given in Table 1. Once again, the isolated

point in the unstable zone corresponds to the system working as given by Table 1.

Figure A.5 shows, in terms of Ewel (The parameter in Welander's paper [1])

and qss the neutral stability curves as a function of HTC. In this figure there are

three differentiated types of curves: a) one corresponding to laminar flow

(HTC=1 000), b) one corresponding to a transitional flow (HTC=5000) and c) those

for fully turbulent flow (HTO=1 0000, 20000 and 50000). The curve corresponding

to laminar flow (HTC=1 000) is coincident with the first part of the transitional curve

(5000 L-T). At Reynolds number 2300 the friction law changes, imposing a

discontinuity in the neutral curve that follows the turbulent pattemn. This situation

can not arise in nature and is a consequence of the standard friction law transition

adopted. Figure A.6 shows a more detailed plot of this situation, considering the

additional case HTC=3000. It is important to point out that the friction crisis exists,

then, intermediate values can not be defined;* except by using a suitable

interpolation (as RELAP5 and every predictive code do). This problem does not

arise when the Reynolds number is used as the dependent variable. What is more

questionable is, perhaps, using a standard steady state law of friction for this

analysis. However, this is usually the standard practice.

30
CONCLUSIONS

The analysis of Reference [1 ] was extended to consider turbulent flow. This

allowed analyzing the behavior of a simple system and to specify a problem used

as a benchmark for a large-scale system thermal-hydraulic code. In this way, the

effects of numerical damping on unstable flows were quantified. While obtaining

the results shown, the influence of the laminar to turbulent flow transition became

evident. Then, considering smooth flow transition criteria and the effects of time

variation on the friction law may be of importance for the complete definition of the

problem.

31
6 6I i I I I I I I I

5 ýSTABLE REGION

4

3L

E

UNSTABLE REGION

1

0 I- I - I I I I I I I

0 50 100 150 200 250 300 350 400 450 500

(X

FIGURE A.1 The neutral stability curve for turbulent flow

32
8 25FADfl

CFAOO

-IFA)0

St~i.

FIGRE .2 ap stbilty ro obaind fom reerecemoal outo

of te flw ReereneA[5 eqatio, 

33
-o- from stability map

-o- calculated

*$ test problem

6

5

4

C 3

100 200 300 400 500

(X

FIGURE A.3 Comparison of the theoretical and the neutral stability limit

as obtained from Figure A.2

34
25

20

15

8Wel 10

I I I I I I

-o-<- Laminar flow

-oý- Turbulent flow

0 Test case

5 1-

OLI I I I I I I

0 100 200 300 400 500 600

cc

FIGURE A.4 Neutral curves for the fully laminar and fully turbulent flows

35
I I I I I

60 -HTC

1.03

0a-- 1.e4

40 - A- 2.e4

5.e4

EWel -- 5.03 L-T

20

0

2 4 6 8 10

FIGURE A.5 Neutral stability curves for laminar, transitional and turbulent flows

36
180

160

IoHTC=5.e31

140 I

H-HTC=3.e3

120 -Laminar Flow

100 -Turbulent Flow

80

S-Wel Not defined

60

40-A0

20

00

-20 -J B

2 4 6 8 10

qss

FIGURE A.6 Neutral stability curves for transitional flows

37
APPENDIX B

MEMORANDUM

Date: July 1991

SUBJECT: LIMITATIONS OF THE USE OF A "HEAT EXCHANGER"

.APPROXIMATION FOR A POINT HEAT SOURCE

Produced by: G.M. Grandi and J.C. Ferreri

Sector: Divisi6n Modelos Fisicos y Numd~ricos, Ger. Seg. Rad. y Nuclear, CNEA,

Av. del Libertador 8250, 1429 Buenos Aires, Argentina

The objective of the following analysis is showing why a one cell

approximation (loosely speaking a lumped parameter approximation) of a heat

exchanger can not be used to get the correct value of the steady state

temperature jump. This elementary condition holds for any representation using a

first order upwind approximation. A correction factor is derived to deal with this

situation. The working hypotheses of Welander's [1]i original formulation are

adopted in order to define a theoretical benchmark.

Let us consider the equation for the temperature approximation in the

tubes:

DT k T. q () T T)(B.1)

where, accordingly with Welander's [1] nomenclature, q is the volumetric flow rate,

T is the fluid temperature (assumed as uniform over each tube cross section),

To(s) is the temperature at the source wall, A is the tube cross section area, k is a

38
proportionality coefficient (assumed as constant) for the heat transfer, t is the time

and s is the space co-ordinate measured along the tube length.

For steady-state (ss) it is:

USS DT = k -(To (s) - T) (B.2)

where:

Uss = qss / A (B.3)

To is set as To ± AT in the source/sink. The rest of the tubes are adiabatic. Then,

dividing the previous equation by To and letting G -it is:

USS.~ k-(1-0) (BA4)

Expression (4) is valid in the source.

Let us consider a cell piecewise continuous variation of T along the source,

then, under an upwind approximation (Figure B.1), it is:

But: Ow =-0p, then:

Op= 1 1

1_+___s 1+2-V(B6

k-As

where V is the non-dimensional velocity.'

If the approximation for the temperature in the cell is defined as the mean

of ep and Gw, then:

The following table shows the values obtained with both approximations as a

function of acc/.

39
a/c V G (Welander) E)- 1 Op=1

___ _ ___ __ _ ___ __ ___ __1+2-V 2-V

0.1 0.1 0.999 0.833 5.0

1.0 0.648 '0.648 0.436 0.772

2.0 0.958 0.479 0.343 0.52

10.0 2.222 0.222 0.184 0.227

0.01 0.01 1.00 0.98 50.0

0.5 0.417 0.834 0.545 2.398

50.0 4.98 0.0096 0.0912 0.010

TABLE B.1: Values of Op computed using different approximations for the

source temperature.

As may be observed from Table B.1, expression (B.7) is a better

approximation for higher flow-rates. The inverse happens with expression (B.6).

Now, let's pose the following question:

Which is the expression to be used for the heat generated

at the source, under a FTUS approximation, to recover the

exact value (i.e. Welander's theoretical value) of Op?

The answer comes considering the following analysis. Referring once again

to Figure BA1 and integrating (B.2), it comes:

T T ASk

where As is the heated length, then:

T = To +(T - TO) -exp(-1 /V)

where V = Uss!/ (k.As).

The mean value of an equivalent heat source on As is:.

(k -(T -To)) = -s (T. - T) - [1 -exp('-1/ V)] (B.8)

Let us verify if this expression gives correct values for the fluid temperature (i.e.

coincident with the values given by Welander's analysis) at the source, when the

latter is computed via an upwind approximation. Using a FTUS discretization, the

previous expression leads to:

Tp-Tw =(TO -Tw).F .(B.9)

where:

40
F =1 -exp(-1 /V)

If we consider that Tp -Tw and introducing non-dimensional variables, then,

replacing in (B.8):

2F

This expression gives values coincident with the ones of Welander's analysis.

In conclusion, using a definition for AT that may be obtained from (13.8)

allows obtaining correct values for the steady state temperature jump at the

source. This applies to a discrete upwind approximation. This, in turn, brings

correct values for the steady state flow rate.

p

~TwJ -4- & ______________ ______________

Figure 1 A sketch of the node definitions for simple upwinding

41

APPENDIX C

RESULTS FROM THE APPLICATION OF RELAP5IMOD3.x CODE

0

S

LIST OF CASES CONSIDERED AND THEIR NODALIZATIONS

PLOTS OF SELECTED RESULTS

in Figures CA to C.10

Fluid Pressure

100010000

130010000

*

*

lower plenum

upper plenum

Fluid Temperatures

100010000

110010000

110050000

110100000

130010000

120100000

120050000

120010000

*

*

*

*

*

*

*

*

lower plenum or source

ascending leg 1

ascending leg 5

ascending leg 10

upper plenum or sink

descending leg 10

descending leg 5

descending leg 1

Mass Flow Rates

100010000

100020000

300000000

*

*

*

lower plenum to ascending leg

lower plenum to descending leg

upper plenum to TDV

43
Test No. 1 (code RELAP5IMOD3.O)

Cross section of components 100, 110, 120, 130 = 7.854x1 o-3 m2

Hydraulic diameter of components 100, 110, 120, 130 =0.1 m

Length of components 100 and 130 = 1 m

No. of nodes in the legs (110 and 120) = 10

Roughness parameter = 10-7 m

Pressure in components 400 = 105 Pa

Initial fluid temperature 293.15 K

Geometry of heat structures :cylindrical

Heat transfer area = 0.628 m2 (internal diameter = 0.2 m, length =1 m)

Thickness = 10~-4 m

Thermal conductivity = 100 W/(m2K)

Volumetric heat capacity = 10~5 j/(m3K)

Initial lower and upper structure temperature =293.15 k

Final upper structure surface temperature = 293.15 K

Final lower structure outer surface temperature = 303.15 K (after 0.1 s)

Imposed heat transfer coeff icient at lower and upper structures = 1000 W/(m2K)

Time step =0.5 s

Initial loop flow rate = 0.1 kg/s

Test No. 2 (code RELAP5IMOD3.0)

As Test No. 1 with:

Length of components 100 and 130 = 0.5 m

Structure heat transfer area = 0.314 m2 (internal diameter = 0.2 m,

length = 0.5 m)

Imposed heat transfer coefficient at lower and upper structures = 2000 W/(m2K)

44
Test No. 3 (code RELAP5IMOD3.O)

As Test No. 1 with:

Length of components 100 and 130 = 0. 1 m

Structure heat transfer area = 0.0628 m2 (internal diameter = 0.2 m,

length= =0.1 m)

Imposed heat transfer coefficient at lower and upper structures =10000 W/(m2K)

Test No. 4 (code RELAP5IMOD3.O)

As Test No. 3 with:

No. of nodes in the legs (110 and 120) = 25

Test No. 5 (code RELAP5/MOD3.O)

As Test No. 3 with:

No. of nodes in the legs (110 and 120) = 50

Test No. 6 (code RELAP5IMOD3.O)

As Test No. 3 with:

No. of nodes in the legs (110 and 120) = 90

Test No. 7 (code RELAP5IMOD3.O)

*As Test No. 6 with:

Initial loop flow rate = 0.001 kg/s

Test No. 8 (code RELAP5IMOD3.O)

As Test No. 7 with:

No. of nodes in the legs (110 and 120) = 40 (searching for neutral stability)

Test No. 9 (code RELAP5/MOD3.O)

As Test No. 7 with:

No. of nodes in the legs (110 and 120) = 38 (searching for neutral stability)

Test No. 10 (code RELAP5/MOD3.O)

As Test No. 7 with:

No. of nodes in the legs (110 and 120) = 36 (searching for neutral stability)

Test No. 11 (code RELAP5/MOD3.O)

As Test No. 7 with:

No. of nodes in the legs (110 and 120) = 34 (searching for neutral stability)

Test No. 12 (code RELAP5IMOD3.1)

As Test No. 8

Test No. 13 (code RELAP5IMOD3.1)

As Test No. 6

45
Test No. 14 (code RELAP5IMOD3.O)

As Test No. 7 with:

No. of nodes in the legs (110 and 120) =37 (searching for neutral stability)

Test No. 15 (code RELAP5/MOD3.O)

As Test No. 14 with:

Time step = 1.0Os

Test No. 16 (code RELAP5IMOD3.O)

As Test No. 14 with:

Time step = 0.25 s

Test No. 19 (code RELAP5IMOD3.O)

As Test No. 9 with:

Pressure control junction 300 closed after 3000 s

ADDITIONAL CASE SERIES

Test No. X1 (code RELAP5IMOD3.O)

As Test No. 1 with:

Structure heat transfer area = 0.0314 m2 (internal diameter =0.1 m,

length = 0.1 m)

Test No. X2 (code RELAP5IMOD3.O)

As Test No. 2 with:

Structure heat transfer area = 0.0314 m2 (internal diameter =0.1 m,

length = 0.1 m)

Test No. X6 (code RELAP5IMOD3.O)

As Test No. 6 with:

Structure heat transfer area = 0.0314 m2 (internal diameter =0.1 m,

length = 0.1 m)

Test No. X14 (code RELAP5/MOD3.O)

As Test No. 14 with:

Structure heat transfer area = 0.0314 m2 (internal diameter = 0.1 m,

length = 0.1 m)

46
Tests performed using RELAP5IMOD2 for the

Benchmark Calculation

Nodalization

400

115

113

110

Test No. F2

Cross section of loop components = 1 .26x1 0 -3m2

Hydraulic diameter of loop components = 0.04 m

Overall length of components 100-102-104 and 120-122-124 =0.5 m

Length of components 100, 104, 120, 124 = 0. 1 m

Length of nodes in pipes 102 and 122 = 0.08 m, 0.05 m 0.04 m, 0.05 m, 0.08 m

Length of pipes 110, 115, 130, 135 =5mr

No. of nodes in components 110, 115, 130, 135 =50

R oughness parameter = 10 m r

P ressure in component 400 = 10 5Pa

Initial fluid temperature 293.15 K

Geometry of heat structures: cylindrical

Length of structures = 0.04 m

Thickness = 5x1 0 5m

Thermal conductivity = 100 W/(m 2K)

Volumetric heat capacity = 105 j/(rn K)

47
Initial lower and upper structure temperature = 293.15 K

Final upper structure surface temperature = 293.15 K

Final lower structure outer surface temperature = 303.15 K (after 100 s) 2

Imposed heat transfer coefficient at lower and upper structures = 2500 W/(m K)

Time step = 0.5 s

Initial loop flow rate = 0.1 kg/s

Steady-State Period = 10 s

Test No. F3

Cross section of loop components = 3.1416x1 0 -4m2

Hydraulic diameter of loop components = 0.02 m

Overall length of components 100-102-104 and 120-122-124 = 0.5 m

Length of components 100, 104, 120, 124 = 0.1 m

Length of nodes in pipes 102 and 122 = 0.08 m, 0.05 m 0.04 m, 0.05 m, 0.08 m

Length of pipes 110, 115, 130, 135 =5 m

No. of nodes in components 110, 115, 130, 135 =50

Roughness parameter = 10 7m5

Pressure in component 400 = 10~ Pa

Initial fluid temperature 293.15 K

Geometry of heat structures :cylindrical

Length of structures = 0.02 m

Thickness = 5x10 5m

Thermal conductivity = 100 W/(m 2K)

Volumetric heat capacity = 1 05 J/(m 3K)

Initial lower and upper structure temperature =293.15 K

Final upper structure surface temperature = 293.15 K

Final lower structure outer surface temperature = 303.15 K (after 100 s)

Imposed heat transfer coefficient at lower and upper structures = 2500 W/(m 2K)

Time step = 0.5 s

Initial loop flow rate = 0.2 kg/s

Steady-State Period = 10 s

Test No. F4

Cross section of loop components = 3.1416x1 0 m2

Hydraulic diameter of loop components = 0.02 m

Overall length of components 100-102-104 and 120-122-124 = 0.5 m

Length of components 100, 104, 120, 124 = 0.1 m

Length of nodes in pipes 102 and 122 = 0.08 m, 0.05 m 0.04 m, 0.05 m, 0.08 m

Length of pipes 110, 115, 130, 135 =5 m

No. of nodes in components 110, 115, 130, 135 = 50

48
-7 Roughness parameter = 10 mn

Pressure in component 400 = 10 Pa

Initial fluid temperature 293.15 K

Geometry of heat structures :cylindrical

Length of structures = 0.04 mn

Thickness = 5x1 0 5m

2 Thermal conductivity = 100 W/(m K)3 Volumetric heat capacity = 105 J/(m K)

-Initial lower and upper structure temperature =293.15 K

Final upper structure surface temperature = 293.15 K

Final lower structure outer surface temperature = 303.15 K (after 100 s) 2

Imposed heat transfer coefficient at lower and upper structures = 2500 W/(m K)

Time step = 0.5 s

Initial loop flow rate = 0.5 kg/s

Steady-State Period = 10 s

Test No. F5

As Test No. F2 with:

Only 10 nodes in the legs (in place of 100 nodes as in test F2)

49
WinGmfl 2.02.24-1997

1.50

1.25

1.00

.75

.50

.25

0

i

99

0

99

99

U

Test #1

I I I IXXX FERI MFLOWJ100010000

YYY FERI MFLOWJ100020000-

ZZZ FERi MFLOWJ300000000

-Y

-.25

-.50

-.75

.1.00 0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Figure C.la TEST #1 Base Case using RELAP5/MOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1im, Heated Length =i1m,

Initial Mass Flow Rate=O.1 Kg/s)

YAkGrt- 2. .022-W927

99

E

U.

27.00

26.00

25.00

24.00

23.00

22.00

21.00

20.00

19.00

0 1000. 2000. 3000. 4000. 5000. 600. 7000. 8000. 9000.

Time (s)

Figure C.1lb TEST #1 Base Case

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure CAla.

50
W14r.t1 .q2-24-1997

1.50

1.25

1.00

.75

U

0

IL

U

U

.50

.25

0

-.25

-.s0

-.75

-1.00

0 1000. 2000. 3000. 4000. 5000.

Time (s)

6000. 7000. 8000. 9000.

Figure C.2a TEST #2 Sensitivity of results to heated length and RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1im, Heated Length = 0.5m,

Initial Mass Flow Rate=0.1 Kg/s)

Wh4.GluLV .02-24-1997

CL

E

ILA

27.00

26.00

25.00

24.00

23.00

22-00

21.00

20.00

19.00

0 1000. 2000.- 3000. 4000. 5000. 6000. 7000. 8000. 900.

Time (s)

Figure C.2b TEST #2 Sensitivity of results to heated length

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.2a.

51
i

a,

0

a,

a,

a,

1.50

1.25

1.00

.75

.50

.25

0

-.25

-.50

-.75

W.G,.f2.1 .02.24-la,?

TEST #3

XXX FER3 MFLOWJlOOOlOOOO

YYY FER3 MFLOWJ100020000

ZZZ FER3 MFLOWJ300000000

Ix x

-1 .00

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Figure C.3a

28.00

27.00

26.00

TEST #3 Sensitivity of results to heated length and RELAP5/MOD3.O

Mass Flow Rates in Legs and TDV junction, 10 nodes per leg,

(Length of Components = 1im, Heated Length = 0.1m,

Initial Mass Flow Rate=O. 1 Kg/s)

9

E

25.00

24.00

23.00

22.00

21.00 - 20.001.)

19.00

0

Figure C.3b

1000. 2000. 3000. 4000. 5000.

Time (S)

6000. 7000. 8000. 9000.

TEST #3 Sensitivity of results to heated length

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.3a.

52
1.50

1.25

1.00

.75

.SO

.25

0

0

LL

WAGml2.1.02-24-1597

Test #4

XXX FER4 MFLOWJIOOO10000

YYY FER4 MFLOWJlD0020000

ZZZ FER4 MFL.OWJ300000000

I AIUVM I 

-.25

.%50

-.75

-1.00

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 900.

Time (s)

Figure C.4a TEST #4 Sensitivity of results to nodes number and RELAP5/MOD3.O

Mass Flow Rates in Legs and TDV junction, 25 nodes per leg,

(Length of Components = 0.4m, Heated Length = 0.1lm,

Initial Mass Flow Rate=O.1 Kg/s)

Wimf 2.1. 02-24-IS?

5.

U

S

a

S

0.

E

0

.5

28.00

27.00

26.00

25.00

24.00

23.00

22.00

21.00

20.00

19.00

0 1000. 2000. 3000. 4000. 5000.

Time (s)

6000. 7000. 8000. 9000.

Figure C.4b TEST #4 Sensitivity of results to nodes number

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.4a.

53
WI.GMM .1 02.24-IN?

2.00

1.50

1.00

M5

.2

W

TEST #5

XX(X FERS MFLOWJ100010000

'YYY FER5 MFLOWJ100020000 -

AfMAA~~A YWWxM ZZZ FER5 MFLOYVJ300000000 w_

f~f\A ~mr\RwvvvPA

1h SA V

..50

-1.00

.1.50

0 1000. 2000. 3000. 4M0. 5000. 6000 7000. 8000. 9M0.

Trime (s)

TEST #5 Sensitivity of results to nodes number and RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 50 nodes per leg,

(Length of Components = 0.2m, Heated Length = 0. 1 m,

Initial Mass Flow Rate=0.1 Kg/s)

Figure C.5a

34.0

32.0

30.0

WrCGUZI .02-24-1997

ca

E

28.0

26.0

24.0

22.0

20.0

18.0

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

rime (s)

Figure C.5b TEST #5 Sensitivity of results to nodes number

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.5a.

54
WMW~t2.1 .02.2-1997

2-00

1.50

1.00

.50

*0

U

0

01

U

U

I I I I I I I a

TEST #6

XXX FER6 MFLOWJ100010000

YYY FER6 MFLOWJ100020000-

ZZZ FER6 MFLOWJ3000000 AV\kWM/!pNN. AAA1f\AM 1A .v

f~/N\ ~ ffYrVL Y~\ I I I II I I

-.50

-1.00

-1.50

Figure C.6a

32.0

30.0

0 1000. 200. 3000. 400. 5000. 600. 7000. 8000. 900.

Time (s)

TEST #6 The recovery of unstable behavior and RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 90 nodes per leg,

(Length of Components =0.l11m, Heated Length = 0.1lm,

Initial Mass Flow Rate=0. 1 Kg/s)

VMGmfZ1 .02-2&-1997

U

S

U

C.

2

U

28.0

26.0

24.0

22.0

20.0

18.0

Figure C.6b

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

TEST #6 The recovery of unstable behavior

Fluid Temperatures at Source, Ascending Leg Cell 5,

Sink and Descending Leg Cell 5, conditions as in Figure C.6a.

55
W,,.Wt2.2.02-2S.1997

34.0

32.0

30.0

2.

S 26.0

E

*~24.0

U.

22.0

20.0

18.0-

-1.25

Figure C.6c

28.00

298.00 Fi

27.00

-1.00 -.75 -.50 -.25 0 M2 ' .50 JS5 1.00 1.25

Mass Flows (Kgfs)

TEST #6 The recovery of unstable behavior

Source Fluid Temperature vs. Mass Flow Rate,

Conditions as in Figure C.6a.

W.C.&t = 2.02-2S-1997

2.

S

S

a

S

S

a.

ES

I.-

b

26.00

25.00

24.00

23.00

22.00

21.00 L

-1.25 -1.00 -.75 -.50 -.25 0 .25

Mass Flows (K91s)

.50 .75 1.00 1.25

Figure C.6d TEST #6 The recovery of unstable behavior

Sink Fluid Temperature vs. Mass Flow Rate,

Conditions as in Figure C.6a.

56
Wr4mf 2.1 - 02-24-1997

2.00

1.50

1.00

-.50

-10

-1.50

Figure C.a

0 1000. 2000. 300. 400. 500. 600. 7000. 8000. 900.

Time (s)

TEST #7 The recovery of unstable behavior starting from (almost) rest

and RELAP5IMOD3.O

Mass Flow Rates in Legs and TDV junction, 90 nodes per leg,

(Length of Components =0.1 1 m, Heated Length = 0.1 m,

Initial Mass Flow Rate = 0.001 Kg/s)

WbrG.12.1-02-24-I99

2.00

1.50

1.00

I

0

C

S

S

U

.50

-.50

-1.00

-1.50

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Figure C.7b TEST #7 The recovery of unstable behavior starting from

Mass Flow Rate in Ascending Leg, 90 nodes per leg,

(Length of Components =0.l11m, Heated Length = 0.1lm,

Initial Mass Flow Rate = 0.001 Kg/s)

(almost) rest

57
W.GMa2II-02.24-1997

34.0

32.0

MX FER7 TEMPF1 00010000

YYY FER7 TEMPF1 30010000

30.0 Z FER7 TEMPF11I00SOOOO

a VW FER7 TEM F12 005O~

S 28.0

S 26.0

:2 240/

22.0 II

20.0

18.0 I

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Figure C.7c TEST #7 The recovery of unstable behavior starting from (almost) rest

Fluid Temperatures at Source, Ascending Leg Cell 5, Sink and

Descending Leg Cell 5, conditions as in Figure C.7a.

W14MMI 2.102.241991

34.0

Test #7

32.0

XXCX FER7 TEMPFIOOOIOOOO

30.0 -YYY FER7 TEMPF130010000

30.0

260

0. 

.2

26.0 Y IV PP

22.0V

20.0

18.0 I

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

Figure C.7d TEST #7 The recovery of unstable behavior starting from (almost) rest.

Fluid Temperatures at Source and Sink, conditions as in Figure C.7a.

58
VWCM? L2 - 2-25-1997

34.0

32.0

XXX FER7 TEMPFlDO010000

30.0

* 28.0

* 26.0

E

:2 24.0

22.0

20.0

18.0 ________________________

-1.25 -1.00 a.75 -.50 -.25 0 .25 .50 .75 1.00 1.25

Mass Flows (Kg/s)

Figure C.7e TEST #7 The recovery of unstable behavior starting from (almost) rest

Source Fluid Temperature vs; Mass Flow Rate,

Conditions as in Figure C.7a.

YA&4f .t2 -02-2S.1907

E

A

32

28.00 1 1 1 1- -

Test #7

27.00

26.00

25.00

24.00

23.00

22.00

21.00

20.00 -

-1.25 -1.00 -.75 -.50 -.25 0 .25

Mass Flows (KqIs)

.50 .75 1.00 1.25

Figure C.7f TEST #7 The recovery of unstable behavior starting from (almost) rest

Sink Fluid Temperature vs. Mass Flow Rate, conditions as in Figure C.7a.

59
Y&4*f 2.2 -02-26,1997

1.20

1.00

.80

.80

.40

a

.2

Ua

a

a

.20

0

-.20

Figure C.8.a

0 1000. 200. 3000. 4M0. 5000. 600. 7000. 8000. 900.

Time (s)

TEST #14 Search for Neutral Stability and RELAP51MOD3.O

Conditions as in TEST #7

Mass Flow Rates in Ascending Leg and toward TDV, 37 nodes per leg,

(Length of Components =0.1l1m, Heated Length = 0.1lm,

Initial Mass Flow Rate = 0.001 Kg/s)

a

a

a

0

0.

2

0

I-

*0

.5

U32.0

30.0

28.0

26.0

24.0

22.0

20.0

18.0

WICWuf22 .02-24W7

S Test #14

)O(X FE14 TEMPF100010000

YYY FE14 TEMPF130010000

]x

ft

-I

0 1000. 2M0. 3000. 400. 5000. 60M. 7000. 8000. 900.

rime (s)

TEST #14 Search for Neutral Stability

Fluid Temperatures at Source and Sink conditions as in Figure C.8a.

Figure C.8b

60
VA,~G~t 2.2.- 02-2S-1"7

34.0 ---- T

Test

32.0

30.0

* 28.0

S26.0

E

2 24.0

22.0

20.0

0 .10

Figure C.Bc

28.00

27.00 - Test A

26.00

S 25.00

24.00

E 23.00

21.00

20.00

19.00 L...

.20 .30 .4O M5 .60 .70 .80 .90 1.00

Mass Flows (Kg1s)

TEST #14 Search for Neutral Stability

Source Fluid Temperature vs. Mass Flow Rate,

Conditions as in Figure C.8a.

W,GCsf 2.2 .02-2S-lgg7

0 .10 .20 .30 .40 M5 .60 .70 .80 .90 1.00

Mass Flows (K91s)

Figure C.8d TEST #14 Search for Neutral Stability.

Sink Fluid Temperature vs. Mass Flow Rate, conditions as in Figure C.8a.

61
WV&,GafZ2.2.O2S.1997

2.00

1.50

1.00

£50

0

01

a

2

-.£0

-1.00

Figure C.9a

0 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000. 9000.

Time (s)

TEST #19 Sensitivity of results of TDV isolation and RELAP5IMOD3.O

Mass Flow Rates in Ascending Leg, 37 nodes per leg,

(Length of Components = 0.263m, Heated Length = 0. 1 m,

Initial Mass Flow Rate=O.001 Kg/s)

X 10 -5m'?~ .22¶9

as I

a

a

.2

LI.

a

a

a

30

25

20

i5

10

5

.0

0 1000. 2000. 3000. 4000. 5000.

Time (s)

6000. 7000. 8000. 9000.

Figure C.9b TEST #19 Sensitivity of results of TDV isolation at 3000s

Flow rates to/from TDV, Test #9 With valve left open up to 9000s

Conditions as in Figure C9.a

62
2.00

1.50

1.00

.50

0

U

0

U

0,

U

W.,G.t±2.202.2S-10q7

Comparison Test #9 & Test #19

XXX FE19 MFLOWJ100010000

'YYY FER9 MFLOWJ100010000

A I

..So

-1.00

0. 1000. 2000. 3000. 4000. 5000. 6000. 7000. 8000.

Time (s)

9000.

Figure C.9c TEST #19 Sensitivity of results of TDV isolation at 3000s

Flow rates in Ascending Leg, Comparison of Tests #9 and #19

Conditions as in Figure C9.a

VAGmutZ.2.02~-?99

U

0

U

U

U

1AO

1.20

1.00

.80

.60

.40

.20

0

-.20

0 1000. 2000. 3000. 4000. 5000.

Time (s)

6000. 7000. 8000. 9000.

Figure C.10 TEST #14 Search for Neutral Stability, conditions as in TEST #7

Mass Flow Rates in Ascending Leg using RELAP5/MOD3.O and MOD2,

37 nodes per leg, Length of Components =0.1llm, Heated Length = 0.1lm,

Initial Mass Flow Rate = 0.00 1 Kg/s

63
APPENDIX D

INPUT DECK FOR WELANDER'S PROBLEM

64
*Natural Circulation Input deck for Welander's Problem

* Case: alpha = 338.94-epsilon =2.3 90 Nodes per Leg

100 new transnt

*time steps

*mint maxt mie mae

201 200. 1.e-7 5.e-2 07003 500 5000 5000

202 1000. 1.e-7 5.6-2 07003 500 5000 5000

203 l.e6 1.e-7 S.e-2 07003 500 5000 5000

* minor edit requests

*pressure

;0 0000 *lwrpeu

301 p 100010000 * lower plenum

: fluid temperatures

303 tempf 100010000 * lower plenum

304 tempf 110010000 * ascending leg 1

305 tempf 110050000 * ascending leg 5

306 tempf 110100000 * ascending leg 10

310 tempf 130010000 * upper plenum

307 tempf 120100000 * descending leg 10

308 tempf 120050000 * descending leg 5

309 tempf 120010000 * descending leg 1

*mass flow rates

31 foj 10100 *lwe lnmt sedn o

311 mflowi 100010000 * lower plenum to aescending leg

----------------------------------------------------

* trips ----------------------------------------------------

*pressure control

527 time 0 ge null 0 -1. 1 *

528 time 0 ge null 0 1.e06 I1*

603 527 xcr 528 n *

*heater temperature table

513 time 0 ge null 0 -1.e06 1 * start of heat, and cool.

: cooler temperature table

514 time 0 go null 0 -l.e06 1 * start of beat, and cool.

*heater htc table

515 time 0 ge null 0 -1.e06 1 * constant htc

*cooler htc table

516 time 0 ge null 0 -l.e06 1 * constant htc

*end trip signal

512 time 0 ge null 0 10000;. 1 * end of problem

600 512 * end of problem

65
*_ hydraulic components -- - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - ---

*lower plenum

1000000 low.ple branch

1000001 2 1

1000101 7.854e-03 0.1000 0. 0.

1000200 000 1.00e5 0.10465.6

1001101 100010000 110000000 0.

1002101 100000000 120000000 0.

1001201 0.99 0. 0.

1002201 -0.99 0. 0.

0. 0.0000 1.e-7 0.100 00000

2.675e6 0.

0. 0. 000000

0. 0. 000000

* ascending leg

1100000 asc.leg pipe

1100001 90

1100101 7.854e-03 90

1100301 0.111111111111 90

1100401 0. 90

1100601 90.0 90

1100801 1.Oe-7 0.0 90

1101001 00000 90

1101101 000000 89

1101201 000 1.0005 0.106743e6 2.675e6 0.0

1101300 1

1101301 0.99 0.000 0. 89

: descending leg

1200000 desc.leg pipe

1200001 90

1200101 7.854e-03 90

1200301 0.111111111111 90

1200401 0. 90

1200601 90.0 90

1200801 1.0.-7 0.0 90

1201001 00000 90

1201101 000000 89

1201201 000 1.00.5 0.102557e6 2.675e6 0.0

1201300 1

1201301 -0.99 0. 0. 89

0. 90

0. 90

*upper plenum

1300000 upp.ple branch

.1300001 2 1

1300101 7.854e-03 0.1000 0. 0.

1300200 000 1.O0eS 0.10465.6

1301101 110010000 130000000 0.

1302101 120010000 130010000 0.

1301201 0.99 0. 0.

1302201 -0.99 0. 0.

0. 0.0000 1.e-7 0.100 00000

2.675e6 0.

0. 0. 000000

0. 0. 000000

* pressure control for steady state

3000000 pc.vlv valve

3000101 130010000 400000000 7.85de-03 0.1 0.1 000100 1. 1. 1.

3000201 1 0. 0. 0.

3000300 mtrvlv

3000301 603 528 0.25 1.

* pressure control volume

4000000 pc.vol tmdpvol

4000101 0. 0.5 10. 0. 0. 0. 4.e-5 0. 00000

4000200 000

4000201 0. 1.0e5 0.084e6 2.675e6 0.

4000202 1. 1.0eS 0.084e6 2.675e6 0.

4000203 1.e6 1.0.5 0.084e6 2.675e6 0.

66
* tructures ------------------------

*lover plenum beaters and upper plenum coolers

11001000 2 11 2 1 0.05

11001100 0 1

11001101 10 0.0501

11001201 1 10

11001301 0. 10

11001400 0

11001401 293.15 11

11001501 -900 0 1900 1 0.10

11001502 -901 0 1901 1 0.10 2

11001601 100010000 00000 3902 1 0.10

11001602 130010000 00000 3903 1 0.10 2

11001701 0 0.00 0.0 0.0 1

11001702 0 0.00 0.0 0.0 2

11001901 0. 100. 100. 0. 0. 0. 0. 1. 2

----------------------------------------------------

*material tables ----------------------------------------------------

20100100 tbl/fctn 1 1 ----------------------------------------------------

b eat structure conductivity ----------------------------------------------------

20100101 1.0 100.

20100102 5000.0 100. ----------------------------------------------------

b eat structure volumetric heat capacity ----------------------------------------------------

20100151 1.0 1.e5

20100152 5000.0 1.eS

----------------------------------------------------

*general tables ----------------------------------------------------

*heater temperature

20290000 tem 513

20290001 -1.0 293.15

20290002 0. 293.15

20290003 0.1 303.15

20290004 1.e6 303.15

*cooler temerature

20290100 te= 514 20290101 -1.0 293.15

20290102 0. 293.15 20290103 100.1 293.15

20290104 1.e6 293.15

: heater htc

20290200 htc-t 515

20290201 -1.0 20000.

20290202 0. 20000.

20290203 0.1 20000.

20290204 10.0 20000.

20290205 1.e6 20000.

67
*cooler htc

20290300 khtc-t 516

20290301 -1.0 20000.

20290302 0. 20000.

20290303 0.1 20000.

20290304 10.0 20000.

20290305 l.e6 20000.

----------------------------------------------------

*control variables ----------------------------------------------------

*none

68
NRC FORM 336 U.S. NUCLEAR REGULATORY COMMISSION 1. REPORT NUMBIER

(2q(Ansignd by NRC. Add Vol.. Supp.. Rev.. NRCR 110UZRPHCDAASHE and Addendum Numbems Ifany.) 320130 311113 HCDT HE

2. TITLE AND SUBTTE(eekh~mn nM rvm NUREG/IA-Ol 51

Verification of RELAP/MOD3 With Theoretical and Numerical Stability3. DTREOTUBSD

Results on single-Phase, Natural Circulation in a Simple Loop 3 AERPR ULSE

February 1999

4. FIN OR GRANT NUMABER

5. AUTHOR(S) 6. ~TYPE OF REPORT

J.C. Ferreri, ARN Technical

W. Ambrosini, USP____ ________

7. PERIOD COVERED llncusive Dows)

8. PERFORMING ORGANIZATION -NAME AND ADDRESS (If ARQ provide D&Ws, Oftie or Region, U.S. Mifw Reguihy yCornmisuionL and meiQV~ eddass:* if oontwkr

prwide nerne and mailfing address.)

Autoridad Regulatoria Nuclear Universita degli studi Pisa, Facolta dit ingegneria

Av. del Ubertador 8250 Dipartimento di Costruzioni Meccaniche e Nucleari

1429 Buenos Aires, Argentina Via DiotisaMv 2, 56126 Pisa, Italy

9. SPONSORING ORGANIZATION - NAME AND ADDRESS (ITPK be 'Sameas above' ,f conlrcbca provide NfRC Di v~ikn cfoe orRegion, U.S. NuclawRegulay Coimmsskiv

and maftn sdvhsss.

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

10. SUPPLEMENTARY NOTES

11. ABSTRACT.p words aliess)

The theoretical results given by Pierre Welander are used to test the capability of the RELAP5 series of codes to predict Instabilities

in single-phase flow. These results are related to the natural circulation in a loop formed by two parallel adiabatic tubes with a point

heat sink at the top and a point heat source at the bottom. A stability curve may be defined for laminar flow and was extended to

consider turbulent flow. By a suitable selection of the ratio of the total buoyancy force in the loop to the friction resistance, the flow

may show instabilities. The solution was useful to test two basic numerical properties of the RELAP5 code, -namely: a)

convergence to steady state flow-rate using a "lumped parameter' approximation to both the heat source and sink and, b) the

effect of nodalization to numerically damp the instabilities. It was shown that, using a simple volume to lump the heat source and

sink, it was not possible to reach convergence to steady state flow rate when the heated (cooled) length was diminished and the

heat transfer coefficient increased to keep constant the total heat transferred to (and removed from) the fluid. An algebraic justification of these results is presented, showing that it is a [imitation inherent to the numerical scheme adopted. Concerning the

effect of nodalization on the dam ping of instabilities, it was shown that a "reasonably fine" discretization led, as expected, to the

damping of the solution. However, the search for convergence of numerical and theoretical results was successful, showing the

expected nearly chaotic behavior. This search lead to very refined nodalizations. The results obtained have also been verified by

the use of simple, ad hoc codes. A procedure to assess the effects of nodalizations on the prediction of instabilities threshold is

outlined in this report. It is based on the experience gained with the aforementioned simpler codes.

12. KEY WORDS/DESCRIPTORS (Lis word a ptvases Um.t willassist resee'ches in, tocef Umg thepodv) 13. AVAILABUIMYSTATIEMEN

RELAP5IMOD3 unlimited

14. SECURi1Y C.AsSSFICATION

unclassified

unclaasified

15. NUMVBER OF PAGES

16. PRICE

NRC FORM 335 (2439)


Federal Recycling Program

NUREG/IA- 0151 'VERIFICATION OF RELAP5/MOD3 WITH THEORETICAL AND NUMERICAL STABILITY

RESULTS ON SINGLE-PHASE, NATURAL CIRCULATION IN A SIMPLE LOOP

FEBRUARY 1999

UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, DC 20555-0001

SPECIAL STANDARD MAIL

POSTAGE AND FEES PAID

USNRC

PERMIT NO. G-67

OFFICIAL BUSINESS

PENALTY FOR PRIVATE USE, $30