UUV

VIKlibrary, B.W. Bldg

FEB 1 B65

technical <7tote Cfto. 301

ON THE FORMULATION AND NUMERICAL EVALUATION

OF A SET OF TWO-PHASE FLOW EQUATIONS

MODELLING THE C00LD0WN PROCESS

Stephen Jarvis, Jr.

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
THE NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards is a principal focal point in the Federal Government for assuring

maximum application of the physical and engineering sciences to the advancement of technology in industry and commerce. Its responsibilities include development and maintenance of the national stand- ards of measurement, and the provisions of means for making measurements consistent with those standards; determination of physical constants and properties of materials; development of methods

for testing materials, mechanisms, and structures, and making such tests as may be necessary, particularly for government agencies; cooperation in the establishment of standard practices for incorporation in codes and specifications; advisory service to government agencies on scientific and technical problems; invention and development of devices to serve special needs of the Government; assistance

to industry, business, and consumers in the development and acceptance of commercial standards and

simplified trade practice recommendations; administration of programs in cooperation with United

States business groups and standards organizations for the development of international standards of practice; and maintenance of a clearinghouse for the collection and dissemination of scientific, technical, and engineering information. The scope of the Bureau's activities is suggested in the following

listing of its four Institutes and their organizational units. Institute for Basic Standards. Electricity. Metrology. Heat. Radiation Physics. Mechanics. Applied Mathematics. Atomic Physics. Physical Chemistry. Laboratory Astrophysics.* Radio Standards Laboratory: Radio Standards Physics; Radio Standards Engineering.** Office of Standard Ref- erence Data. Institute for Materials Research. Analytical Chemistry. Polymers. Metallurgy. Inorganic Materials. Reactor Radiations. Cryogenics.** Office of Standard Reference Materials. Central Radio Propagation Laboratory.** Ionosphere Research and Propagation. Troposphere

and Space Telecommunications. Radio Systems. Upper Atmosphere and Space Physics.

Institute for Applied Technology. Textiles and Apparel Technology Center. Building Research.

Industrial Equipment. Information Technology. Performance Test Development. Instrumentation.

Transport Systems. Office of Technical Services. Office of Weights and Measures. Office of Engineering Standards. Office of Industrial Services.

* NBS Group, Joint Institute for Laboratory Astrophysics at the University of Colorado. ** Located at Boulder, Colorado.
NATIONAL BUREAU OF STANDARDS

technical ^©te 301

Issued January 4, 1965

ON THE FORMULATION AND NUMERICAL EVALUATION OF A SET OF

TWO-PHASE FLOW EQUATIONS MODELLING THE C00LD0WN PROCESS

Stephen Jarvis, Jr. National Bureau of Standards

Boulder, Colorado

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U. S. Government Printing Office Washington. D.C. 20402

Price: 5 5 cents

CONTENTS

Page

Abstract 1

1. Introduction 1

2. Basic Equations (2V2T) k

3- The Reduced System (1V2T) 10

k. Auxiliary Routines 23

5- Controls and Blocks 29

6. Computation Summary 31

7- Glossary 36

8. References kO

Appendix - The Finite Difference Method of Lax 4l

Figure 1---- U2

Figure 2 k-3 Figure 3-------- hk

Figure k— V?

Figure 5

^-6 Figure 6 k'J Figure 7------ k8

111

On the Formulation and Numerical Evaluation of a Set of

Two-Phase Flow Equations Modelling the Cooldown Process

Stephen Jarvis, Jr.

A model of transient two -phase flow in a pipe is con- structed in Eulerian coordinates assuming a single velocity

but independent temperatures in the two phases . Experiments

on the numerical integration of the system for Cooldown

problem by both lax and Courant-Isaacson-Rees methods indicate that a very fine spatial difference net must be used to

compensate for the numerical diffusion essential to computational stability if a second surge is to be realized.

1. Introduction

This report is intended to summarize progress to date on the formulation and numerical integration of systems of equations modelling general

two-phase flows in a pipe. Principal application would be to the study

of the development of choked flows in tubes of varying cross -section, and

to the Cooldown problem. The latter problem, the cooling of warm transfer

lines resulting from the flow of cryogenic liquids, has been the target

of the present study.

In a private communication, a set of Cooldown equations was proposed

based on conservation laws for a single fluid, and three distinct phases

were admitted: all gas, all (compressible) liquid, and a two -phase medium

based on a quality weighted average of gas and liquid at a single velocity

(V), pressure (p), and temperature (t) . State changes in the two-phase

domain were confined to the vapor -pressure curve with quality (y) satisfying < y < 1. The complicated thermodynamic structure of the fluid

(N_) was approximated by least squares fit to empirical curves. The

numerical approach taken, Picone's method, involved polynomial approximation in x and t using the governing equations to determine coefficients.

This method is not well suited, however, to hyperbolic systems of partial

differential equations (PDE^) whose solutions may have discontinuous

derivatives at interior points (and in the case of non-linear equations, even discontinuous variables, or shocks), although it is applicable to

elliptic and parabolic systems for which a high degree of smoothness in

the solutions is the rule . Attempts to force polynomial fits to
discontinuities results in highly unstable "behavior, which was, indeed,

observed. Finite difference schemes are generally more successful with

hyperbolic systems

In addition to the restrictive assumption of homogeneity on the twophase medium in the above model the location of phase boundaries is a serious problem. If Ax is the spatial difference interval, it is useless

to attempt to resolve phase domains on a scale even as small as Ax. A

reasonable value of Ax, based on computing time, is 10 feet for the

Cooldown problem, so that it is clear that if the flow is not insensitive

to the structure' of phase domains of this scale, the problem is computationally impossible. Clearly the small scale spontaneous generation

and degeneration of two -phase pockets is beyond the scope of the numerical

approach in general, and broad scale averages are at all times under

consideration

A fundamental property of continuous hyperbolic systems is that the

flow at a point (x,t) is completely determined by the flow at time t - At

on the interval (x - XAt, x + XAt) (where At is small), the so-called

"domain of dependence". Here X is the highest signal speed in the fluid

at (x,t) . For given Ax then, the smaller is X, the shorter is the

computation time required to integrate over a given real time interval.

