NUREG/JA-0073

-. ABEEW-M2590

International

Agreement: Report

Time Step and Mesh Size

Dependencies in the Heat

Conduction Solution of a

Semi-Implicit, Finite Difference

Scheme for Transient Two-Phase Flow

Prepared by

R. O'Mahoney

Winfrith Technology Centre

United Kingdom Atomic Energy Authority

Dorchester, Dorset, DT2 8DH

United Kingdom

Office of Nuclear Regulatory Research

" U.S. Nuclear Regulatory Commission

Washington, DC 20555

April 1992

Prepared as part of

The Agreement on Research Participation and Technical Exchange

under the International Thermal-Hydraulic Code Assessment

and Application Program ([CAP)

Published by

U.S. Nuclear Regulatory Commission
NOTICE

This report was prepared under an international cooperative

agreement for the exchange of technical information.ý Neither

the United States Government nor any agency thereof, or 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.

Available from

Superintendent of Documents

U.S. Government Printing Office

P.O. Box 37082

Washington, D.C. 20013-7082

and

National Technical Information Service

Springfield, VA 22161
NUREG/IA-0073

AEEW-M2590

International

Agreement Report

Time Step and Mesh Size

Dependencies in the Heat

Conduction Solution of a

Semi-Implicit, Finite Difference

Scheme for Transient Two-Phase Flow

Prepared by

R. O'Mahoney

Winfrith Technology Centre

United Kingdom Atomic Energy Authority

Dorchester, Dorset, DT2 8DH

United Kingdom

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555

April 1992

Prepared as part of

Te Ageemtent on Resac Partciaton and Te ,chca Exchange

une th ntrational Thra-yrulicCoeAss en

and Application Program (ICAP)

Published by

U.S. Nuclear Regulatory Commission
NOTICE

This report is based on work performed under the sponsorship of the

United Kingdom Atomic Energy Authority. The information in this

report has been provided to the USNRC under the terms of the

International Code Assessment and Application Program (ICAP)

between the United States and the United Kingdom (Administrative

Agreement - WH 36047 between the United States Nuclear Regulatory

Commission and the United Kingdom Atomic Energy Authority Relating

to Collaboration in the Field of Modelling of Loss of Coolant

Accidents, February 1985). The United Kingdom has consented to the

publication of this report as a USNRC document in order to allow

the widest possible circulation among the reactor safety community.

Neither the United States Government nor the United Kingdom or any

agency thereof, or any of their employees, makes any warranty,

expressed or implied, or assumes any legal liability of

responsibility for any third party's use, or the results of such

use, or 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.
AEEW - M 2590

TIME STEP AND MESH SIZE DEPENDENCIES IN THE HEAT CONDUCTION

SOLUTION OF A SEMI-IMPLICIT, FINITE DIFFERENCE SCHEME FOR

TRANSIENT TWO-PHASE FLOW

R O'Mahoney

Summary

This report examines, and establishes the causes of, previously

identified time step and mesh size dependencies. These

dependencies were observed in the solution of a coupled system of

heat conduction and fluid flow equations as used in the

TRAC-PFl/MODl computer code.

The TRAC-PF1/MODl computer code employs a semi-implicit, finite

difference solution scheme to solve the differential equations

describing heat transfer and two-phase fluid flow; it is commonly

used to analyse loss-of-coolant accidents in Pressurised Water

Reactors.

The report shows that a significant time step size dependency can

arise in calculations of the quenching of a previously unwetted

surface. The cause of this dependency is shown to be the

explicit evaluation, and subsequent smoothing, of the term which

couples the heat transfer and fluid flow equations. An axial

mesh size dependency is also identified, but this is very much

smaller than the time step size dependency.

The report concludes that the time step size dependency

represents a potential limitation on the use of large time step

sizes for the types of calculation discussed. This limitation

affects the present TRAC-PFl/MODl computer code and may similarly

affect other semi-implicit finite difference codes that employ

similar techniques. It is likely to be of greatest significance

in codes where multi-step techniques are used to allow the use of

large time steps.

Safety and Engineering Science Division

Winfrith Technology Centre

July 1989

AEEW - M. 2590 i

CONTENTS

PAGE

SECTION 1 INTRODUCTION 1

2 DESCRIPTION OF TRAC-PFI/MODl 1

3 TRAC-PFI/MOD1 QUENCHI NG RESULTS 2

4 DETAILED EXAMINATION OF CONDUCTION TERMS 4

4.1 Finite Difference Equation 4

4.2 Surface-to Fluid Effects 4

4.2.1 Explicit Evaluation and 5

Smoothing Effects

4.2.2 Time Step Size Effects 6

4.3 Axial Conduction Effects 6

4.4 Quench Front Profiles 7

4.4.1 Calculation with 0.25 mm 7

Minimum Axial Mesh

4.4.2 Calculation with 0.25 mm 8

Mesh and 0.3 ms Time Step

4.4.3 AXIAL Term Time Step8

Dependency

4.4.4 Calculations with no Axial 9

Conduction

4.4.5 Axial Mesh Size Effects 10

5 SUMMARY AND CONCLUSIONS 10

6 REFERENCES 12

AEEW - M 2590V

FIGURES

PAGE

1 Rod surface temperatures at 5 elevations, for 4 different 13

mesh sizes. Min axial meshs are: CONT=2.5 mm,

SHORT=0.25 mm, LONG=0.1 mm, THICK=0.05 mm.

2 Axial profile of rod surface temperature at 4 different 14

times. Profiles at 0, 10, 20 and 30 seconds, for 0.1 mm

axial mesh.

3 Rod surface temperatures at 5 elevations, for 3 15

calculations. Calculations are: CONT=0.25 mm,

SHORT=0.05 mm, LONG=O.25 mm + 0.3 ms step.

4 Surface-to-fluid heat transfer coefficient vs temperature, 16

at 20 seconds. TRAC calculation with 0.25 mm min mesh +

theoretical heat transfer.

5 Surface-to-fluid heat flux vs temperature, at 20 seconds. 17

TRAC calculation with 0.25 mm min mesh + theoretical heat

flux.

6 Surface-to-fluid heat transfer coefficent vs temperature, 18

at 20 seconds. TRAC calculation with 0.25 mm min mesh, no

smoothing + theoretical heat transfer.

7 Surface-to-fluid heat tranfer coefficent vs temperature, at 19

20 seconds. TRAC calculation with 0.25 mm min mesh, 0.3 ms

time step + theoretical heat transfer.

8 Surface-to-fluid heat flux vs temperature, at 20 seconds. 20

TRAC calculation with 0.25 mm min mesh, 0.3 ms time step +

theoretical heat transfer.

9 Rod surface temperatures at 7 elevations, for 3 no-axial 21

calculations. Calculations are: CONTO0.25 mm,

SHORTO0.25 mm. + 0.3 ins, LONG=0.05 mm + 0.3 ins.

10 Surface-to-fluid heat flux vs temperature, at 20 seconds. 22

TRAC calculation, no axial conduction, 0.25 mm minioum mesh

+ theoretical heat flux.

11 Surface-to-fluid heat flux vs temperature, at 20 seconds. 23

TRAC calculation, no axial conduction, 0.25 mm min mesh,

0.3 ms time step + theoretical heat flux.

12 Heat conduction equation: quench front profile at 24

20 seconds. TRAC calculation with 0.25 mm, min mesh.

13 Heat conduction equation: quench front profile at 25

20 seconds. TRAC calculation with 0.25 mm min mesh

(exploded view).

AEEW - M 2590 i Vii
FIGURES (Continued)

PAGE

14 Heat conduction equation: quench front profile at 26

20 seconds. TRAC calculation with 0.25 mm min mesh, 0.3 ms

time step.

15 Heat conduction equation: quench front profile at 27

20 seconds. TRAC calculation with 0.25 mm min mesh, 0.3 ms

time step (exploded view).

16 Surface-to-fluid heat flux vs temperature, at 20 seconds. 28

TRAC calculation with 0.25 mm min mesh, 0.3 ms step,

reduced CHF + theoretical heat flux.

17 Heat conduction equation: quench front profile at 29

20 seconds. TRAC calculation with 0.25 mm min mesh, 0.3 ms

step, reduced heat flux.

18 Heat conduction equation: quench front profile at 30

20 seconds. TRAC calculation with no axial conduction,

0.25 mm min mesh.

19 Heat conduction equation: quench front profile at 31

20 seconds. TRAC calculation with no axial conduction,

0.25 mm min mesh (exploded view).

20 Heat conduction equation: quench front profile at 32

20 seconds. TRAC calculation with no axial conduction,

0.25 mm min mesh, 0.3 ins, time step.

21 Heat conduction equation: quench front profile at 33

20 seconds. TRAC with no axial conduction, 0.25 mm min

mesh, 0.3 ms step (exploded view).

22 Surface-to-fluid heat flux vs, temperature, at 20 seconds. 34

TRAC calculation with 0.1 mm min mesh, 0.3 mns time step +

theoretical heat transfer.

23 Heat conduction equation: quench front profile at 35

20 seconds. TRAC calculation with 0.1 mm min mesh, 0.3 ms

time step.

24 Heat conduction equation: quench front profile at 36

20 seconds. TRAC calculation with 0.1 mm min mesh, 0.3 ms

time step (exploded view).

25 Rod surface temperatures at 5 elevations, for 3 37

calculations. Calculations are: CONTO0.25 mm,

SHORT=0.25 mm + 0.3 mns, LONG=0.l mm + 0.3 ins.

AEEW - M 2590 vi Viii
NOMENCLATURE

A Surface area

Cp Specific heat at constant pressure

h Heat transfer coefficient

K Thermal conductivity

qK Volumetric heat generation rate

Q Heat transfer (energy)

r Radial cylindrical coordinate

p Density

T Temperature

t Time

z Axial cylindrical coordinate

(SI units)

AEEW - M 2590 i

i X

1 INTRODUCTION

A previous study, [1], examined certain axial effects in the heat

conduction solution of the transient, two-phase flow computer

code TRAC-PFl/MODl [2]. Calculations which simulated the

quenching of the surface of a nuclear fuel rod were seen to have

time step size and, to a lesser extent, axial mesh size

dependencies. The purpose of the present paper is to examine and

explain these dependencies. Similar dependencies may well arise

in other computer codes which employ semi-implicit, finite

difference solution schemes.

Section 2 of this paper gives a brief description of the

TRAC-PFl/MODl computer code. This section concentrates on the

particular aspects of the code that are relevant to this study.

Section 3 presents some results from the TRAC-PFl/MODl

calculations which demonstrate the time step size and axial mesh

size dependencies.

In Section 4 a more detailed examination is made of the

individual terms that contribute to the heat conduction equation.

Various graphical surfaces are generated by over-plotting the

results from several successive time steps.

Finally, Section 5 presents the overall conclusions of this

study.

2 DESCRIPTION OF TRAC-PFl/140D1

TRAC-PFl/MODl is used to perform analyses of Loss-of-Coolant

accidents and other transient's in Pressurised Water Reactors

(PWR's). It is also used to analyse a wide range of related

thermal-hydraulic experiments.

The basic operation of the code is to solve the time-dependent

partial differential equations describing two-phase flow (water

and steam) and heat transfer, by finite difference methods. The

heat transfer equations are treated by using a semi-implicit

differencing technique. The fluid dynamics equations are solved

for one-dimensional components, such as pipes, using a mnultistep

procedure that allows the material Courant condition to be

violated. For a three-dimensional component, such as the reactor

vessel, a semi-implicit differencing scheme is used. The

combined finite-difference equations form a system of coupled,

non-linear equations. They are solved by a Newton iteration

procedure for each time step.

One aspect of the numerical scheme that is relevant to the

subsequent discussion in this paper relates to the coupling

betwen the beat transfer equations and the hydrodynamic

equations. The heat transfer equations might, for example, be

used to model the two-dimensional heat conduction within a heated

cylindrical rod. The coupling with the hydrodynamics equations

AEEW - M 25901 1
takes place via the surface heat transfer between the rod and the

surrounding fluid. This surface heat transfer will be dependent

on the rod surface temperature and several of the fluid's

properties; it provides a surface boundary condition for the heat

conduction equation and contributes to the energy and mass

conservation equations for the fluid. The surface boundary

condition for the heat conduction equation, at time step (n4-l),

is of the form:-

hT n +1 - Tn+l) (1) K br h sTurface - fluid)

