AlllDE EM7T17

PUBLICATIONS

NATL INST OF STANDARDS & TECH R.I.C. NBSIR i A11102247917

Yeh, T. T/A computer code for gas-llquld QC100 .U56 NO.86-3414 1986 C.2 NBS-PUB-C A Computer Code for Gas-Liquid

Two-Phase Vortex Motions:

GLVM

All 10224791

T. T. Yeh

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Gaithersburg, MD 20899

July 1 986

QC

100

• U56 Veronautics and Space Administration

8R-^U1U Space Center

5 g Center, FL 32899

1986

C. 2

RESEARCH

INFORMATIOM CENTER NBSIR 86-3414

A COMPUTER CODE FOR GAS-LIQUID

TWO-PHASE VORTEX MOTIONS:

GLVM

T. T. Yeh

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Gaithersburg, MD 20899

July 1986

Prepared for: National Aeronautics and Space Administration

John F. Kennedy Space Center

Kennedy Space Center, FL 32899

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler, Director
er»»

, HDRA3^3fl ‘ rtorrAfrWCRw^

ft*'*'

' _

'

• >

' ,^q}i/oo,«Ao «o=i aaao flTnjqr/io^A'.it ^

;P MOfTOM X3TWOV 3eAH«W3V/T

•V JO 3

ft*T

.2r,y

•l.'n. .f3M0 ui^ ii9 '

|.'1i » '1 iUSitftS Ajt

^ J

-Qi) iniiiin^'

tt ? icq? ybtmifil #jeiw3-

_)» i«rn«iD

d I'

' U

m. » »

‘

i#^*iva v#Ajr..A iMmi .»o«A<iw<iTi ^0 y^ii'fUf
Table of Contents

Abstract ii

List of Notations iii

List of Figures v

I.- Introduction 1

II. Model Equations 2

III. Numerical Method 5

IV. Program Details 10

V. Summary 14

VI. Acknowledgments 15

VII. References 15

Appendix A. Code Listing 16

Exhibit A. A Sample Input 32

Exhibit B. A Sample Output 33

Abstract

A computer program aimed at the phase separation between gas and liquid at

zero gravity, induced by vortex motion, is developed. It utilizes an explicit

solution method for a set of equations describing rotating gas-liquid flows.

The vortex motion is established by a tangential fluid injection. A LaxWendroff two-step (McCormack's) numerical scheme is used. This program can be

used to study the fluid dynamical behavior of the rotational two-phase fluids in

a cylindrical tank. It provides a quick/easy sensitivity test on various

parameters and thus provides the guidance for the design and use of actual

physical systems for handling two-phase fluids.

Key Words: computer code; gas-liquid separation; numerical modeling; two-phase

vortex motions

ii
' j^J ' Ji' e ‘ r i b i«'n i

ttiStiKyir: c.;rr , p&eL ‘j^i^liillj^i

‘•L

') ^

’ >1

I

imjirr.i .Y -I i 4

bn*

^

’^.’4$ii«liili#‘4i

fii^

X-K
List of Notations

A , A , a ak

Added mass coefficients

*d’ *dk Drag coefficients

Body force density

‘^Ij’ Sk’ C^, Generalized coefficients

dk

Bubble diameter

^2 Liquid (droplet) diameter

Effective interfacial force density

n

a Exponent used for w

S K

Exponent used for

P Pressure

R Tank radius

Re P-V R/u^» Reynolds number

Jet opening, R^^ < r < R^^

Minimum radius considered in the numerical analysis

r Radial coordinate

t Time

V . J Averaged jet velocity

"n Gas radial velocity

(\l > Liquid radial velocity

V

s Velocity scale

1

<D 1 > Gas tangential velocity

iii
w

w

a

a

ak

dk

1

2

T

0

U

1

<P

oa

Liquid tangential velocity

Gas velocity

Liquid velocity

Weighting function for added mass coefficients

Weighting function for drag coefficients

Gas volume fraction

(1-a^), liquid volume fraction

Exponent for diameter variation

Circumferential coordinate

Gas dynamic viscosity

Liquid dynamic viscosity

u, + u, ^ total effective viscosity

k k

Turbulence or eddy viscosity

u /p , kinematic viscosity

i\ K

Gas density

Liquid density

Averaged density of k-phase

*

''a * “2‘’2^

Angular velocity

iv
List of Figures

Page

Figure 1. Two Step Difference Scheme 3^

(Backward Predictor-Forward Corrector Version)

Figure 2. Velocity Vector Distributions 35

a) Gas phase, b) Liquid phase

The annular region between two dashed circles

is the region of injection

Figure 3. Gas Volume Fraction Distributions 36

V

I. Introduction

Mechanical systems have been devised for producing artificial gravity

fields to spin-up liquids in containers. These involve rotating mechanisms that

are cumbersome and, more importantly, have moving parts that can wear out.

Here, liquid rotation created by fluid injection is considered. The detail

analysis of the two-phase vortex model can be found elsewhere [1]. In this

report, the details of the computer code are described.

The computer program was developed to study the fluid dynamical behavior of

two-phase fluids in a tank at zero gravity. The phase separation between gas

and liquid, induced by vortex motions, is of primary interest. The program

utilizes an explicit solution method for a set of equations describing rotating

gas-liquid flows. The vortex motion is established by a tangential fluid

injection. A Lax-Wendroff two-step (McCormack's) numerical scheme is used in

the computer program. This scheme uses a conservation form of a system of

equations together with an auto time step feature.

The program was developed and tested on an HP-1000 minicomputer. The HP1000 's FORTRAN 77 is based on the American National Standards Institute (ANSI)

77 standard programming language FORTRAN (ANSI X3. 9-1978). The HP FORTRAN 77

has extensions to provide a more structured approach to program development and

more flexibility in computing for scientific applications. It fully implements

the Military Standard Definition (MIL-STD-1 753) of extensions to the ANSI 77

standard. In order to make the computer code more useful for other computer

systems, modifications have been made so that the code is closer to the ANSI 77

standard and thus less system dependent. Some limited extensions are still kept

in order to produce the code in the HP-1000. Since the graphic routines are

system dependent and must be modified with their equivalents at each computing

facility, the original graphic code has not been included in this report. All
lines preceded by ”*V'’ are originally adopted to use the vector operation

package supplied by Hewlett Packard. The speed of the code can be increased by

replacing many ”do-loop" operations in the code with high speed vector

operations. Little effort is required to incorporate the vector operation into

the code if the vector operation package is in the system.

Thus with limited effort, the program can be adapted easily to most systems

accepting the ANSI 77 standard FORTRAN. For example, the EMA (Extended Memory

Area) statements may have to be removed from the code for some computers. Also

double precision real numbers (Real * 8) could be replaced by single precision

real numbers.

II . Model Equations

The vortex induced model is based on a two-phase, two-fluid continuum [2].

It incorporates several interactions between phases; namely fluid drag and

virtual mass effects and it can be modified to include additional interaction

effects. Detailed analysis of the model has been reported in Ref. 1. A brief

summary of the system of equations is given below.

The equations for the conservation of mass and momentum for the two fluid

two-phase model in an one-dimensional, axisymmetric case

1

9ra, 9ra, V

k k rk

3t 9r

0

2
2 ^ „

2 _ . ^ IH

3t 3r '^k 9k ^k pk

^ 3r

\ ^kl ^ 3r “il^092, ^

+ a, C,, r (V - V ^) k dk rl r2

3t 3r “k^rk^9k “k “2.^r9Z

^ ®9H

+ a C r (V - V ) k dk ^91 92^

for k = 1 and 2 and with

Si ° '“1V2 " ^

S2 ' '“l“2’’l

^ ^

Cfi = ia^p^ + A^)/<p >

S2 Si - ^a'<^ >

S2 '“iS " >

3
^dl

^

^d2

= >

and

<p > = a^a^p^p2 + A^ ^ <^2^2^

The effective stresses are modeled as

T , = 2y, 3V , /3r rrk k rk

a, ^ ’^Q 1 =1^1 ^ /r)/3r

r0k 9rk k 0k

X = 2 u V /r 90k ^k rk

with

e t ^k

= "k " \

and the interfacial forces are modeled in the form of

Ml - A^(?2 -

Vi) * ^

(V^ - ?i)

is the force density acting on the phase 1 by the phase 2. and are the

added mass and drag coefficients, respectively.

The incompressibility condition is reduced to * ^2^r2

”

Q is the net radial outflow,

r In the program = 0 is assumed, since the mixture pumped out Is Injected

immediately back into the tank at the nearby location. The net volume or mass

in the system is effectively unchanged except for the net change on the angular

momentum. Thus, the pump system (withdrawal and injection) acts as a body for-'e

4
on the mixture at the nozzle location. The net momentum gain is thus the

momentum introduced into the system minus the local momentum pumped out.

Therefore, we will model this pumping dynamic by body forces without considering

the mass transfer. That is, the body force density will be replaced by

a. p, V

.

k j — —

the net momentum gain, —

r

—

(V.n - V)

2ir j

injection speed.

at the nozzle location. where V

.

J is the

III . Numerical Method

The complete solution of the complicated system of equations can only be

obtained through numerical methods. An improved Lax-Wendroff , two-step scheme,

(also referred to as MacCormack’s method) [3, 4] is adopted for solving this

time-dependent problem. This non-centere<f differencing scheme, using a full

step backward prediction and forward correction version, requires no explicit

artificial viscosity if a proper stability condition is satisfied. Using this

technique for solving fluid flow problems is very efficient and has been in

widespread and successful use for some time. It is good both for the timeaccurate computation of steady and unsteady flow problems. The general features

of the scheme are: i) its explicitly conservative form, ii) it is a two-step

predictor-correction type, iii) it is three point, two level - that is, the

n+1 n solution of f^ at level n+1 depends only on three values of f^ at level n, and

iv) it is second-order accurate in time and in space.

For using the MacCormack's numerical technique, the system of equations can

be expressed in the conservative form as:

W = F + P + gG + S

t r r r Here the subscripts (t and r) denote partial differentiation with respect to t

and r, respectively, and W, F, P^, gG^ and S are column matrices with five

5
elements. All the components of F, P , gG and S can be regarded as functions

r r of the components of W which are the independent variables. The fundamental

theory of the MacCormack’s scheme is briefly given below.

For second order accuracy, the solution could be written as

= W® + At W

° W

°

,,o At ,, o ^ At 0 ^ o.

= W + — W. + :r- (W. + At W. . ) 2 t 2 t tt

where

and

= I

(W° + At W^°) + I