Since liquids in general have very high signal speeds, it is obvious

that compressibility will be a major factor in the computing economics

Furthermore, few two-phase problems of interests will involve "hammer"

(liquid-shock), so that considering the liquid to be incompressible

would be both reasonable and economical. A further simplification that can be made which will probably not be more significant than those already

discussed is to assume the gas to be a perfect gas and approximate the

vapor-pressure curve by a straight line to reduce computation.

The flow of incompressible liquid (having infinite signal velocity)

can only be handled in the finite difference scheme by considering the

dynamics of complete liquid columns bounded at each end by compressible

media. Thus the problem of scale length and phase boundaries again appears
Since we are in any case restricted to the consideration of twophase flows which do not depend critically on "phase resolution, " it

appears that the whole matter could "be greatly simplified by allowing

only two -phase flow everywhere, 0< Ym^Y ^ Y* < / , placing then a minimum and a maximum on the quality y(x,t) . A small volume of gas

then, while storing insignificant quantities of energy and momentum,

adds compressibility to the "all-liquid" flow, Y - Yi • A. suffi - ciently small quantity of liquid can similarly have no effect on the

"all-gas " flow, y

- /u . Taking this approach, we must generalize our models now to admit

both a liquid and a gas temperature to avoid being confined to the vapor

pressure curve. This extension appears to be a worthwhile generalization

in that heat transfer between gas and liquid is controllable, while supercooling and super-heating are admissible. A system of five first-order

partial differential equations (system "1V2T") is needed to describe the

flow, involving a single velocity V"(x,t), pressure p(x,t), gas density

r(x,t), gas temperatures ><, » *, */

, and liquid temperatures /^ ( *, v. (Longitudinal heat conduction in the wall is ignored also, so that the

wall temperature is controlled by a single first -order equation in time

coupled to the five mentioned above only through source terms and is of no mathematical consequence)

For further generality in a one -dimensional model a 6th equation

was added (forming system "2V2T") to allow both a gas velocity V"(x,t)

and a liquid velocity U(x,t), and requiring a description of momentum

coupling between the two fluids. But unless the difference I ^ ' W is

adequately restricted (by friction, for example), the initial value

problem for this system is not "well-posed", i.e., there are fundamental

hydrodynamic stability problems involved. (This is evidenced by complex

characteristics). This 6th order system was dropped, then, in favor of

the 5th order (one velocity) system pending a study of the meaning of

this instability.

In any fluid system, two natural coordinate systems are available:

the system fixed in space (fixed in the pipe) is called the Eulerian
system, while a system fixed in the fluid is called Lagrangian. Only in one -dimensional transient flows does the Lagrangian approach offer some advantages over the Eulerian, although the fixed boundary conditions (at

the pipe ends) of the Eulerian scheme go over into variable boundary con- ditions . From the coding point of view, the finite difference method in

Lagrangian coordinates has the great disadvantage that the difference net

shifts in and out of the pipe and may get too dense or too sparse to admit

good approximation to the dynamical equations. There is involved frequent

interpolation and renumbering of net points . Furthermore, in the (2V2T)

system, the occurrence of two velocities complicates the Lagrangian system,

(and might eliminate the natural advantage of the Lagrangian system)

On the other hand, the Lagrangian schemes would have minimized the

difficulty caused by the great variations in density and momentum density

that the two phase system presents, and might have described the movement

of "liquid points" without the large mass diffusion which turns out to be

the most serious problem of the Eulerian approach.

Experimentation with a Lagrangian formulation would be very desirable. It would involve no new principles, but some tedious coding

difficulties

2. Basic Equations (2V2T)

We shall assume the pipe has a variable cross-sectional area A(x)

shared between a gas stream of cross -sectional area zA and an incompressible liquid stream of cross -sectional area (l-z)A, so that z(x,t) is a volumetric quality.

Let s be the constant liquid density, and E and F be specific energies of gas and liquid.

In terms of source strength/unit length, S., the conservation laws

J are readily written for the two streams

Liquid Continuity,

(A(i-b)s) + (A(i-i)s u) ~ S a)
Gas Continuity,

(Air) f (AzrV) x - Sa

Liquid Momentum,

(A0-*)sU) t +(A0-i)sU> + AO-'Jyfl,

Gas Momentum,

(AirV)jt