The surface to fluid heat transfer contribution to the energy

equation, for time step (n+l), is of the form:-

0 sufae hn (surface - Tftliid) At (2)

to fluid

The point of particular significance in this heat transfer

coupling is that the surface heat transfer coefficient is

evaluated explicitly; it is calculated using rod and fluid

conditions from the previous time step. In later sections of

this paper it is shown that this explicit evaluation, taken

together with the smoothing that is applied to the heat transfer

coefficient, can significantly affect the calculated surface heat

transfer.

3 TRAC-PFl/MODI QUENCHING RESULTS

The calculations originally reported in [1] were hypothetical

simulations of a 1 m, vertical, length of nuclear fuel rod inside

a cylindrical pipe. The calculations were initialised with the.

rod temperatures sufficiently high that the surface, for

elevations above the very bottom, could not be wetted. A

constant flow of water was introduced at th *e bottom of the pipe;

the resulting cooling and ultimate quenching of the rod surface

by the fluid, was then calculated.

Some typical results from the TRAC-PFl/MODl calculation are

presented in Figure 1. This Figure shows rod surface

temperatures, at five elevations, plotted against time, for four

separate calculations. The differences between the four

calculations lie in the size of the smallest axial mesh used in

the finite difference representation of the fuel rod. This mesh

is separate from the mesh used to solve the fluid flow equations,

which was unchanged. It can been seen from Figure 1 that there

is a wide variation in the times at which the rod surface

temperature, for any particular elevation, quenches (ie falls

rapidly to the fluid saturation temperature). It is not

immediately apparent why changing the axial mesh size should have

this effect.

AEEW - M 25902 2
The reason for wanting to change the axial mesh can best be

explained by reference to Figure 2. This Figure shows axial

profiles of the rod surface temperature at successive times, for

one of the calculations represented in Figure 1. It can be seen

from Figure 2 that shortly after the start of the calculated

transient a sharp, or steep, temperature gradient develops; this

gradient, or quench front, effectively separates the hot

unquenched region from the cooler quenched region. As the

transient continues this quench front progress along the rod.

The reason for changing the axial mesh size in the original

TRAC-PFl/MODl calculation was to identify and examine the effects

it might have on the quench front progression.

The quench front region itself is typically only a few

millimetres wide. The TRAC-PFl/MODl solution scheme attempts to

resolve this very steep temperature gradient by inserting an

extra row of heat conduction mesh points, wherever the

temperature difference between adjacent surface nodes exceeds a

user-input value. This value is typically 3*K for mesh points in

the vicinity of the quench front. In order to prevent an

excessively large number of mesh points being used the user also

specifies a lower bound on the axial mesh spacing that can 'have

an extra row of mesh points inserted. The four calculations

represented in Figure I used differing values of this lower

bound; the effective minimum mesh sizes were 2.5 mm, 0.25 mm,

0.1 mm and 0.05 mm. Figure 1 shows that reducing the lower bound

causes the quench front to progress more quickly; it also causes

the quenching to occur at slighly higher surface temperatures.

The semi-implicit nature of the heat conduction solution in

TRAC-PFl/MODl leads to additional complications in trying to

understand the apparent mesh size dependency.

TRAC-PFI/MODl uses a two-dimensional (r,z) cylindrical heat

conduction equation. Azimuthal symmetry is assumed. The

differential equation can be written as

6T .4SS 1 a T a 6

The finite difference form of equation (3), implemented in

TRAC-PFlIMOD1, has implicit differencing in the radial (r)

direction and explicit differencing in the axial (z) direction.

The explicit differencing used for the axial term in Equation (3)

leads to a stability restriction on the maximum time step size

(Atmax) for a particular minimum axial mesh size (Azmin). This

restriction is of the form -

Atmax '" AZ2min

AEEW - M 25903 3
Thus, for the four calculations represented in Figure 1, changing

the lower bound on the mesh size has also changed the time step

size in the calculations. The effect of reducing the time step

alone can be judged from Figure 3. This Figure shows results

from three calculations; a large mesh size case, a small mesh

size case, and a case with a large mesh size but a time step

restricted to 0.3 millisecs. (0.3 millisecs was the average time

step size of the small mesh size calculation). Figure 3 shows

that most of the effect seen in reducing the mesh size is in fact