(W° + At

= I

(wP + W°)

+ At is the predicted value,

= W° + At W ^ is the corrected value.

The superscripts denote the time-level of the information and subscripts denote

the partial derivative with respect to either time t or space r. Specifically,

superscripts 0 and 1 are the initial and the completely advanced time (here two

steps) plane; p and c are the predicted (1st step) and corrected (2nd step) time

plane. Thus, W ^ is the time derivative of W evaluated at the initial time, and

is time derivative of W evaluated at the predicted time.

Fig. 1 shows the diagram of the two step difference scheme used in the

computer program. Due to the difference scheme, the spatial location after each

step in time is a half grid off from the original one. Thus, the spatial offset

6
which resulted from a backward predicting step will cancel with those of the

forward correcting step.

Numerically, the predicted values are

= i.

) ”i 2 ^”i-l/2 '^i+1/2^

where

= i (W.°^ . W.°) . [(F.° - F.°^) .

. At (3O

, 3^0^, ,

and the corrected value is evaluated at the predicted time place, that is at

“in/2-

w“

L W.° + At1 t

= W

° + — [(F ^

i Ar i+1/2 F ^ ) + i-1/2^

(g ^ + g

^ ) ^®i + l/2 *^i-l/2^

<h*i/2

-

+ — (S^ +3^^ )+AtP^ 2 ^^i+1/2 ^i-1/2'' ^i

P c Here pressure correction terms at half and full time

steps respectively. Thus, for each time step, the advance is carried out in two

7
steps: a full step backward predictor, and then a forward corrector. As

indicated in the diagram, the subscript i is the regular mesh spatial location

at which solution is to be advanced, i _+ 1 is the spatial location of regular

mesh points immediately to the right and left of the location i, i + 1/2 is the

location midway between i and i + 1 or between i - 1 and i at the predictor

plane. Thus, for each time step as the procedure advanced, the outermost data

points at the boundary are not updated through the numerical scheme. The values

at the boundary are to be given through some suitable boundary conditions. The

numerical procedure utilizes a uniformly preselected spatial mesh and variable

time increment. To avoid a singularity at the center of the core region, a

finite radius is used for the inner boundary. The tank radius R is the outer

boundary. The time step is determined at each time step to ensure numerical

stability [5]. For a finite grid size Ar, the maximum time step At is given by

At, = 1/[|C

dk rk /Ar + Ar

“2*'2 =. 2

>’

Where k = 1 and 2. The minimum At^^ (with some rounding off) is used for the

time step. Normally, the technique with the time step condition gives fairly

good numerical stability. However, in critical conditions numerical damping can

be added either for damping oscillations due to large gradients or for

accelerating the calculation by increasing the time step. A damping factor, D

thus was added in the program as

W.

1 ID

(1-D) + (W^_^ + D

8
1_ XD

where is the value obtained based on the two-step scheme, and W is the

value after the damping factor D is added. A typical value of D = 0.2 can be

used for debugging the program. If no damping factor is desired, D = 0 should

be used.

The computer program was written in a Fortran 77 based computer code. The

code will permit evaluation of the effects of various parameters which control

the fluid dynamical behavior. These include tank size, fluid properties, such

as density and viscosity, etc., characteristic gas bubble and liquid drops

sizes, and relative location of injection nozzles.

A sample input and its output are shown on Exhibits A and 3, respectively.

The initial conditions for the gas and liquid volume fractions are taken to be

25J gas and 75J liquid. These fractions are uniformly distributed over the

circular cross-section of the cylindrical tank. Initially both fluids are at

rest. Other parameters can be found in the sample input in Exhibit A. The

resulting velocity distributions and gas volume fraction as function of time for

the sample run are given in Figs. (2) and (3) respectively. The velocity

distributions are displayed along equally spaced rays at different times to

enable clear observation. These velocity vector fields indicate all flows are

primarily in angular rotation with gas phase tending to move inward and liquid

phase trying to move outward, as expected. As the result of these radial

movements, the volume fraction distribution is also changed with time. And as

expected, the gas volume fraction is increasing at the inner region and

decreasing at the outer region as shown in Figure 3. More detailed results have

been reported in Ref. 1 .

9
IV. Program Details

The complete computer code is listed in Appendix A. The code consists of a

main program, GLVM and several subroutines. It is written in subroutine form

such that each subroutine performs an individual task. Each logical part is

clearly isolated and it can be easily modified to reflect different modelings

for the interfacial forces. The interactive input mode with self-instruction is

used for easy parameter insertion. Many instructive internal documentations are

included in the program. In the following, each subroutine is listed with a

brief description of its major function.

1) GLVM, the Main Program.

*To initialize data and start the program: Logical unit to save data (LUS),

(Logical Unit is 1 for terminal, and 6 for printer), job identification

notes (NOTES) i data file name for saving data (NAMR), initial time (TO),

final time (TMAX) , time interval for data output (DTPRT), etc.

*To control the calling sequences to the other subroutines.

*To check the time step.

*To save, print (and plot) the output data.

*To obtain the predicted and corrected values in the two step, numerical

scheme.

*To impose boundary condition.

*To update the data, time, and step number for the new time step.

*To provide a shutdown procedure either in normal (e.g., t > t ) or max abnormal (e.g.. At is too small) conditions.

2) INIT, Initialization.

*To input the test parameters, initial conditions and set-up the initial

column matrix W.

10
Default values are provided for most of the parameters. The default values

are listed at each interactive input step. If the default value is acceptable,

a comma is inputted.

Some of the relevant symbols used are listed below:

ALMT Limit values of a^, ALMT(l) < < ALMT(2).

DAMP Numerical damping factor. Normally set to 0.

DEN1D2 Density ratio, P^_/p2

DO

Y

Base diameter, i.e., d, = d , a, k ok k

DS Density scale =

EVF effective eddy viscosity factor.

GAMMA Diameter exponent Y

IVTX Type of simple initial flow: 0 = at rest, 1 = simple rotation, 2 =

H.amrael-Oseen Vortex, 3 = G.I. Taylor Vortex. This is needed only when

there is no data file (NAMR) given for an initial condition.

IW Boundary condition at wall (for tangential component). 0 = free-slip

(no-skin friction), 1 = non-slip. Also, when I = 1, a factor of (1-

w

r)*^’^ was included on IVTX flow to simulate an initial power law

boundary layer.

MM Size of data arrays, MM > NG. MM appears in many subroutines.

MU dynamic viscosity for phase k.

MUEF (1 + EVF) * MU. Effective viscosity.

MU1D2

^l'^^2’ Viscosity ratio.

NA, ND n ,n , Exponents for weighting function for drag and added mass

3. D

coefficients, A and A,.

a d
NAMR Data file name for initial condition, if any.

NG Number of grid points used.

OMEGA 0), initial rotation speed, V. = or, when IVTX > 0.

y

PS

2

Pressure scale, DS * VS

QJ,VJ Injection flow rate and speed.

RE Reynolds number, p V R/uw 3 ^

RJ1,RJ2 Jet opening, RJl < r < RJ2

RPEAK Location of the peak speed of the initial vortex, if IVTX > 1

RTANK Tank radius = Length scale LS.

VPEAK Peak speed of the initial vortex, if IVTX > 1.

TJ1,TJ2 Injection time, TJl < t < TJ2.

TS Time scale = LS/VS

VS Velocity scale.

The format of the data file for the initial condition (if any) is a six

column and NG row data file, where NG is the number of grid points. The column

sequence is K, a^(K), ^91^^^’ '^92^^^' where K is the grid point

number, and the rest of the terms are the gas volume fraction, radial gas

velocity, radial liquid velocity, tangential gas velocity and tangential liquid

velocity respectively. The data format is free.

3) DELA, AA

To determine the fraction of grid size in which the injection is made.

0 £

AA _< 1. This is used to define the location of jet. The region of

injection could cover several full or fractions of grid sizes.

4) DERIVl, Derivative

To get the first derivative of a data array using a center difference sch- me

except the two end points in which three points near the boundary are used.

12
5) DGCOEF, Generalized Coefficients

To calculate the added mass, drag and all the generalized coefficients (A ,

di,

A, and C. .). This is the heart of the modeling,

d iJ

The effective coefficients are modeled as:

A = a

"d

= * ''d2“d2

^1 == + 202/(1 + 3a^))

= a^02Pj_/(ct2 + 2o^/(l + 302 ))

^dl

2

= 18

^d2

= 18 - 02 / 0 . 8 )^82 ^)

al

na,/ na na^ * = OL^ /(o^ + O

2

)

^a2

=

•

= 1 - w

al

'^dl

= nd-, nd nd^

= O

2

/(a^ + QI 2

)

^^d 2

' 1 “di

6 ) DWPDE, Partial Differential Equations.

To evaluate the values of the increments on the column matrix W from the

partial differential equations. This is the major part of the McCormack’s

scheme. In each complete time step this routine will have to be called twice

7) FNDDT, At

To determine the suitable time-step size.

8) FSOFW, Column matrices F and S.

To determine the convective matrix F and the source matrix S.

13
9) JET, Injection.

To determine the momentum source due to the jet injection.

10) SIZES

To determine the gas bubble and liquid droplet sizes. In the model the

sizes were modeled to be functions only of the volume fraction, i.e.

d, = d , a, . k ok k

Different models for size distributions could be easily adopted here.

1 1 ) TAUOFW

To determine the stress tensor t and its derivative.

12) UOFW

To convert the column matrix W into the physical independent variables,

such as a, V , V^.

r 0

V. Summary

A computer program aimed at the phase separation between gas and liquid at

zero gravity, induced by vortex motion, is developed. The vortex motion is

created by fluid injections. The computer program uses a FORTRAN 77 based code

and HP-1000 minicomputer. It is flexible and accepts various input parameters

for different flow conditions. Other interaction effects can also be added or

modified easily. This program can be used to study the fluid dynamical behavior

of the rotational two-phase fluids in a cylindrical tank. It provides a

quick/easy sensitivity test on various parameters and thus provides the guidance

for the design and use of actual physical systems for handling two-phase fluids.

1 4
VI. Acknowledgments

We would gratefully acknowledge the support received from NASA's Kennedy

Space Center to carry out the model-development work described herein. Specific

thanks go to Mr. Frank Howard whose involvements and inputs on this project have

been most helpful.

VII. Refer ences

1. Yeh, T. T. (1986) "A Vortex-Induced, Gas-Liquid Separation in a Cylindrical

Tank at Zero Gravity", NBSIR 86-3322.

2. Drew, D. A. (1983) "Mathematical Modeling of Two-Phase Flow", Annual Review

of Fluid Mechanics, Vol. 15.

3. MacCormack, R. W. (1969) "The Effect of Viscosity in Hypervelocity Impact

Cratering", AIAA Paper 69-354.

4. Chima, R'. V. and Johnson, G. M. (1985) "Ef/icient Solution of the Euler and

Navier-Stokes Equations with a Vectorized Multiple-Grid Algorithm", AIAA

Vol. 23, No. 1.