* (AzrV* + Azfl K

Liquid Energy,

[A('-*)*(iU\F)]

+[A0-*UU(iut *F)+A(i-JOf,] ~ Sr

Gas Energy,

[A*r(J V

3 tE)}

+ [AitV(i V\E) + AzV+]- S<

Wall Temperature,

(2)

(3)

00

(5)

(6)

T - S (7)

y, t l

We assume a state equation for the gas,
Calling _e^ and ^c- the mass rates of evaporation and condensation per

unit length of pipe at (x,t), we have

o = ^ - *> and (9) ^ a. Let Jc £e Jc he perimeters of contact "between gas -wall, liquidwall, and gas -liquid, respectively, at the section (x,t), and -? t^ "X the associated shear stresses (wall on gas, wall on liquid, liquid on

gas). Assume gravitational acceleration Q(x) in positive - x direction given. The scl. and -<z, are Von Neuman - Richtmyer non-linear

pseudo -viscous forces introduced to permit shocks, and are described

later. Then we set:

S 3 , - A s (i- z) y

- r£ 2£ + r i.

(io)

and

S,r &<-*%- ^

£

i - r< £< (id

f.* V - /* U - A 2 -f^t

Let Q Q o be heat transfer rate from wall-into-gas and liquid-intogas per unit area, and Q the wall-liquid rate. Then:

- j£

(12)

and

V *" Vf A +fA - 1,-*,
We have taken the following point of view in formulating these sources

When a quantity of liquid JL A*- ±s evaporated, its specific momentum

and specific enthalpy are "brought to those values in the gas stream

at the expense of the liquid stream, and contrariwise for condensation.

For example, the gas momentum source term is -£ V- ^c LJ rather than

jt-U-ycj V ' C^11 a continuous three-dimensional model, the phase change

takes place on an interface on which U - V >a ~ ' £ . The option

here arises from the fact that in the one -dimensional model, the dependent variables are averages over portions of stream cross -sections .) This

has worked out better than the alternative because the mass -£- A t is

normally taken to be due to some part of the energy transferred into the

liquid in time At. Hence, in a small finite interval At, there is less

danger of removing more energy from a stream than is available. Thus,

S--h„(tt

*.+ 1,',)/A (ik)

A considerable effort was expended in integrating this system by

[3] the Lax method (Appendix) although it fails to fit the model of

"Conservation Form" by the presence of the term -p (Az), in (k) and

the analogous term in (5). The instability which plagued this approach

was probably not due to this fact, however, but was due to the fact the

system may become non-hyperbolic, i.e., it may have complex -valued

signal speeds, as is shown below. Since the Lax finite difference form

is subject to the "domain of dependence" restrictions associated with

the maximum signal speed in the medium, a calculation of signal speeds was made. By expanding the (2V2T) system to a quasi-linear system of

the form

a . u + 6, a. =0

\j tit <</ h*

[2] the signal speeds X are found as solutions of the system:
Following is a 6th order algebraic equation in \ whose coefficients are functions of the flow variables (but not their derivatives):

o = (A - u) (\ - v)\(X-of[^ - v - (A - v)

3

]

-^(x-v)[^v + »Qi- v)]

(15)

where

CO =

S 2

*>*<

dr

(cO > c)

(~l > o)

(v < o)

(16)

In general, £-<-> will be small, since F> .OS will be taken, and will

vanish in the case of all gas where 7 - / . When to is small, and

O i=- V } we find the roots

A =

where

94.

\

V

u

Uf U) v*

<f> + O^)

(U- v)((~.-v)V + * °)

'/?

(17)

(18)

8
Then as U -> V

4

,^-c- v v (u- v)

which is imaginary for

V(u- v) < o

For large <-0 ,

\ - < o

V

V (i- "Vj

(19)

V _ y-

^K/ all real#

On the other hand, if there is sufficient friction between the streams

so that | V - U| is not large at interior points of the (x,t) domain (as

must be expected for small <-o ), the above limits do not apply. Let us

assume the friction is large enough that

U ~ V + * CO (20)

where 0~ is bounded as <-*} -* ° . We obtain the roots

f V

u

A - (21)

\

u + nJ

7Ul-») (i±\f^7 sU-~ \J &(*$

which are all real for small co at least . No attempt was made to

determine how stringent is the restriction on the interface friction

and initial values of the functions imposed by the requirement of real

roots, and a study of this (2V2T) -system was postponed in favor of the
(1V2T) -system. (improperly posed initial value problems, as would "be one with complex roots, occur in unstable problems, as for example the

flow following the collision of two liquid masses).

3 . The Reduced System (1V2T)

Henceforth, we shall discuss only the system with a single velocity

V obtained by adding (3) and (k) and setting (J= V throughout. We

obtain an equation for the conservation of total momentum which is

(A

{

V) t (A(V\ A4) r

4Ajt - S, (22)

f

5 I- i + f (l- g) (23)

^3 " A

f f

- TJ 4 - *j £j

AQ- *)*.„- A *~ - S * S

(sk)

Tv~ n

*?i 3( 3i

The resulting system (1V2T)., is entirely in "Conservation Form" (A, a given function: see Appendix) and was programmed in the lax scheme.

Difficulties are encountered when area variation occfurs due to the

structure of the Lax Derivatives . An alternate (lV2T)_. system, also in

"Conservation Form" is obtained by removing A(x) from the time derivative,

Let B(x) = /A QO; then

ZGi>- A* /a - -B\/B (25)

z

t -((/-») y) -o-*jv2 -- S,b/s (26)

(fi) + (rzV) + nVZ -<$ B (27)

J t J X

10
(28)

(r Z (^V\E))

t+

(^V(iV',E,^)l ^

\i-di Vr F))

t +(b-i)V(j v'*f+ i^))

(30)

This was programmed in the lax scheme (using ^ = ^-= °) • The result was entirely unsatisfactory due to the essential averaging process implicit

in the system. It appeared that the CIR (Courant-Isaacson-Rees) method

would avoid much of this spreading, and could yet handle shocks by inclusion

of the concentrated pseudo -viscous diffusion {^i-c, ^p) • These terms

must he considered as source terms, however, and not a part of the signalgenerating structure basic to the CIR methods . The presence of the diffusion terms required the inclusion of an upper bound on the quantity

^/fcj *) based on stability in parabolic systems. Only when

strong shocks are present is this term comparable to the maximum signal

speed limit. The term used was f • / I ] /shock diffusion!

'?!]'[32t*IV(x*Ax)- V(x)\) I condltion

J

A t s< { (31)

t'l M \ M/

/domain of dependence!

"| condition

11
The equations treated by the CIR method are the quasi -linear forms

derived from (lV2T)>_, and called (lV2T)_:

F f Vz x - ((-*) 1/ = D (32)

** * ^. * f Vx - D, (33)

Vt * 1/1/ + i(PE M

+ Pr )- b (3*0

(35)

f

jt > x ^f l>£ JX^ >x/

where for convenience we have "broken the sources into the following

forms:

(36)

D. - 5

u

* D.

a kffJ* Pt3 (efiJ (37)

12
D„ - V2 0-*)

D «- b

(3 - i

K-~ rV 7

»„-- 33 2

* X A3/ " 7" KVj f T*0* 4

/)33 K^3

&,- £(?.", - ?,*,-^-<^*2

* 7 (n * , * *, x* + ^

J

+i

- [z-.f- i+)/ r

D+ 3 r

a * 4" ff^^^-I.O-yr-^^

r V

+ Y (Cr*** * **«* A)

D

4l

« <rY/-*J

-

/

v.

vT3 J^/-*J (r- E - ^

(38)

13
where

v- Bi.

Sjl.

(39)

(ho)

T

*'"(**)'(&)' (&2

< O

> o

(M)

"ft*

afso

<^m "2V8>j

(te)

where a is a number between 1 and 2 (in the case of air at standard

conditions, this has the effect of spreading steady shocks of moderate

strength over three to four net intervals)

14
This system is of the form

« b

L t U* - D. 6 - /,•-,

^

(^3)

(summed on repeated index) with

fe

l > -

/- V o c'

3 V c"

o c" V

/c*' c

41 C

HS

c

fl

V. C n

c

r?

o

\

o

o o

^- 3 + C o

C*t

o

c

rt

\/

j

m

r

c'

3 - - 0-*)

c"- r/i

C* a - -

v/

C"- A>

C*'- -•%

C* J = ^-/>ca/r/- ;J

c

f/ . - \A£ /f»C»-

^

c"< = -C^fjV,c"

c"- M

^r4 /. /" 1 I / ^- s? f

c

Ct5)

V.

15
We shall assume the perfect gas thermodynamics, which gives

where \ is the ratio of specific heats . The signal velocities now

satisfy the 5th order system:

f(A- vT - o

7

V- P

J P*

While it has not been proved that the cubic has three real roots for

all realizable values of the flow variables, these roots were computed

from the explicit form for cubics, and no complex roots were encountered,

(the program stops in the event of complex roots). By Descartes' Rule of

Signs, it is readily proved that, if all 5 roots are real, there are at

least k roots (X , ..., X^ ) having the same sign as V, and at most 1,

(5) X , with the same sign as (-V) . In the case of "supersonic" flow, all

X have the same sign as V. As y —» <*>, only the first two of the cubic

terms persist:

(K-V)+-o , ana (A-vJ- "&-?(<-*} v

The cases Z -* / lead to the pure gas forms,

A

Cjj = V

and (48)

C*J i <rj A

w y ~ y ± fc+M

16
while the near-liquid case H--> o gives

A ® = —c WW 1/

and (49)

giving infinite signal speeds as 2 -=» ° . However, we shall in practice

take z > z min = . 05 • Then

< 3 6 yj (^+- Ivl)

giving a lower signal speed than pure gas alone (since *p- > 20, normally)

Thus with a reasonable restriction on z, z > 0.05, the limiting signal

velocity will normally be that in the gas

If all the eigenvalues (signal speeds) A of

l\

(W - b

1 *

I - o (50)

are real, we assume there exists a non-singular transformation satisfying

(p)

?

1 (^'-x (e) r*)-° (5i)

This has not been proved rigorously, but the following argument makes

it a reasonable assumption for this physical system. If the h x k

matrix obtained from (58) by deleting row 2 and column 5 is non-singular,

then of course the 5 x 5 matrix (58) will also be non-singular. This is

related to the fact that equations (32-35) of (1V2T) do not depend on F.

For the k x k system, however, a non-singular 2 * exists in general when

the eigenvalues are distinct, which is "almost always." A small

17
independent change in V could get one out of the special case, however,

yet cannot he expected to change the mathematical structure of the system.

Equation 51 defines a set of left eigenvectors £ corresponding to

X , (p = 1, . . .,5) . Then applying this transformation to (28), we get

or (52)

[3] This is the so-called canonical form for hyperbolic systems 3 in which

each equation (index p) has derivatives in a single direction in the

(x, t ) plane defined by the "characteristic lines" / :

r p

where V is a parameter along / Then

u - r \

<-J . + A <u t

>*

so that

JV

P

z

pt

t

J f

u =2

^p

1 & (p-i.-,<r)

^

is, for each p, an ordinary differential equation.

It is in the canonical form (^3) (using (V?)) that the CIE finite

difference scheme is analyzed. When the coefficients and initial data

satisfy certain smoothness conditions, then as A x -=» o f the finite difference scheme is

z z

pL \-^-

A t

+ £zpi A®[± u>Hx+A\t)~ ^(^))_ 7 /y

(5*0

id
. ( AJL1

where A t £ /v*^*\> \ [\\ j and. one uses ( + ) according to

Ik, A ]

l -

is stable and converges to a solution of (V3) in

some finite, Ax- independent time interval AT > of the initial line.

While theoretically valuable, the analysis does not give useful estimates

of AT or the error for Ax. Certainly the procedures could not be continued beyond the formation of the first discontinuity, and might fail before

that. The addition of the non-linear pseudo-viscous terms (Ul-^2) has in

practice permitted continuation beyond the formation of shocks, but

estimates of the error must be inferred from knowledge of the correct

solution or experimentation with different values of Ax. Because of the

large number of arithmetic operations involved, a minimum Ax(AT) exists

for a given digital word length below which accumulated round-off errors exceed the truncation error. For smaller values of Ax, all significance

in the result will be lost. In the problem at hand, however, it appears

that, with the 8-digit word available, the computing time involved for

small Ax is a more stringent restriction by far than the round-off error.

IZ

^ i j Taking Z7 to be the inverse of 2 * ,

*''< 2'' -

*'''' Z*> - f' (55)

and we return the system (52) to a form where each equation gives a single time derivative:

u' Z 2"V'Aw

'i?'' -D l 66)

where the summation is explicitly shown because of the multiple occurrence of the index p. The index p on U is to indicate that, in

a finite difference scheme according to (5*0 > one should use right- or

\(P) left-hand difference approximations according as A is negative or positive

19
Since

i~P

2

r /*"(V"- A

w f, )J? -

by (51); using (55) we obtain

ft""- \

v Dzu

- o .

C

(57)

so that Z

*Q)

is a right eigenvector corresponding to the eigenvalue

\W,(*) s= (1,...,5).

(3)

As V-» <?, there is exactly one root, X u , of the cubic (V7)

which vanishes . Write

then explicitly, we may choose:

/"

/ o-*J

-\

z^

•*-t/"

?z \j(A^U3lii M0-?)] -i(A

[3)-i)

r

A - (A^ A

(r)

J ^

a

V

(58)

-^

corresponding to eigenvalues

A

v " - V. I/.

>> > °V A

(4) A w (59)

20
The first approach to the finite difference numerical integration

was to use (58) and the computed inverse Z. in the finite difference

form of (56) . Difficulties were encountered due to widely disparate

orders of magnitude and "ill-conditioning" in (58) . (The computed product

generated diagonal elements near unity, but with some off-diagonal elements of the order 10 .) Rather than go to double precision, a satisfactory alternative was derived "based on a fact earlier mentioned so that

we may put

£-Ouvo (Af») - S-eyY^/ (V) (J-', ;*)

J

We may then rewrite (56) in the following form:

u

>t *

-h K

t3

(60)

-D c

where

oC =

P

and ( 2

to x<5)

U)

I

-/

Cv> o)

tV<c)

(61)

/ a

U)

v < o)

Or)

(62)

n l 21 ) are left-and-right eigenvectors corresponding

Specifically, we take

21
r

H - V

3+ ,. 4-1 11 , +* . (s) ? '*~~ ^c^l c

"

c ~- o"(c" *

S3 _ (c^-x^x^-v)

c

"hi

ss

Z = O

\

and

I

ir c

43. C

ti

ir

3S \

{s> - V

c

IS

i r

13 s~ 3f C" C

I

s~sr i

ir

(63)

(00

C C 1- c c

V.

c

,2 (^-v)

22
For the 'subsonic case , Q - ' , X* ' has multiplicity 1, so 2 and

-7 i(r) ^ •?/(*}

zt are orthogonal- Orthonormality gives the unknown scale jr : z

/(f) ~ I / £ E

W*I^, I (65)

^-. Auxiliary Routines

Primary emphasis so far has been on constructing a stable, controllable and flexible scheme "FLOV for handling two -phase flows based on

the hydrodynamic model (1V2T) . Auxiliary factors, such as boundary

conditions, initial conditions and sources terms, have not been intensively studied or varied except in so far as they have directly interacted

with our primary purpose. To the extent that they have not interacted,

simple representative types were taken for the auxiliary factors, and

elaborated only to a degree that seemed necessary.

As far as was possible, the auxiliary factors have been isolated

in subroutines to improve organization and facilitate changes

In addition to the basic program (FLO^) using the CIR method with

equations (60-64), three subroutines are used: BNDRY, INVC and KAPA.

Subroutine BNDRY is used to set values of the dependent variables

on the end points, so that only interior points are treated by FIDk.

In cases I and II, BNDRY has been based on the assumption of an infinite

reservoir at the left end, with source free, adiabatic flow from the .right

end at atmospheric pressure. Some studies were made with a finite tank

BNDRY to determine tank pressure - surge pressure coupling, but this

effect was determined not to be basic to the second surge formation.

Subroutine HWC computes and sets the initial values of the dependent variable (including A(x) determined from a functional form) at both

interior and end points

The subroutine KAPA computes for interior points the source strength

jz^ and ^o- and the source coefficients U' of (38) .> including the

friction terms ( T^ ~g ), the heat transfer rates ( *? %, « %i)}

the wetted perimeter per unit area ( K e K» K^ ) and the non-linear

23
pseudo -viscous stresses and work terms ( /*0, /y /(] ) via (kO-k-2)

The ( K e K a l<i ) describe the character of the flow, K indicating

the degree of dispersion of the media. They are prescribed functions of

the flow variables and geometry. In cases (i and II), we have taken

K£ = (j-z) P (M)

where P is the pipe perimeter,

K

- i P (67)

9

and

k - 3 «i (68)

This gives volume -proportionate wall wetting in gas and liquid streams.

Numerical experiments with annular flow, Kg = P^ K-~ &

}

K = {j~2) P , showed sharply curtailed At as a —=> O because of large heat transfer

into small quantities of liquid. Flow experiments indicate annular flow

is not the common mode in Gooldown, but that °'^l should be much

larger

Tke ( T a T. ) are taken as normal friction stresses:

r = r Vlvl f /?

7

(69)

t, - s V/v/ f

t

/t

where J - t, - GZ3 for Reynolds number vd. *d, z greater than

2000, and / /, « -^- . for Reynolds number below 2000. Here

214-
(70)

£ " * /w

^ ^-c, 1/ fi (x)

For the equation of state in the gas, we take

f

= A, r T (t\ r

o.3(lo)

7)

(71)

A

with gas and liquid energies,

(72)

F - Cvl (Tt - Tj

c =- o. fi (i°y

We take the following linear approximations to the pressures and latent

heat from data for N^ on the vapor pressure curve:

and

^a." K, '

/ ^^ 7"-

-4./3*(/0)'t &ZIUQ*)' T

(73)

and

\, 2, +£j T = zj?(fo) f - i^o(io) 7 T

Since A - (£"- FJ , we find £*^ « O. Q 1 1 (f o)

(Ik)

s,t.

and / c = 5 3

*

25
2

We take the heat transfer rates per cm per degree Kelvin from wall * * * to gas Q , from wall to liquid 9 £

, and from liquid to gas q to ^" h

fix 1 '"* "be a simple form of the Dittus-Boelter correlation:

* * * '

i — °-?

1]- %" ^JBC*) «.

%] - %V ^

t iBU) &I

(75)

1*. - <" d,i/86tl ^ (76)

and put

•/ A = o.i(io)

a (3^3

Then for heat sources we have

f"j*7- 5 Vj'fc'" V

<i %e - K* 1* (t„- 7-J

U„?. - *. %* Cr* - ^

(77)

Since we are using only "bulk heat transfers and mean temperatures,

no precise estimates of evaporation or condensation at the interface can

"be made. Such an estimate would depend on temperature gradients at the

interface, and the relative values of interface temperature and interface

saturation temperature. We distinguish three sources of evaporation and

condensation: ( jz, ^c ) at interface, ( ^a ^j ) at wall, ( -e^ •«•, ) spontaneous for supersaturated cases . Finally,

(*,, ^>) - ( **>,; **, ) + U^ ^a J * (^ <**) (78)

26
Then we have used the following forms

i O

/ K *^ '-fc-zj

{T<Tftt (f))

(79)

^-^A (Tj>T.tW)

I Q

- T.

(

z * Z«,J

A [Ts , t >rt

(80)

where

% t - srv^c. [ e- F j o.Qrhof] (81)

X.

<

o

<2 <t ll(Tl- 77,

J

/I

ft;* 77, J

(82)

fc > 7;,,

and

27
r o

^ u

v. V~ K

3 ^g l_Jjl

'^'

A

(83)

The 0~) are "efficiencies" of conversion of transferred heat to evaporation and condensation. Also,

r

X

3 =

o (* > 2 t?Ax)

r< * Z.t

r,o-*)h-vjA b:>r,j

m

and

^r

*i*(r„r -TjA

(85)

We have used CT = 0, \T f ^ - f

, ^ *tj S~) , and <r •= 0. 00/' for

(^ = 3

j 6). This latter was chosen simply because it gave quantities for

( ^ /cz

) of the same order as had been obtained in typical cases for ( -£

a ^c a ).

2 V

In the calculation of (4l, ^2) the calculation of <=r-" was modified

28
to

f

- I.

the minuend being chosen to remove the part of the derivative due to

variable geometry, as computed from the steady adiabatic flow equations

(A(i-*)vl x -a

(An- v) ,- o

J X

In the equation (7) for wall heat transfer, the heat transfer rate was taken to be

<A ^- -77^^ = ""^

(87)

where P„ is the perimeter of the uniform pipe, S its thictaess, ^

its mass density, and dpu_,its heat capacity.

5 . Controls and Blocks

This finite difference scheme, called FIDk, is finally then in the

form

m

l (*;tfAt)= <^

l (\t) * At(D

l (x,t)- E'fx, t)J (88)

At least three complete passes through the net of points in x are involved

in FL04 for a single time step. On the first pass, boundary values are set,

and (At) , satisfying (3l)> is obtained. We also require At to be less

than twice the value of At used on the previous time line . D and E are also computed. On the second pass, (u'J s. u ' (

j t f A t) is found at each

29
x, and tested to be certain that the following limits are satisfied:

( l2

r -2l « (OLl)- 1*1

\r+ - r / t (GLR)-lrl

[\J

+ -y\ ,< ^w((uuJy (uLuJ.|v/J (89)

W-Tj < (ullJ-/t-J

Tested for positivity is the argument of the logarithm in (92). If

any of these tests fail, we reset At = 2"At, and repeat the second pass

In the event that the computed value z falls outside the range

7 x 2: $ z? _ . . , we set z + to the nearer of these two limits, and

use equation (32) in reverse to solve for the --o or -£. (positivevalued) required. This latter value is used in the remaining equations.

In case I and II, we have used

U L2 - OL&- UUO - (jLt- ULL-- G. I ULV= <3. / Uo)

-to (90)

In the third pass, the new values of the dependent variable

(z+ (t,x), etc.) are stored as old values (z(t,x), etc.) in preparation

for computing a new line. BWDRY is used here.

In view of the fact that 1 - z must be permitted to become very

we replaced

- <r small, we have used /— £ = C in all algebraic expressions, and

A 2 - F At (91)

by its exact equivalent

A c^= -J^ ( I - cf At F\ (92)

and similarly for x - derivatives

30
6. Computation Summary

Computations were carried out on an IBM 7090 with an 8-digit floating

point mantissa. The finite difference scheme "based on (60-6k) appears to be

stable and to give solutions of the proper character during the growth and

early stages of decay of the first surge, a real time of about 2 seconds

In Case I, a peak pressure of 17 atmospheres was calculated. For these

computations, Ax = 290 cm was taken (20 stations in the pipe length),

yielding a computing rate of 0.6^ flow seconds per minute of computing.

Since flow times of the order of 10 seconds or more are desired (l6 minutes

computing), it does not seem reasonable to reduce Ax significantly, especially since cutting Ax by two quadruples the computing time for a given

real time interval . After the first surge, liquid was flowing back into the tank, and gas was flowing rapidly from the down stream end of the pipe. The pressure

fell rapidly until it began to near the tank pressure. At this level,

the pressure no longer fell, but was sustained partly by boiling of a non-diminishing liquid residue in the pipe and by non-diminishing gas

energy and gas density. A quasi-steady state thus resulted, and no second

surge evidenced itself.

The reason that energy, pressure, and liquid phase are maintained

in spite of the rapid transport of these quantities out of the pipe at

each end can be seen by writing the equations (60) in the form

U

>t *

^

l ^

As-) = - x(3 Z A (93)

31
and

P

L ^ J (Case I) (9^a)

where we have approximated the difference in right- and left-hand derivatives by a second derivative within the accuracy of the finite difference

approach. The equations (93) a**e diffusion equations,

v tr) ^iW is) i -o<"AZ_ £ > O , the diffusion terms in the right-hand side

tending to smooth out irregularities in the functions and their derivatives

Under the smoothness conditions of the CIR proof, these terms vanish as

/S. X —> O , but remain large enough to stabilize the numerical method.

In the computation described above, however, these terms dominate the

convective terms on the left-hand side when the pressure has fallen to

nearly tank pressure. This diffusion also probably accounts in part for

the earlier surge peak times encountered. This numerical diffusion is,

of course, a physical interpretation of the effect of truncation involved

in the finite difference equations

A few attempts have been made to circumvent these difficulties, and

allow for the formation of multiple surges without reducing the stability

of the system. The only moderately successful one was to introduce

multipliers which tended to reduce the diffusion when it was already

small

/ fu:.).-(«U(

/ c \ I [f i \ I ( Case I3: F ) (9^)

l =Ft i * / , and t~ - I if the right- and left-hand

derivatives have opposite signs. Using this form (Case II) a small

second surge was observed, while the peak pressures in the first surge

reached only 11 atmospheres

Because so severe a flow change appears to result from adjustment of

purely computing parameters, it might on first sight seem to invalidate

32
"both results, "but the diffusion is unquestionably overvalued in each case. Any scheme which retains stability while reducing the diffusion should

give improved results . It might well be worth testing other ideas in

this regard. Two alternatives are (l) greatly reducing Ax, giving a practical computing method only for problems in a very small geometry

or a very short time, or (2) attempting to program a Iagrangian scheme,

in which the diffusion required for stability might be much less, and

better liquid control might be maintained. It does seem that because

of the high mass, energy and momentum density of the liquid, a good

representation of liquid movement is essential to a meaningful computing

scheme

It should he noted that in (93 )> *, the factor in the liquid

continuity (z) equation, can be set to zero, or even made negative without interfering with the stability of the scheme, so long as o < 2- < I is retained. In some tests here, further improvements have been made in

avoiding liquid diffusion, but results are not conclusive. Making F

dependent on flow-direction and o~x~' > has also been tried inconclusively.

Some realism might be gained here, in reducing the requirement of continuity in the liquid, and emphasizing transport. (indeed, for small z with droplets, z - waves cannot propagate.)

Cases I and II were computed with the initial values: (CGS system)

Ax <x $ J~ A x J~A x < x * Jo A X

V o a

c.S'oGs-(io) o. i o i J (l o) 1

a. «W / H (t o)~' a. (i 3. & ( i o)~*

c ot~ o. ?? ??r

7c* K. 3 oc ° K.

T£ 70°K. t,f

c K.

X 7 e ' K. 3oc* K.

33
with the single addition / u ( 4 A X^ oJ= ~J O K. to admit heat leakage.

Another variation was made on FI04 (using F = l) to test a model of

valve opening time T and to attempt to match computationally the strong

dependence of peak pressure on valve-opening time observed experimentally.

While not thoroughly understood, this dependence may "be due to liquid jet

formation at the valve . If this were so, the effect would he due then to

spatial variation which cannot be expected to come out of a one -dimensional

analysis . We used

B(*,t). BM =p C° « * < tJ

B(x,t)'~ S(x) (t ±T) (95)

An increase of peak pressure did occur as opening time increased,

hut this trend did not reverse itself for large opening times as is

observed experimentally. Again the liquid diffusion dominated transport

when the valve opened slowly; the liquid that was thus introduced downstream of the valve boiled, and the valve restrictions prevented the

required back flow which would have prevented the sky-rocketing pressure

at 10 seconds. Furthermore, the valve was represented by a transient

area variation which was not small and stationary as was assumed in the

derivation of the equations

Again, finer difference nets would have permitted better valve

representation without the large area variations introduced here:

Valve Opening Time (Sec) Peak Pressure (Atmosphere) Time

o.i 13.8 0.30

1 15.9 0.U9

3 19.3 0.59

5 39.0 0.43

10 ? (See figures 3> ^> and 5«)

In Figures 1 through 5 > "the pressure -ft , volumetric quality z, and

velocity V are plotted at each of the twenty pipe stations for times given

by the baselines (arrows on time scale) . Above each baseline, h is

measured in atmospheres (absolute )/o.l inch, z in units/inch, V in 1000

3^
cm/sec/O.l inch.

Figure 6 shows the pressure in atmospheres for Case I (lower) and

Case II (upper) as a function of time for station 7 (x = 7 A x), the

location of the pressure peak in each case. The circles are at the points

in time for which results were printed, less frequent in Case I. Figure 7 is an experimental curve given in a private communication.

Note the spread in experimental data.

Clearly, the worst defects of the FLOU model are in time of first

peak and peak decay rates . While perhaps these could be improved "by

modifying the geometry B(x) and improving the empirical heat and momentum

transfer formulae, it is certain that the numerical diffusion is the major

culprit. With the present costs of high speed computing, it is not practical to use Ax small enough to evaluate the model's potential effectiveness

in predicting the Cooldown process.

ACKNOWLEDGEMENT The author wishes to express his gratitude to Messrs.

R. V. Smith and W. G. Steward for their cooperation and advice, and to

Mr. J. H. Wilson for many useful discussions and invaluable aid in

programming

35
7. Glossary

Physical quantities are measured in CGS units. Subscripts g,£,o

refer respectively to gas (or wall-to-gas), liquid (or vail-to-liquid)

and interface (or liquid-to-gas), and are replaced "below by the generic

subscript V . Units are shown in square brackets

a. .: coefficient matrix of time derivatives for quasilinear hyperbolic

system

A : gas constant in state equation for perfect gas. (eq.7l)

A(x) : cross -sectional area of pipe (A(x,t) used in study of effect of

valve -opening time) . [cm ]

b : coefficient matrix of x-derivatives for quasilinear hyperbolic system.

B(x): l/A(x). [cm

-2

] /c } ^o . : mass rates of condensation per unit length of pipe, (eq.78)

cm. fee.

[orgs ]

C : specific heat at constant volume. [' ~ -1 ^ j

I e^l^ 1 C : specific heat of pipe wall. I v c^'

D., D. . : source terms for (1V2T) . with partition, (eqs.37>38)

-C ,-C>.: mass rates of evaporation per unit length of pipe, (eq.78) I'^Zri^'c I

E: specific energy of gas. (eq.72) /

' *j J

1 F: specific energy of liquid, (eq.72) I *—~

j F : a non-linear diffusion coefficient, (eqs .9^-a,9^b) [l]

F*: typical time derivative, (eqs. 91*92 only)

f : friction factor, [l]

Q (x) : component of gravitational acceleration along pipe . / -

a I /ON l^JJL- 1 h : wall heat parameter, (eq.07) / —t~ / WX r ., <-~i<3 -*

h^, : heat conductivity. [ ~ ;t

k^: i,B. (eq.39) [**»'']

rt.Je.c-. <\. J

36
Jt^ ' perimeter of contact on pipe cross-section, (eqs.10-13) [cm]

/jo. : pressure. J-£_r J

/p. (t) : vapor pressure, (eq.73) I I

J

P(r,E): general equation of state for vapor, (eq.8)

P (x) , r : pipe perimeter, and its value on uniform section, [cm]

J

/ e^x ] q* : heat transfer rate per area per degree, (eqs -7^,76)1 Cma j cc

o^- /

q*-*: constants in Dittus-Boelter correlation. (eq.75,76) [l]

q : heat transfer rate per unit area, (eq.77) I ffr~:

T J ^n : pseudo -viscous stress. (eqs Al,42) / —+z—3

r: gas density. [ ~T 3 /

ftcL : Reynolds Number, (eq.70)

s: liquid density. / CfA *

J

S. : sources, (eqs. 9-1*0

Jy . : sources. (eqAo)

1

t: time, [sec]

T: valve -opening time, (eq.95) [sec]

T , (fa): saturation temperature, inverse of /h (T) . [ °K]

T : fluid temperature . [ °K]

T : wall temperature. [ °K]

w

T : "base temperature for eqs . 72. [ °K]

u : dependent variable vector: (z,r,V, E, F)

(u , )_, (u , )_: right- and left-hand x-derivatives

x 'x R' v 'x L

few. / U: liquid velocity. [

" Xe J

ULZ, UIR, ULU, ULV, ULG, ULL: limits, (eq.89)

V: gas velocity. [fzZ.i

37
x: coordinate along pipe axis, [cm]

y: quality (mass of gas per mass of fluid mixture), [l]

z: fraction of pipe cross-section containing gas. [l]

Z: pipe expansion rate, (eq.25) [cm ]

ZMAX'

ZMHT

: lDOUnds on admissible z. (eq. 90)

z , Z : matrices of left and right eigenvectors. (eqs .51, 57>63, ^4)'

Y : ratio of specific heats in gas phase, [l]

d : pipe wall thickness, [cm]

<3 : Kronecker delta

Ax: fixed space increment, [cm]

At: variable time increment. (eq.3l)[sec]

A : fluid signal speeds . / J^ J

\*: latent heat of vaporization. (eq.74) I ~ ?

' / ^r : viscosity. I'ZZiZc. J

^ : (eqs. 16,46) ["^.J

v : (eqs. 16,46) / e-^j

p : mean density, (eq. 3l) / "^—i

/

p : pipe wall density. /-5—S /

\ V I C m . J sr = - £^ (i-z) [1]

CT". : coefficients in formulae (79-83). [l]

T : shear stresses, (eq.69) / ^vT*~ /

to : (eq.l6) [l]

(2V2T) : system of conservation equations (l-8) allowing each fluid a

velocity and temperature

(1V2T). : system of conservation equations with one velocity obtained by

replacing (3,4) in (2V2T) by (23,24)

38
(1V2T) : system of modified conservation equations (26-30) with one B

velocity.

(lV2T) p

: system of quasilinear equations (32-36) with one velocity used in

FL(A numerical evaluations

39
8. References

1. R. Courant, E. Isaacson, M. Rees, "On the Solution of Non-linear

Hyperbolic Differential Equations by Finite Differences," CPAM V,

pp. 2^3-255 (1952).

2. R. D. Richtmyer, "Difference Methods for Initial Value Problems,"

Interscience, 1957*

3« G. Forsyth and W. Wasow, "Finite Difference Methods for Partial

Differential Equations," Wiley, i960.

k. F. W. Dittus, and L. N. K. Boelter, Univ. of Cal. Pubs., Eng. Vol. 2,

p. kk3, 1930.

kO
Appendix

The Finite Difference Method of lax.

Lax's method assumes a system of hyperbolic partial differential

equations in "Conservation law" form:

where, for our cases, (r y

<? J h ) are algebraic functions of the flow

variables (

2- f"j V^ T~ 7~^ ) . The computation is advanced in time

from a given line of values given at equally spaced (2 A x) values of x on a staggered net by substituting in (A-l) the following expressions:

A , f

\ ^ -f

l (*J iAt)-i[f

l (x+A *A^f

L(x-A*J)]

(A-2)

ji At

and

g

l"(x,tj- ?

L'f"* x,^- f (*-**- ti

(A-3)

2 A X

After computing a line, the flow variables are computed from the

f (t + At,x) . On theoretical grounds, we must have that At/ be less

than the maximum signal speed occurring in the flow at time t. The

averaging resulting from the average used in (A-2) is necessary to the

stability of the method, but causes too great a diffusion of flow

quantities in the Cooldown problem to be useful. There does not appear

to be a rigorous theory of convergence and stability for the lax method

comparable to that given for the CIR method (39)

kl



Figure I

to

1.2

FLOW STRUCTURE IN X AT

A SEQUENCE OF TIMES

CASE II: F' *l (SEE Eq. 94b)

2.8

2.6

2.4

2.2

2.0

1.6

1.4

0.8

0.7

0.6

0.5

0.4

0.3

0.2

Figure 2

m

<T> a

CD 'sT Z UJ

rr 0> 2

UJO

Z>

hCOcr O3

UJ

III cr UJ

-> h- UJ

_J co co < — >

u_O o

_i

u_ o

to

UJ< ) 2 d z O II UJ J- 3 UJ

L 1

toc r o

1- > s> O_• o N ^ o2 o

Q. -

C3 > c- C3

44

en

1.2

INFLUENCE OF VALVE OPENING

TIME ON FLOW STRUCTURE

F

1 = I, T = O.I sec (SEE Eqs. 94a, 95)

Figure 3

o

UJ

CLO

UJ

Z>O

uj a: > h-co 5 o

UJ o

a>

tri aUJ

UJ

UJ

COo

cu

to

J- i1 \

IT)c t o > > | O_ O M V o2 o

Q.

*.">

<3-

z>en

1.2 ^

INFLUENCE OF VALVE OPENING

TIME ON FLOW STRUCTURE

F

1 = I, T= 1.0 sec. (SEE Eqs. 94a, 95)

! rv

Figure 4

m

o>

o"

UJ Z cr. . IxJ D </) CL 1- CT O 3

UJ

UJ or UJ > H UJ

_J co U) <

*"' > £ o

lj_ o <l> O

_lu.

(/)

III o ( ) ^ CO z o II UJ

Z)

_J

UJ2

1-

M

m

en

1.2

1+6

INFLUENCE OF VALVE OPENING

TIME ON FLOW STRUCTURE

F

1 I, T=3.0 sec (SEE Eqs. 94a, 95)

\ . J

r

Figure 5

PRESSURE VERSUS TIME AT A

POINT 20 FEET FROM THE VALVE

case n F'gfe i (SEE Eq. 94b)

15

p 10

CASE I F

1 = I (SEE Eq. 94Q)

Figure 6

^7

TYPICAL EXPERIMENTAL CURVES OF PRESSURE

AND LIQUID VOLUME FLOW-RATE DURING FIRST SURGE

CYCLE FOR ONE DRIVING PRESSURE

c

*E\

o

en

ul

<

o

_1

LlUJ

O

>

fluid: liquid nitrogen

TRANSFER LINE: 200 ft. X 0.75 in.

0. D. x 0.67 in I.D. COPPER,

VACUUM INSULATED.

INLET DRIVING PRESSURE 62.5 psig,

TEMPERATURE 75.7° K, VALVE OPENING TIME: APPROX. 1/2 Sec

EXPERIMENTAL

COMPUTED

DRIVING PRESSURE

(62.5 psig)

Figure 7

1+8
\