due to the resultant reduction in time step. This time step

size, and to a lesser extent mesh size, dependency is further

examined and explained in Section 4.

4 DETAILED EXAMINATION OF CONDUCTION TERMS

The previous''section highl~ighte'd the fact that reducing' the time

step size us'ed in the quenching calculatio~ns had changed the

results. In particular, it had caused the rod surface to quench

at a faster rate and from a higher temperature. To a lesser

extent, reducing the axial mesh size had a similar effect. This

behaviour is now examined in more detail by considering the

individual terms of the heat conduction solution.

4.1 Finite Difference Equation

TRAC-PFl/MODl solves a finite difference form of equation (3);

this is obtained by applying an integral method to an appropriate

differential volume. If the resulting finite difference equation

for each node is divided by pCp and by the node volume, then an

equation of the form:

TOTAL = GENERATION-+ RADIAL + AXIAL (3a)

(where each term is in *K/sec)

can be written for each node in turn. The heat generation occurs

internally within the rod so that for the surface nodes the

GENERATION term in Equation (3a) will be zero. For the nodes of

intere-st iEn this section, ie close to the quench front, the

automatic mesh refinement will cause all the node sizes to be

close to the minimum allowed.

4.2 Surface-to-Fluid. Effects

Plots presented later in this section show the individual terms

of Equation (3a), for the surface nodes, drawn as a function of

the wall temperature. First, however, it is useful to examine

one component of the RADIAL term, namely the surface heat

transfer between the rod and the coolant.

Figure 4 shows a plot of surface heat transfer coefficient versus

surface temperature. Results from the TRAC-PFl/MODl calculation

with a 0.25 mm. minimum axial mesh are displayed as a sequence of

points, drawn as numbers. The results are taken from each

surface node for 11 consecutive time steps at approximately

20 seconds into the calculation. The fluid conditions will

AEEW - M 25904 4
normally only change slightly during 11 time steps so it is

reasonable to expect that the points representing heat transfer

coefficient versus wall temperature will lie on a curve. In

Figure 4 the points labelled "1" are from time step 1 of the

sequence and so on. Points labelled ""~ and "A" are for time

stops 10 and 11 respectively. The curve traced out by the points

labelled "I" to "A" is the effective heat transfer curve for this

particular calculation,' at 20 seconds. Figure 4 also shows the

theoretical heat transfer curve derived for the particular fluid

conditions present in the TRAC-PFl/MODl calculation. This curve

was calculated by evaluating the TRAC-PFI/MODl heat transfer

correlation separately, in a stand-alone manner, for the range of

surface temperatures of interest.

Figure 4 shows that once a surface node is cooled below

approximately 600*K its surface heat transfer coefficient

increases sharply. A theoretical maximum is shown to be reached

at appoximately 470*K; this corresponds to the point of critical

heat flux. However, the most striking feature of Figure 4 is the

fact that the achieved, or effective, heat transfer curve is

significantly below the theoretical curve. Many values are

40-50% below the theoretical values and the critical heat flux

temperature appears to be 20*K lower. These differences are

further highlighted in Figure 5 which shows the surface heat flux

values corresponding to the coefficients given in Figure 4.

4.2.1 Explicit Evaluation and Smoothing Effects

The differences observed between the effective and theoretical

heat transfer curves arise from two separate aspects of the

TRAC-PFl/MOD1 solution scheme. Firstly the explicit evaluation

of the surface heat transfer coefficients; this means, for

example, that the surface temperature from the previous time step

is used to evaluate the new coefficient. Secondly, the smoothing

and limiting techniques applied to the calculated heat transfer

coefficient; 55% under-relaxation is used (55% old-time value +-

45% new-time value), followed by the restriction that,

essentially, the resulting new value is no more than twice the

old-time value. These techniques are applied on a per-time step

basis and not on a per-unit time basis; thus, for example, during

the rapid increase in coefficient shown in Figure 4 some time

step size dependency will occur.

The TRAC-PFl/MODl results shown in Figure 6 will allow these two

aspects of the solution scheme to be considered separately. The

results shown in Figure 6 are from a calculation in which the

surface heat transfer smoothing and limiting have been removed.

The theoretical heat transfer curve has been derived for the

fluid conditions present at the end of the time step sequence.

The effect of the explicit evaluation of the surface heat

transfer coefficient can be clearly seen in Figure 6 for time

steps 3 onwards (ie points numbered 3-9, * and A). For example,

the point marked 'W4', at approximately 485*K, has a heat transfer

value that corresponds to the theoretical curve evaluated at the

temperature of the point marked "3", close to 500*K. Similarly

the point "5' value corresponds to the point "4" temperatures and

AEEW - M 25905 5
so on. (This correspondence does not work in Figure 6 for the

points marked '3", "2" and "1" because the fluid conditions at

those time steps were slightly different to those used to derive

the theoretical curve). In other words the surface temperature

T, is at time step (n+l), but the heat transfer coefficient h, is

at time step (n). This point is confirmed in the formulation of

equation (1).

Figure 6 demonstrates that, in a region where the heat transfer

coefficient is changing rapidly, the explicit evaluation of the

coefficient can lead to a significant deviation of the effective

heat transfer curve from the theoretical one. In Figure 4 the

deviation also includes the under-relaxation and limiting

effects; the difference between the effective and theoretical

heat transfer is greater, particularly with regard to the peak

value.

4.2.2 Time Step Size Effects

The calculation for which results were* presented in Figure 4 used

time steps that were in the range of 5-10 milliseconds. Figure 7

shows results from an equivalent calculation in which the time

step was constrained to be no greater than 0.3 milliseconds. The

TRAC-PFl/MODl calculated values have been drawn every 24 time

steps, ie every 7.2 milliseconds, as this corresponds to the

average time step size of the earlier calculation. Figure 7

shows that reducing the time step size has caused the effective

heat transfer curve to follow closely the theoretical curve.

Figure 8 shows the surface heat flux values corresponding to the

coefficients given in Figure 7. A comparison with Figure 5

emphasises the effect of reducing the time step size.

Clearly, reducing the time step size has led to an increase in

the effective surface heat flux for surface temperatures between

approximately 450*K and 620*K. This is likely to be a

significant factor in explaining the time step size effect seen

in Figure 2, for example. However, as the next subsection shows,

the presence of axial effects must also be taken into account.

4.3 Axial Conduction Effects

In an attempt to isolate the separate contributions of the RADIAL

and AXIAL terms of equation (3a) several TRAC-PFl/MODl

calculations were carried out with the AXIAL term artificially

set to zero. This prevents any axial conduction of heat within

the rod. The results, shown in Figure 9, are somewhat

surprising, With the AXIAL term removed the calculations show

virtually no sensitivity to either time step or axial mesh size.

Figure 10 shows the surface heat flux, plotted as a function of

surface temperature, for the first NO-AXIAL conduction

calculation. The TRAC-PFl/MODl results are similar to those

shown in Figure 5 for the standard calculation; the effective

heat flux curve is again significantly below the theoretical

curve. Figure 11 shows the equivalent results from the NO-AXIAL

conduction calculation with the time step restricted to.

AEEW - M 25906 6
'0.3 milliseconds. The TRAC-PFl/MODl results now closely follow

the theoretical curve in a similar manner to the standard

calculation results presented in Figure 8.

This shows that reducing the time step size in a calculation

without axial conduction causes the effective surface heat flux

curve to follow closely the theoretical curve. However, this

does not affect the overall quenching behaviour to any

significant effect. The time step size effect seen in the

standard calculations must, therefore, depend on more than just

the change in the effective surface heat transfer.

4.4 Quench Front Profiles

In previous sections effective surface heat transfer curves have

been generated by over-plotting heat transfer values from a

sequence of consecutive time steps. A similar technique can be

used to generate a quench front profile of the individual terms

of equation (3a), for nodes at the rod surface.

4.4.1 Calculation With 0.25 mm Minimum Axial Mesh

Figure 12 shows the quench front profile at 20 seconds for the

standard TRAC-PFlIMOD1 calculation with a 0.25 mm minimum axial

mesh. Points labelled "A" represent the magnitude of the AXIAL

term of equation (3a), points labelled "R" represent the RADIAL

term and points labelled "T" represent the TOTAL term, ie-tF-e

sum of the AILand RADIAL term. For the sequence of 11 time

steps plott-e-d in Figu-re 12 the points representing the separate