5. Peyret, R. and Taylor, T. D. (1983) "Computational Methods for Fluid Flow"

Springer-Verlag, New York.

15
Appendix A

Code Listing

iGLOM T=00004 IS ON CR T4 USING 00136 BLKS R=0000

0001

0002

0D03

0004

0005

0006

0007

0 0 08

0009

0010

0011

0012

0013

0014

0015

0016

0017

0018

0019

0020

0021

0022

0 0 23

0024

0025

0026

0 0 27

0028

0 0 29

0030

0031

0032

0 0 33

0034

0035

0036

0037

0038

0039

0040

0041

0042

0043

0044

0045

0046

0 0 47

0048

0049

0050

0051

0 0 52

0053

0054

0 0 55

0056

0 0 57

0058

FTN77

$EMA /DATA/ , /UWUU/ , /COEFF/ , /SOURCE/ , /FANDS/ , /TAU/

$FILES 1 .2 PROGRAM GLUM< ,99), <860425.1537)

CC

C

Cc

c

cc

c

c

c

c

cc

c

c

c

c

c

THIS, PROGRAM WAS DEVELOPED TO STUDY THE FLUID DYNAMICAL BEHAVIOR

OF A ROTATIONAL TWO-PHASE FLUIDS ( GAS/L lOUID ) IN A CYLINDRICAL TANK.

THE VORTEX MOTIONS ARE ESTABLISHED BY TANGENTIAL FLUID INJECTION.

THE PROGRAM WAS DEVELOPED ORIGINALLY BY T.T. YEH OF NBS

IT WAS BASED ON HP'S FORTRAN 77(ANSI 77+MIL-STD-l 753

)

WHEN WHO WHAT

8502XX TTY ZERO-G FUEL TRANSFER, START-UP STAGE.

LAX-WENDROFF 2-STEP SCHEME (FULL STEP PREDICTION+CORRECTION)

WITH NUMERICAL DAMPING FACTOR ( NORMALLY SET TO ZERO)

CONSERVATION FORM, VAR I ABLE < AUTO ) TIME STEP

REAL»8

INTERFACIAL FORCES: DRAG, ADDED MASS

PUMP CONDITION: MOMENTUM SOURCE BUT NO MASS SOURCE

850715 TTY GENERALIZED EQUATIONS AND COEFF. Cij

851018 TTY IN ANSI 77 STANDARD< WITH A LITTLE EXCEPTION FOR

TESTING IN HP-1000)

C »»«* INTERNAL SUBROUTINES *****

C DELA., DERIVl, DGCOEF, DWPDE , FNDDT , FSOF-W

,

C INIT, JET, SIZES, TAUOFW and UOFW

C MOST OF THE LIST OF NOTATIONS ARE GIVEN IN SUBROUTINE .. INI

T

C CHARACTER NAMR*16, N0TES»72

INTEGER I,IOS,J,JTIME(5) , K , MM , NPRT , NT

Z ,IW,LUP ,LUS,NG,NGM1 ,NGM2

PARAMETER <MM=101)

REALi*8 BA(MM,2) , DDT , DT , DTMAX , DTMIN , DTPRT , P ZERO

X ,DW1, T,TMAX,TPRT, VJT, VDR<2)

Y ,RJ1 ,RJ2,TJ1 ,TJ2,QJ,VJ

1 , ALMT,U, V, ALP ,P ,R , W , UP , WN , DW , RDP , RH , F,S

4 ,BR,BRH, RH0,MUEF,V18,NA,ND

5 , DO, GAMMA, DAMP, DR

6 ,TRR,TRA,TAA,RTRR,RTRA, C, CPA, CD

COMMON

Y /JETS/ RJl ,RJ2,TJ1 ,TJ2,QJ,VJ

Z /CONTP/ IU,LUP,LUS,NG,NGM1 ,DAMP,DR

1 /ALPLMT/ ALMT(2)

2 /COEFF/ C<MM,2,2) ,CPA(MM,2) ,CD(MM,2)

3 /DATA/ U(MM,2) ,V(MM,2) ,ALP(MM,2) ,P(MM) ,R(MM)

4 /DRAGl/ RH0(4) ,MUEF(2) ,V18(2) ,NA,ND

6 /FANDS/ F(MM,5) ,S<MM,5)

7 /SOURCE/ BR(MM) ,BRH(MM)

8 /TAU/ TRR(MM,2) ,TRA(MM,2) ,TAA<MM,2) ,RTRR(MM,2) ,RTRA(MM,2)

9 /WWUW/ W<MM,5) ,UP(MM,5) ,UN(MM,5) ,DU(MM,5) ,RDP(MM) ,RH(MM>

EQUIVALENCE (WN,BA>

C ****** RHO(l) < RH0(2) (i.e. PHASE-1=GAS, PHASE-2=L IQU

I

D

)

C DTPRT TIME STEP FOR PRINTOUTCAND PLOT)

C

16
0 0 59

0060

0061

0062

0063

0064

0065

0 0 66

0067

0068

0 0 69

0070

0071

0072

0073

0074

0075

0076

0 0 77

0078

0079

0080

0081

0082

0083

0084

0085

0086

0087

0 0 88

0089

0090

0091

0092

0093

0094

0095

0 0 96

0097

0098

0099

0100

0101

0102

0103

0104

0105

0106

0107

0108

0109

0110

0111

0112

01 13

0114

01 15

0116

01 17

0118

LUP=1 ! LU FOR PRINTING DEBUG DATA(=1 TERMINAL)

LUS=6 ! LU FOR STORING DATA<*6 PRINTER)

7 F0RMAT<2X,A,3I5)

8 F0RMAT(2X,A,3(1PE12.4)

)

URITE(1,7) 'Enter lu for sawing data. D.F.=',LUS

READd,*) LUS

C KEET JOB TIME FOR FUTURE REFFERENCE

C CALL EXEC(11,JTIME,JTIME(D)

IF (LUS .NE. 1 .AND. LUS . NE . 6) THEN

C ******** Define a file nane for stroing output *»»****»*«*»*«»»

URITEC 1

, ' (2A) ' ) 'Enter FILE NAME for sawing data.'

READ( 1

,

' (A) ' ) NAMR

LUS=90

OPEN(LUS,FILE=NAMR ,IOSTAT=IOS,STATUS='NEW' ,ERR=999)

ENDIF

URITE< 1

, ' (A) ' )

' Enter N0TES(<73 CHAR,) for the job'

READ( 1

, ' (A) ' ) NOTES

URITE(LUS, ' (3H1 ,A)') NOTES

C URITE(LUS, ' (514) ' ) JTIME

C To set-up the initial condition.

CALL INIT

NGM2=NG-2

VDR ( 1 >=2 . »MUEF ( 1 ) /DR*«2 ! For deter«ining tine step

UDR (2)*2 . »MUEF(2)/DR»»2

T = 0 . ! TMAX=5,

DTPRT=0 .2 URITE(1,S) 'Enter INITIAL and

READ(1,») T,TMAX

URITE(1,8) 'Enter TIME STEP f READ(1,«) DTPRT

DTMIN=1.0D-6 ! DTMAX=DTPRT

NT=0 ! TPRT=T

NPRT=0

PZERO=0.D0 ! DT=DTMIN

INITIAL TIME

FINAL TIMES. D . F . =

'

, T , TMAX

r output. D.F.=',DTPRT

SET MINIMUM TIME STEP

Ti«e step nuMber

Pressure at center core

10 NT=NT+1

CALL UOFW(W,R,NG)

PRINT OUT AT SELECTED TIME

IF(T .GE. TPRT .OR. T . GT . TMAX) THEN

IF(NT .GT. 1) THEN

NPRT=NPRT-M

TPRT*TPRT+DTPRT»( 1 . +DNINT( ( T-TPRT ) /DTPRT )

)

FORMAT( IHl ,2( A5, 14, A5, 1PE9 . 3)

)

URITE(LUP ,21 ) 'NP=' ,NPRT , 'T=' ,T, 'NT=' ,NT, '0T=' ,DT

WRITE (LUS, 21 ) 'NP=' , NPRT , ' T=

'

, T

, 'NT=' ,NT, 'DT*' , DT

WRITE (LUS, ' (A5,A6,5A1 1

)

' ) ' J' , 'ALPl U1 'U2' , 'VI V2'

17
0119

0120

0121

0122

0123

0124

0125

0126

0127

0128

0129

0130

0131

0132

0133

0134

0135

01 36

0137

0138

0139

0140

0141

0142

0143

0144

0145

0146

0147

0148

0149

0150

0151

0152

0153

0154

0155

0156

0157

0158

0159

0160

0161

0162

0163

0164

0165

0166

0167

0168

0169

0170

0171

0172

0173

0174

0175

0176

0177

0178

DO 30 J»1 ,NG

30 URITE(LUS, MI5,F6. 4,5(1 PEI 1 .3))')

1 J , ALP (J,l),U(J,l),U<J,2),y(J,l),V(J,2),P<J)

ENDIF

ENDIF

IF( T .GT. TMAX) GOTO 9999

C TO SOLVE THE DIFFERENTIAL EQUATIONS

C USING 2 STEP LAX-WENDROFF SCHEME

C MacCornack's « 0 thod.. BACKWARD PREDICTOR, FORWARD CORRECT

C CENTER DIFFERENCED ON TAU

C FIRST TO GET THE SPECIAL VARIABLES AND THEIR SPACIAL DERIVATIVES

CALL DGCOEF(ALP ,NG) ! DRAG AND GENERIZED COEFF

.

CALL TAUOFW(MUEF,DR,R,NG) ! STRESS

VJT=VJ

IF( T .LT. TJl .OR. T . GT . TJ2) VJT=0 . ! NO INJETCTION

CALL JET(BA,RHO,BR,V,VJT,NG)

CALL FSOFW(W,BA,R,NG)

C DETERMINE THE TIME-STEP SIZE

IF( NT .GT. 1) THEN

CALL FNDDT(DT,DR ,VDR,NG)

DT=DT+5 . *DTMIN/1 0 . **NT

IF(DT .GE. DTMIN) THEN

I=»DLOGl 0 (DT/DTMIN)

DDT=DTMIN*1 0 . **I

DT=DINT(DT/DDT+0 .001 )«DDT

IF( DT .GT. DTMAX) DT=DTMAX

ELSE

WRITE(LUP, ' <5X,A,1PE13.3) '

)

1 'STOP DUE TO TOO SMALL TIME STEP. DT-',DT

GOTO 9999

ENDIF

ENDIF

C BACKWARD PREDICTOR