terms trace out an effective quench front profile.

The role of the AXIAL term can be readily seen from Figure 12.

At the high temperature end of the region the AXIAL term is

negative ie tending to cool the rod surface. ThIn-ct for

temperatures above 5500K the AXIAL term makes up almost all of

the TOTAL term. (The RADIAL It-e-r-m as positive values above

approximately 585*K because the heat being transferred from

inside the rod to the surface exceeds that being transferred from

the surface to the fluid). At the low temperature end of the

region the AXIAL term is positive ie it is opposing the cooling

rate genera-ted by the larger negative RADIAL term. Thus, the

overall effect of the AXIAL term is to transfer heat from the

high temperature end to -the low temperature end where the RADIAL

term, largely governed by the surface-to-fluid heat flux, is

large and negative.

In Figure 12 the magnitude of the TOTAL (WT/bt) term becomes

small, for temperatures above approx=imately 655*K. This

corresponds to the temperature of the "knee" in the temperature

versus time plot for the 2.5 millimetre minimum mesh calculation,

shown in Figure 1. For temperatures above this value the rod

surface is cooled comparatively slowly. However, for

temperatures below this value the rate of temperature fall

increases very rapidly, until the surface is quenched. It can be

seen from Figure 12 that, at least for this calculation, the

AEEW - M 25907 7
-temperature at which this knee occurs is governed by the onset of

the large negative AXIAL values.

The movement of the quench front region along the rod can

conveniently be characterised by the movement along the rod of

the knee in the temperature profile. The actual temperature at

which the knee is maintained will be dependent on the details of

the heat conduction solution within the quench front region

itself. Figure 13 presents an exploded view of the AXIAL and

RADIAL terms taken from Figure 12 in the region of the knee. For

temperatures above 660*K it can be seen that the AXIAL term is

essentially zero and the RADIAL term is negative (jie cooling the

surface) and increasing i7n -magnitude with increasing surface

temperature. For temperatures below 660*K the RADIAL term

becomes negligible and the AXIAL term very rapidl-ybecomes large

and negative. This large.neg7t"Ive AXIAL term rapidly cools the

cladding surface and allows the que-nc-h front or temperature knee

to move forward.

4.4.2 Calculation With 0.25 mm Minimum Mesh and 0.3 ms Time Step

Figure 14 presents the quench front profile at 20 seconds f or the

standard calculation with the reduced time step size. Comparison

with Figure 12 shows that the magnitude of the peak negative

RADIAL term has increased significantly; this is in line with the

increaed surface heat flux seen by comparing Figure 8 with

Figure S. The magnitudes of the peak AXIAL terms (positive at

low temperatures, negative at high temperatures) have also

increased significantly, leading to increased magnitude TOTAL

term values. In particular, the increased magnitude nega~tiv-e

AXIAL terms at high temperatures have moved the temperature of

the knee up by approximatley 20*K. This is borne out by

comparing the temperature versus time profiles shown in Figure 3.

An exploded view of the AXIAL and RADIAL terms close to the knee

is given in Figure 15. C-omparison -wit-hFigure 13 shows that the

RADIAL terms ahead of the knee, from the two calculations, lie

approximately on the same curve.

Figure 14 shows that with the reduced time step the magnitude of

the AXIAL (and hence TOTAL) term increases more rapidly, as the

surf-ace temperature f-alls below the knee temperature, than for

the standard calculation shown in Figure 12. This is consistent

with the observed faster progression of the knee in the small

time step calculation.

Thus the observed time step size dependency appears to be related

to the increased magnitude AXIAL terms at high temperatures. Two

questions remain unresolved however: why are the AXIAL terms

increased in magnitude, and why does the calculatio-nwith no

axial conduction show no time step size dependency. These two

questions are now addressed in turn.

4.4.3 AXIAL Term Time Step Dependency

The significant increase in the AXIAL term magnitude shown in

Figure 14 could be due to two possirble effects. Firstly, the

AEEW - M 25908 8
large increase in the peak RADIAL term magnitude will have

changed the axial temperature profile in the quench front region.

This is likely to change the AXIAL term values as they are,

essentially, derived from the axial temperature profile.

Secondly, reducing the time step size may in itself have changed

the AXIAL term values as they are evaluated explicitly. To

resolve this issue a calculation has been performed using the

reduced time step size but with the surface-to-fluid heat flux

modified so that it remains at the level shown in Figure 5 rather

than the increased level shown in Figure S. This was achieved by

reducing the critical heat flux value (CHF) used by TRAC-PFl/MODl

in evaluating the heat transfer coefficients.

Figure 16 shows the effective surface heat flux curve from this

new reduced-time step calculation. It is in fact quite close to

the effective curve presented in Figure 5 for the original

calculation. Figure 17 shows the quench front profile for the

new calculation at 20 seconds. Both the RADIAL and AXIAL term

curves are very similar to the correspond'in-gcurves sh-own in

Figure 12 for the original calculation. Thus the AXIAL term

values have no time step size dependency of their -own (Twithin

the time step range considered) but rather they reflect the time

step size dependency of the RADIAL term. This in turn reflects

the time step size dependency of-the surface-to-fluid heat flux;

as previously shown this is due to the explicit heat transfer

evaluation and smoothing techniques inherent in the solution

scheme.

4.4.4 Calculations With No Axial Conduction

Figures 10 and 11 showed the effective surface-to-fluid heat flux

curves for two calculations with no axial conduction. Reducing

the time step size to 0.3 milliseconds caused the effective curve

to follow the theoretical curve (Figure 11) but did not, however,

change the overall quenching behaviour (Figure 9).

Figure 18 shows the quench front profile for the large time step

calculation. As the AXIAL term is zero the TOTAL term is simply

equal to the RADIAL term. The effective RADIAL term curve in