CALL DWPDE(DW,RDP, DR , DT , RH , NGMl > ' INCREAMENT

DO 40 J=1 ,NGM1

DO 40 1=1 ,5

40 WP( J, I >=0 . 5»(W< J+1 , I)+W( J, I ) )+DW( J, I ) ! BASE+INCREAM

C PREDICTION DATA COMPLETED, CONTINUE FOR CORRECTION

CALL UOFW(WP,RH,NGMl

)

CALL DGCOEF(ALP,NGMl

)

CALL TAUOFW(MUEF,DR,RH,NGMl

)

CALL JET{BA,RHO,BRH,V, VJT,NG)

CALL FSOFW(WP ,BA,RH,NGM1

)

C FORWARD CORRECTION

CALL DWPDE(DW,P(2) , DR , DT , R ( 2 ), NGM2

)

P(1 )=PZERO

DO 50 J=1,NGM2

! MOMENTUM SOURCE

! CONVECTIVE -F AND SOURCE-S

I INITIAL INPLUSE TREAMENT

! ROUND OFF TIME STEP

18
0179

0180

0181

0182

0183

0184

0185

0186

0187

0188

0189

0190

0191

0192

0193

0194

0195

0196

0197

0198

0199

0200

0201

0202

0203

0204

0205

0206

0207

0208

0209

0210

021

1

0212

0213

0214

0215

0216

0217

0218

0219

0220

0221

0222

0223

0224

0225

0226

0227

0228

0229

0230

0231

0232

0233

0234

0235

0236

0237

0238

50

P < J+1 )=P< J) + <0 .25»<RDP( J)+RDP ( J + 1 ) )+0 ,5*P( J+1 ) )/R ( J)

DO 50 1=1,5

UN<J+1,I)=.25»(WP<J+1,I)+UP<J,I)) + .5»<W<J+1 ,I)+DW<J,I) ) P(NG)=P(NGM1

)

C 2ND STEP (PREDICTION+CORRECTION) COMPLETED

IF< NT .EQ, 1 ) GOTO 10

C Estirtation oP initial P coMpleted, return to the initial condition

C and start to advance the program in time.

»«»«««««««««««««««»«*« M »*««««««

C DATA W AT THE NEW TIME STEP COMPLETED

C IMPOSED B.C. *6.4

DW1=.5»(WN(2, 1 )+U<2, 1 ) )-0

.

ALPd ,1 ) = < ,5»<W<1 >5)+U(2,5))-DUl*DT/DR)/RH<l ) IF(ALP(1,1) .LT. ALMT(D) ALP ( 1 , 1 ) =ALMT ( 1

)

IF<ALP<1,1) ,GT. ALMT<2)) ALP ( 1

, 1

) »ALMT ( 2

)

WN(1 ,5)=ALP( 1

, 1 )*R(1

)

DUl»0.-0 ,5»(UN<NGM1 .1)+U<NGM1 .1))

ALP<NG,1 )=( ,5»<U(NGM1 ,5)+U

IF<ALP<NG,1) .LT. ALMT(1)>

IF(ALP<NG,1) .GT. ALMT(2))

WN < NG , 5 ) *ALP ( NG , 1 > »R ( NG

)

NG,5) )-DUl»DT/DR )/RH(NGMl

)

ALP(NG,1 )=ALMT(1

)

ALP<NG, 1 )=ALMT(2)

UN(1,1)=0. ! NO RADIAL VEL

.

WN(1 ,2)=0

.

UN(NG, 1 )=0

.

UN<NG,2)=0

.

UN(1,3)=UN<2,3)*WN(1 ,5)/UN(2,5)»R<l)/R(2) ! OMEGA=CON.

WN < 1, 4 ) =WN < 2 , 4 ) * ( R ( 1

) -WN <1,5))/(R(2) -WN (2,5))*R(1)/R(2)

IF<IW .EQ. 0) THEN

WN(NG,3)=WN(NGM1 , 3 ) *WN< NG , 5 ) /WN ( NGMl , 5 ) *R < NG ) /R < NGMl

)

WN(NG,4)=WN(NGM1 ,4)*<R<NG)-WN(NG,5) )/(R(NGMl )-WN(NGMl ,5)

)

1 *R<NG)/R<NGM1

)

ELSE

WN(NG,3)»0. ! NON-SLIP AT WALL

WN(NG,4)=0

.

ENDIF

C ARTIFICIAL DAMPING

DO 60 1*1,5

W(1 ,I) = <1 ,-DAMP)*WN(l , I )+DAMP*WN(2, 1

)

W(NG,I)=(1 .-DAMP)»WN<NG,I)+DAMP»WN(NGM1 ,1)

DO 60 J*2,NGM1

60 W(J,I) = <1 , -DAMP ) *WN< J , I ) +DAMP* < WN( J-1 ,I)+WN(J+1,I)-WN(J,I))

C SOLUTION FOR THIS TIME STEP COMPLETED

90 T-T+DT

GOTO 10

! UPDATE TIME AND CONTINUE TO THE NEXT

999 WRITE(LUP,7)

WRITE(LUP,7)

WRITE<LUP,7)

'OPEN FILE FAILED ON

NAMR

' lOSTAT*' , lOS

FILE: /

9999 CONTINUE

19
0239

0240

0241

0242

0243

0244

0245

0246

0247

0248

0249

0250

0251

0232

0253

0254

0255

0256

0257

0250

0259

0260

0261

0262

0263

0264

0265

0266

0267

0268

0269

0270

0271

0272

0273

0274

0275

0276

0277

0278

0279

0280

0281

0282

0283

0284

0285

0286

0287

0288

0289

0290

0291

0292

0293

0294

0295

0296

0297

0298

C CALL EXECCll ,JTIME,JTIME<1 )

)

C URITE(LUS, ' (35X,5I4) ' ) JTIME

CLOSE (LUS)

END

»«»««*««««**»«»*«»«»» »»«••«»

$EMA /DATA/, /WUWU/, /SOURCE/

SUBROUTINE INIT ,< 860425 . 1 537

>

C TO SET-UP THE INITIAL CONDITIONS

INTEGER I,IOS,ITLOC,IVTX,

Z ,ITIME,IW,LUP,LUS,NG,NGM1

PARAMETER (MM*101)

CHARACTER NAMR*16

REAL»8 ALMT, W , WP , WN , DU , RDP , RH

1 ,U,V,ALP,P,R, F,S

4 ,RHO,MUEF,V18,NA,ND, DO, GAMMA

7 ,BR,BRH, DAMP, DR

Y ,RJ1 ,RJ2,TJ1 ,TJ2,QJ,VJ

X ,A2,DELA,DEN1D2,D1 ,D2,PI ,RE,RMIN,RTANK

X ,DS,L3,VS,TS,PS

X ,ALP10, OMEGA, RPEAK,RJB,VPEAK

X ,EUF<2) ,VT,MU<2) ,MU1D2

COMMON

1 /ALPLMT/ ALMT(2)

Y /JETS/ RJl ,RJ2-,TJ1 ,TJ2,QJ,UJ

Z /CONTP/ IW,LUP,LUS,NG,NGM1 ,DAMP,DR

3 /DATA/ U(MM,2) ,V(MM,2) ,ALP<MM,2) ,P(MM) ,R(MM)

4 /DRAGl/ RH0<4) ,MUEF(2> ,yi8(2) ,NA,ND

3 /DSIZE/ D0(2) ,GAMMA(2)

7 ./SOURCE/ BR (MM) ,BRH(MM)

9 /UUWU/ U(MM,5) ,UP(MM,5) ,UN(MM,5) ,DW(MM,5) ,RDP(MM) ,RH(MM)

DATA PI/3. 141596D0/

C »*«*** NOTES:

C RHO

C ALMT

C DAMP

C DS,LS,yS,TS

C RE

C RTANK

C IVTX

C

C RPEAK,VPEAK

C OMEGA

C EVF,MUEF

C DO, GAMMA

C RJ1,RJ2,TJ1,

C TJ2,QJ,VJ

C

PHASE-1=»GAS, PHASE-2»LIQUID

DENSITY, RHO(l) < RH0<2>

LIMIT VALUES FOR ALPl, ALM ( 1 ) < ALP 1 < ALM (

2

) NUMERICAL DAMPING FACTOR, NORMALLY =0.

DENSITY, LENGTH, VELOCITY AND TIME SCALES

REYNOLDS *=Vs»R TANK *»RHO ( 2 ) /MU < 2 ) TANK RADIUS

TYPE OF INITIAL FLOW.

( 0=Rest , 1 *P ur« r o Ta T i on , 2»H . 0 , 3“GIT Vortex)

VORTEX PARAMETERS

PURE ROTATION. V=»OMEGA*r

EFFECTIVE VISCOSI TY , MUEF=* ( 1 t-EVP ) »MU

DIA. PARAMETERS: D=DO*ALP »»GAMMA

TO DEFINE JET SIZE, PUMPING TIME,

VOLUME FLOW RATE AND INJECTION MEAN SPEED

5 F0RMAT(2X,A,2(X,A)

)

7 F0RMAT(2X,A,3I5)

8 F0RMAT<2X,A,3(1PE12.4)

)

9 F0RMAT(A25,2(1PE15.4)

)

92 FORMAT(X,7( IPEl 1 . 4)

)

20
0299

0300 C DEFINE THE PARAMETERS FOR THE PROBLEM.

0301

0302 RTANK=1.D0 ! COULD BE SET TO 1 (M)

0303 RH0(2)=1 . OOOD+3 ! " (KG/M»»3)

0304 MU(2)=1 .514D-3 ! “ (KG/M-S)

0305

0306 C WRITE(1,8) 'ENTER RTANK<M) OR DEFAULT ',RTANK

0307 C READ(1,») RTANK

0308

0309 C The values of RTANK , RHO ( 2 )

, and MU<2) could all be set to 1, since

0310 C the length and density scales are based on RTANK, and RH0<2) and the

0311 C value of the viscosity MU(2) can be coHbined into and specified by the

0312 C Reynolds nunber RE. Thus all character istic sea les ( LS , VS , TS,and DS ) 0313 C are fixed after RE is given,

0314

0315 RE=1.0D5 ! LS««2/ < NU»TS ) =LS*VS/NU

0316 WRITE(1,8) 'enter Reynolds no., RE. D.F.=',RE

0317 READ(1,») RE

0318 LS=RTANK ! LENGTH SCALE (M)

0319 DS=RH0<2) ! DENSITY SCALE (KG/M»»3)

0320 VS=RE»MU(2)/RH0(2)/LS ! VELOCITY SCALECM/S)

0321 TS=LS/VS ! TIME SCALE (S)

0322 PS=DS»VS**2 ! PRESSURE SCALE

0323

0324

0325 C AFTER THIS POINT ALL VARIABLES ARE BAESED ON THE CHARA. SCALES

0326 C i.e. ALL VARIABLES ARE DIMENSIONLESS

0327

0328

0329 DEN1D2=1 .293D0/1 .000D3 ! D1/D2

0330

0331

0332

MU1D2=1 .71D-5/1 .514D-3