Figure 18 is _sir`=har to the RADIAL term curve shown in Figure 12

for the standard calculation. However, the lack of an AXIAL term

means that the TOTAL term becomes small at a lower temperature

than in Figure flT=,e the temperature knee is maintained at a

lower temperature. This is confirmed by the temperature versus

time profiles shown in Figure 9. Figure 18 also shows that the

magnitude of the TOTAL term increases slightly less rapidly, as

the surface temperature falls below the knee temperature, than

.for the standard calculation. This is consistent with the slower

progression of the quench front in the calculation with no axial

conduction. Figure 19 shows an exploded view of the TOTAL term

in the region of the temperature knee. A comparison wi~ththe

standard calculation results, given in Figure 13, shows that the

RADIAL term values lie essentially on the same effective curve.

Howevr, in the no axial conduction calculation the knee is

maintained at a lower temperature.

AEEW - M 25909 9
'The quench front profile for the no axial conduction calculation

with the reduced time step size is presented in Figure 20. The

peak magnitude of the RADIAL term has increased, compared to

Figure 18, in line with -the increase in the surface-to-fluid heat

flux shown in Figure 11. However, at the high temperature end of

the region the values are unchanged. This is further borne out

in the exploded view shown in Figure 21.

Thus in a calculation with no axial conduction, although the peak

TOTAL term magnitude is increased, the TOTAL term values at the

high temperature end of the quench front region are unchanged

when a small time step size is used. This is consistent with the

observation that the overall quench front movement is unchanged

when a small time step size is used.

4.4.5 Axial Mesh Size Effects

Having established and examined the time step size dependency it

is now worthwhile examining any mesh size effects. TRAC-PFl/MoDl

will normally automatically reduce the time step size when small

axial mesh sizes are used, because of the explicit evaluation of

the axial terms. Therefore, to establish any genuine mesh size

effects a comparison has to be made with a calculation that

already uses a sufficiently small time step size.

A calculation has been performed using a 0.1 mm. minimum axial

mesh and a 0.3 millisecond time step size. Figure 22 shows the

effective surface-to-fluid heat flux curve for this calculation;

this can be compared to the curve in Figure 8, which used a

0.25 mm minimum mesh. The two effective curves are very similar;

the smaller mesh curve lies slightly closer to the theoretical

curve at the peak value. The smaller mesh size gives more nodes,

and hence a better resolution, in the region of peak surface-tofluid heat transfer values.

Figure 23 shows the quench front profile for the new calculation;

this can be compared to Figure 14 for the larger mesh size

results. The results are again very similar apart from the peak

AXIAL and RADIAL values at the low temperature end of the region.

An exploded vie-w of the AXIAL and RADIAL terms at the high

temperature end of the region is shown in Figure 24.Threut

are very similar to those shown in Figure 15 for the larger mesh

size calculation. This suggests that the overall quench front

progress should be very similar for the two calculations.

.The surface temperatures versus time for the new calculation are

shown in Figure 25. A comparison is made with the larger mesh

size calculation and also the original larger time step size

calculation. The new calculation shows that there is a small

axial mesh size dependency, but that it is very small compared to

the time step size dependency.

5 SUMMARY AND CONCLUSIONS

The purpose of this paper is to examine and explain the time step

and mesh size dependencies observed in calculations of the

ABEW - M 2590 110
quenching of a nuclear fuel rod. Both effects have been shown to

arise from an underlying dependency in the surface-to-fluid heat

transfer. The time step dependency occurs because the heat

transfer coefficient is evaluated explicitly, ie using values

from the previous time step, and because under-relaxation is

applied to the newly calculated coefficient. This dependency

will be particularly noticeable whenever the heat transfer

coefficient is changing significantly from one time step to the

next, such as occurs during quenching. The smaller mesh size

dependency appears to arise from changes in the spatial

resolution at the calculated heat transfer coefficient close to

its peak value.

The paper has shown that changes in the surface-to-fluid heat

transfer affect the overall quenching behaviour by virtue of

changing the axial temperature profile; this changes the axial

conduction terms in the overall rod conduction equation. It is

changes in the axial conduction terms, at the high temperature

end of the quench front region, that alter the overall quench

front progression. In calculations where the axial conduction

term was artificially removed, changes to the surface-to-fluid

heat transfer did not affect the overall quenching behaviour.

The findings can be summarised in the following conclusions:-

5.1 The studies described in this report have identified a

significant time step size dependency in the solution

obtained from a coupled system of heat transfer and twophase flow partial differential equations.

5.2 The time step size dependency of the solution arises from

the time step size dependency of the surface-to-fluid heat

flux; this flux is the coupling between the heat transfer

equations and the fluid flow equations. The dependency

occurs as a result of the explicit evaluation of the

surface-to-fluid heat transfer coefficient, and as a result

of the time 6tep-to-time step smoothing techniques applied

to the coefficient.

5.3 For the TRAC-PFl/MODl quenching calculations described in

the report the time step size dependency of the solution

dissappears if the axial conduction term of the heat

conduction equation is removed. This is because the

surface-to-fluid heat flux time step dependency affects the

overall solution only by changing the size of the axial

conduction terms.

5.4 The studies described in this report have also identified a

small axial mesh size dependency; this is, however, much

smaller than the time step size dependency. This

dependency again appears to arise from a small mesh-size

dependency of the surface-to-fluid heat transfer, 'mainly in

the region of high and rapidly changing heat transfer

values.

AEEW - M 2590 111
5.*5 The time step size dependency represents a potential

problem in the use of the TRAC-PFl/MODl code, with regard

to running times. The numerical solution scheme for onedimensional components employs a multistep procedure that

allows the material Courant condition to be violated. This

ability to use large time step sizes will be restricted if

small time steps are needed for the heat transfer

evaluation part of the scheme. Further work is needed to

improve or replace the explicit heat transfer evaluation

and to remove the time step size dependency from the heat

transfer smoothing techniques.

6 REFERENCES

(1) O'Mahoney, R. A Study of Axial Effects in the

TRAC-PFl/MOD1 Heat Conduction Solution During Quenching.

AEEW - M 2552, PWR/HTWG/P(89)686, June 1989.

(2) TRAC-PFl/MODl. An Advanced Best-Estimate Computer Program

for Pressurised Water Reactor Analysis. Los Alamos

National Laboratory Report. LA-10157-MS, NUREG/CR-3858.

AEEW - M 2590 112
TRAC-PF1/MOD1 013.0 (bO3e)

0

900-

8501

800.

750-

700-

650

S600

550-

500-

450-

4001.

FIGURE 1

HIGH TEMPS

LOW FLOWS

5 10 15 20 25

REACTOR TIME ,SECONDS

30

lWinfr it h

XI

ROD SURFACE TEMPERATURES AT 5 E~LEVATIONS, FOR 4 DIFFR NT ESH SIZES

MIN AXIAL MESMS ARE: CONT=2. 5MM, SHORT=O. 25MM, LONG-O. 1MM, THICKO0.05MM
TRAC-PFI/MODI 013.0 (bO3e)

I,

rn

1100

1000

900

800

700

600

500

Annf

. ......... --

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

-- - ---- - -- - - -- - - -- -- - - --- -- - - - - -t . . . .. . - - - - - - -

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

T- T-0T3 ,01 

-- m -- - - - ------r -- - -- - -- - - - -- - - - -- - - - - -- -- -- -- --

...........

HIGH TEMPS

LOW FLOWS

-u

0 20 40

DISTANCE ALONG IROD

60 80 100

FIGURE 2 WinEri

-11

AXIAL PROFILE OF ROD SURFACE TEMPERATURE AT 4 DIFFERENT TIMES

PROFILES AT 0, 10, 20 & 30 SECS, FOR 0.1M MAIN AXIAL MESH
TRAC-PF1/MODI 013.0 (bO3e)

0

lN,

a

-4

N

(A

900 -

850-

800-

750-

700

s1 650

S600-

550-

500-

450-

400 -

0

FIGURE 3

HIGH TEMPS

LOW FLOWS

W in Fr It h

5 10 15 20 25 30

REACTOR TIME ,SECONDS

-q,

m

L4

ROD SURFACE TEMPERATURES AT 5 IELEVATIONS, FOR 3 CALCULATIONS

CALCULATIONS ARE: CONT=O.25MM, SHORT=O0.O5NM, LONG=O0.25MM + O.3MS DT
TRAC-PFl/MOD1 013.0 (bQ3e)

120000 I

to

1000001

800001

C14 60000

-- -- -- - - \- --- -- - - - - ---- U--------------A- - - - -- --------

< THEORETICAL HEAT TRANSFER

-- --- - - - - - .4. ----- -- - - - - - - - - - ---. . . . .

:1 SYMBOLS:1 - 91 *p A ARE TRAC CALCUJýATION VALUES

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

-- - -- - I - - - - - - -- - - - - - - - - - - - -

----- ----I

N,---- ------- -------------- - -----------

40000~

20000

OL

400 450 500 550 600

DEG. K

650 700

SURFACE TEMPERATURE

FIGURE 4 Winfrithl

-r,

SURFACE-TO-FLUID HEAT TRANSFER COEFF. VS TEMPERATURE, AT 20 SECONDS

TRAC CALCULATION WITH 0.25MM MIN MESH + THEORETICAL HEAT TRANSFER
TRAC-PF1/MOD1 013.0 (bO3e)

mq

I.'

('I to

0

7000000

6000000

5000000

4000000

4c

I-I

-.4

-o

3~

-4

0

~0

0,

to

-.4

I~3

U'

3000000

2000000 -

1000000 --

0 - 400 450 500 550 600 650 700

SURFACE TEMPERATURE IDEG. K

FIGURE 5 Win fri thl

0n

Lc

SURFACE-TO-FLUID HEAT FLUX VS TEMPERATURE , AT 20 SECONDS

TRAC CALCULATION WITH 0.25MM MIN MESH + THEORETICAL HEAT FLUX
TRAC-PF1/MOD1 013.0 (bO3e), NO HEAT TRANSFER SMOOTHING

120000

mn

rn

IQ

00

-o

C)

co

to

100000k --------

:6-

, I

-I. I'

'I

"8' * - . - I. - -

'I

47

I,

- - -I --------- - - - -

I THERETICAL HEAT TRASFER

-- - - - - - - - - - - - - - - - - - - ---------- 80000k --------

04' A

IIliI~bz I - ýl I F iia TKRJ4L i ~JAX4ALUN VAWtkS

- -- - - - - - - - - - - -.- - ---.- - - - - - --.- - - - - - - - - -

60000k ---------7/

40000

9

20000k1---- -- ---- - - - - - - - - - - - - - - - - - - -

*A-4

* ~~~7 Q D7 It A~f~ - 01

40

I

i

0 450 500 550 600

DEG. K

650 700

SURFACE TEM4PERATURE

FIGURE 6 lWinfrithl

c-)

SURFACE-TO-FLUID HEAT TRANSFER COEFF. VS TEMPERATURE, AT 20 SECONDS

TRAC CALC. WITH 0.251MM MIN MESH, N NO SMOOTHING MOHN + THEORETICAL HOEIA H..TT.
TRAC-PF1/MOD1 013.0 (bO3e)

rn

f-n

to

0

120000

100000 -----------

i

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

SYMBOLS;1 - 9p *1 A ARE TRAC CMOWUI~fON VALUES

--------

80000

clqI 60000k --------- I,

40000 ..--

-o

-4

L"

5\

20000 k-.7 --------------- -------------------------

01

40

I Ion NX M&LOANtahrmMiLAL-AL-m-m

0 450 500 550 600

,DEG. K

650 700

SURFACE TEMP2ERATURE

FIGURE 7 Wtinfrit

-n SURFACE-TO-FLUID HEAT TRANSFER COEFF. VS TEMPERATURE, AT 20 SECONDS

TRAC CALC. WITH 0.25NM MIN MESH,03M 0.3MS TISEP+ TIMESTEP + THEORETICAL HORICLH. H.T.
TRAC-PF1/MOD1 013.0 (bO3e)

0

8'

7000000

6000000

5000000

4000000

3000000

00

-4

2000000-

1000000 -

0-

400 450 500 550 600 650 700

SURFACE TEM4PERATUJRE I DEG. K

FIGURE -8 lWinfritdi

-Ti

C)

C

P1

02

SURFACE-TO-FLUID HEAT FLUX VS TEMPERATURE , AT 20 SECONDS

TRAC CALC. WITH 0.25MM. MIN MESH, 0.3MS TIMESTEP + THEORETICAL H.T.
TRAC-PFl/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

mn

coI

to

900-

850-

800-

750

700

650

600-

550-

500-

450-

400L

0

HIGH TEMPS

LOW FLOWS

5 10 15 20 25 30

REACTOR TIME ,SECONDS

FIGURE" 9 Iintril

ROD SURFACE TEMPERATURES AT 7 ELEVATIONS, FOR 3 NO-AXI.AL CAMCS

CALCULATIONS ARE: CONT=O .25MM, SHORT=O .25NM+ . 3M5, LONG-. O5NM+ . 3MS

Nd

-In

m

to
TRAC-PF1/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

7000000

'Ii

~1~

01

to

0

cq

4'

-u

0-%

CA

01 OEII-I -wm

400 450 500 550 600 650 700

SURFACE TEMPERATURE r DEG. K

FIGURE-10 1WInfitjh

C-n

C)

SURFACE-TO-FLUID HEAT FLUX VS TEMPERATURE, AT 20 SECS

TRAC CALC, NO AXIAL CONDN, TRA CAC, O MIN AIALCONN,0.25MM MESH + THEORETICAL HEAT FLUX
TRAC-PF1/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

rn

to

0

7000000 -

6000000 -

5000000 -

4000000 -

3000000 -

2000000 -

1000000 -

0-

400 700

SURFACE TEMPERATURE I DEG. K

FIGURE 11 Winfr'it hý

-I,

C

I,,

SURFACE-TO--FLUID HEAT FLUJX VS TEMPERATURE,

TRAC, NO AXIAL CONDN, 0.25MM MIN MESH, 0.3

AT 20 SECONDS

,MS STEP + THEORETICAL H.T.
TRAC-PF1/MOD1 013.0 (bO3e)

rn

#0

to

-4

00

to

14)

U)

CDrz~Q

6000

4000

2000

0

-2000

-4000

-6000

-8000

-10000

44

* SYMBOLS As R, T ARE TRAC CALCULAT~ION VALUES

A < AXIAL TERM (A).

*R e2RADIAL TERM (R)

-- - - -- - - - - - - - -- - - I - - - - - - - - - - - - - - - - - - - - - -

00 450 500 550 600

DEG. K

650 700

SURFACE TEPERATURE

FIGURE 12

HEAT CONDUCTION EQUATION: QUENCH FRO]

TRAC CALCULATION WITH O.25NvM MIN MESH

Winfrithi

-n,

0

NTPROFILE AT 20 SECONDS
TRAC-PF1/MOD1 013.0 (bO3e)

5 It

20 -- - - - - - - - - - I -- - - -A1- - - - -1 - -- - - - - - - -. . . . .

0a -25'9

-10

SMC TEPRr 9 G 9

-IUE13 1W 9r

HEA CODCTO EQAIN QUNC FRN PRFL AT2EOD

TRA CACLTO WIT 0M MI9EHEPLDDVE

- r9
TRAC-PF1/MOD1 013.O (bO3e)

6000 P1

P1

I.3

Ut to

0

4000

2000

-----------

A

---------- :4--

SYMBOLSA, R, T AR~E TRAC CALMULTI101 VALUJES

A 1 AXIAL TERM (A)

A

-A...A.---------------------------------------.... ...- .... ........

A n ca R:R RR RRW R; I U ~ ~ J I Ii - I 1=~~~~~ -

0'4

r-4

0

U)

rzl

-2000

-4000

-6000

-8000

-10000

40

.. .. 4 . ---- -... 

0 P

A R

ARA

R~ TT

T I IT T I TOTAL DTDT (T)

T T T1 PT

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

R < RADIAL 'TERM (R)

-- - - - - - - - - - - ----9 - --- - - - - -- - - - - - - - - - - - - -- - - - - - - - - -----------

0 450 500 550

SURFACE TEM~'PERATURPE

600

DEG. K

650 700

FIGURE 14 lwinfith

~1~

C

I,,

HEAT CONDUCTION EQUATION: QUENCH FRONT PROFILE AT 20 SECONDS

TRAC CALCULATION WITH 0.25I1M4 MIN MESH, 0.3MS TIMESTEP
TRAC-PF1/MOD1 013.0 (bO3e)

[I.

IT'

mn

R

R

R

R AXIA EWM (A)

A A 0 I I II I

A

A'

A - a

-51----- ------------- --- ------ * --- --- --- --- --- --

if Urn * if lyw a * *-* "if

U I)J .LhLL

9 1

C') .... ... .. . .... .I -10 1 - -- - - -- - - - - -- - - - --- - - - -- - - -

-15 k------ 1--------- I--------

-o

-4

-20 .. - - ... --..-S- J --- --- ---

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

A .

-25

60'0 620 640 660 680 700

DEG. K

720 740

SURFACE TEM4PERATURE

FIGURE 15 1WInfrit

~'1

5

c

m

-A

U,

HEAT CONDUCTION EQUATION: QUENCH FRONT PROFILE AT 20 SECONDS

TRAC, CALCULATION WITH 0.25MM MIN MESH, 0. 3MS TIMESTEP [EXPLODED VIEW]
TRAC-PF1/MOD1 013.0 (bO3e), REDUCED HEAT FLUX

7000000

:SYMBLS 1 - %~ *,, A ARE PEDUCBD-CH MW: CALCATICU MMIE

'ii

m

r

(J~

~0

0

60000001-.

5000000

4000000 L-

(~~l

30000001--

* I THE~OPMEICAL WWEVA!FD (UMODIFIEP)

- - - - --- 

-N - - - - - - - -- - -- -- ---.- -- --- --.-.- - -- --

---- -...---- 

------- ------ ----. . .. . ... .. . ..--- ---.--- --

9

* ~7 *

--- I ..-------

------.

5

2000000W00

10000001--

0L

40

* 4 A 3 ~*-.

0 450 500 550

SURFACE TEMIPERATURE

600

DEG. K

650 700

FIGURE 16 1Winfrithl

-I,

SURFACE-TO-FLUID HEAT. FLUX VS TEMvPERATURE, AT 20 SECONDS

TRAC CALC WITH 0.25MM MIN MESH, 0. 3MS STEP, REDUCED CHF + THEORETICAL
TRAC-PF1/MOD1 013.0 (bO3e), REDUCED HEAT FLUX

I

rv,

to

0

6000

4000

SYMOLS: At R, T AME TRAC C LATIO VALUE

~AA A <AXIAL AA TERM :(A)

A A'

20001

A A:

0 n n9I n n

m I nh1 m r----------.1 -- WEB.........

C,,

-2000 k----

R~

A AA

ARR

It R K K ' IAm~

Av~~ Av~ At' T T~ T ~ ' T~i LJLLU IIUL 9.L_

R *

RR RADIAL TERM (R) -4000ý- --

-u

N

t

-4

C-,

N

-u

to

S.-.

N

U,

-6000

-8000 -----------

I . . . I.. . . . . . . . . . . . . . .

-100001U

400 450 500 550

SURFACE TEMPRATURE

600

DEG. K

650 700

1WInfr'ithi FIGURE 17

C)

-4

HEAT CONDUCTION EQUATION: QUENCH FRONT PROFILE AT 20 SECONDS

TPAC CALCULATION WITH 0.2 5M! MIN MESH, . MS ST-EP, EDCE HFAT FLUX
TRAC-PF1/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

0

6000

4000

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

2000 F ---------------- --------------- - - - - - - - - - - . . . . . . . .

0 isissal i -. ITI [a T ST ITT

U)

T:

: T T

-2000 k- - < TOTM 6Mln (T) = RAatAL TERM .T T. TIý- - - - -- - -- - - - - - -- - - - -

-4000 k---

0

0-%

N

(if

-6000k --------

-8000 -- - -- -- - -- -- * I I

* I I

* I I

* I

* I I

* I I

-100001

40 0 450 500 550

SURFACE TEMPERATURE

600

DEG. K

650 700

FIGURE 18 1WInfrithl

cxo

HEAT CONDUCTION EQUATION: QUENCH FRC

ThAC CALCULAION WITH NO AXIAL CONDN,

)TPROFILE AT 20 SECONDS

0.25MM MINMESH
TRAC-PF1/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

2~).

m

P1

I~J

U'

~0

0

(.A~

I-.

0

-5

TI T I I I

T I I1 f Co

C,

rz~

-101-----------...

TOTAL DT~ir (T) = RADEAL TEPM

T

f V

-151 - -- ------ -T--- - ---- ------ - --- -- .... .. :-- -- -

-U

=

-g

C)

-U

£0

U'

-20 -- - - - - - - - -

ITI

- - - - - - -

-251

6(0 620 640 660 680 700

DEG. K

720 740

SUJRFACE TEMAPERATURE

FIGURE, 19 lWinfrdith

C

HEAT CONDUCTION EQUATION:QU

TRAC CALC WITH NO AXIAL CONDN,

,NCH FRONT PROFILE AT 20 SECONDS

0.25MM MIN MESH 0.21~'4 t EXPLODED MN MSH VIEW ]
TRAC-PF1/MOD1 Y13.0 (bO3e), NO AXIAL CONDUCTION

i

6000

fyi

r

'-U~0

0

(J~

-U

N

z

-g

C,

N

-u

w

(0

-.1

01

SYMBOLS: T AR~E TRAC CAICULAMICN VALUES

4000 k- - - - - ---- - - - - - .

2000

0

.*.. . .. . . . .. .. . . - - - - - - - - - - - - - - - - - - - - - - - -

49MM" i i iI1 muTTT=T9X TIT N :: :

U) aT

-2000 ------------

*~ 1T T T

T

TTOTAL DTDT (T) = RWDIAL TERM

-40001ý - - -- --- a - - - - - -I- - - - - - -- - - - - - - - - - - -:- - - - - -

T s

-6000 I- - - - - - .. . . .. . .- - - - --- - - - -- - - - -

-8000

-i Aflflf

-- - - - - - . . . . . . . . . : - - - - - - - - - - -

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

4--~00 450 500 550 600

,DEG. K

650 700

SUR~FACE TEMPERATURE

FIGURE 20 1WInfri th1

-yi

C)

C

P1

0

HEAT CONDUCTION EQUATION:QU

TRAC CALC WITH NO AXIAL CONDN,

NCH FRONT PROFILE AT 20 SECONDS

0.25MM MIN MESH, 0. 3MS TIMESTEP
TRAC-PF1/MOD1 013.0 (bO3e), NO AXIAL CONDUCTION

rn 0

-5 -- - - - - - - - -

----------

'T Ti .......... ------ 'I

U)

-10 T .. . .

TOTAL DTdT (T) =RADALTIM 4j T y

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

-15

--4

14

N

CA

-201 --------

T

--- -- --- ---- - ----- - ---- --- - - - - - - --..- -.-..

-251

600 620 640

SURFACE

660 680 700 720 740

TEMPERATURE I DEG. K 1WInfrih

FIGURE 21

In

t'5

HEAT CONDUCTION EQUATION:

TRAC WITH NO AXIAL COND,

QUENCH FRONT PROFILE AT 20 SECONDS

0.25MM MIN MESHFO.3MS STEP [EXPLODED VIEW]
TRAC-PF1/MOD1 013.0 (bO3e)

7000000

6000000

I'3

01

0

50000001

40000001

C4J

3000000

- - - - ---------------- ~---- ------------- ------------ ------------ ......

6 SIMB0LS 1 - 9, *1 A jRE TRA CMOJAIION OE

*~ V-

-- ----------- I---6x--------L-------------

------------.-------------------- A...........

- - -1 -

2000000t

-u

N

-I

CI)

N

-u

w

10

U'

1000000

0t

400 450 500 550 600

DEG. K

650 700

SURFACE TEMPERATUJRE

FIGURE 22 1WInfrithl

-I, SURFACE-TO-FLUID HEAT FLUX VS

ThAC CALC WITH 0.1MM MIN MESH,

ITEM1PERATURE,

0. 3MS TINESTEAP

AT 20 SECONDS

+ THEORETICAL H.T.
TRAC-PF1/MOD1 013.0 (bO3e)

*A

6000

m

~0

A

*A SYMBOLS:At R, T ARE TRAC CALCULATION VALUES

4000 k----------- IA i

2000

0

-- - - - - - - - - -

( i---- AXIA TEM (A)

A~~~~ ~ ~ : I

A

AA

A m S al

cl)

--- --- --- 7 T -T. rT -

A A

-2000 A - - - - - - - - -

* P t~lT

yT T TTT rTOTAL DTDT (T)

-40001k------

* R

R- - - - - - - -- - - - - - -- - - - -

* R A

c"

to

*-4

0,

-60001. ------

-8000 -- - - - -- -- - -- --

~< : ADIAL:TEFNM(R)

R

-100001L

40

I - - a I a

0 450 500 550 600

,DEG. K

650 700

SURFACE TEPERATURE

FIGURE 23 1WInfrith1

c)

;o

HEAT CONDUCTION EQUATION: QUENCH FRONT PROFILE AT 20 SECONDS

TRAC CALCULATION WITH 0.1WM MIN MESH,0.3STESE 0.3MS TIMESTEP
TRAC-PF1/MOD1 013.0 (bO3e)

ON

(~)

-u

-4

-u M4 660 680 700

SURFACE TEMPERATURJE ,DEG.IK

FIGURE 24

HEAT CONDUCTION EQUATION:QU

TRAC CALC WITH 0. 114 MIN MESH,

1Winfrithl

-I,

'NCH FRONT PROFILE AT 20 SECONDS

0. 3MS TRIESTEP [ EXPLODED VIEW
TRAC-PF1/MOD1 013.0 (bO3e)

800

r~i

m

N

U' to

0

Li.)

-U

N

I

-4

C,

N

(0

N

C,'

750 HIGH TEMPS

700 LOW FLOWS

650

600

550

500

450

4001L

0 2 4 6 8

REACTOR TIME

10 12 14

rSECONDS

16 18 20

FIGURE 25

-I,

;o

mn

ROD SU1WACE TBP4ERATURES AT 5 IELEVATIONS, FOR 3 CALUT.ATIONS

CALCS ARBE: CONT=O.25NMISHORT=O.25MM + O.3M3 DT,LONG-O.1MM + O.3MS DT

DISTRIBLITIC0

PWR HEAT TRANSFER AND HYDRAULICS 'VrRKING GROUJP

VxM.MWE Coney

Dr C A Cooper

Dr G RKiuber

Mr I Brittain

I& K G Pearson

Dr D B Litton

Dr P AW Batby

bt D Cbucili

Mr P C Hall

Dr P R Farmer

Dr L Daniels

Mr D K 7Ing

Mr B Chojnowski

Dr M El-Shanawany

Mr K Meyer

Mr P Light foot

W1T2

W1C

wre

N&C

NNC

cEXB

CEGB

CEGB

TBC ILeatherhead

342/B41

233 /A32

201/A32

338/B41

Booths Hall, OCielford Road, I'nuts ford, Cheshire

WA16 9OZ

Booths Hall, Qielford Poad, Knutsford, Cheshire

Springfields Works, Saiwick, Preston, Lancs

TDV, Barnett Way, Barnwood, Gloucestershire

National Power, Romn C1016, Wourtenay House, 18

Warwick Lane london BZ4P 3EB

Bldg 392, AERE Harwell

CLF 125, SRD Cuicheth, Warrington

Marchwood Enigineering Labs, Marchwood, Southanipton

National Power, Nuclear Safety Branch

Oourtenay House, 18 Warwick Lane, London EC4P 3EB

PPG, Booths Hall, Qielford Ibad, Knutsford, Cheshire

PPG, Booths Hall, Chelford Road, Knutsford, Cheshire

WV=, DORCHESTER, DORSE~r, fIT2 8DH

Dr I

Dr C

Dr A

Dr M

Dr A

Mr R

Dr W

Dr B

Mr G

Mr P

H Gibson

Richards

T D Butlarx5

K Denham

J Wickett

0 'Mahoney

M.Bryce

J Holmes

Ward

I R. Jones

236/A32

260/A32

343 /B41

336/B41

266/A32

263/A32

262 /A32

264/A32

239/A32

231 /A32

(Suimmary Only)

(10 co~pies)

39 AEEW -M 2590
MC BOO"H'S FALL, CHELFORD ROAD, KNUTFIORD, CH1ESHIRE

Library

Mr G FA Bland

Wr C J Drayner

Mr R W Clarke

Wr K T Routledge

Dr D B Newland

Mr N Hobson

CEGB, ITD, 18 1kRWICK LANE, LONDONT, EC4P 3EB

Dr G Ahmed

CEGB, PPG, BOOMS HAIM, CHELFORD ROAD, ioanUsFt)RD, CHESHIRE

Dr N Buttery

UKAEA, RISLEY NPUE, WARING~TtO, 1VA3 6AT

Reports Library

UKAEA, AERE HARUMEL, DID=C~, OXON, OX11 ORA

Dr D Hicks G35, Bldg 77

Dr G F Hewitt TH 392.7

CEGB, TEC, KELVIN AVENUE, LEATHERHEAD, SURREY

Dr G C Gardner (Suimffary Only)

Dr S Ridzinski

CEGB, BERKELEY NUCLEAR LABORATORIES, BERKELEY, GLWUCES7TRSHIRE, GL13 9PB

Dr D Withringtcn (Suminary Only)

CEGB, NATIONAL POWER, NUCLEAR SAFErY BRANCH, 18 MMRWCK LANE

LONDON, EC4 4EB

Mr J R. Harrison

UNIVERSITIY OF STRATHCLYDE, JAMES VEIR BUILDING, 75 MDMT40SE STREE, GLASGOW

Gl lXJ

Mr T Callander

Mr W Denpster

RR&A, PO BOX 31, DERBY, DE2 8SW

Mr M Mclntyre

AEEW - M 2590 440
NUCLEAR. INSTALIATIOtNS INSPECM~RATE, HEALTH AND SAFETY EDCEOTIVE, ST PETERS

HOUSE, STANLEY PRECINCT, EALLIOL ROAD, BCXYIL, L20 3LZ

Mr C Potter

Mr J F Campbell

Dr S A Harbi scn

(summary onl~y)

(Sununary only)

IN'IERNTIOAL CODE ASSESSHE~T AND APPLICATIWNS PRO)GRAMMETh

USNRC, 5650 NICO1SWt LANE, ROCKVIIL, WMURUtAN 20814, USA

Wr G Rhee

Wr D Bessette

(2 copies)

E3GW (IIHO) INC, UNE, P0 BOX 1625, IDI&D FAILS, IDAHO 83415, USA

mr

Mr

Dr

R Shultz

R Hanson

G Wilson

(2 copies + single sided copy)

LOS ALAMDS NXTIOtNAL LABORATORY, PO BOX 1663, LOS AIAWOS, NEW MEXIC 87545, USA

Mr

Dr

Dr

Mr

N Scbnurr

R A Nelson

M Capiello

R Jenks

(2 copies)

PAUL sCHERRE INsmrIu, CH-5 303, W:RW2LINGEN~, SWITZERLAND

Dr S N Aksan

Dr P Coddington

JAERI, TOKAI-MJRA, NAKA-GUN, IBARAKI-KEN, 31911, JAPAN

Dr Y Murao,

Dr H Akixtoto

DERS, CEA, BP NO 6, F-92260, FONTENA-AUX-ROSES, FRANCE

Dr R. Pochard

LIBRARIES

Dou nrea y

Harwell

RislIey

Windscale

Winfrith

-2

-2

-4

-2

-4

41

AEEW -M 2590

NRC FORM 335 U.S. NUCLEAR REGULATORY COMMISSION 1. REPORT NUMBER

(2-89) lAsiged by NRC. Add Vol., Sig. Rv..

NRCM 1102.,p ddnu NVbr . n.

3201.3202 BIBLIOGRAPHIC DATA SHEET NUREGICA-007umes.If3nr

(See inlstrutitons on the feverse) AEREW-M25900

2. TITLE AND SUBTITLE I AE-29

Time Step and Mesh Size Dependencies in the Heat Conduction 3. DATE REPORT PUBLISHED

Solution of a Semi-Implicit, Finite Difference Scheme for MONTH YEAR

Transient Two-Phase Flow April 1992 4. FIN OR GRANT NUMBER

____ ___ ___ ___ ___ ___ ___ ____ ____ ___ ___ ___ ____ ___ ___ ___ A4682

5. AUTHORIS) 6. TYPE OF REPORT

R. O'Mahoney Technical

7. PERIOD COVERED (minciese Dares)

0 PER~FORMING ORGANAIZATIONJ - O.5AUC Atdf ArlflDC0C .... , .,. - ..

nsme end mailing address.) I - r, U - ,ufeo ein . ulwRgltr omsin n ~~n ~e;i ojw rvd

Winfrith Technology Centre

United Kingdom Atomic Energy Authority

Dorchester, Dorset, DT2 8DH

United Kingdom 9. SPONSORING ORGANIZATION - NAME AND ADDRESS (if NRC. type me as. e it contrac vo,~ provide NRC Division~ Offsce or Region. U.~ Nuclear Regulatory Commission.

and mailing address.)

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555

10. SUPPLEMENTAIRPY NOTES

11. ABSTRACT 1200 words or less)

This report examines, and establishes the causes of, previously identified time step

and mesh size dependencies. These dependencies were observed in the solution of a

coupled system of heat conductionand fluid flow equations as used in the TRAC-PF1/MQD1

computer code. The report shows that a significant time step size dependency can

arise in calculations of the quenching of a previously unwetted surface. The cause of

this dependency is shown to be the explicit evaluation, and subsequent smoothing, of

the term which couples the heat transfer and fluid flow equations. An axial mesh size

dependency is also identified, but this is very much smaller than the time step size

dependency. The report concludes that the time step size dependency represents a

potential limitation on the use of large time step sizes for the types of calculation

discussed. This limitation affects the present TRAC-PF1/MOD1 computer code and may

similarly affect other semi-implicit finite difference codes that employ similar

techniques. It is likely to be of greatest significance in codes where multi-step

techniques are used to allow the use of large time steps.

12. KEY WORDS/DESCR:PTORS (List wordsoriphrases that will wasist esearchers In locating the report.) 13. AVAILABILITY STATEMENT

Unl imi ted

Time Step, Mesh Size, Dependencies, Heat Conduction Solution, 14. SECURITY CLASSIFICATION

Transient Two-Phase Flow (This Pag)

Uncl assi fied

(This Report)

Uncl assi fied 15. NUMBER OF PAGES

16. PRICE

NRC FORM 3351249)


THlIS DOCUM2ENT WAS PRINTED USING RECYCLED PAPER

NUREG/IA-0073 TIME STEP AND MESH SIZE DEPENDENCIES IN THE HEAT CONDUCTION SOLUTION OF

A SEMI-IMPLICIT, FINlTE DIFFERENCE SCHEME FOR TRANSIENT TWO-PHASE FLOW

APRIL 199

UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, D.C. 20555

FIRST CLASS MAIL

POSTAGE AND FEES PAID

USNRC

PERMIT NO. G-67

OFFICIAL BUSINESS

PENALTY FOR PRIVATE USE, $300