RH0(2)=RH0(2)/DS

1 MU1/MU2

0333

0334

0335

0336

MU<2)=MU(2)/(DS»LS»VS)

RHO( 1 )=DEN1D2*RH0(2)

MU< 1 )=MU1D2»MU(2)

! = 1/RE

0337 ALMT(1)=0. 00010 1 MIN. OF ALPl

0338 ALMT(2)=0 . 99999 1 MAX. OF ALPl

0339 DO( 1 )=1 .D-2/LS ! GAS DIAMETER at ALP1=»1

0340 D0(2)=l .D-2/LS 1 LIQUID DIAMETER at ALP2=1

0341

0342

GAMMAd )=2.D-1

GAMMA<2)=2. D-1

j 0343 EVF(2)=1 . D3 1 TURB.+PHASE-DISPERSION EFFECTS

0344 EVF( 1 )=DEN1D2/MU1D2»EVF(2 ) ! MODEL

0 345 DAMP=0 .

0

I NUMERICAL DAMPING FACTOR < e . g .=. 2

)

0346 NA=4. ! WEIGHTING EXP. FOR ADM

0347

0348

ND=4. ! WEIGHTING EXP. FOR DRAG

0349 WRITE(1,8) 'Enter DENSITY and VISCOSITY ratios.

0350

0351

WRITE(1,8) 'D.F.*' ,DEN1D2

READ<1,») DEN1D2,MU1D2

,MU1D2

0352

0353 URITE(1,8) 'Enter BASE DIAMETERS: D01,D02'

0354 WRITE(1,8) 'D.F.=',DO

0355 READ<1,») DO

0356

0357 URITE(1,8) 'Enter SIZE EXPONENT: GAMMAl , GAMMA2'

0358 WRITE(1,8) 'D. F. =»', GAMMA

21
0359

0360

0361

036x2 0363

0364

0365

0366

0367

0368

0369

0370

0371

0372

0373

0374

0375

0376

0377

0378

0379

0380

0381

0382

0383

0384

0385

0386

0387

0388

0389

0390

0391

0392

0393

0394

0395

0396

0397

0398

0399

0400

0401

0402

0403

0404

0405

0406

0407

0408

0409

0410

041

1

0412

0413

0414

0415

0416

0417

0418

READ<1,») GAMMA

WRITE(1,8) 'Enter weighting exponent: NA,ND. D.F,=',NA,ND

READU,*) NA,ND

WRITE(1,8) 'Enter GAS VOLUME FRACTION 1 ini ts : ALMT 1, ALMT2 .

WRITE(1,8) 'D.F,*',ALMT

READ<1,») ALMT

URITE(1,8) 'Enter eddy wisicosity factor. D.F.=',EVF

READd,*) EVF

IW=1

WRITE(1,7) 'Enter wall condition , l=nonslip , 0=slip . D.F.',IW

READd, ») IW WRITE(1,8) 'Enter nunerical damping factor. D.F.=',DAMP

READd,*) DAMP

DO 10 K=l,2

V18(K)=18.*MU(K)

10 MUEF(K)=MU<K)»d .+EVF(K) ) ! EFFECTIVE VISCOSITY FOR STRESS

RH0(3)=RH0d )»RH0<2)

RH0(4)=RH0d )-RH0(2)

RMIN=0.1 ! MINIMUM FLOW RADIUS IN THE TANK

NG=101 ! OF GRID POINTS USED

NGM1=NG-1

DR = d . -RMIN)/NGM1

C Initial c leann ing-up

.

DO 15 J=1,MM

DO 15 K=1 ,6

W( J,K )=0 . DO

15 U<J,K)=0.D0

C MOMENTUM SOURCE, JET CONDITIONS

RJ1=8.5D-1

RJ2*9.5D-1

VJ=10. ! TANGENTIAL INJECTION SPEED

TJ1=0

.

TJ2=«10 . WRITEd,8) 'Enter JET SIZE defined by RJ1,RJ2. D . F . - '

, R J1 , R J2

READd,*) RJ1,RJ2

WRITE(1,8) 'Enter INJECTION SPEED AND TIME RANGE, VJ,T1,T2'

URITEd,8) 'D.F.=' ,VJ,TJ1 ,TJ2

READd,*) VJ,TJ1,TJ2

QJ=<RJ2-RJ1 )*VJ ! JET VOLUME FLOU RATE

DO 20 J=t ,NG

R( J)=RMIN+< J-1 )*DR

RH< J)=R(J)t-0 .5*DR

BR < J)=DELA( (R < J) ) ,DR ,R J1 ,RJ2)/(2 . *PI ) ' JET DISTRIBUTION

20 BRH(J)=DELA( (RH(J) ) ,DR,RJ1 ,RJ2)/(2. »PI) I PER RADIAN

C SETUP INITIAL CONDITIONS

IVTX=0

22
0419

0420

0421

0422

0423

0424

0425

0426

0427

0428

0429

0430

0431

0432

0433

0434

0435

0436

0437

0438

0439

0440

0441

0442

0443

0444

0445

0446

0447

0448

0449

0450

0451

0452

0453

0454

0455

0456

0457

0458

0459

0460

0461

0462

0463

0464

0465

0466

0467

0468

0469

0470

0471

0472

0473

0474

0475

0476

0477

0478

OMEGA=0

.

yPEAK=0

.

RPEAKsRMIN

NAMR='SiHple uortex'

URITE(1,5) 'Enter data FILE NAME for initial cond., if any;'

WRITE(1,5) 'D.F,=',NAMR

READ(1,'(A)') NAMR

IF< NAMR .NE. ',' .AND. NAMR . NE . 'Simple vortex') THEN

C INITIAL CONDITION FROM A GIVEN FILE NAMR.

OPEN ( 99 , F ILE=NAMR , IOSTAT=IOS , STATUS= ' OLD '

, ERR=299

)

DO 25 J»1,NG ! INITIAL VALUES FROM FILE NAMR

25 READ<99,*) K , ALP < J , 1 ), U < J , 1 )

, U ( J , 2 )

, V ( J , 1 ), V C J , 2

)

CLOSE (99)

ELSE

C TO DEFINE INITIAL CONDITION.

ALP10=2.5D-1 ! INITIAL GAS VOL. FRACTION

WRITE(1,8) 'Enter initial value of alpl. D.F.='.ALP10

READ(1,») ALPIO

URITE(1,7) 'Enter type of vortex: 0=At rest,l=pure rotation'

URITE(1,7) '2=H.O. ,3=GIT. D.F.=',IVTX

READ(1,») IVTX

IFdVTX .GT. 0) THEN

IF (IVTX .GT. 1) THEN

URITE(1,3) 'Enter PEAK VEL . and LOCATION for classic vortex'

WRITE(1,8) 'D.F.=' ,VPEAK,RPEAK

READ(1,») VPEAK,RPEAK

IF(RPEAK .LE. 0.) RPEAK=1 . ! SINGULAR AT ZERO

ELSE

WRITE(1,8) 'Enter CIRCULAR SPEED(rad . /unit time).D.F.='

1 , OMEGA

READ(1,») OMEGA

ENDIF

ENDIF

DO 30 J=1 ,NG

ALP( J, 1 )=ALP1

0

VT=OMEGA»R( J) ! PURE ROTATION

IFdVTX .GT. 0) THEN

RJB=R( J)/RPEAK

IF( IVTX .EQ. 1) THEN ! H.O. VORTEX

VT=VT+1 . 398»VPEAK/RJB*( 1 ,-DEXP(-l . 25643»R JB*»2 )

)

ELSE I G.I.T. VORTEX

VT=VT+VPEAK»RJB»DEXP ( ( 1

. -RJB«*2)/2.

)

ENDIF

ENDIF

DO 30 K=*l ,2

IFdU .EO. 0) THEN

V(J,K)=»VT ! NO WALL

ELSE

V( J,K )=VT»( 1

. -R d) )»»0 . 1 ! BOUNDARY LAYER

ENDIF

30 U(J,K)=0. ! NO RADIAL VEL.

ENDIF

DO 40 J*1,NG ! FORM U FOR NUMERICAL CAL.

23
0479

0480

0481

0482

0483

0484

0485

0486

0487

0488

0489

0490

0491

0492

0493

0494

0495

0496

0497

0498

0499

0500

0501

0502

0503

0504

0505

0506

0507

0508

0509

051 0 051

1

0512

0513

0514

0515

0516

0517

0518

0519

0520

0521

0522

0523

0524

0525

0526

0527

0528

0529

0530

0531

0532

0533

0534

0535

0536

0537

0538

ALP<J,2)=1 .-ALP(J,1

>

P(J)=O.DO

W(J,5)=R(J)»ALP(J, 1

)

DO 40 K = 1 ,2

W ( J , K ) =ALP < J , K ) *U < J , K ) »R < J

)

40 W(J,K+2)=ALP(J,K)»V(J,K)»R(J)

PRINTOUT PARAMETERS

URITE

WRITE

WRITE

WRITE

WRITE

WRITE

WRITE

WRITE

URITE

WRITE

WRITE

URITE

(LUS,5)

<LUS,5)

<LUS,9)

(LUS,9)

<LUS,9)

(LUS,9)

(LUS,9)

(LUS,*)

<LUS,9)

(LUS, 9)

(LUS, 9)

(LUS, 9)

'INITIAL CONDITION FILE:',NAMR

'••DIMENSION UNITS ARE IN MKS^»'

'DENSITY SCALE(kg/M3) ' ,DS

'LENGTH SCALE=RTANK, (m)' ,LS

'VELOCITY SCALE(m/5) ' ,VS

'TIME SACLE(s) ' ,TS

'PRESSURE SCALE(Pa)' ,PS

'Reynolds nuMber, Re', RE

'Jet size, RJl ,RJ2' ,RJ1 ,RJ2

'Tangential jet, QJ,VJ',QJ,VJ

'Injection t i«e , T J1 , TJ2 '

, TJl ,TJ2

URITE(LUS, ' (/33X, "PHASE-l 8X , “PHASE-2 ") ' ) WRITE (LUS, 9) 'Density' ,RHO( 1

) ,RH0(2)

URITE(LUS,9) 'Viscosity ', MU

WR ITE( LUS , 9 ) 'Eddy viscosity factor',EVF

URITE(LUS,9) 'Base dia . ' ,DO

WRITE(LUS,9) 'Size exp.', GAMMA

WRITE(LUS,9) 'Phase 1 ini t s

'

, ALMT ( 1 )

, 1

. -ALMT ( 2

)

WRITE(LUS,^)

WRITE(LUS,*)

WRITE(LUS, ' ( “ OTHER CONSTANS : lU , I VTX , NA , ND , DAMP , VPEAK , R PEAK

1 *,0MEGA,D1/D2,MU1/MU2") ' ) WRITE(LUS, ' (2I3,2F4. 0 ,F5. 2) ' ) lU , I VTX , NA , ND , DAMP

WRITE(LUS,92) VPEAK , RPEAK , OMEGA , DENI D2 , MUl D2

URITE(LUS,»)

WRITE(LUS, ' (II 0,F10 .2) ' ) NG,RMIN

RETURN

299 WRITE(LUP,5) 'OPEN FILE FAILED ON INPUT FILEi',NAMR

WRITE(LUP,7> 'IOSTAT=' , lOS

STOP 111

END

••»»»»*»**»**»»••••••»•••••*•••*»«•••*••*•*••*•••*•*****»***»***•

REALMS FUNCTION DELA ( R , DR , RJ 1

, RJ2 ),( 860425 . 1 537

>

C TO DETERMINE THE EFFECTIVE NOZZLE SIZE AT EACH GRID LOCATION

C THE SIZE IS IN THE FRACTION OF GRID SIZE DR (i.e. 0<DELA<1)

REAL*8 R ,DR,RJ1 ,RJ2, R1,R2

R1=R-0.5^DR

R2=R1+DR

DELA=0 .DO IF(R1 .GE. RJ2 .OR. R2 . LE . RJl) RETURN

IF(R1 ,LT. RJl) R1=RJ1

IF(R2 ,GT. RJ2) R2=RJ2

DELA=(R2-R1 )/DR

RETURN

24
0539

0540

0541

0542

0543

0544

0545

0546

0547

0548

0549

0550

0551

0552

0553

0554

0555

0556

0557

0558

0559

0560

0561

0562

0563

0564

0565

0566

0567

0568

0569

0570

0571

0572

0573

0574

0575

0576

0577

0578

0579

0580

0581

0582

0583

0584

0585

0586

0587

0588

0589

0590

0591

0592

0593

0594

0595

0596

0597

0598

END

SUBROUTINE DERI VI (Y,DY,DX,N2) , <860425. 1537>

C GET 1ST DERIVATIVE, USING CENTERED DIFFERENCE

REAL»8 DX,Y(1 ),DY(1) ,C

EMA Y,DY

C=5 .D-l/DX

DO 10 J=2,N2-1

10 DY< J)=C»<Y( J+1 )-Y< J-1 )

)

*V CALL DWSUB(Y(3) ,1 , Y, 1 ,DY<2) , 1 ,N2-2)

DY< 1 )=<Y<2)-Y( 1 ) )/DX ! BASED ON 3-END PTS

DY(N2)=(Y(N2)-Y(N2-1 > )/DX

DY(1 )=2.*DY(1 )-DY<2)

DY(N2)=2. »DY(N2)-DY(N2-1 )

*V CALL DUSMY<5. D-l/DX, DY, 1 ,DY, 1 ,N2)

RETURN

END

EMA /COEFF/ . SUBROUTINE DGCOEF < ALP , N2 ),< 860425 . 1 537

>

C CALCULATE THE DRAG, ADDED MASS AND GENERIZED COEFF,

INTEGER J,MM,N2

PARAMETER (MM=«101)

REAL»8 ALP(MM,2)

EMA ALP

REAL»8 C, CPA, CD, RHO , MUEF , V 1 8 , NA , ND

X ,AA,AA1 ,AA2,AD,AD1 ,AD2,A1 ,A2,A12,DA12,DB2, D 1

, D2 , WTl , UT2 ,

X

COMMON

1 /COEFF/ C(MM,2,2) ,CPA(MM,2) ,CD(MM,2)

4 /DRAGl/ RH0(4) ,MUEF(2) ,V18(2) ,NA,ND

DO 50 J=1 ,N2

A1=ALP(J, 1

)

A2=ALP( J,2)

A12=A1»A2

C TO GET DRAG COEFF. AD

CALL SIZESCDl ,D2,A1 ) AD=V18( 2)»Al/( A2*D1*D1 ) ! =AD1 IF A2>.78

IF( A2 .LT. ,78) THEN

AD2=»V18(1 )»A2/< (1 ,-A2/.8)*D2)**2

IF(AD2 .LT, AD) THEN

WT1=A2*»ND

WT2=A1»»ND

AD=<AD*UT1+AD2*UT2)/(WT1+WT2)

ENDIF

ENDIF

C ADDED MASS COEFF. AA

25
0599

0600

0601

0602

0603

0604

0605

0606

0607

0608

0609

0610

061

1

0612

0613

0614

0615

0616

0617

0618

0619

0620

0621

0622

0623

0624

0625

0626

0627

0628

0629

0630

0631

0632

0633

0634

0635

0636

0637

0638

0639

0640

0641

0642

0643

0644

0645

0646

0647

0648

0649

0650

0651

0652

0653

0654

0655

0656

0657

0658

AAl=A12»RH0(2)/(Al+A2/( .5+1 ,5*A1 )

)

AA2=A12*RH0(1 )/(Al/< .5+1 .5*A2)+A2)

WT1=A2»*NA

WT2=A1»»NA

AA= ( AA 1 «WT1 +AA2»UT2 ) / ( UT 1 +UT2

)

C THE GENERIZED COEFF . CPA, C, AND CD

DB2=A12»RH0(3)+AA»(RH0( 1 ) »A 1 +RHO < 2 ) »A2

)

CPA< J, 1 )=A1«<A12*RH0(2)+AA)/DB2

CPA< J,2)=A2»(A12*RH0(1 )+AA)/DB2

C( J,1 ,2)*AA/DB2

C< J,2, 1 )=C<J,1 ,2)

C(J,1 ,1 )=A2*RH0(2)/DB2+C<J,1 ,2)

C( J,2,2)=A1»RHQ(1 )/DB2+C( J,2, 1

)

CD< J, 1 )=-A2*RHQ(2)»AD/DB2

CD< J,2)=A1»RH0< 1)»AD/DB2

50 CONTINUE

RETURN

END

$EMA /DATA/ , /FANDS/ , /TAU/ , /COEFF/

SUBROUTINE DUPDE ( DU , RDP , DR , DT , RR , N2 ),< 860425 . 1 537 > C

C TO GET DU OF THE PDEs

C INTEGER J, JPl ,K,KP2,MM,N2

PARAMETER (MM=»101)

REAL»8 DW<MM,5) ,RDP(MM) , DR,DT,RR(MM)

EMA DW,RDP,RR

C NOTES 1 COEFF, C =»C*ALP WHEN THIS IS CALLED

REAL*8 C, CPA, CD, U,V,ALP,P,R, F,S

8 ,TRR,TRA,TAA,RTRR,RTRA

4 ,RH0,MUEF,U18,NA,ND

X ,ALP1 ,ALP2,CP1 ,CP2,DTDR,DU1 ,DU1DT,G1 , G2 , HDT , UT

,

W J 1 ,UJ3,UJ4

COMMON

1 /COEFF/ C(MM,2,2) ,CPA(MM,2) ,CD(MM,2)

3 /DATA/ U(MM,2) ,U(MM,2) ,ALP <MM,2) ,P(MM) ,R (MM)

4 /DRAGl/ RH0(4) ,MUEF(2) , V18<2) ,NA,ND

6 /FANDS/ F(MM,5) ,S(MM,5)

3 /TAU/ TRR(MM,2) ,TRA(MM,2) ,TAA(MM,2) ,RTRR(MM,2) ,RTRA(MM,2)

DTDR=DT/DR

HDT=0 . 5*DT

DO 10 J=l,N2+l ! CHANGE C TO ALP*C

C( J, 1

, 1 )=ALP ( J, 1 )«C( J, 1 , 1 ) C< J,2, 1 )=ALP( J,2)»C( J,2, 1

)

C( J, 1 ,2)=ALP ( J, 1 )»C< J, 1 ,2)

10 C(J,2,2)=ALP(J,2)*C< J,2,2)

DO 20 J»1 ,N2

JP1=J+1

DU< J,5)=DTDR)*(-F( JPl ,5)+F(J,5) )+HDT»(S( JPl ,5)+S( J,5) > DO 25 K=1 ,2

26
0659

0660

0661

0662

0663

0664

0665

0666

0667

0668

0669

0670

0671

0672

0673

0674

0675

0676

0677

0678

0679

0680

0681

0682

0683

0684

0685

0686

0687

0688

0689

0690

0691

0692

0693

0694

0695

0696

0697

0698

0699

0700

0701

0702

0703

0704

0705

0706

0707

0708

0709

071 0 0711

0712

0713

0714

0715

0716

0717

0718

G1=.5*(C(JP1 ,K,1)+C(J,K,1) ) G2= . 5* ( C ( JP 1 , K , 2 ) +C ( J , K , 2 ) ) DU( J,K)=DTDR»(-F< JPl , K ) +F ( J , K ) +G1 » < RTRR < j P 1 , 1) -R TR R < J , 1 )

)

1 +G2»(RTRR(JP1 ,2)-RTRR( J,2) )

)

2 +HDT»(S< JPl ,K)+S( J,K)

)

KP2=K+2

25 DW< J,KP2)=DTDR*(-F( JPl ,KP2)+F(J,KP2)

1 +G1»(RTRA< JPl , 1 )-RTRA( J, 1)

)

2 +G2*(RTRA( JPl ,2)-RTRA( J,2) )

)

3 +HDT*(S( JPl ,KP2)+S( J,KP2)

)

C DP FOR PRESSURE CORRECTION

CP1=0 .5*(CPA( J, 1 )+CPA( JPl , 1 )

)

CP2=0 .5*(CPA< J,2)+CPA< JPl ,2)

)

IF< -DU<J,1) .GT. DW(J,2)) DW < J , 1

) =-DU < J , 2 ) ! DP>=0

RDP ( J)=(DU( J, 1 )+DU( J,2) )/(CPl+CP2)

DU( J,1 )=DW< J,1 )-CPl»RDP( J)

DU( J,2)=-DW< J,1

)

20 RDP ( J)=RDP ( J)/DT

RETURN

END

*EMA /COEFF/,/DATA/,/UWWU/

SUBROUTINE FNDDT ( DT , DR , UDR , NG ),< 860425 , 1 537

>

C DETERMINE THE TIME-STEP SIZE

INTEGER I , J,LUP,MM,NG

PARAMETER (MM=101)

REAL»8 DT,DR,VDR(2)

REAL»8 C, CPA, CD, RHO , MUEF , V 1 8 , NA , ND

9 ,y,WP,UN,DW,RDP,RH, U,V,ALP,P,R

X ,DUM1,DUM2

COMMON

2 /COEFF/ C(MM,2

3 /DATA/ U(MM,2)

4 /DRAGl/ RH0<4)

9 /UUUU/ U<MM,5

DATA LUP/1/

,2) ,CPA(MM,2) ,CD(MM,2)

,U(MM,2) , ALP (MM, 2) ,P(MM) ,R(MM)

,MUEF(2) ,V18<2) ,NA,ND

) ,WP(MM,5) ,WN(MM,5) ,DU(MM,5) ,RDP(MM) ,RH(MM)

DUM1=0

.

DO 10 J=1,NG

DUM2=-CD( J, 1

)

1 t-DAFS(U( J, 1 ) )/DR

2 +UDR( 1 )«ALP( J, 1 )*C( J, 1

, 1 ) 3 +UDR(2)*ALP( J,2)*C( J, 1 ,2)

IFCDUMl .LT. DUM2) DUM1=DUM2

DUM2=CD( J,2)+DABS(U( J,2) )/DR

1 +UDR<2)*ALP(J,2)»C(J,2,2)

2 >UDR(1 )»ALP(J,1)»C(J,2,1 ) IF(DUM1 .LT, DUM2) DUM1=DUM2

10 CONTINUE

' 1/St

! CONVECTIVE

' VDR=2»MUEF/DR**2

! CROSS VISICOSITY

C FIND THE MAXIMUN OF DU

27
0719

0720

0721

0722

0723

0724

0725

0726

0727

0728

0729

0730

0731

0732

0733

0734

0735

0736

0737

0738

0739

0740

0741

0742

0743

0744

0745

0746

0747

0748

0749

0750

0751

0752

0753

0754

0755

0756

0757

0758

0759

0760

0761

0762

0763

0764

0765

0766

0767

0768

0769

0770

0771

0772

0773

0774

0775

0776

0777

0778

CALL DUMAXC I ,DU, 1 ,NG) ! VECTOR OPERATION

»V DUM1=DABS(DW(I, 1 )

)

»V CALL DUMAX<J,DW(1 ,2) ,1 ,NG)

*V DUM2=DABS<DW< J,2)

)

DT=1 .DO/DUMl

RETURN

END «»««««««« *««»««««*»«««««»»«*«*«»«««»«««»«««««««««««»**««««« «««»»««»*

$EMA /DATA/ , /FANDS/ , /TAU/ , /COEFF/

SUBROUTINE FSOFU ( U , BA , RR , N2 ) ,<860425. 1537>

CC CALCULATE THE CONVECTIVE-F AND SOURCE-S TERMS

C INTEGER J,K ,KP2,L,MM,N2

PARAMETER (MM=101)

REAL*8 U,U,ALP,P,R, F,S, RHO , MUEF , V 1 8 , NA , ND

2 ,TRR,TRA,TAA,RTRR,RTRA, C, CPA, CD

X ,ALPD(2) ,RDU,RDV

REALMS U(MM,5) ,BA(MM,2) ,RR(1

)

EMA U,BA,RR

COMMON

1 /COEFF/ C(MM,2,2) ,CPA(MM,2) ,CD<MM,2)

3 /DATA/ U<MM,2) ,V(MM,2) ,ALP<MM,2) ,P(MM) ,R(MM)

4 /DRAGl/ RH0<4) ,MUEF<2) , V18(2) ,NA,ND

6 /FANDS/ F(MM,5> ,S(MM,5)

8 /TAU/ TRR<MM,2) ,TRA(MM,2) ,TAA(MM,2) ,RTRR (MM,2) ,RTRA(MM,2)

»V CALL DWMPY<W, 1 ,U, 1 ,F, 1 ,2*MM)

*V CALL DWMPY(U(1 ,3) ,

1 ,U, 1 ,F( 1 ,3) , 1 ,2«MM)

*V CALL DWMOV(U, 1 ,F< 1 ,5) , 1 ,MM)

DO 20 J=1,N2

F< J,5)=W< J, 1

)

S( J,5)=0

.

RDU=RR ( J)*<U< J, 1 )-U< J,2)

)

RDV=RR( J)»(V( J, 1 )-V< J,2)

)

DO 20 K=1 ,2

KP2=K+2

F(J,K)=U(J,K)»U( J,K)

F( J,KP2)=W( J,KP2)*U( J,K

)

S( J,K )=ALP< J,K )»< V( J,K )**2+CD( J,K )«RDU-C( J ,K , 1 )»TAA( J,

1

)

1 -C( J,K ,2)»TAA( J,2)

)

20 S< J,KP2)=ALP ( J,K )»<-U< J,K )»V< J,K )+CD( J ,K )»RDU+C( J ,K , 1 )*< BA( J,

1

)

1 +TRA( J, 1 ) )+C( J,K ,2)»(BA( J,2)+TRA( J,2) )

)

RETURN

END

SUBROUTINE JET < BA , RHO , BR , V , UJ , NG

)

C

C INJECTION MOMENTUM SOURCE

C

28
0779

0780

0781

0782

0783

0784

0785

0786

0787

0788

0789

0790

0791

0792

0793

0794

0795

0796

0797

0798

0799

0800

0801

0802

0803

0804

0805

0806

0807

0808

0809

081 0 0811

0812

0813

0814

0815

0816

0817

0818

0819

0820

0821

0822

0823

0824

0825

0826

0827

0828

0829

0830

0831

0832

0833

0834

0835

0836

0837

0838

PARAMETER (MM=101)

REAL»8 BA (MM, 2) ,RHO< 1 )

, BR ( 1 )

, V < MM , 2 >

, VJ

EMA BA,BR,V

DO 10 J=1,NG

IF< BR(J) .GT. 0.) THEN

Q=BR(J)»yj ! Flow rare/radian

BA( J, 1 )=Q»(VJ-V( J, 1 ) )»RHO< 1 > ! NET MOMENTUM GAIN

BA( J,2)=Q»<VJ-V( J,2) )*RH0<2)

ELSE

BA( J,1 )=0 .DO BA< J,2)=0 . DO

ENDIF

10 CONTINUE

RETURN

END

*»**»«**»»»»*»»»*»*»***«*»***»»**»»»*»»**«»»»*»»*»**«»*****»**»»**»**»

SUBROUTINE SIZES(D1 ,D2,ALP1 >,<860425. 1537)

C TO DETERMINE THE PARTICLE DIAMETERS

C REAL»8 D1 ,D2,ALP1 , DO, GAMMA

COMMON /DSIZE/ D0(2) ,GAMMA(2)

Dl=DO(l >»ALP1**GAMMA(1

)

D2=D0(2)*(1 . ODO-ALPl )»*GAMMA(2)

RETURN

END

«*«*««»««*««««««»««»««« »«««««»»»»*««»»*« »»««*«« »*»««««»*«*

SEMA /DATA/,/TAU/

SUBROUTINE TAUOFU < MU , DR , RR , N2 )

, <360425. 1537)

C STRESSES AND THEIR DERIVATIVES

C REAL*8 MU<2) ,DR,RR(1)

EMA RR PARAMETER (MM=101)

REAL»8 U,V,ALP,P,R

2 ,TRR,TRA,TAA,RTRR,RTRA

X ,TAMU,TMU

COMMON

3 /DATA/ U(MM,2) ,V(MM,2) ,ALP(MM,2) ,P(MM) ,R(MM)

8 /TAU/ TRR (MM,2) , TRA(MM,2) ,TAA<MM,2) ,RTRR (MM,2) ,RTRA(MM, 2)

DO 50 K=1 ,2

CALL DERI VI <U(1 ,K ) ,TRR ( 1 ,K ) ,DR ,N2)

»V CALL DWMPY(ALP( 1 ,K) , 1 ,TRR ( 1 ,K )

, 1 ,TRR( 1 ,K) , 1 ,N2)

*V CALL DWSMY(2 . »MU(K ) ,TRR ( 1 ,K )

, 1 ,TRR ( 1 ,K )

, 1 ,N2)

»V CALL DWDIV(V<1 ,K) ,1 ,RR,1 ,TAA(1 ,K) ,1 ,N2) ! V/R

DO 1 0 J=1 ,N2

10 TAA( J,K )=V< J,K )/RR < J)

CALL DERI VI (TAA< 1

, K > , TR A ( 1, K )

, DR , N2

)

I dU/dR

! »ALP

! »2*MU=TRR

! dO/dR

TMU=2 . *MU<K

)

DO 20 J=1 ,N2

TAMU=TMU*ALP( J,K)

TRR ( J,K )=TAMU*TRR( J,K)

29
0839

0840

0841

0842

0843

0844

0845

0846

0847

0848

0849

0850

0851

0 852

0853

0854

0855

0856

0857

0858

0859

0860

0861

0862

0863

0864

0365

0866

0867

0868

0869

0870

0871

0872

0873

0874

0875

0876

0877

0878

0879

0880

0881

0882

0883

0884

0885

0886

0887

0888

0889

0890

0891

0892

0893

0894

0895

0896

0897

0898

TR A ( J , K ) =MU ( K ) »ALP ( J , K ) *RR ( J ) *TR A ( J , K ) TAA ( J , K ) =TAMU»U ( J , K ) /R R ( J

)

RTRR<J,K )=RR< J)*TRR< J,K

)

20 RTRA<J,K)=RR(J)»TRA(J,K)

*v

»v*v

CALL

CALL

CALL

DWMPY(RR, 1

, TRAC 1 , K )

, 1

, TR A ( 1

, K )

, 1

, N2

)

DUMPY (ALP (1 ,K) ,1 ,TRA(1,K) ,1 , TRACI ,K),1 ,N2)

DWSMYCMUCK) , TRACI , K )

, 1, TR A ( 1

, K )

, 1, N2

)

! «R

! *ALP.

' «MU=TRA

»v

*v»v

CALL

CALL

CALL

DUDIVCUC 1 ,K )

, 1 ,RR , 1 ,TAA( 1 ,K ), 1 ,N2) ! U/R

DWMPYCALPC 1 ,K) , 1 ,TAA< 1 ,K )

, 1 ,TAA( 1 ,K ), 1 ,N2)

DUSMYC2 . »MU(K ) ,TAA( 1 ,K )

, 1 ,TAA( 1 ,K ), 1 ,N2)

I «ALP

! »2*MU=TAA

»v

»v

50

CALL DWMPYCRR , 1 ,TRR(1 ,K )

, 1 ,RTRR( 1 ,K ), 1 ,N2)

CALL DWMPYCRR , 1, TR A < 1

, K )

, 1

, RTR A ( 1, K ), 1, N2

)

CONTINUE

RETURN

END

! R»TRR

' R»TRA

$EMA /DATA/

SUBROUTINE UOFU < W , RR , N2 ) ,<860425.1537)

C CONVERTS U TO THE INDEPENDENT VAR TABLES ( U , V , ALP

)

PARAMETER (MM=101)

REAL»8 U(MM,5)*RR(MM)

EMA W ,RR

.

REAL*8 ALMT, U,V,ALP,P,R

X ,UJ6

COMMON

X /ALPLMT/ ALMT(2)

3 /DATA/ U(MM,2) , V(MM,2) , ALP (MM,2> ,P <MM) ,R (MM)

»V CALL DUDIV(U( 1 ,5) , 1 ,RR , 1

, ALP , 1 ,N2)

C CHECK VOLUME FRACTION 6 FLOW DIRECTIONS

DO 50 J=1 ,N2

ALP( J, 1 )=W( J,5)/RR < J)

IF(ALP<J,1) .LE. ALMT(l) .OR. ALP(J,1) . GE . ALMT(2>) THEN

IF(ALP(J,1) .LE. ALMT(l)) THEN

ALP( J, 1 )=ALMT< 1 ) W< J,3)=U( J,4)*ALMT(1 )/( 1 . D0-ALMT< 1 )

)

ELSE

ALP( J, 1 )=ALMT(2)

W<J,4)=U(J,3)»(1 .D0-ALMT<2) ) /ALMT (2)

ENDIF

W( J, 1 )»0 .DO W<J,2)»0 .DO

W< J,5)=ALP( J, 1 )*RR< J)

ENDIF

ALP( J,2) = l .D0-ALP( J, 1 )

IF( W(J,2) .LT. 0.) THEN ' PHASE-2 DOES NOT MOVE IN

W< J,2)=0

.

W( J, 1 ) = 0 . ENDIF

30
0899

0900

0901

0902

0903

0904

0905

0906

0907

0908

0909

0910

0911

0912

0913

0914

U< J,1 )=W<J,1 )/U<J,5)

V( J, 1 )=W< J,3)/U< J,5)

WJ6=RR( J)»ALP<J,2)

U( J,2)=W( J,2)/WJ6

J,2)=W< J,4)/UJ6

CONTINUE

CALL DUDIU(W,1 ,U( 1 ,5) , 1 ,U, 1 ,N2)

CALL DUDIV<W(1 , 3 ), 1

, W (1 , 5 )

, 1

, V , 1

, N2

)

CALL DWHPY(ALP( 1 ,2),1,RR,1,U<1,6),1,

CALL DWDIU(U(1 ,2) ,1 ,U<1,6),1,U(1,2),

CALL DUDIV(W<1 ,4),1,U(1,6),1,V(1,2),

RETURN

END

N2)

1 ,N2)

1 ,N2)

31
Exhibit A

:GLVM

Enter

90

A Sample Input

lu. for saving data. D . F . = b

Enter hILE NAME tor saving data.

TS153r :LB

Enter N0TES<<73 CHAR.) for the job

SAMPLE RUN OF TEST 153

enter Reynolds no., RE. D.F.= 1.0000E+05

« Enter

D,F.=

DENSITY and VISCOSITY ratios. ' .

1 .2930E-03 1 . 1295E-02

Enter

D.F.=

BASE DIAMETERS: D01 ,D02

1 .OOOOE-02 1 .OOOOE-02

f

* Enter

D. F.

SIZE EXPONENT: GAMMA 1 ,GAMMA2

= 2.0000E-01 2.0000E-01

* t Enter ueighting exponent: NA,ND. D.F.= *4 . OOOOE t-OO ^.'JOOOE + OU

« « Enter

D.F.=

GAS VOLUME FRACTION 1 imi ts : ALMT1 . ALMT2 .

1 .OOOOE-04 9.9999E-01

Enter eddy visicosity factor. D.F.= 1 . 1 4'^8E'-‘J2 I.OuQOEt-03

1000,1000

Enter wall condi t ion , i. = nonsl ip . 0 = si IP . D.t. ’

Enter nujTiericai damping factor. D . F . =• 0.0000E-*-00

Enter JET SIZE defined by RJ1.RJ2. D.F.= 3.5000E-01 9.5000E-0

« * Enter

D.F,=

1 ,0,1

Enter

D.F.=

INJECTION SPEED AND TIME RANGE. VJ.T1,T2

1.0000E+01 O.OOOOE+00 1.0000E+01

data FIl.E iNAME for initial cond . i^ any:

Simple vortex

Enter initial value of alpl. D.F.= 2.5000E-t.M

Enter type of vortex: 0=At rest.1=pure rotation

2=H.O. ,3=GIT. D.F.= 0

Enter

0.^01

Enter

.01

1 NP =

1 NP =

INITIAL and FINAL TIMEs. D.F.= O.UOOOE^OO 5.00UOE+00

TIME STEP for output. D.F.= 2.0000E-01

1 T=0.00GE+00 NT- 2 DT-i.OOOE-06

2 T = l.200E-02 NT= b DT-3.000F.-n3

32
Exhibit B

A Sample Output

TS153 T--00004 13 ON CR LB USING 00024 PLKS R = 0000

0001 1 SAMPLE RUN OF TEST 153

0002 INITIAL CONDITION FILE: . 0003 ^DIMENSION UNITS ARE IN MKS**

0004 DENSITY SCALE(kg/fn3) ! . OOnOE + 03

0005 LENGTH SCALE=RTANK,(m) 1 .OOOOE+00

0006 VELOCITY SCALE<m/s) 1 .5I40E-01

0007 TIME SACLE(s) 6.6050E+00

0008

0009

PRESSURE SCALE(Pa) 2.2922E>01

0010 Reynolds number, Re 1 .0000E^05

001 1 Jet Size, RJ1 , RJ2 8.5000E-01 9.5000E-01

0012 Tangential jet, QJ , V

J

1 .OOOOE-01 1 .OOOOE+00

0013

0014

Injection time,TJ1,TJ2 O.OOOOE+00 1 .OOOOE+00

0015 PHASE1

PHASE-2

0016 Dens 1 ty 1 .2930E-03 1 . OOOOE + OO

0017 U iscos i ty 1 . 1295E-07 1 .OOOOE-05

0018 Eddy viscosity factor 1 .OOOOE+03 1 .0000E + 03

0019 Base dia. 1 .OOOOE-02 1 .OOOOE-02

0020 S 1 ze exp

.

2.0000E-01 2.0000E-01

0021 Phase limits 1 .OOOOE-04 T .0014E-05

0022

0023

0024 OTHER CONSTANS: IN , IVTX . NA , ND . DAMP . VPEAK , RPEAK , OMEGA , D 1 /D2 . MU I /Mi

0025 1 0 4. 4. 0.00

0026 0. OOOOE+OO 1 .OOOOE-01 0. OOOOE+OO 1 .2930E-03 1 . 1295E-02

0027

0028 101 . 10

0029 1 NP = 1 T=0.000E+00 NT= 2 DT=1 .OOOE- 06

0030 J ALPi U1 U2 VI V2

0031 1 .2500 O.OOOE+00 O.OOOE^OO O.OOOE+00 O.OOOE+00 0 .OOOE + 0

0032 2 .2500 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 o.oooE+n

QUTP UT IN THE BETWEEN OMITTED

0230 97 .2500 -8.685E- 10 2.895E-10 3.513E-04 3.474E-04 3.969E-0

0231 98 .2500 -1 .775E-10 5.918E- 1 I 1 .268E-04 1 .250E-04 3.972E-0

0232 99 .2500 -2.785E-1

1

9.284E-12 4.319E-05 4.256E-05 3.973E-0

0233 100 .2500 -5.225E-12 1 .742E-12 1 .553E-05 1 .531E-05 3.373E-0

0234 101 .2500 O.OOOE+00 O.OOOE+00 O.OOOE+00 O.OOOE+00 3.973E-n

THE REST OF THE OUTPUT IS OMITTED

33
Complete etep

Correction step

Prediction step

Initial value

Figure 1. Two Step Difference Scheme (Backward Predictor - Forward Corrector Version)

34
35

Figure

2.

Velocity

Vector

Distributions

The

annular

region

between

two

dashed

circles

Is

the

region

of

Injection.
lx? ‘uoijODjj aujn|OA sdq

36

Kluure

CaM

Vultaoe

Fraction

Dial

r

ibutiona
U.S. DEPT. OF COMM. 1. PUBLICATION OR 2. Performing Organ. Report No. 3. Publication Dace BIBLIOGRAPHIC DATA

SHEET (See instructions)

REPORT NO.

NBSIR 86-3414 JULY 1986

4. TITLE AND SUBTITLE

A Computer Code for Gas-Liquid Tx^70-Phase Vortex Motions: GLVM

5. AUTHOR(S)

T. T. Yeh

6. PERFORMING ORGANIZATION (If joint or other than NBS, see in struction sj 7. Contract/Grant No. national bureau of standards

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234

9.

SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)

National Aeronautics and Space Administration

John F. Kennedy Space Center

Kennedy Space Center, FL 32899

10.

SUPPLEMENTARY NOTES

[~ ! Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11.

ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bi bliography or literature survey, mention it here)

A computer program aim*ed at the phase separation between gas and liquid at

zero gravity, induced by vortex motion, is developed. It utilize’s an-explicit

solution method for a set of equations describing rotating gas-liquid flows.

The vortex motion is established by a tangential fluid injection. A Lax-Wendroff

two-step (McCormack's) numerical scheme is used. This program can be used to

study the fluid dynamical behavior of the rotational two-phase fluids in a cylindrical tank. It provides a quick/easy sensitivity test on various parameters

and thus provides the guidance for the design and use of actual physical systems

for handling two-phase fluids.

12.

KEY WORDS (S/x to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) computer code; gas-liquid separation; numerical modeling; two-phase vortex motions

13.

AVAILABILITY

Unlimited

I I For Official Distribution. Do Not Release to NTIS

I I Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

14. NO. OF

PRINTED PAGES

44

15. Price

Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $9.95

USCOMM-DC 6043-P80
fe«<>JSS^»> - =-^- ' I' • ii»— M •! I

; _ i n-...jii t [ » [i*^i-ViT^ri i I II Iiiii r

"

''' jij«-»i— " _ - '•.'•

• 1^

W‘A?f& f «,»< ' J

- T11 H£v

7jt0u» «ma°3jfF^

-S' tolS j^boD i654y.<j«roO:'Ji

T

' % .156 W p*"i ' ' 4

? aJlAflA*': Lr^"

»:;-^ • m%r\ % ^ 7 irj . fn J fl $ ‘upr ‘ 4ja •« A»

rt yu4;A>y

.-vr

**- J‘ 5

*<19 JlpBW^X j

» 3 '

T*71^; ,iT

“’*'*'51 '1® i> . -Si 4

^ --"I T~-- y

rj’

S'^ *1 ‘^ft'»'“*'» *> ' ' '**j[^

iJJJSjte- *^*‘ ^ ^ • ‘’ ^ ^ •• njnpo# ^

_ ..A *'%'% ‘A ^ C ^ ^ ^ ate. ^ ^.te^.A^L.«teS

Li<il5 ixiJasgtwu :.%t*^4if tiT* i
a

ot

C'- .