UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Predicting the effect of vertical pipe orientation by the method of characteristics during sudden discharge… Mazaffari-Nejad, Hooshang 1978

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1979_A7 M69.pdf [ 4.12MB ]
Metadata
JSON: 831-1.0094637.json
JSON-LD: 831-1.0094637-ld.json
RDF/XML (Pretty): 831-1.0094637-rdf.xml
RDF/JSON: 831-1.0094637-rdf.json
Turtle: 831-1.0094637-turtle.txt
N-Triples: 831-1.0094637-rdf-ntriples.txt
Original Record: 831-1.0094637-source.json
Full Text
831-1.0094637-fulltext.txt
Citation
831-1.0094637.ris

Full Text

i  PREDICTING THE EFFECT OF VERTICAL PIPE ORIENTATION  BY THE METHOD  OF CHARACTER ISTICS DURING SODDEN DISCHARGE OF A TWO-PHASE FLUID  BY HOOSHANG |MOZAFFARI-NEJAD B. S. , North C a r o l i n a State U n i v e r s i t y , 1971 H. S., Oregon S t a t e U n i v e r s i t y , 1974  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES i n t h e Department of Mechanical He Accept  Engineering  T h i s T h e s i s as Conforming t o the Required  THE UNIVERSITY OF BRITISH COLUMBIA December, (c)  1978  Hooshang Mozaffari-Nejad, 1978  Standard  In p r e s e n t i n g  t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r  an advanced d e g r e e a t the U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I further agree that permission  for extensive  copying o f t h i s thesis  f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s .  I t i s understood that copying o r p u b l i c a t i o n  o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my written  permission.  Department o f  Mechanical Engineering  The U n i v e r s i t y o f B r i t i s h Columbia 2075 Wesbrook P l a c e V a n c o u v e r , Canada V6T 1W5 Date  A p r f T ?, 1979  ABSTRACT  Present  research  was  undertaken to study  o r i e n t a t i o n d u r i n g discharge opening. ranges  pipe  of s a t u r a t e d water through a sudden  A n a l y t i c a l i n v e s t i g a t i o n s have been made f o r v a r i o u s of i n i t i a l  pressure,  s t u d i e d by others. used e n t i r e l y with  the e f f e c t of  i n c l u d i n g high  A U m long 32 mm  pressures  i n diameter  previously  pipe  f o r t h i s i n v e s t i g a t i o n to f a c i l i t a t e  has  been  comparison  a v a i l a b l e experiments. When the i n i t i a l  pressure  i n the pipe was  just s l i g h t l y  atmospheric, body f o r c e s were  obviously  important  became e s p e c i a l l y pronounced when  role.  The  d i r a c t l y comparing pressures,  it  effect  upflow  has  versus  expected  downflow.  to  Even  play  above  at  higher  been seen i n flow v i s u a l i z a t i o n experiments  that the flow becomas s t r a t i f i e d  in a h o r i z o n t a l pipe.  When the  pipa i s o r i e n t e d v e r t i c a l l y  (up or down), t h i s w i l l no longer  t r u e s i n c e t u r b u l e n t mixing  occurs.  The r e s e a r c h to  was  undertaken because of i t s d i r e c t  the l o s s - o f - c o o l a n t a c c i d e n t  licenses be met. operation  for The  an  (L3CA).  In  order  be  relevance to  obtain  o p e r a t i n g r e a c t o r s , c e r t a i n s a f e t y s t a n d a r i s oust reactor  must  be  shown  to  be  safe  in  normal  as w e l l as i n the case of a h y p o t h e t i c a l a c c i d e n t .  Different  models  have  s i m p l e s t o f which i s the equal  been  used  to  analyze  tha LOCA,  v e l o c i t y - e q u a l temperature  (EVET)  model, which c o n s i d e r s both phases i n e q u i l i b r i u m a t a l l  times.  the EVET model good model  (used e n t i r e l y  for vertical  void f r a c t i o n  were  study  histories  made  the e f f e c t of g r a v i t y on pressure at  the  c l o s e d end.  i n v e s t i g a t e d because of i t s d i r e c t speed of sound.  Void f r a c t i o n , s e c t i o n , was  r e l e v a n c e to mixture  velocity,  were made f o r water a t i n i t i a l  35 atm  saturation m  (at the top of the p i p e ) .  was  density  c a l c u l a t e d at the  pressures of 1.2,  A l l cases s t a r t e d j u s t  with a s m a l l void f r a c t i o n of .001  diameter  column  pipe  exerting  (Figure an  1).  The  approximate  i n a 4 m loog  latter  of  35  atm.  For  1.2  pressure  g r a v i t y f o r c e to pressure f o r c e times  Euler  the .18  number) MPa  pressure was The "3P" gravity  3. 5  of  H  only  or .12  (reciprocal  defined as 3P'"  case, and  at .032  40000  Pa  was  was  MPa, of  the r a t i o of Froude  about 33%;  1.1 Jt f o r the  case  the  number  i t was  when  25%  initial  MPa.  r a t i o was  effect.  atm.,  1.8,  h y d r o s t a t i c head of water  comparable to the former pressures, but q u i t e d i f f e r e n t from  for  the  end. Runs  and  and  A l s o , t o t a l discharge or mass flow r a t e , a  combined measure of d e n s i t y and open  to be a  using the method of  r a t i o of vapor area to pipe area at any c r o s s  and  thought  flow.  S e v e r a l , computations c h a r a c t e r i s t i c s to  i n t h i s work) was  When  shown to be an important GP=1.1%, there was  measure of  the  neligible difference  between the pressures along the pipe f o r the upflow and  downflow  iv  case-  In both lower  initial  pressure  cases,  however,  there  seemed t o be a c l e a r d i f f e r e n c e between the p r e s s u r e s f o r upflow and downflow at t h e c l o s e d end.  V  TABLE OF CONTENTS Page ABSTRACT -.  i i  LIST OF FIGURES  vii  LIST OF TABLES  viii  LIST OF SYMBOLS  ix  ACKNOWLEDGEMENTS  xi  I.  -  1  1.1  General ....................................  1  1.2  P r e v i o u s Work  5  INTRODUCTION  1.3 II.  -  1.2.1  General  5  1.2.2  P r e v i o u s experiments ..............  6  1.2.3  EVET Model  7  1.2.4  EVUT Model  8  1.2.5  UVOT Model  9  Scope o f Present Work ......................  10  GOVERNING EQUATIONS:  12  2.1  C o n s e r v a t i o n Equations .....................  12  2.2  C o m p a t i b i l i t y R e l a t i o n s ....................  13  I I I . NUMERICAL PROCEDURES:  15  3.1  General ...........  3-2  Two-Phase Algorithm  3.3  Open-End -  3.4  Closed-End  3.5  Expansion Fans  IV. RESULTS AND DISCUSSIONS: 4.1  General .......  4.2  Results f o r A i r  4.3  R e s u l t s f o r Helium  4.4  R e s u l t s f o r Water  V. CONCLUSIONS: ..... VI. REFERENCES: .. VII. APPENDICIES: 7.1  Appendix A  7.2  Appendix B  7.3  Appendix C  7.4  Appendix D  7.5  Appendix E  7.6  Appendix F  LIST OF FIGURES Figure 1 2  Page  A Pipe at Quiescent S t a t e with Some I n i t i a l Void O r i e n t e d V e r t i c a l l y f o r Upflow Blow-down  ..---......105  Pressure H i s t o r y a t Closed End S h o r t l y a f t e r Pipe Burst ................................ 106  3  C h a r a c t e r i s t i c s Diagram i n z - t Plane  4  G e n e r a l Two-phase Algorithm  ....................107  to Advance  S o l u t i o n along C h a r a c t e r i s t i c ........................... 108 5  Have I n t e r a c t i o n s at the Boundaries  6  Expansion Fan a t the End of Duct Suddenly Opened t o Atmosphere Pressure H i s t o r y a t C l o s e d End of  7  Duct f o r A i r Problem  .....................109 .110  ....................................111  8  Expansion  fans f o r Helium Gun  112  9  Pressure a t C l o s e d End (P0=3.5 MPa)  10  Choking E f f e c t a t Open End  11  Void F r a c t i o n a t Closed End (P0=0. 12 MPa)  12  Void F r a c t i o n a t C l o s e d End (P0=0.18 MPa)  13  Pressure H i s t o r y a t Closed  14  Mass Flow Rate a t the Open End  15  Computer Program ........................................119  113 ....114  End (Low Pressure Case)  115 ......116 117 118  viii  i  LIST OF TABLES Table  Page  1  Parameters f o r Gas Dynamic Problems  122  2  R e s u l t s f o r Helium  123  3  Parameters f o r High-pressure Water Problem .............. 124  NOMENCLATURE A  =  duct c r o s s s e c t i o n flow area  (m )  a  =  sound speed (m/s)  aO  =  i n i t i a l sound speed (m/s)  Al  =  v o l u m e t r i c c o n c e n t r a t i o n of vapor a t a s e c t i o n  d  =  pipe diameter  I  =  frictional  g  =  g r a v i t a t i o n a l a c c e l e r a t i o n (m/s )  h  =  mixture  1  =  p i p e l e n g t h (m)  P  =  p r e s s u r e (Pa)  q  =  r a t e of heat added per u n i t mass of f l u i d  R  =  mixture  s  =  the known d i r e c t i o n  S  =  entropy  t  =  time (s)  T  =  temperature  u  =  mixture  x  =  guality  y  =  characteristic  Z  =  e l e v a t i o n above an a r b i t r a r y datum  z  =  space c o o r d i n a t e (m)  2  (m)  f o r c e per u n i t mass o f f l u i d  (N/kg)  2  enthalpy = xhg * (1-x)hf  d e n s i t y = AlRg • (1-Al)Rf  (J/kg)  (H/kg)  (kg/m ) 3  (J/kg K)  (K)  v e l o c i t y (m/s)  direction (m)  Subscripts g  =  vapor phase  f  =  liquid  i  =  interface  w  =  wall  phase  ACKNOWLEDGEMENTS  Author wishes t o express h i s deepest s u p e r v i s o r . Dr.  E.  G.  P.  essential  Hill  at the  and  is  McDonald,  initial  D.  Kawa,  and  Many  thanks  Dr.  W.  ( W h i t e s h e l l , Manitoba) f o r t h e i r many understanding  throughout the r e s e a r c h .  proof-reading checking  useful  p a r t of the t h e s i s and  some  of  o t h e r s f o r having  the  mathematical  Dr.  B.  of  AECL  Mathers  explanations  derivations;  given t h e i r help and  be  to  Steeves  with  proved t o  Thanks to Gary Ford  Alan  his  advice  Useful d i s c u s s i o n s  part of the r e s e a r c h  acknowledged.  Dr.  to  Hauptmann without whose h e l p and  t h i s work would not have been p o s s i b l e . Dr.  appreciation  for  and for  double-  and to many  support  whose  the  National  names  are  too numerous to mention h e r e . Support Council  of  facilities  for Canada  this is  were provided  research gratefully by  from  acknowledged.  Research Computing  the U n i v e r s i t y Computing Center.  .1  I.  INTRODUCTION  1. 1 General Because o f i t s d i r e c t r e l a t i o n s h i p t o the s o - c a l l e d l o s s o f c o o l a n t a c c i d e n t , (LOC&), the a n a l y s i s of t r a n s i e n t flow  and  associated  heat t r a n s f e r phenomena ( f l o w - b o i l i n g ) i s  f o c u s of ouch c u r r e n t r e s e a r c h .  I n t e r e s t i s c e n t e r e d around the  development of more g e n e r a l c o n s e r v a t i o n governing  equations  steam-water  for  laws  from  which  f l o w - b o i l i n g are d e r i v e d and improved  numerical techniques f o r s o l u t i o n of these governing Presence  of  time  flow  explore  than  of  problem  the  not  that flow  an  even  more  steady-state  usually  i t is  uniformly  equations.  as an a d d i t i o n a l v a r i a b l e i n the system makes  t r a n s i e n t two-phase  equiliorium,  the  can  not f u l l y  distributed.  challenging  flow. not  developed Thus,  area  to  For the t r a n s i e n t  reach  thermodynamic  and t h e two phases are interfacial  transport  processes have to be c o n s i d e r e d . Unsteady t h r e e - d i m e n s i o n a l of  the  complete  flow-boiling  problems a r e most r e p r e s e n t a t i v e situations  s t u d i e s , however these are very d i f f i c u l t suitable  in  reactor  t o handle.  safety  By use  of  c r o s s - s e c t i o n a l area averaging t e c h n i q u e s , a range o f  2  unsteady one-dimensional  analyses  considered.  the  Models  of  of  two-phase  actual f l u i d  flow  can  be  flow regime are then  formulated and a p p l i e d t o the problem t o - o b t a i n s o l u t i o n s . simplest  model  used  to  handle  v e l o c i t y - e q u a l temperature medium  problems i s  the  equal  (EVET) model, which assumes the  flow  as being a homogeneous, bubbly mixture (pseudo-gas) .  a r e s u l t o f t h i s assumption the assumed  such  The  to  have  equal  vapor  velocities  and and  liquid  phases  temperatures.  are  A more  complicated model i s the equal v e l o c i t y - u n e q u a l temperature model,  which  velocities  assumes but  sophisticated  the  1  liquid  different  (EVUT)  and vapor phases have e q u a l  temperatures.  model (which has met  As  A  relatively  with l i t t l e s u c c e s s to-date)  i s the unequal v e l o c i t y - u n e q u a l temperature  (UVOT) model,  which  assumes v e l o c i t i e s and temperatures are d i f f e r e n t f o r l i g u i d  and  vapor phases. Flashing,  which  is  conversion  of  part of the l i q u i d  vapor due to d e p r e s s u r i z a t i o n , o c c u r s when l o c a l p r e s s u r e to  o r below the s a t u r a t i o n v a l u e f o r any given temperature.  latent  heat  needed  to  f u l l n u c l e a t i o n to occur from the  state  depends  on  the  homogeneous e q u i l i b r i u m  drops The  vaporize the l i q u i d i s s u p p l i e d t o the  two-phase i n t e r f a c e from the bulk of the l i q u i d . for  to  model  initial  chosen.  The time taken depressurization  Models  such  (EVET) model do not allow any  as  the  time  for  3  departure  from  thermal  EVOT) models, however,  equilibrium. allow  a  few  Nonegulibrium milliseconds  (such  as  "relaxation  time". The  present  wave-tracinq analyze  technique,  gas  phase flow  work  dynamic  numerical  the  which  method of c h a r a c t e r i s t i c s o r  has  been  used  extensively  problems, to analyze one-dimensional  problems.  characteristics  uses  This  method  dz/dt=u*a,  difficulties  is  u-a.  present  based  It  in  on  greatly  other  a  mesh  The  grid  characteristic  lines.  grid  These  but p o i n t s of  pressure d i s t u r b a n c e s propagate;  tracing.  The  main  advantage  paths  propagation o f e r r o r s . must  propagate  along  are  of  a  of  intersection  followed  thus the i d e a of wave  technique  in  throughout  Since d i s c o n t i n u i t i e s  which  sound  in  two-phase  such  the  i s i t s lack of  in  the  solution  the c h a r a c t e r i s t i c s , methods which do  t r a c e the c h a r a c t e r i s t i c s w i l l d i f f u s e any S i n c e speed  model  are l i n e s along which s m a l l -  amplitude  characteristic  the  o f p o i n t s a t which the s o l u t i o n i s c a r r i e d  out i s not the usual space-time of  of  difference  c a l c u l a t i o n s , making i t u s e f u l f o r benchmark s o l u t i o n s and testing.  two-  reduces  finite  to  not  discontinuities.  mixture  is  a  very  s e n s i t i v e f u n c t i o n of void f r a c t i o n , s t a y i n g w i t h i n the scope of the  above  algorithm  is  not  always  possible.  c h a r a c t e r i s t i c changes i t s s l o p e ( i n the z - t plane)  In a d d i t i o n a abruptly  in  4  moving  from the l i q u i d t o the two-phase r e g i o n .  The reason f o r  such an abrupt change i n s l o p e i s the sudden drop i n  the  sound  In r e g i o n s  speed  from  the  liquid  to two-phase v a l u e .  where such r a p i d but continuous changes occur change  in  speed  of  as  a  local  result  sound, one needs t o r e l o c a t e the  advanced  point on the s a t u r a t i o n l i n e  ( i . e . enthalpy  liquid  value) or to use a very f i n e mesh.  value  to  Interpolations  two-phase  are  performed  along  the  changes  of  from  characteristics  the  to  s u b s t a n t i a l l y reduce the s p u r i o u s n u m e r i c a l d i f f i c u l t i e s . Hancox less  (1) has shown that a wave-tracing technique has much  numerical  diffusion  than  finite  However i t i s much slower and more addition  because  difference  difficult  to  program.  is  extermely  sensitive  to  irreqularities  p r o p e r t i e s such as £, SR/SP, 3R/8h, h, or quantities  In  of the nature of i t e r a t i v e technique (Newton-  Baphson method) r e q u i r e d to o b t a i n the s o l u t i o n at it  methods.  each in  the  externally  such as dA/dz, dZ/dz, d, and K / l .  point, fluid  specified  I f such terms are  d i s c o n t i n u o u s , d i f f i c u l t i e s can be expected. Initially, applying  the  the r e s e a r c h s e t out to study the p o s s i b i l i t y of method  of  hodograph  transformations  (used  s u c c e s s f u l l y to analyze i s e n t r o p i c gas-dynamic problems) t o  two-  y  phase  problems  (appendix  C  and  D).  However  d i f f e r e n t i a l e q u a t i o n s proved t o have undetermined  the  partial  coefficients.  5  Hence, closed-form The  a n a l y t i c a l s o l u t i o n s could n o t be o b t a i n e d .  present  research  has  focussed  o r i e n t a t i o n o f the body f o r c e ( g r a v i t y ) at a  on  the e f f e c t  range  of  of  initial  pressures.  1.2  P r e v i o u s work  1. 2. 1 General Solving  a complete unsteady  problem i s a f o r m i d a b l e t a s k . attempted  by  Panton (2),  three-dimensional f l o w - b o i l i n g  Thus averaging methods have  Delhaye(3),  V e r n i e r (4) , B a n e r j e e ( 5 ) ,  Mathers (6), A g e e ( 7 ) , Lyczkowski (8) t o make the problem These  averaging  methods  change  t h r e e - d i m e n s i o n a l t o an unsteady  and  terms  between  one-dimensional formed,  problem. which  A set include  (mass, momentum, and energy) between the phases  each  complementary  simpler.  the problem from an unsteady  of c o n s t i t u a t i v e r e l a t i o n s h i p s a r e then transfer  been  sets  phase of  and  the  equations  boundaries. to  ( c o n t i n u i t y , momentum, and energy) , and  These  the governing are called  are  equations "external*  1  interface constitutative relationships. To  handle t h e flow problem, one needs t o s o l v e the d e r i v e d  one-dimensional results  from  system  of  equations  numerically.  Predicted  numerical c a l c u l a t i o n s a r e then compared with t h e  6  experiments t o check accuracy of the o b t a i n e d s o l u t i o n .  1.2.2  P r e v i o u s Experiments F r o n t i e r work i n  Edwards  and  0 Brien(9) ,  of subcooled l i q u i d range  for  flow-boiling f  experiments  who s t u d i e d the r a p i d  were 4 m long  a l l runs, the pressure gage a  at  pressure  A l l pipes  chosen  for  (with v a r i o u s d i a m e t e r s ) . the  be  closed  end  of  the  reason  For  the  for  pipe  Thermal  the  dip.  and Hather(10) d i d another set of blow-down experiments  with 4.0 meter long p i p e s (with their  The  s m a l l dip 1 ms a f t e r r u p t u r e had o c c u r r e d .  n o n e g u i l i b r i u m was e x p l a i n e d to Edwards  by  these experiments was from 500 to 2500 p s i g , and the  experiments  recorded  done  depressurization  from long, h o r i z o n t a l p i p e s .  temperatures v a r i e d from 467 F t o 636 F. these  were  tests  saturation  started  pressure.  Baner jee (11) ,  which  diameters).  Most  of  from a s u b c o o l e d c o n d i t i o n , 35 bars above Another experiment was done by Uancox was  essentially  Edwards and Mather f o r 32 mm The  various  experiments  comparable  to  that  and of  diameter p i p e .  mentioned  above  d e p r e s s u r i z a t i o n t r e n d at t h e c l o s e d end.  showed  a  similar  The blow-down time i n  a l l cases was l e s s than a second. Numerical  models  temperature model, or  such the  as  equal  the  equal  velocity-equal  velocity-unequal  temperature  7  model  have  predicted  lower  measured e x p e r i m e n t a l l y . been  attributed  d e p r e s s u r i z a t i o n times than those  The reason f o r  to the assumption  this  discrepancy  has  of e q u a l v e l o c i t i e s f o r both  phases.  1.2.3  The  Equal V e l o c i t y - E q u a l Temperature  T h i s i s the s i m p l e s t model reactor  blow-down  used  problems.  to  (EVET)  predict  solution  B r i t t a i n (12) ,  velocities  and  temperatures  to  M a r t i n i (13),  T u r n e r ( 1 4 ) , Mathers(15), Karwat(16), and Moore(17) t h i s method i n t h e i r s o l u t i o n s .  Model:  incorporated  In t h i s model phases have e q u a l (homogeneous, bubbly  pseudo-gas).  The only a d d i t i o n a l r e l a t i o n s needed, were those f o r t r a n s f e r o f momentum and heat from the wall to the f l u i d . Edwards and Mather (10) compared their  own  the two was closed  end  numerically temperatures  experiments. not good.  at  pressure  experimentally)  phases  a l l times.  were  assumed  Second  the  d e p r e s s u r i z a t i o n o c c u r r e d sooner t h a t i n since the  phase actual  "EVET"  model  G e n e r a l l y speaking comparison  F i r s t the i n i t i a l  (measured because  their  the  was to  dip not  with  between at  the  predicted  have  equal  p r e d i c t e d long-term actual  situation  v e l o c i t i e s were c o n s i d e r e d equal i n the model. blow-down  experiments  (oecause  of  flow  In  being  s t r a t i f i e d ) the vapor has a higher v e l o c i t y than the l i q u i d .  8  Experiments  indicated  that  the  EVET model p r e d i c t e d t h e  r i g h t t r e n d , but improvements were needed more  accurate.  to  make  this  model  I n p a r t i c u l a r , extensions to allow f o r thermal  n o n e q u i l i b r i u m and unequal v e l o c i t i e s of phases vere needed. was of  c l e a r that improvements could n o t be made by simply changing constituative  1.2-4  relations.  The E g u a l V e l o c i t y - U n e q u a l Temperature  (EVUT) Model  T h i s model assumes t h a t vapor and l i q u i d same  velocities  situations  example. initially  but d i f f e r e n t  where  depressurization  this  of  temperatures.  model  seems  have  similar  necessary.  Early  showed t h a t d u r i n g r a p i d pressure  undershoots  as an  velocities  b u t the l i q u i d may be extremely  water,  the  There are c e r t a i n  During tiie blow-down, vapor and l i q u i d  high e n t h a l p y  values.  phases  of blow-down phenomenon may be mentioned  Edwards and O'Brien (9) of  It  are  superheated.  depressurization the  Banerjee(18) and F e r c h ( 1 9 ) , with a simple  saturation  relationship  i n t e r f a c i a l heat t r a n s f e r , have shown c l o s e agreement between  t h e i r numerical code and the e x p e r i m e n t a l l y measured (Figure  2) .  They  t r a n s f e r between  undershoot  chose two values of time c o n s t a n t s f o r heat  phases,  the  agreement with t h e experiments.  smaller  of  which  gave  better  9  So  f a r i t has become c l e a r t h a t s i m p l e r r e p r e s e n t a t i o n s o f  transient flow-boiling, trend  (EVET and EVDT models)  do  predict  the  o f the experiments g u i t e well and could be s u f f i c i e n t f o r  many c a s e s .  Another l e v e l of  more complex  predictions.  sophistication  This l e v e l i s  is  met  reguired  for  with unequal the  v e l o c i t y - u n e q u a l temperature (UVUT) model, which produces  other  difficulties.  1.2.5  The Unequal V e l o c i t y - U n e q u a l Temperature Sets  of  equations  of  this  type  (separate  momentum, and energy e q u a t i o n s f o r each phase) characteristics.  step  continuity,  produce imaginary  I t was demonstrated by Bryce(20) that  the s c h o o l of thought, the s i m p l e s t "slip")  (UVUT) Model  whatever  UVUT model (using a form  of  d i d not g i v e converging s o l u t i o n s as mesh s i z e and time  were  refined.  Thus  although  numerically  dispersed  s o l u t i o n s are a v a i l a b l e , t h e i r s i g n i f i c a n c e i s unknown. Banerjee(5)  and  Mathers(6)  recently  proposed an unequal  v e l o c i t y , unequal temperature, unequal pressure that  has  real  characteristics  speeds f o r s e p a r a t e somewhat  flows.  (UVUTUP)  model  and gave r e a l i s t i c propagation  Stuhmiller(21)  also  suggested  a  d i f f e r e n t approach with i n t e r f a c i a l p r e s s u r e d i f f e r e n t  from phases.  Under some c o n d i t i o n s t h i s model  system o f e q u a t i o n s with r e a l c h a r a c t e r i s t i c s .  also  led  to  a  10  None o f these unequal systems with r e a l c h a r a c t e r i s t i c s are entirely  convincing  d i f f i c u l t i e s show up. the  model,  but  since  there  are  conditions  More work i s needed not only  a l s o to d e r i v e a r e l i a b l e s e t of  to  where advance  constituative  relations.  1.3  Scope of Present Work The p r e s e n t r e s e a r c h i n t e r e s t has been to see the e f f e c t o f  g r a v i t y i n the blow-down problem, using EVET model and that a l l problems s t a r t e d o f f with a s m a l l void  fraction  (void  t o t a l area at any use  a  more  fraction  and  ( f i g u r e 1), say  .001,  i s r a t i o of the a r e a of void  pipe c r o s s - s e c t i o n ) .  sophisticated  assuming  I t would seem proper  complicated  to to  model such as EVUT  (which a l l o w s departure form thermal  eguilibrium).  However from  the r e s u l t s obtained by Hancox, and  Mathers  the  e f f e c t s p r e v a i l e d around 200-300 ms were unimportant At  high  to b o i l i n g e f f e c t s at the c l o s e d  pressure  thus,  t o p r e v a i l even a t high  encompassed  and  end.  regions, h o r i z o n t a l flow i s s t r a t i f i e d .  upward or downward f l o w s .  were expected  gravity  from the time of rupture  S t r a t i f i c a t i o n , however, i s not expected vertically  (1),  solving  a  range  to  be  significant  in  Hence the g r a v i t y e f f e c t s pressures.  The  research,  of problems i n c l u d i n g  the  11  reactor f a i l u r e  problems.  Tnermodynamic d e r i v a t i o n s have  been  sound speed as a f u n c t i o n of pressure pressure  and  enthalpy  properties  Keenan and Keyes (22).  to  obtain  the  and temperature i n s t e a d o f  (since the sound speed i s a f u n c t i o n o f  d e n s i t y and d e n s i t y a f u n c t i o n of pressure steam-water  made  used  are  and  enthalpy).  The  from continuous f u n c t i o n s i n  12  IX.  2. 1  Conservation  GOVERNING  Equations:  The homogeneous e q u i l i b r i u m to used  illustrate  averaged  (EVET) model has  field  in  gas-dynamic  equations  problems.  result  when  over the c r o s s - s e c t i o n , and  direction  of  flow.  hence  cross-section.  technique,  One-dimensional  flow  variables  vary  only  2  The  resultant  i n t o the f o l l o w i n g : *  (du/dz) = a  conservation  equations  i n the  2  may  be  (from Appendix A)  (2.1)  [ (q*uF) dR/dP| h - (u/A) (dA/dz) J = C2, (2.2)  DP/Dt*a R(iu/8z) = a t - (q+uF) 2>R/dh| p- (uR/A) (dA/dz) j=C3, 2  are  are t h e same a t any  Du/Dt+(aP/a»z)/R = -F - g(dZ/dz) = C l , Dh/Dt+ a  adopted  In t h i s model i t i s a l s o assumed that the  l i q u i d and vapor v e l o c i t i e s and temperatures  rearranged  been  the method of c h a r a c t e r i s t i c s s o l u t i o n  extensively  transient  EQUATIONS  2  (2.3)  where D/Dt a  2  = 8/2>t • u(3/tfz) ,  =1./ £ SR/dP|h • (dB/dh | P)/R J,  dfi/ap|h = p a r t i a l d e r i v a t i v e pressure a t c o n s t a n t enthalpy  1  of  density  with r e s p e c t to  13  Equation the  (2.1)  energy  i s the momentum, (2.2)  equation.  In  the c o n t i n u i t y , and  (2.3)  d e r i v i n g these equations d e n s i t y i s  used i n the form B = fi (P,h). Equations differential  (2.1)  throuqh  (2.3)  are  quasi-linear  e q u a t i o n s of the h y p e r b o l i c type.  of h y p e r b o l i c  equations  (for  example  see  partial  From the  Appendix  theory  B) ,  the  c h a r a c t e r i s t i c d i r e c t i o n s can be d e r i v e d : dz/dt = u, u*a,  u-a.  These r e p r e s e n t r e p e c t i v e l y the p a r t i c l e path, the disturbance  propagation  in  the  pressure d i s t u r b a n c e propagation  2.2  flow,  and  the  o p p o s i t e to d i r e c t i o n of f l o w .  ordinary d i f f e r e n t i a l equations,  equations, (Appendix  which  apply  along  these  called  dz/dt =  u+a,  dz/dt =  (C3 • B a d  ) dt = K d t ,  u-a, dP - Ba du =  compatibility  characteristics  B)  dP + fia du =  along  of  Compatibility Relations The  along  direction  pressure  (C3-BaC1)dt = L d t .  are:  14  along  dz/dt = u: dP - Bdh = (C3 - BC2)dt = M dt  where K, L , and H are c o n s t a n t s . P h y s i c a l l y , the f i r s t in  pressure  as  characteristic  related line  (a  equation r e l a t e s a s t e p enthalpy  change  energy e g u a t i o n ) .  in  two e g u a t i o n s represent a step change to  velocity  change  in  form of momentum e q u a t i o n ) . change  crossing  in  pressure  as  a  a characteristic line  crossing The  a  final  result  of  (a form of  15  III-  3.1  NUMERICAL PfiOCEDUEES  General A f i n i t e d i f f e r e n c e a l g o r i t h m was  along the c h a r a c t e r i s t i c s i n t h e (Figure  3).  z-t  Characteristics  Geometry long  and  chosen  .032  experiments. with a .001 suddenly 1).  m  for  in  (along  left  finite  were the open end and  to  which  small  the  right  running,  compare  At time  t=0* one  a pipe 4 a  with  available  difference  algorithms  task  of  the c l o s e d  advancing  water  end o f the pipe  the subsequent blow-down o c c u r r e d  was  (Figure  were a l s o used to  The boundaries  considered  end.  the  s o l u t i o n was  because of s t a r t i n g the blow-down with s m a l l  g r e a t l y reduced  initial  Al  (void  As a r e s u l t of t h i s s t a r t i n g two-pnase c o n d i t i o n , a  s p e c i a l t r a n s i t i o n algorithm  3.2  nodes  a l l the computer runs was  advance s o l u t i o n s on the boundaries.  mixture) was  known  (or t r a c e of the p a r t i c l e s ) .  diameter  void f r a c t i o n -  Special  traction).  from  I n i t i a l l y , the p i p e c o n t a i n e d high pressure  r u p t u r e d and  The  plane  lines  d i s t u r b a n c e s t r a v e l ) chosen were the running, and the p a r t i c l e path  used to advance s o l u t i o n s  not needed.  Two-phase Algorithm  (from one-phase l i q u i d  t o two-phase  16  The a l g o r i t h m opposite-running  used f o r s o l v i n g the point of i n t e r a c t i o n waves  inside  the  pipe  of  was c a l l e d "IIPHAS".  There were 10 unknowns z,t,u,P, and h r e l a t e d t o p o i n t s 3 and (Figure  4) .  Ten  nonlinear  equations  were used t o s o l v e f o r  these unknowns: s i x o f these e q u a t i o n s were directions other and  the c h a r a c t e r i s t i c  and t h e i r corresponding c o m p a t i b i l i t y r e l a t i o n s ; t h e  4 were o b t a i n e d 2-  The  using  i n t e r p o l a t i o n s between known nodes 1  Newton-fiaphson  method  was  used  t o solve  nonlinear  equations.  point  were kept; s o l u t i o n s f o r p o i n t 0 were d i s c a r d e d  the  3  3-3  i n z - t plane  Open-End  Algorithm:  For  open  the  end  boundary M  a choking c r i t e r i o n replaced  the r o u t i n e  c o m p a t i b i l i t y eguation.  M  Closed-End  CL0SED  w  as  the  "OPEN" was used, that  the l e f t - r u n n i n g c h a r a c t e r i s t i c and This required  that P3=Pfi i f u3<a3, (Figure 5-b).  Algorithm  For t h e c l o s e d end of the 4 m long M  (i.e.,  with the d i f f e r e n c e  otherwise u3=a3 and P3 was found i t e r a t i v e l y  3.4  at  progressed).  which was e x a c t l y the same as I I P H A S  its  these  Only the s o l u t i o n s f o r f i v e p r o p e r t i e s  p a r t i c l e s were not mapped c o n t i n u o u s l y  algorithm  0  s e t up t h e c o e f f i c i e n t  t e s t tube,  the  routine  matrix t o solve 5 equations and  17  5 unknowns f o r p o i n t solution  (Figure 5.a).  characteristics (meaning  3.5  3,  were  using  the  Expansion  of  8  helium) emanating analyze  Here o n l y l e f t - r u n n i n g and p a r t i c l e path used  to  determine the advanced  the  pressure c a s e s , an i n i t i a l  waves  for  water  problem  and  original  value  expansion  11 waves f o r a i r and  were advanced (Figure 3 ) .  dropped  atmospheric found  at  from  iterated  until  convergence  on  To s e t i n i t i a l from  equal increments on 8  dP*Ba du=0, 2  using  was  accounted  important), for  u=a  because  appropriately.  occurred  the  formulae (within a  results.  to  p r e s s u r e s (when the g r a v i t y  hydrostatic  head  was  not  As a r e s u l t o f t h i s shortcoming  c h a r a c t e r i s t i c s c r o s s e d one another and caused  producing meaningless  The  The expansion f a n method d i d not seem  s a t i s f a c t o r i l y a t low s t a r t i n g  effect  to  at g r i d s of  isentropically  Enthalpy was found from dP-Bdh=0 (Figure 6).  work  fan  drop dp from one fan t o t h e one immediately succeeding  t o l e r a n c e o f .5 m/s).  the  to  The v e l o c i t y u was  pressure  and  solutions  c o n d i t i o n s on the waves, pressure  were  solution  from the r u p t u r e end a t time z e r o was used  these waves i n t e r s e c t i n g one another  it.  of  Fans  high i n i t i a l  (consisting  waves.  method  t h e p r o p e r t i e s at p o i n t 3 ) .  For  its  Hewton-Baphson  difficulties,  18  Since i t was d i f f i c u l t  a t these  pressures t o i n i t i a l i z e the  problem p r o p e r l y , the idea o f u s i n g an expansion and  8  wavelets  initially  along  (gradual the  t=0  pressure axis.  reduction) This  method  f a n was dropped were  chosen  was s u c c e s s f u l  because i t enabled posing the h y d r o s t a t i c e f f e c t s c o r r e c t l y .  19  IV  4.1  RESULTS AND DISCUSSIONS  General Tiie method o f c h a r a c t e r i s t i c s a l g o r i t h m was  computer  runs f o r a i r , helium, and water.  at an i n i t i a l l y This  s a t u r a t e d c o n d i t i o n with  eliminated  two-phases Accuracy  the  of  and  to  0.001  difficulty  void  fraction.  i n going from one to saturation  line) .  t h e program was checked by making p r e l i m i n a r y  helium  Eleven f a n s were  problems,  pressure water problem.  make  Hater was c o n s i d e r e d  ( i . e . r e l o c a t i n g advance point on  with a i r and helium. air  spurious  used  chosen  to  analyze  runs the  and e i g h t waves t o analyze t h e high  These expansion waves emanated from the  open end a t t=0«- due t o a sudden  rupture.  The  expansion  fan  method worked w e l l f o r high p r e s s u r e s s i n c e g r a v i t y d i d not play an  important  than 1.5%).  role  ( r a t i o of g r a v i t y t o pressure ,GP, was l e s s  Using more than 8 waves f o r water d i d not seem t o  make much d i f f e r e n c e i n the accuracy o f r e s u l t s .  For the water  runs a t lower p r e s s u r e s , however, t h e expansion f a n method d i d not prove to be the c o r r e c t way o f a n a l y z i n g the problem. The r a t i o GP was 33% f o r P0=1.2 atm, and 25* f o r P0=1.8 atm (meaning measurable cases  the  r a t i o of gravity to the  pressure  force).  s t r a t e g y was t o i n i t i a l i z e the problem  In  both  along the t=0  20  a x i s to account end  f o r h y d r o s t a t i c head.  The  pressure a t the  was allowed t o drop g r a d u a l l y , so t h a t a c e n t e r e d  fan d i d not  exist.  tried.  was  It  Several  rates  of  conditions  f a n method was  prevailed  because  much l e s s than the l o c a l speed not  seem  organized.  c a u s i n g computational  were  case l e s s than 5  d e p r e s s u r i z a t i o n r a t e did not a f f e c t the numerical when the expansion  expansion  depressurization  found t h a t f o r the .12 MPa  open  calculations,  used f o r P0=1. 2 atm,  subsonic  the v e l o c i t y a t the open end of sound.  The expansion  Initializing  the  was  fan  Some of the f a n s c r o s s e d over one difficulties.  ms  did  another problem  from t=0 a x i s , however, seemed to work g u i t e w e l l . DZ/dz=1  defined  upflow,  and  dz/dz=-1 d e f i n e d  Each run was  c a r r i e d f o r 30 time advances of about 140  4.2  for Air  fiesults  An i s e n t r o p i c problem with a i r accuracy  of  diameter  was  7.  MPa  subsonic. used  to  the  algorithm.  A  4  was  chosen  to  of i t s o r i g i n a l value  (Table 1).  Isentropic relations including state initialize  the  problem  isentropic.  check  the  along  mm  pressure  of  The flow  was  eguations  were  the f a n s a t t=0 + .  r e t a i n the p a r t i c l e path c h a r a c t e r i s t i c s , the flow completely  ms.  meter long p i p e with 32  d e p r e s s u r i z e d r a p i d l y from an i n i t i a l  to .738  downflow.  can  not  To be  Hence a very s m a l l amount of f r i c t i o n  21  was  i n t r o d u c e d as an approximation  f r i c t i o n f a c t o r chosen was  of i s e n t r o p i c flow.  f=1.E-6  The  (making f r i c t i o n f o r c e  F=(2f/d)u^) very s m a l l to model i s e n t r o p i c  situation.  The a l g o r i t h m progressed i n z - t plane with i n t e r a c t i o n s reflected  waves  emanating from  the  the open end.  those of Owczarek  4.3  from  end  with  oncoming  R e s u l t s compared q u i t e  (23) f o r i s e n t r o p i c flow  waves  well  with  (Figure 7).  R e s u l t s f o r Helium I s e n t r o p i c problems may  hodograph  transformations  s o l u t i o n s are hypergeometric case of helium, tube  g e n e r a l l y be s o l v e d a n a l y t i c a l l y (using  Biemann  functions  Invariants).  (Appendix  C).  In  by The the  with k=5/3, the a n a l y t i c s o l u t i o n f o r a ruptured  i s q u i t e simple.  time and  closed  of  Parkinson  (z and t) can be expressed  (24)  showed t h a t the space  as v e l o c i t y and  sound  speed  and (u  a) as f o l l o w s : z(u,a) = 2tu(9a2-3u2*27a02 J/(27a ) 3  t ( u , a ) = 6 (9a -u2*9a02)/(27a3) 2  where  aO  and  respectively  a  are  the  initial  and  l o c a l speed  of sound,  (Figure 8)-  The same i n i t i a l c o n d i t i o n as f o r the a i r problem was to  solve  the  helium  problem  compared with a n a l y t i c a l v a l u e s .  numerically Accuracy  was  and to  used  results  were  within  0.1%  22  (Table  2).  The  friction  r e t a i n the p a r t i c l e path  4.4  again chosen as 1.E-6  f a c t o r was  characteristic.  R e s u l t s f o r Sater Three runs  condition  were made with water.  (P0=3.5  HPa,  and  emanating from the open end were  Values  for  For  T0=243 were  the  C),  Initial  (Table 3).  (vertical line i n z-t plane),  intermediate  waves  were  Hence i n i t i a l c o n d i t i o n s were s e t  high  On  on  the  fans  conditions the  trailing  choking  iterated  up  pressure  e i g h t expansion  chosen.  s e t on l e a d i n g edge of the fan  edge of the f a n  occurred.  isentropically. expansion  fans.  Boundary c o n d i t i o n s were c o n s i d e r e d at c l o s e d and open ends. the  closed  pressure was choking  end  the  velocity  was  zero while at the open  s e t equal to ambient i f the flow was  occurred  iteratively.  (u=a),  local  pressure  A two-phase m u l t i p l i e r between  chosen to determine f r i c t i o n  factor  w  subsonic. was  0  n  solved  and  M  1"  downflow Hence  in  even  horizontal  the though  pressure flows  history were  little  at  c l o s e d end  fundamentally  to v e r t i c a l o r i e n t a t i o n s (because  being s t r a t i f i e d , but v e r t i c a l g r a v i t y e f f e c t was  flow  At end If for was  (25).  There d i d not seem to be much d i f f e r e n c e between upflow  very  to  being  sensed  and  (Figure 9).  different  from  of h o r i z o n t a l flow turbulent  i n the pressure  mixing), history.  23  In t h i s case GP at  about 90 ms  was was  small felt  (1.1%).  The hump i n the  pressure  to be due to the i n t e r a c t i o n of  curve  compression  waves i n the medium. The speed of sound was fraction was  and  depeneded  about 110 m/s  a very  sensitive  l i t t l e on p r e s s u r e .  when v o i d f r a c t i o n  was  .5  function  I t s minimum  and  density.  value  Sound speed i n i t i a l l y  enthalpy,  decreased  each s u c c e s s i v e fan  (as v o i d f r a c t i o n i n c r e a s e d t o .5) and  i n c r e a s e d to 177 m/s  (when void f r a c t i o n  1) .  Quality  throughout  the  run  increased  was  less  f r a c t i o n , though, i n c r e a s e d q u i c k l y to 0.9 For P0= 1.2 e i g h t wavelets The  atm  ( s l i g h t l y more than  than  then  -5  0.3.  atmospheric  to Void  pressure),  s e t the i n i t i a l c o n d i t i o n s from the t=0  axis.  h y d r o s t a t i c head of water column exerted a pressure 40000 Pa, comparable t o i n i t i a l  top  pipe  the  in  some  cases.  g r a v i t y f o r c e t o pressure f o r c e , GP, MPa  from  on  and h i g h e r v a l u e s .  approximately of  void  (Figure 10),  Pressure decreased on each s u c c e s s i v e fan as d i d temperature,  of  For .12 was  33%;  pressures MPa,  at  of the  the r a t i o o f  i t was  25% f o r  .18  case. The «*GP  M  r a t i o showed to be an  gravity effect.  In both lower  MPa,  fraction  the  c l o s e d end  void  f o r both  up  measure  of  pressure cases, i . e . , .12 and  increased and  important  down  the .18  g r a d u a l l y with time a t the flow  (Figure  11  and  12).  24  Pressure  at  the  c l o s e d end f o r downflow i n c r e a s e d s l i g h t l y t o  above i t s o r i g i n a l v a l u e at around 90 os f o r both MPa  pressure cases ( F i g u r e 13) .  t o wave i n t e r a c t i o n of compression pressure r e g i o n  (bottom o f p i p e ) .  .12  and  The hump was thought t o be waves coming  from  the  .18 due high  25  V CONCLUSIONS I t was found d u r i n g t h e course of the present r e s e a r c h t h a t the  method  of  characteristics  is  a  very a c c u r a t e n u m e r i c a l  technigue f o r a n a l y z i n g two-phase problems . comparing  numerical  results  obtained  f o r helium  a v a i l a b l e c l o s e d form a n a l y t i c a l s o l u t i o n solution). technigue  T h i s was shown with  (hypergeometric  Tne method of c h a r a c t e r i s t i c s , however, i s a since  by the  series costly  convergence f o r each advanced p o i n t i n t h e z - t  plane i s r e q u i r e d .  Hence i t i s  a  comparisons with other numerical  good  method  for  benchmark  techniques.  Reynolds number f o r the blow-down water problems was o f the order  of  10 . 6  Hence  as  many  other  authors have done, t h e  v i s c o s i t y e f f e c t s were c o n s i d e r e d n e g l i g i b l e i n the flow and  were  prevailed,  neglected  from  the  analysis.  however, a t the wall i n the form  Comparison  between  experimental  Viscosity  medium effects  of f r i c t i o n .  results  of  vertically  o r i e n t e d p i p e s and t h e p r e d i c t e d numerical values should be made as  a  next  s t e p t o see i f the EVET model g i v e s f a i r  f o r such flow regimes mixing  and  are  not  (since v e r t i c a l stratified  g e n e r a l body f o r c e s were found especially differences  a t lower in  initial  flows  involve  prediction turbulent  as a r e h o r i z o n t a l f l o w s ) .  to s i g n i f i c a n t l y  pressures.  There  In  a f f e c t the f l o w , were  pronounced  v o i d f r a c t i o n and pressure h i s t o r i e s f o r upflow  26  and downflow case at t h e c l o s e d end the mass flow r a t e s f o r upflow the  pipe  of the pipe.  were q u i t e d i f f e r e n t  higher).  f o r the case of .12  this  work  available  HPa,  of  whereas  (when  pressure  Accurate property e g u a t i o n s are needed  otherwise  the r e s u l t s f o r sound speed in  to g r a v i t y  and downflow, at the open end  they were r a t h e r s i m i l a r f o r the case of .18 HPa was  Due  while  tabulated  are meaningless  trying  to  ( t h i s was  calculate  thermodynamic  sonic  property  underlined  speeds  from  values f o r Freon-  114). Generally expansion at  for  Hence at  any  the  found  orientation  along  the  pipe  pressure  the  pipes.  initial  was  not  initial  Expansion  pressure  simulating for  the  to  gradually  drop.  a r a p i d opening. opening  the case of .12  unimportant. was  because  This  I t was  assumed  f a n s c o u l d not be  pressure c a s e s , s i n c e i n i t i a l uniform  of  was  the  T h i s i s because  was r a t h e r  Instead of an expansion f a n , the pressure at allowed  conditions,  t o be adeguate.  pressures g r a v i t y e f f e c t  however, f o r low in  initial  fan method was  high i n i t i a l  uniform  high  used,  pressure  h y d r o s t a t i c head. the  found  found t h a t l e s s  open  end  adeguate than  was for  5  ms  d i d not make a d i f f e r e n c e i n f i n a l r e s u l t s i n HPa.  For f u t u r e work, f u r t h e r r e s e a r c h of the u n f i n i s h e d work i n hodograph  transformation  for  non-isentropic  problems  is  27  suggested  (Appendix  experimental c l o s e d form  and  D).  This  Bust  be done by mixed  p e r t u r b a t i o n techniques g i v i n g r i s e to semi-  solutions.  I t i s evident from the l i t e r a t u r e survey t h a t so f a r i n a l l p r e v i o u s work, the two-phase and well studied.  s u p e r c r i t i c a l r e g i o n s have been  The n e a r - c r i t i c a l r e g i o n could be an  interesting  area f o r f u t u r e r e s e a r c h . S i m i l a r to the work done by Hancox(1), i n the present the f i n a l form of the governing of d e n s i t y i n the form R=R (P,h). because i n the  e q u a t i o n s are d e r i v e d making The  compressed l i q u i d s t a t e , d e n s i t y i s a r a t h e r  errors resulted I t may  use  reason f o r such a c h o i c e i s  i n s e n s i t i v e f u n c t i o n of pressure and enthalpy and  small.  research,  from forming prove,  hence the  diffusion  the d e n s i t y d e r i v a t i v e s are q u i t e  with f u r t h e r s t u d i e s , to be  advantageous  to use the d e n s i t y i n a l t e r n a t e forms such as R=R(?,u) or R=R(P,T), where u i s the i n t e r n a l energy. Work may  a l s o be pursued  s o p h i s t i c a t e d technique  i n the area of developing a more  f o r handling complicated  with a s s o c i a t e d changes of area and  pipe net-works  with the shock-waves p r e s e n t .  28  VI- REFERENCES 1  Hancox, W. T., Bathers, W. G., and Kawa, D. , A n a l y s i s o f Transient Flow-Boiling; Application of the Method of Characteristics, AIChE Paper 42, 15th N a t i o n a l Heat T r a n s f e r Conference, August 1975.  2  Panton, B., "Flow P r o p e r t i e s f o r the Continuum Viewpoint o f a Non-Equilibrium G a s - P a r t i c l e M i x t u r e , J . F l u i d Mech., 31, 273-303, 1968.  3  Delhaye, J.M. , and Achard, J . L. ,On the Averaging O p e r a t o r s Introduced i n Two-Phase Flow Modelling, Invited Paper presented a t 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t TwoPhase Flow, T o r o n t o , August 1976.  4  V e r n i e r , P. , and Delhaye, J . H. ,General Two-Phase Flow Equations A p p l i e d to the Thermo-Dynamics of B o i l i n g Water R e a c t o r s , Energie P r i m a i r e 4, No1-2, 1968-  5  Banerjee, S., Ferch, R., Mathers, W.G., McDonald, B.H., Dynamics of Two-Phase Flow i n a Duct, Proc. Sixth I n t e r n a t i o n a l Heat Transfer Conference, Toronto, August 1978.  6  Mathers, S.G., F e r c h , R.L., Hancox, W.T., and McDonald, B-H., Equations f o r T r a n s i e n t Flow-Boiling i n a Duct, Invited Paper Presented a t 2nd CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow, P a r i s , June 1978.  7  Agee, L.J. , Banerjee, S., Duffey, R. B., and Hughes, E. D., Some Aspects of Two-Fluid Models and t h e i r Numerical Calculations, Invited Paper Presented at 2nd CSNI Specialists Meeting on T r a n s i e n t Two-Phase Flow, P a r i s , June 1978.  8  Lyczkowski, R.H., T h e o r e t i c a l Bases o f the D r i f t - F l u x Eguations and Vapor Drift-Velocity, Proc. Of  Field Sixth  29  International 1978-  Heat  Transfer  Conference,  Toronto, August  9  Edwards, A. B. , and O'Brien, T. P., S t u d i e s of Phenomena Connected with D e p r e s s u r i z a t i o n of Mater B e a c t o r s , J o u r n a l BNES V o l . 9, pp. 125-135, A p r i l 1970.  10  Edwards, A. B., and Mather, D. J . , Some UK S t u d i e s B e l a t e d t o L o s s of C o o l a n t Accident, Topical Meeting on Mater Beactor Safety, Salt Lake C i t y OSAEC Report CONF-730304 (1973.)  11  B a n e r j e e , S., and Hancox, M.T., On the Development Methods f o r A n a l y z i n g T r a n s i e n t Two-Phase Flow, I n t . Multiphase Flow, i n press, 1978-  12  B r i t t a i n , I . , and F a y e r s , F . J . , A Beview of Becent Mork i n the M o d e l l i n g and S o l u t i o n of t h e Hydrodynamic E g u a t i o n s i n the LOCA Code, EELAP-OK, I n v i t e d Paper presented at t h e 1st CSNI Specialists Meeting on Transient Two-Phase Flow,Toronto, August 1976.  13  Martini, B., P i e r n i n i , G.C., and Sandri, C, A Oned i m e n s i o n a l T r a n s i e n t Two-Phase Flow Model and I t s I m p l i c i t Finite D i f f e r e n c e S o l u t i o n , I n v i t e d Paper Presented a t the 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow, Toronto, August 1976.  14  T u r n e r , H.J., and Trimble, G.D., C a l c u l a t i o n of T r a n s i e n t Two-Phase Flow, I n v i t e d Paper presented a t the 1st CSNI Specialists Meeting on T r a n s i e n t Two-Phase Flow, Toronto, :.' August 1976.  15  Bathers, M.G., Zuzak, H.M., McDonald, B.H., and Hancox, W.T., On F i n i t e D i f f e r e n c e S o l u t i o n s t o the T r a n s i e n t Flow. Boiling Eguations, I n v i t e d Paper presented a t the 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow, Toronto, August 1976.  of J.  30  16  Karwat, H., BBOCH-S, A Code to I n v e s t i g a t e Blowdown of B o i l i n g Hater Beactor Systems, Hucl. Eng. Design 9, 363381, 1968.  17  Moore, K.V., and B e t t i g , W.H. , RELAP4-A Computer Program f o r T r a n s i e n t Thermal-Hydraulic A n a l y s i s , A e r o j e t Nuclear Company Beport ANCB-1127, 197 3.  18  Banerjee, S., F e r c h , B . L., and Mathers, I . G. , the Dynamics o f Two-phase Flow i n a Duct, Paper Submitted to 6th I n t e r n a t i o n a l Heat T r a n s f e r Conference, August 1978.  19  Ferch, B.L., Method of Characteristics Solutions for Noneguilibrium T r a n s i e n t Flow B o i l i n g , Submitted to I n t . J. Multiphase Flow, 1978-  20  Bryce, W.H., Determination of the H y p e r b o l i c Begime of the Hydrodynamic Equations Modelled i n the LOCA Code EELAP-OK, I n v i t e d Paper Presented at the 1st CSNI S p e c i a l i s t s Meeting On T r a n s i e n t Two-Phase Flow, Toronto, August 1976.  21  S t u h m i l l e r , J.H., The I n f l u e n c e of I n t e r f a c i a l Pressure F o r c e s on the C h a r a c t e r of Two-Phase Flow Model E g u a t i o n s , Int. J . Multiphase Flow, 3, 551-560, 1977.  22  Keenan, J.H., and Keyes, F.G., Steam, Wiley, N.T., 1958.  23  Owczarek,J.A., Textbook Co.,  Fundamentals p 330,  of  Thermodynamic P r o p e r t i e s of  Gas  Dynamics, I n t e r n a t i o n a l  1964.  24  Parkinson,G.V., The Short-Chambered Helium Besearch Board of Canada, T e c h n i c a l Memo.  25  Hancox, 8. T., and Nicoll, W. Westinghouse Canada L t d . , Beport, 1972.  Gas Gun, Defence TMAB-60, 1960.  B. ,  Unpublished  31  26  Carnahan, L.W.,  A p p l i e d Numerical Method,  Wiley, p 321,  1969.  27  F r i t z , John, "A Note on •Improper' Problems in Partial Differential E q u a t i o n s , " Communication on Pure and A p p l i e d Math., 7, 591-594, 1955.  28  fiichtmeyer, fi. D. And Morton, K. W. , " D i f f e r e n c e Methods f o r I n i t i a l Value Problems," I n t e r s c i e n c e , New York (1967).  29  Payne, L. E., Bounds i n the Cauchy Problem f o r the L a p l a c e Equation , Arch, f o r Bat. Mech. and Anal., 5, 3545 (1960).  30  S c h a e f e r , P. 8 . , On Uniqueness, S t a b i l i t y and Pointwise Estimates i n the Cauchy Problem f o r Coupled Elliptic Equations, Quart. Of Appl. Math., 321-328, 1973.  31  Lax, P. D., Differential and M a t r i x Theory, Comm. 194, 1958.  32  J a r v i s , Stephen J r . , "On the Formulation and Numerical E v a l u a t i o n o f a Set of Two-Phase Flow Equations Modelling the Cooldown P r o c e s s , " BBS Tech. Note No. 301, January 4, 1965.  33  Siegmann, E. B., Theoretical Study of T r a n s i e n t Sodium Boiling in Beactor Coolant Channels Utilizing a Compressible Flow Models" Argonne Nat. Lab., ANL-7842, J u l y 1971.  34  Boure, J . A., "GEVATBAN: A computer Program t o Study Two Phase Flow Dynamics, European Two Phase Flow Group Meeting, Home, CONF-720686-3, June 1972.  35  Boure, J . A., F r i t t e , A. , G i o t , H. And Reocreux, M. , Choking Flows and Propagation of S m a l l D i s t u r b a n c e s ,  Equations, Difference Equations Pure and Appl. Math., 11, p 175-  32  European 1973.  Two  Phase Flow Group Meeting, B r u s s e l s , Paper FL,  36  Gidaspow, D i m i t r i , Lyczkowski, B. W. , S o l b r i g , C. W. , Hughes, E. D. And Mortensen, G. A., C h a r a c t e r i s t i c s of Unsteady One-Dimensional Two-Phase Flow. American Nuclear S o c i e t y T r a n s a c t i o n s , 17, p 249-250, November 1973.  37  Boure, J . A., Dynamics Propagation de Petites 4456, October 1973.  38  In h i s writeup P r o f e s s o r Gidaspow shows t h i s analytically f o r the i n c o m p r e s s i l b l e l i m i t . Numerical r e s u l t s show t h a t the same i s t r u e i n the c o m p r e s s i b l e case; see Gidaspow e t . Al., Trans. ANS Winter Meeting (San F r a n c i s c o , November 1973) , p 249.  39  Hancox, W.T., and Banerjee, S-, Numerical Standards i n Flow B o i l i n g A n a l y s i s , Nuclear S c i . and Eng. 64, p 106-124, 1977.  40  Wallis, G. B., One-dimensional McGraw B i l l , *969.  41  S h a p i r o , A. H., The Dynamic And Compressible F l u i d Flow, V o l . 2, Ch.  42  Necmi, S. , and Hancox, W. T. , An Experimental and T h e o r i t i c a l I n v e s t i g a t i o n o f Blow-Down from a Horizontal Pipe, 6th I n t . Heat T r a n s f e r Conference, Toronto, 1978.  43  Bousseau, J . C , Blowdown Experiments and I n t e r p r e t a t i o n , European Two-phase Flow Group Meeting, H a i f a , I s r a e l , 1975.  44  Mathers, W. G., Hancox, W. T.,  on  des Ecoulements Diphasiques : Perturbations, Grenoble CEA-R-  Zuzak, 8 . Finite  -  Two-phase  Flow, Ch.  9,  Thermodynamics Of 24, Bonald Press Co.  W., McDonald, B. G., Difference Solutions to  and the  33  T r a n s i e n t F l o w - B o i l i n g E q u a t i o n s , Proceedings of the NEA Specialits Meeting on T r a n s i e n t Two-phase Toronto, August. 1976 0., the O i l and Gas J o u r n a l , V o l .  53, p.  OECD Flow  45  Baker, 1954.  175-195,  46  Banerjee, S., Hancox, H.T., J e f f r i e s , fi.B., and Sulatisky, M-I.,. T r a n s i e n t Two-Phase Flow d u r i n g Blowdown with Heat A d d i t i o n , AIChE Symposium S e r i e s , Heat T r a n s f e r V o l . 174, 1978.  47  Premoli, A., and Hancox, V.T., An Experimental Investigation of Subcooled Blow-down with Heat A d d i t i o n , I n v i t e d Paper presented a t the 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow.  48  Arrison,N.L., Hancox,H.T-, Sulatisky,M.T., and Banerjee,S. , Blowdown of a R e c i r c u l a t i n g Loop with Heat A d d i t i o n , Heat and F l u i d FlOw i n Water Beactor S a f e t y , I n s t i t u t i o n o f Mech. Engs. London, p 77-83, 1977.  49  Hadamard,J., L e c t u r e s on Cauchy's Problem i n L i n e a r D i f f e r e n t i a l Eguations, Yale D n i v e r s i t y Press, New Conn. , 19 23.  50  Boure,J.A.,and L a t r o b e , A., On well-Posedness of Two-Phase Flow Problems, Paper No. 76-CSME/ CSCHE-9 presented at the 16th National Heat T r a n s f e r Conference, St. Louis, M i s s o u r i , August 1978.  51  Bound Table Discussion on Two-Phase Flow, I n t e r n a t i o n a l Heat T r a n s f e r Conference, Tokyo, 1974.  52  Soo, S.L. , Multiphase Mechnics of S i n g l e Phase Flow, P h y s i c s F l u i d s , 20, p 568-570,  Partial Haven,  Component 1977.  Fifth  Two-  34  53  Chao,B.T., Sha,H.T., and Soo,S.L., On I n e r t i a l Coupling i n Dynamic E q u a t i o n s of Components i n a Mixture, Int. J . Multiphase Flow, i n press.  54  Boure, J.A., F r i t t e , A . A . , Giot,M.H., and Beocreux,M.L., H i g h l i g h t s of Two-phase C r i t i c a l Flow, I n t . J . M u l t i p h a s e Flow, 3, p 1-22, 1976.  55  Boure,J.A., On a U n i f i e d Presentation of the Equilibrium Two-Phase Flow Models, Non-Equilibrium Phase Flows ASME, N.Y., p 1-10, 1975.  56  T r u e s d e l l , C . , and Toupin,B., The C l a s s i c a l F i e l d Theories, Handbuch der P h y s i k , V o l . 3/1, Springer V e r l a g , B e r l i n , 1960.  57  Standard,G., The mass, momentum, and energy heterogeneous flow systems, Chem. Engng. 1964.  58  S o l b r i g , C . H., McFadden, J.H., Lyczkowski,B. V. , and Hughes,E.D., Heat Transfer and F r i c t i o n Correlations Beguired to D e s c r i b e Steam-water Behaviour i n Nuclear Safety Studies, Paper Presented a t 15th N a t i o n a l Heat T r a n s f e r c o n f e r e n c e , San F r a n c i s c o , 1975.  59  Banerjee,S., A S u r f a c e Benewal Model f o r I n t e r f a c i a l Heat and Mass T r a n s f e r f o r Two-Phase Flow, I n t . J . M u l t i p h a s e Flow, i n p r e s s .  60  Yanenka, N.M., The Method o f F r a c t i o n a l Time Steps: The Solution of Problems of Mathematical P h y s i c s i n S e v e r a l V a r i a b l e s , ed. M- H o l t , S p r i n g e r - V e r l a g , N. Y., 1971.  61  Marchuk,G.I., A T h e o r e t i c a l Heather F o r e c a s t i n g Model, Dok. Akad, S c i . USSB, 155, p 10-12, 1965.  NonTwo-  equations f o r S c i . , 19, p 227,  35  62  Douglas, J . , A l t e r n a t i n g D i r e c t i o n Methods f o r Three-Space V a r i a b l e s , Numerische Hathematik, 4, p 41-63, 1962.  63  Lax,P.D., Beak S o l u t i o n s of Nonlinear H y p e r b o l i c Eguations and T h e i r Numerical Computation, Comm. Pure A p p l . Math., 7, p 159-193, 1954.  64  Lax,P.D. , and Hendroff,B., Difference Order of Accuracy f o r S o l v i n g H y p e r b o l i c Pure ApplMath., 17, p 381, 1964.  65  Evans, M. E. , and Harlow, F. H. , The P a r t i c l e - i n - C e l l Method f o r Hydrodynamic C a l c u l a t i o n s , LASL Beport No. LA-2139, N. Mexico, 1957.  66  Gentry,B.A., Martin,B.E., and Daly, B.J., An E u l e r i a n Differencing Method f o r Unsteady Compressible Flow Problems, J . Coop. P h s i c s , 1 , 87-118, 1966.  67  HcGormack,B.W., Numerical S o l u t i o n of the I n t e r a c t i o n o f a Shock . Have with a Laminar Boundary Layer, i n Proc. 2nd I n t . Conf. on Num. Methods i n F l u i d Dynamics, ed. By M. H o l t , V o l . 8, S p r i n g e r - V e r l a g , N.Y., 1971.  68  Porching, T.A., Murphy,J.H., and B e d f i e l d , J.A., S t a b l e Numerical Integration of Conservation Eguations f o r H y d r a u l i c Networks Nuc. S c i . Engng. 43, p 218-225, 1971.  69  Mulpuru, S.B., and Banerjee, S., An I m p l i c i t Method f o r Shock-Have C a l c u l a t i o n s , submitted t o AIAA J o u r n a l , 1978.  70  Beichenback,V.H-, and D r e i z l e r , H., Ober den Querschmitt Sanderungen and G i t t e r i n Kanalen anf S t o s s w e l l e n , Z e i t . Agewandte P h y s i k , 12, p 62-7 1, 1960.  71  Von Neumann, J . , and Bichtmeyer, B.D., A Mehtod f o r the Numerical Calculation of Hydrodynamic Shocks, J . Appl.  Schemes with High E q u a t i o n s , Comm.  36  P h y s i c s , 21, p 232-257, 1950. 72  McDonald, B.H., Hancox,8. T. , and Mathers, H.G., Numerical S o l u t i o n s of the T r a n s i e n t F l o w - B o i l i n g Equations, Invited Paper Presented a t 2nd CSNI S p e c i a l i s t s Meeting o T r a n s i e n t Two-phase Flow. P a r i s , June 1978.  73  Courant,B., I s a a c s o n , E., and Bees, H., On the S o l u t i o n of N o n l i n e a r H y p e r b o l i c Equations by F i n i t e D i f f e r e n c e s , Comm. Pure Appl. Math., 5, P 243-255, 1952.  74  Curtis, A.C., and B e i d , J.S., FOBTRAN Subroutines f o r the S o l u t i o n of Sparse S e t s of L i n e a r Equations, AERE-B-6844, U.K. Atomic Energy Research E s t a b l i s h m e n t , Harwell, 1971.  75  Harlow,F.H., and Amsden, A.A., Numerical C a l c u l a t i o n of Almost Incompressible Flow, J . Comp. P h y s i c s , 3, p 80-93, 1968.  76  Hirt,C.W., and Oliphant,T.A., S0LA-PLO0P: A Non-Equilibrium Drift F l u x Code f o r Two-Phase Flow i n Networks, I n v i t e d Paper Presented a t t h e 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow, Toronto, August 1976.  77  Graf, 0., Bomstedt, P., and Merner,H., A p p l i c a t i o n o f t h e ASHE Method t o Two-Phase Flow Problems, I n v i t e d Paper Presented a t 2nd CSNI S p e c i a l i s t s Meeting on T r a n s i e n t TwoPhase Flow, P a r i s , June 1978.  78  Bomstedt, P., and ierner,W.E., E f f i c i e n t High-Order Method f o r t h e S o l u t i o n of F l u i d E g u a t i o n s , Nuclear S c i . and Eng. 64, p 208-218, 1977.  79  Biegel, B., Marechal, A., and Bousseau, J.C., Depressurization d'une capacite en double phase i n s t a l l a t i o n CANON, Center D'Etudes N u c l e a i r e s de Grenoble, Note TT-490, 1975.  37  80  Borgartz,B.0., Goodman, B.M.E., O'Brien, T.P., Bawlingson, H., and Edwards A.R., D e p r e s s u r i z a t i o n S t u d i e s , Phase 2: B e s u l t s of T e s t s 115 and 130, UKAEA Beport SBD-B-115, 1978.  81  I b i d , D e p r e s s u r i z a t i o n S t u d i e s , Phase 3: B e s u l t s 144 and 145, OKAEA Beport SBD-B-77, 1977.  82  Ibid, Depressurization S t u d i e s , Phase 3: B e s u l t s of T e s t s 142 and 143, OKAEA r e p o r t SBD-B-76, 1977.  83  White,E.P., and D u f f e y , B.B., A Study o f t h e Unsteady Flow and Beat T r a n s f e r i n t h e B e f l o o d i n g of Water Beactor C o r e s , Central E l e c t r i c i t y Generating Board Beport BD/B/N 3134. September 1974.  84  Reocreux, H.L., C o n t r i b u t i o n a 1'etude des d e b i t s c r i t i q u e s en ecoulement diphasique eauvapeur, Thesis, Univ. S c i e n t i f i q u e et Medicale de Grenoble, 1974.  85  Reocreux,M.L., Experimental Study of Steam Water Choke Flow, i n v i t e d Paper Presented a t the 1st CSNI S p e c i a l i s t s Meeting on T r a n s i e n t Two-Phase Flow, Toronto, August 1976.  of  Tests  APPENDIX A  DERIVATIONS  39  DERIVATION OF EQUATIONS USEFUL FOB  THE CHABACTEBISTICS:  I t i s d e s i r e d t o t r a n s c r i b e the f o l l o w i n g governing e g u a t i o n s i n t o a new  set of eguations useful f o r  a p p l i c a t i o n of the c h a r a c t e r i s t i c dB/dt+(B/A) (4Au/3z)  uni-directional  directions:  = 0,  (Continuity)  du/dt • (3P/tfz)/R = — F - g(dZ/dz),  (Momentum)  dh/dt - d(P/R)/dt • (P/BA) (3Au)/3z) = g • uF, where  q  is  J/Kg/Sec, and  relation B = B(P,h), P,h  is  N/Kg  (because d e n s i t y i s  i n the compressed  liguid  the speed of sound a s above s e t may  F  1  or M/Sec*. fairly  Using the  independent  of  s t a t e ) , and using an e x p e r e s s i o n f o r  a* = 1./  be developed  (Energy)  ( 3E/3P|h + (1. /B) 2>B/dh | P ) , the  i n t o a d i f f e r e n t s e t of e g u a t i o n s .  A p p l i c a t i o n of the c h a i n r u l e g i v e s : dB/dt = 3B/ap|h  (dP/dt) • 3B/ah|P  where from the speed of sound (dB/iP)Jh = l . / a z -  (dh/dt),  (A.1)  relation,  (1./B) (ifi/3h) |P  (A.2)  and (3B/aP)|h  =  (1./a ) z  (1./B) (SB/ah) |P ) .  (A.3)  Then dB/dt =  |1./a  2  - (1./B) (»B/Sh) |P) (dP/dt) • (3R/»h)|P  . dh/dt.  Upon s u b s t i t u t i o n f o r dh/dt from the energy e g u a t i o n ,  1 3E/aP|h = p a r t i a l of d e n s i t y with r e s p e c t to pressure a t constant enthalpy  40  we get: dB/dt = ( 1 - / a * 2  (1./B) ( a f i / d h ) | P)  (dP/dt) + dB/dh| P  •[ q*uF -  (P/EA) ) (3Au/dz) «• (1./B) (dP/dt) - (P/B ) (dB/dt)] . Now  collecting  the  coefficient  terms f o r dB/dt, and dP/dt, i t  follows: (dB/dt) ( H (P/B ) (SB/ah) |P) = (1./a2) (dP/dt) • (dB/3h) | P 2  [q+uF-  (P/EA) 3Au/3z) 3Now  s u b s t i t u t i o n f o r dB/dt from the c o n t i n u i t y  eguation  results:  -(B/A) (3Au/az) (1.*(P/B ) ( a f i / a h ) | P ) = 2  (1/a*) (dP/dt) «• (d£/3h)|P  r L  g*uF-(P/BA)) (3Au/az) .  Opon c a n c e l l i n g the e g u a l -(P/BA) (3Au/2z) (dB/3h)P terms from both s i d e s , we a r r i v e a t t h e f o l l o w i n g :  — (B/A) (aAu/Sz) = (1-/a?) (dP/dt) + (3B/3h) |P[q + uF], or  - (Bu/A) ( d A / a z ) - B 3 u / a z or  finally:  =  (1./az) (dP/dt)  •  (3B/ah) | P  Qq*uF],  1  41  Equation  (A.4) i s t h e new form of the energy e q u a t i o n .  Hext t o o b t a i n the new form of the c o n t i n u i t y e q u a t i o n , we s u b s t i t u t e f o r dP/dt from the energy e q u a t i o n : d B / d t = (SR/aP) |h dP/dt + 3fl/ah|P dh/dt = ( 3 B / a P ) I h dP/dt • B (1-/a?-aB/aP|h) dh/dt , or dB/dt = ( a f i / a P J h { B (dh/dt)-B£q«-uF-(p/fi2) (dB/dt)-(P/BA) ( 3 A u / a z ) ] | +B(1-/a --dB/aP|h) dh/dt. 2  Collecting  c o e f f i c i e n t s of dfi/dt and dh/dt terms g i v e s :  d E / d t J.1. + (P/B) ( a f i / a P | h ] = B ( 1 . / a 2 - a B / S P Jh+aa/dP|h)dh/dt B(aB/aP|h) £q + uF-(P/RA) (aAu/2z)]Now s u b s t i t u t i n g from c o n t i n u i t y e g u a t i o n f o r dfi/dt g i v e s : -(B/A) (3Au/az) (1. • ( P / E 2 )  a£/ah|P)  = B / a 2 . dh/dt -R (3R  / a P | h ) [ q + uF-(P/RA) (3Au/az)3# or  (1./a ) dh/dt • 1/A 2  (cJAu/az)  = afi/3P|h [q + uFJ  or dh/dt • a. a u / a z = a  2  (3R/aP|h) [q+uF] -u/A a A / S z " ^ .  (A. 5)  (A.5) i s t h e new form of the c o n t i n u i t y e q u a t i o n . Momentum  equation  remains  the  same  as  before.  Thus  c o n t i n u i t y , momentum, and energy e q u a t i o n s a r e r e s p e c t i v e l y :  42  dh/dt+a  2  du/dz = a  2  (A.5)  dE/dP|h [g+uF]-u/A «A/dzJ = C2  (A.6)  du/dt*1./B aP/Sz = -F-gdZ/dt= C l dP/dt * Ra  A  du/az = -a I>B/ah| p[g+up3+ (Bu/A) a v 4  a z j = C3  C1, C2, and C3 are c a l l e d the c o n s t a n t s of i n t e g r a t i o n .  (A.4)  The  l e f t hand s i d e of the above h y p e r b o l i c eguations are the p a r t i a l derivatives  of u, p, and h with r e s p e c t t o space and time  These equations may directions.  be r e a d i l y t r a n s l a t e d  into  (z,t).  characteristic  APPENDIX B  DERIVATION OF THE OBDINABY DIFFERENTIAL EQUATIONS  4U  DERIVATION OF THE ORDINARY  DIFFERENTIAL EQUATIONS:  In Appendix A the f o l l o w i n g s e t of h y p e r b o l i c equations were derived: dh/dt • a 3u/3z = C2  (Continuity)  du/dt+ (3P/az)/fl = C1  (Momentum)  2  dP/dt • Ra* 3u/3z = C3  (Energy)  From here one may pursue t o f i n d  the o r d i n a r y  differential  eguations along the c h a r a c t e r i s t i c s . The above can be w r i t t e n a s : h',t  • uh«,z • a*u»,z = C2  u»,t  • uu«,z + P',z/R = C1  P«,t  + uP«,z tRa - 3u/3z = C3. 2  Now t r a n s f o r m a t i o n coordinates  s (along  o f above t o another  s e t of  orthogonal  which the p r o p e r t i e s and t h e i r d e r i v a t i v e s  are known) and y, (the c h a r a c t e r i s t i c s d i r e c t i o n ) , r e s u l t s the followings: k'#y y'#t • u h , y ,  y',z • a u*,y y',z • h*,s s»,t * uh*,s s',z 2  • a*u" s s #  1 #  z = C2  u*»y y'#t • uu',y y',z • P',y y ,z/R • u»,s s , t • uu',s s',z #  (  • p #y V,z/z = c i f  P'#y y*»t • uP«,y y»,z • Ba*u«,y y*,z • P« ,s s»,t • uP»,s s«,2 • fia^.u s s*,z = C3. f  f  i h» , t = d h / a t , e t c .  45  These e q u a t i o n s  can be w r i t t e n i n matrix  !"aV#z  form as f o l l o w i n g :  f*,t*uy«,z  u' ,y  y,t*uy,z  (y',z)/fi  o  P',y  Ba^y •, z  y' t*uy',z  0  h«,y  #  Since u>,y  , P',y  , h ,y ,  = Known.  (the d i s t u r b a n c e s across the  c h a r a c t e r i s t i c s ) are i n d e t e r m i n a t e , then i n order t h a t the above is  true,  the  the determinant  determinant  of matrix  must vanish. Expanding  gives :  ( y ' , t + u y , z ) [(ySt+uy^z)*-a -(y«,z) ] 2  From dy = 3y/3z dz+dy/dt dt = 0 ,  2  =  0.  i t can be w r i t t e n :  y'#t/y«,z = -dz/dt = -u,  (3.1)  or (B.2)  y'#t/y« z = -dz/dt = -(u±a), #  Eguations characteristic  (B. 1)  and  (B-2)  are  the  S u b s t i t u t i o n from  of c o e f f i c i e n t s r e s u l t s the r e s p e c t i v e l e f t way:  eigenvalues  d i r e c t i o n s . L e f t e i g e n v e c t o r s can be  using these e i g e n v a l u e s .  following  set  or  obtained  (B. 1) i n the  matrix  e i g e n v e c t o r i n the  46  £61  62  63] a*  1-/B  Ba  = 0,  0_  2  or expanding, 61 a« • Ba« 63 = 0 62  (1./H) - 0 0=0,  and the e i g e n v e c t o r can he chosen as Siailarly [6 1  62  [-B  0  1],  f o r dz/dt = - (y» , t) / (y ,z) = u+a: 1  63]  -a  -a  1-/B  Ba*  -a  = 0,  or expanding the above g i v e s :  61 a* - a6 2 • Ba - 63 = 0, 2  62 / f i - a 63 = 0,  -a 51 = 0.  and  The r e s p e c t i v e l e f t e i g e n v e c t o r c o u l d be Co Similarly the  Ba  is  i.0 -Ha  1],  I t i s d e s i r e d now t o transform the governing  continuity  (B.3)  f o r dz/dt = - (y* , t ) / ( y ' ,z) = u-a,  eigenvector  i n t o these  l"3«  directions.  For  (B.4) eguations  the dz/dt = u*a, m u l t i p l y the  equation by *0*, the momentum equation by  •Ba•,  and the energy equation by * 1•, we g e t : Ba Iau/3t *udu/dz +(SP/dz)/S J •  aP/az*u(dP/az)•  B a ^ a u / a z = C3*BaC1, or ap/at+(u*a) a p / a z  * a a £ a u / a t + (u+a>au/az ] = C 3 * a a d  or dF/dt+Ba du/dt = C3*BaCl. Similarly  (B. 5)  along dz/dt = u-a:  dP/dt - Ba du/dt = C3-BaC1  (B.6)  For the dz/dt = u d i r e c t i o n  -B Qh/dt*a 3u/az3I 2  , we s i m i l a r l y w i l l have:  • d P / d t * B a d u / a z = -BC2+C3, 2  48  -B ah/at-Buah/az-Ba^au/az+ap/at-map/az  • B a a u / a z = -BC2+C3, 2  or - B ( a h / a t + u a h / a z ) * a p / a t + u a p / a z = -BC2*C3,  or  finally  dP/dt-Bdh/dt = C3-BC2. The equations equations  (B.5)-(B.7)  (B-7) a r e s e t of o r d i n a r y  (called the c o m p a t i b i l i t y  relations).  differential These  equations can be w r i t t e n i n i n c r e m e n t a l form u s e f u l f o r f i n i t e difference  numerical  computations.  APPENDIX C  SOLUTION TO ISENTBOPIC PROBLEM  50  ISENTROPIC SOLUTION: S i m i l a r but  s i m p l e r a n a l y s i s than t o r a n o n i s e n t r o p i c  (Appendix £) can for the  be done to a r r i v e a t c o m p a t i b i l i t y  an i s e n t r o p i c problem. For  the  problem  relations  case of an i s e n t r o p i c problem,  governing e q u a t i o n s a r e :  ( a £ / a t ) / f l • aR/az = 0 du/dt •  (Continuity)  (dP/fiz)/R = 0  The  equation  sound and  (Momentum).  o f s t a t e i s : ^ dP/aRJS=a ' 2  (C.2)  (where a i s speed of  S, the entropy) .  Thus ap/az =  (ap/afl) (aR/az) = a  Then, the e q u a t i o n s  (C.1)  and  #  R  aR/az.  2  (C.2)  R » t • u R«,z  It follows  ( C 1)  be  • R R',z  (u» t • u u»,z) #  may  • a  2  =  written  0,  fi«,z =  0.  from s i m i l a r a n a l y s i s as p r e v i o u s l y  1 aP/aR|S = p a r t i a l of pressure constant entropy 2 R»,t = d f i / a t , e t c .  as:?  with r e s p e c t t o  (Appendix fi) t h a t  density  at  51  (a/B)£R' t *  (u±a)  #  R\z2 i{,a*rt  «• (u±a)  = 0,  a*rz]  or  (a/B) (dB/dt) Hence (a dB)/B  and  du/dt = 0  along  • u = 2r along dz/dt =  (a dfi)/B -u = 21 along dz/dt = Por p o l y t r o p i c gases  s p e c i f i c heat  dz/dt = u i a .  a dB/B  u+a,  u-a. = 2a/(k-1),  where  is  the  1 are c a l l e d the Biemann I n v a r i a n t s and s i n c e  they are c o n s t a n t along each c h a r a c t e r i s t i c , one  may  c o o r d i n a t e s and i n v e r t the space and  to express as z = z ( r , l ) and t = t ( r , l ) . dz = 3 z / 3 l dl= (u+a) 3t/31 Similarly Then  k  ratio.  The terms r and  to be the new  (C.3)  dl 3z/3r =  assuming  or  time c o o r d i n a t e s  Then from dz =  3z/3l =  (u-a)  c o n s i d e r them  (u+a)  (u*a) dt,  St/31.  (C.4)  3t/3r  (C.5)  t h a t both z, and t are continuous f u n c t i o n s  i n r and 1, the eguations  (C.4)-(C.5)  may  be  combined  follow:  z«,lr =  (u#a)«,l t»,l  • (u*a) t ' , l r  z»,rl =  (u-a)«,l t»,r  • (u-a)  t«,rl.  as  52  S e t t i n g the above e q u a l , we g e t : 2a t»,rl • (u+a)«,l t«,l or  since  A l f a • u = 2r,  Alfa = r * l  where  and  ( u - a ) * , l t»,r = 0,  Alfa-u=21,  A l f a = 2a/(k-1) ;  then u = r - l , and  u = r-1  and  a= (k-  1) ( r * l ) / 2 . Then  2a t • , r l [ - 1* (k- 1) /2] t« ,1 - [ l - (k-1) /2]t« ,r = 0,  or 2a t»,rl * (t«,l • t * , r ) (k-3)/2 = 0, or t«,rl • L ( t * , l where L * For equation  • t«,r)/(r+l) = 0, (C.6)  (k + 1)/2(k-1) and k = (2L+1)/(2L-1) . a i r , L=3, and f o r He, L=2.  The s o l u t i o n t o the  ( C . 6 ) i s hyper-geometric f u n c t i o n s .  The hyper-geometric  s o l u t i o n f o r t h e helium-gas gun i s used to check the numerical r o u t i n e used f o r the present work.  APPENDIX D  ANALYTIC SOLUTION TO NONISENTHOPIC PBOBL  54  HONISENTBOPIC PROBLEM: Initially  , i t was  desired  to  the n o n i s e n t r o p i c problem along problem. as b e f o r e  the  O b t a i n exact r e s u l t s f o r same  To do t h i s , the c o m p a t i b i l i t y  l i n e s as i s e n t r o p i c  relations  were w r i t t e n  (as d e r i v e d i n Appendix B) :  dP*Hadu = (C3+BaC1)ds  along dz/dt=u+a  dP-Badu = (C3-BaC1)ds dP-Bdh = (C3-BC2)ds  along dz/dt=u-a along  I n t e g r a t i o n was then performed  dz/dt=dt=u.  f o r both  (D-1) (B.2) (D.3)  s i d e s of t h e above equations.  Renaming A l f a = ^ R a d u , Beta =^Bdh, 2n = P+Alfa, 21 = P-Alfa,and m = P-Beta, we f i n d t h a t P - n+1,  A l f a = n - l and Beta - n«-l-m.  From there s o l v i n g f o r n, n, and 1, we g e t :  2n =  #  55  Now t r a n s f o r m i n g z and t t o a,n, and 1 d i r e c t i o n s ,  gives:*  z',1 =  (u+a) t«,l  or  z* ,ln =  z»,n =  (u-a) t«,n  or  z*,nl = (u-a)»,l t»,n • (u-a) t * , l n  z',m  = u t*,m  or  Next t a k i n g respect  Z',DD  derivative  t o a,m,l, r e s p e c t i v e l y ,  z*,lnm = (u*a)  z' ,nlm =  the  '  f  n i B  (u-a) • ,1m  (u+a)' ,n t»,l • (u*a) t« , l n  of  the  = U',n  t'a • u  t',mn  above eguations with  we g e t :  t» ,1* (u*a) •,n t ' , l B t ( u t a )  t • , l n * (u*a) t',lnm  t , n • (u-a) »,1 t«,nm • (u-a) «,m t» , l n • (u,  a) t*,lnm  z*,nnl = u»,nl t « , a • u«,n t , " ! • u',1 1  t',»a •• u t«,anl  S e t t i n g once the r i g h t hand s i d e of the f i r s t the  above e g u a t i o n s e g u a l  and next the f i r s t  and second o f  and t h i r d , we g e t :  2a t«,anl • 2 a»," t»,ln + (u*a) »,n t ' , l a - (ua) »,1 t*,mn + (u*a)',an t«,1 -  i  z ' , l n = a*z/aian , e t c .  (u-a) *.,•! t * , n  = 0,  56  or  c o l l e c t i n g terms, i t f o l l o w s : a t',mnl • a»,n t , l m • (uia)',niii t * ,1 f  u * , n l t»,»  +  (u*a)«,m t ' , l n - u«,l t ,mn ,  = 0.  Form there c a n c e l l i n g the t ' , o n l term g i v e s :  (u-a)»,n t* ,1m  • 2 u*,m  t ' . l n * (u+a)«,l t » , m n -  (u+a)*,an t',1  - ( u - a j ' ^ m l t , n •-2 u« , n l t*,m = 0. #  There i s no hope cf s o l v i n g the above as long as having undefined terms l i k e  (u-a) n, etc. ,  #  APPENDIX E  SPEED OF SOU  58  DERIVATION OF THE USEFUL FORM OF SOUND SPEED: I t i s more u s e f u l t o d e f i n e a=a (P,T) i n s t e a d o f a=a(P,h), because, t h e n , t h e property chain  t a b l e s can r e a d i l y be a p p l i e d . Using  r u l e r e l a t i o n s , the f o l l o w i n g s can be w r i t t e n :  da=3h/aP|T dP*dh/dT|P dT  (E.1) (E.2)  dv=3v/dP|T dP+av/aTJP dTThen i r o n  (E.2) ,  av/apjh=av/ap|T*av/aT| p a i / a p j h . But  from  (E.1)  and a l s o from  3T/aP|h=-(dh/aPlT)/(3h/aT|P) ,  (E.2)  Terms such as dB/dP)h meaningful way.  (E.3)  av/3h|P = (2v/aT 1 P) / (3h/3TJ P) .  He may w r i t e : (E.6  2  3R/3hJP = d(1/v)/ah|P = - (Sv/dh | P)/v : ,  (E-7)  2  since  (E.5)  and dB/dhjP need t o be d e r i v e d i n a more  3R/3PJh = 3(1/v)/aP|h = - ( a v / a P | h ) / v , and  and  (E.4)  v=x vg • (1-x) v f ,  (E.3), (E.4), (E.5),  then  (E.6), and (E-7),  aR/ap|h = -{ • x avg/ap|T+(i-x)  combining  equations  we g e t :  avf/ap|i  j - i ax  a v g / a u p • (1-  x) a v f / d T J P ] [x S h g / a p j h * (1-x) a h f / 3 P i h ] / i x a h g / a T | P + ( 1 x) 3 h f / 3 T | P ] j / v  2  and 3R/3h|P = - { I x a v g / S T I P+(1-x) 3 v f / 3 T | P ] / [ x x)  ahf/arj P3]/V .  where the q u a l i t y i n terms or v o i d  2  fraction i s :  x = A l / [ A 1 * {1-A1) v f / v f ] , and the void  f r a c t i o n i n terms of q u a l i t y i s :  A l = x/£ x* (1-x) vf/vg J . \ •  3hg/dTJ P  APPENDIX F  THE COMPUTER PROG BA M  61  The proolem  program was  470/V6-II  IIPHAS  written i n computer.  for  handling  Fortran  two-phase  flow  Language and run  Double p r e c i s i o n  was  used  on for  water Amdahl better  accuracy. A t y p i c a l run complete  for  this  depressurization  program  of  pressure  of  1.8  debugged  with help of I n t e r a c t i v e Program ( I F ) .  used s,  u n i t s MKS  atm.  velocity  temperature  in  m/s,  in  C,  ( i . e . meter,  and  viscosity,  J/kg B  and  water  saturation  the program  kilogram,  in  density in  at  for  was  mainly  second)  were  That means, length i n m, time i n  pressure  p r o p e r t i e s were a l s o used absolute  d i a m e t e r ) , tube  Initially,  throughout the a n a l y s i s .  seconds  deg. C high e n t h a l p y  (.032  System  in  400  f i l l e d , 4 meter long about  m  120  took  in  metric for  pa,  enthalpy kg/m .  J/kg,  Thermophysical  3  units  in  (i.e.,  kg m/s  for  c o e f f i c i e n t of conductivity ,  etc.) A t h r e e d i m e n s i o n a l a r r a y Z(I,K1,K2) s t o r e d a l l c a l c u l a t e d n o d a l v a l u e s throughout the run and a r r a y values were saved i n a file  for  l a t e r use.  1=1,  ... ,10 i n d i c a t e d v a r i a b l e s z, t , u,  P, h, T, B, a , x, and void f r a c t i o n . reflective  K1  was  waves from c l o s e d end of the tube.  time  for  K2 was the space  index r e p r e s e n t i n g the encounter of r e f l e c t e d wave end with on-coming subsequent waves.  index  from  closed  62  DESCRIPTION OF BOOTINES:  HAIN  Contained body of the program. Routine "FAN" c o n t a i n e d in this r o u t i n e i n i t i a l i z e d the expansion waves emanating from open end a t time zero. Other r o u t i n e s c o n t a i n e d i n MAIN were START and EVAL.  START  Read file.  EVAL  N-R i t e r a t i v e r o u t i n e used to evaluate p r o p e r t i e s using SIHUL, the Nonlinear Equation S o l v e r .  IIPHAS  Set up c o e f f i c i e n t matrix f o r i n t e r m i t t a n t waves.  OPEN  Set up c o e f f i c i e n t matrix f o r open end c o n d i t i o n . This d i f f e r e d from the IIPHAS one only i n c o n t a i n i n g the choking c o n d i t i o n f o r the open end.  CLOSED  Set up c o e f f i c i e n t  SIMOL  Newton-fiaphson N o n l i n e a r Equation S o l v e r . Having r e c e i v e d the coefficient matrix A, i t r e t u r n e d the s o l u t i o n matrix X.  RE  Evaluated Reynolds number.  TASTAR  Evaluated force.  FL  E v a l u a t e d the r i g h t hand s i d e constant of the equation of left-running c h a r a c t e r i s t i c .  initial  values  generated  by  Routine  FAN from a  matrix f o r c l o s e d end c o n d i t i o n .  two-phase  multiplier  for  the  frictional  Evaluated the right hand e g u a t i o n of p a r t i c l e path.  side  constant  E v a l u a t e d the r i g h t hand s i d e constant of riqht-running c h a r a c t e r i s t i c . E v a l u a t e d temperature from the enthalpy l i q u i d region.  f o r the  f o r the  equation  in  compressed  64  C C C C C C  C C 50 0 C C  C  8 4 101 9  C  6 100 C  7  IMPLICIT HEAL*8(A-H,0-Z) DIMENSION A (30, 30} ,Z(10, 8,30) ,XOLD (30) ,XINC (30) THIS PBOGBAM TO ADVANCE THE SOLUTION OF TRAVELLING* PBESSURE MV1S F.BOH THE RUPTURE END OF HIGHLY PRESSURIZED 4 H LONG TUBE BY THE METHOD OF CHARACTERISTICS. THE N-RAPSON NONLINBAB EQUATION SOLVES IS USED TO OBTAIN XO THE SOLUTIONS. THE ROUTINE EVAL HANDLES THIS. K5=1, CALL IIPHAS, =2, THE IIPHAS, =3 THE CLOSED, AND = 4 , THE OPEN. ITIME=30 IVAR=10 ISP=8 ITI=ITIME ND=30 NN=ISP NS=ISP €1=1643.60+0 X=4.D+0 M=2 N=NS-t CALL STABT|Z,XOLDvIVAl,ISP,ITI,ND,X,NH) K1 INCREMENT IS IN TIME AND K2 ALONG THE PIPE. CONTINUE DO 50 K1=2,ITIME ...K3 IS A SWITCH. , K3=0 MEANS EVEN TIME ADVANCE AND K3=1 O ... CALLING THE CLOSED END ROUTINE. CALL EVAL{A,Z,X0LD,XINC,K1,1,K3,3,IVAR,ISP,ITI,ND) K2= 1 K5=3 DO 8 1=1,5 Z (I, 1,K1) = XOLD(I) WRITE {6, 101) ( Z (I1,1,K1) ,11=1, 1X>) yKt,K 2, K5 FORM AT (6G1 4 . 7 / 4 G 1 4 . 7 , 313) .... USING IIPHAS 10UTINE TO CALCULATE^ THE INTERMITTANCE POI DO 100 K2=M,N ....CALLING THE IIPHAS-*. CALL EVAL(A,Z,XOLD,XINC,K1,K2,K3,1,1VAB,ISP,ITI,ND) DO 6 1=1,5 Z (I,K2,K1) =XOLD fI+5) K5= 1 WRITE(6, 101) ( Z (I1,K2,Kl)iIT=1, 10) ,K1,K2,K5 CONTINUE K2=NN -...ON THE BOUNDARY THE OPEN ROUTINE IS CALLED.... CALL EVAL f A,Z,XOLD,XING, K1,K2,K3,4 ,IVAR,ISP,ITI,ND) K5=4 DO 7 1=1,5 Z(I,K2,K1)=XOLD(I+5) W RIT E ( 6 , 10 1) (Z (11, K 2 , K1) , 11 = 1, 10) , K1, i K 2 , K 5 CONTINUE STOP END :  50  65  C C C C  SUBROUTINE START (Z,XOLD,IVAR,ISP,ITI,ND,X,NN) IMPLICIT '.REAL*8 (A-H,0-Z) ......INITIALIZES THE XOLD AND Z (I,K2,K1) ARRAYS. K1=1. 1=1 IS Z, 1=2 TIME, 1 = 3 U, 4 P, AND 5 H,6 TEMP, 7 DENSITY,8 SOUND SPEED,: 9 THE QUALITY, AND 10 THE VOID FRACTION. K2 IS THE SPACE NODE ALONG THE TUBE. DIMENSION Z {IVAfi,ISP,ITI) ,XOLD{ND) READ(5,101) { (Z(11,KK2,1) ,11 = 1,IVAR) ,KK2=1,ISP) 101 FORMAT (6G14.7,/,4G14.7) WRITE (6 , 102) { {Z ( 1 1 , KK2, 1) , 11=1, IVAR) ,KK2 = 1,ISP) 102 FORMAT (6G14.7,/, 4GT4. 7) DO 6 1=1,25 6 XOLD{I) =0.D0 RETURN END  66  C C C C  C  SUBROUTINE EV1L (1,2,XOLD,.XING,KT, K 2 ,K3, K5,1?HE,ISP,ITI, ND) IMPLICIT RE1L*8(A-H,0-Z) .....K5=1 CALLS THE IIPHAS ROUTINE, =2 THE IIPHAS, =3 THE -.CLOSED, AND =4 THE OPEN ROUTINE. IT AUTOMATICALLY STORES THE NEW VALUES OF TEMP AND VOID IN Z(6,K2,K1) AND Z (10, K2 ,K1) RESPECTIVELY. DIMENSION A (ND, ND) ,Z (IVAR ,ISP,ITI) , XOLD (ND) , X I H C (80) EPS 1=1-00-10 EPS 2= 1.0 DO ITMAX=6 - . ... .THE ITERATION OF CONVERGENCE. DO 9 I=1,ITMAX IF (K5- NE. 1) GO TO 1 1 CALL IIPHAS<A,Z',X0LB,Kl,K2,O,I,I?AR,ISP,ITI,ND,TEaP3 1>BH03,A3,X3,V0ID3)' N=10 GO TO 14 11 IF(K5.NE.2)GO TO 12 CONTINUE 12 IF (K5. NE. 3) GO TO 13 CALL CLOSED(A,Z,XOLD,K 1, K2,K3,I,IVAR,ISP,ITI,ND,TEMP3 1,RH03, A3,X3,VOID3) N=5 GO TO 14 13 IF (K5.NE.4) GO TO 15 CALL 0PEN(A,Z,X0LD,K1,K2,K3,I,IVAR,ISP,ITI,ND, 1TEMP3, RH03, A3,X3, V0ID3) N=10 GO TO 14 : 15 CONTINUE CALL OPEN1(A,Z,XOLD,K1,K2,K 3,1,IVAR,ISP,ITI,ND, iTEMP3,Rfl03,A3,X3,VOID3) N=5 14 CONTINUE — ..SIHUL IS THE G-S NONLINEAR SOLVER-,-. ' DETER=SIMUL(N,A,XINC,EPS1,1,ND) IF (DETER.NE.O.D+0) GO TO 3 GO TO 1 3 ICONT=1 DO 5 J=1,N . . v . .THE TEST OF CONVERGENCE... IF (DABS (XINC (J) ) -GT.EPS2) ICONT=0 5 XOLD(J)=XOLD(J)+XINC(xJ) IF (K5. EQ.3) GO TO 16 16 CONTINUE IF (ICONT-EQ. 0) GO TO 9 GO TO 1 9 CONTINUE 1 IF (R5.NE.2)GO TO 2 —-.ASSIGNING COMPUTED VALUES TO THE ARRAY Z. 2 Z(6,K2,K1)=TEMP3 Z (7,K2,K1) =RH03 Z (8,K2,K 1)=A3 5  C  C  C  2 (9,K2,K1) =X3 Z (10 K2,KT)=VOID3 IF(K5„NE.3) GO TO 1 GO TO 6 10 CONTINUE 6 BETUBN END r  68  C C C C C  SUBROUTINE IIPHAS<A,Z,XOLD,K1,K2,K3,1,17ftR,ISP,ITI,ND, 1TEMP3, RH03,A3,X3,VOID3) IMPLICIT REAL*8 (A-H,0-Z) ..THIS SOLUTION IS THE APEX Of Ai: TBIINGLE, PT 3, WHOSE BASE POINTS 1 AND 2 ARE KNOWNf WE KNOl Z,TyU,P,H FOR THEM.) THE AFTER THE FIRST TIME CALCULATION IN AN ITERATION ALL CALCULATIONS CONSISTED OF FIXED VALUES FOB 1, AND 2 ABE BYPASSED. DIMENSION A (ND,ND) >Z {IVAE,I S P 1 TI) , X OL D (ND) PR=T.D5 ...FOR ANY ITERATION ARRAY ELEMENTS A ARE RESET TO ZESO.. DO Ii .11=1, 10 DO H JJ=1,11 4 A (II,JJ)=0.D+0 .... AFTER THE FIRST ITERATION MANY STEPS ARE BYPASSED. IF (I. HE. 1) GO TO 5 KK2=K2M .....FOR THE 1ST REFLECTIVE WAVE. .... IF(K1.EQ.2) KK2=K2 Z1=Z (1 ,K2-1,K1) T1=Z(2,K2- 1,K1) U1=Z (3,K2-1,K1) P1=Z(4,K2-1,K1) H1 = Z (5,K2-1,K1) TEMP1=Z (6,K2-1,K1) RH01=Z (7,K2-T,K1) A1=Z (8,K2-1,K1) X1=Z(9,K2-1,K 1) VD1=Z (10,K2- 1,K1) Z2-Z (1,KK2,K1-1) T2=ZC2,KK2,K1-1) U2=Z (3,KK2,K1-1) P2=Z(4,KK2,K1-1) H2=Z (5,KK2,K1-1) DELZ21=Z2-Z1 DELT21=T2-T1 DEL021=U2-U1 DELE21=P2-P1 DELH21=H2-H1 CALL SS (P1,TEMP1,X1,VD1,C,RPH,RHP,V) ; A1=C E1HP=RHP R1PH=RPH RHO1=1.D0/V ALF1=ALFA(P1,TEMP1 , X1 , 01) CALL "SS{P2,TEMP2,X2,VD2,C,RPH,BHP,V) A2=C R2PH=RPH R2HP=RHP RH02=1.DO/V ALF2=aLFACP2,TEMP2,X2,U2) FK1=FK (Pl,TEMPT,X1,U1,At,R1HP,aLF1,RH01) FL2=FLfP2,TEMP2,X2,U2,A2iB2HP,ALF2,SH02) Z3=Z2 1  C  C C  X3=X2 V0ID3=VD2 Z0= (Z1 + Z2) /2.D+0 T0= (T1*T2)/2.D0 00=01 P0=P1 H0=H1 TEHPG=TEHP1 X0=X1 U3=-(-P3>P1) / (A1*RHO 1) P3=P2 H3=H0+ (P3-P0)/RH01 TEMP3=TSAf (P3) . . s . ;aiE I E IN THE COMPRESSED LIQUID STATE? IFfX3.NE.0.D0)GO TO 7 TEHP3=T{H3) GO TO 8 7 X3=QUTY(P3,TEMP3,H3) 8 CALL SS{P3*TEMP3,X3,VD3,C,RPH,RHP,V) A3=C OA31= ((B3+A3) + (U1+A1J)/2.D+0 T3= (Z3-Z1) /UA31+T1 .....INITIALIZING XOLD ARRAY... XOLD (1)=Z0 XOLD (2) =T0 XOLD(3)=U0 XOLD(4)=P0 XOLD (5)=H0 XOLDf6)=Z3 XOLD (7) =T3 XOLD (8) =03 XOLD<9) = P3 XOLD (10) =H3 5 CONTINUE ....THE ARRAY ELEMENTS ABE COMPUTED^.!. A (1, 1)=DELP21 A (1,4) =-DELZ21 A (2, 1)=DELU21 A(2,3)=-DELZ21 A (3 , 1) = DEL H2 1 A (3,5) =-DELZ21 A (4,2) =DELP21 A(4,4)=-DELT21 A (5,1) = -1. D + 0 IF (I-EQ. 1) GO TO 6 Z3=X0LD(6) T3=XOLD (7) U3=XOLD (8) P3=XOLD (9) H3=XOLD (10) TE8P3=TSAT (P3) IF (X3..NE.0. DO) GO TO 17 TEMP3=T(H3) GO TO 18 :  17 X3=QUTY(P3,TEMP3,H3) 18 Z0=XOLD(1) TEMPO=TSAT (PO) IF(XO. NE.Q-DO) GO TO 27 TEMPO=T{HO) GO TO 6 27 XQ=QUTY (PO,TEMPO, HO) 6 ALFO=ALFA(PO,TEMPO,X0,U0) CALL SS(PO,TEMP0,XO,VDO,C,RPH,1HP,V) A0=C R0PH=RPH R0HP=RHP RHO0=1.DO/7 FMO=FM (PO,TEMPO,XO,UO,AO,ROHP,ALFO,RHOO) ....-F*S COMPRISE THE AUGMENTED COLUMN FOR THE MATRIX. F1=- ( (ZO-Z 1) *DELP2 1- (DELZ21) * (P0-P1) ) F2=- { (ZO-Z 1) *DELU2 1- (DELZ 21) * (UO-U1) ) F3=-{(Z0-Z1) *DELH2 1 - (DELZ21) * (H0-H1) ) F4-- ( (TO-T1) *DELP21- (DELT21) * (PO-P 1) ) A(1,11)=F1 A(2,11)=F2 A(3,11)=F3 A(<*,11)=F«I ALF3=ALFA(P3,TEMP3^X3,U3) CALL SS (P3,TEMP3,X3,VD3,C,RPH,RHP,V) ft3=C R3PH=BPH R3HP=RHP RHO3=1iD0/V VOID3=VD3 FK3=FK (P3,TEMP3,X3,U3,A3, R3HP, ALF3, RH03) FL3=FL{P3,TEMP3,X3,U3,A3,R3HP,ALF3,RH03) F«3=FM(P3,TEMP3,X3,U3,A3,R3HP,ALF3,RH03) U30= (U3+UO)/2.D+0 IF(I.EQ. 1) U30=U0 RHO30=(RHG3 + RHO0) /2.D+0 IF {I-EQ.1) RHO30=RHO0 ..;.AT FIRST M30 IS SET EQUAL MO, ETC. FM30= (FM3+FM0) /2. D*0 IF (I-EQ. 1) FM30=FM0 UA31= \ (U3+A3) + (U1+A1) J/2.D+0 I F (I. EQ. 1) U A31=0 1+A 1 ABH031=(A3*RH03+A1*RH01)/2.D+0 IF (I-EQ. 1) ARH03 1=A 1*RH0 1 FK3 1= (FK3+FK 1) / 2. D + 0 IF(I.EQ. 1) FK31=FK1 UB32= {(.03-A3)-* (U2-A2) ) /2. D+0 IFfI.EQ.1) UB32=U2-A2 ARH032= (A3*RH03*A2*SHO2)/2. D+0 IF (I.EQ. 1) ARH032=A2*RH02 FL32= (FL3+FL2) /2.D+0 IF (IiEQ. 1)FL32=FL2 BETA63=(3. D + 0/2.D+0)*ALF0*RHO0*U0**2*(T3-T0)/2-D+O BETA68=(3.D+0/2.D+0)*ALF3*RH03*03**2*(T3-T©)/2.D+0 BETA88 = ARH031+( (3. B+0/2.B+0} *A3**2*ALF3*U3**2*R3HP  1+SHG3*a3*AlF3*03j * (T3-T 1fc/2. D+0 BET 108 =— A8H032 * ( (3.D+9/2.D+0)*A3**2*ALF3*03**2*B3HP 1-BH03*A3*ALF3*03) * (T3-T2)/2.D+0 F5=-{Z3-Z0-tJ30*|T3-T0) ) F6=- (P3-P0-BHO30* (H3-H0) -M30* (T3-T0)) F7=- (Z3-Z1-UA31* (T3-T1) ) F8=- (P3-P1+ARH031* (03-01)-FK3 1* (T3-T1) ) F9=-(Z3-Z2-0B32* (T3-T2) } F10=-(P3-P2-A8H032*(03-02)-FL32*(T3-T2) ) A (5, 2) =030 A (5>3)=- (T3-T0) /2.D+0 A (5,6) = 1.D + 0 A (5 ,7) =-030 A (5, 8) =-(T3-T0)/2-D+0 A(5,11)=F5 A(6,2)=FH30 A (6,3)=BETA63 A ( 6 , 4) =-1.D+0 A (6,5}=BHO30 A(6,7)=-FH3Q A (6,8) =BETA68 A (6, 9) =1.D + 0 A (6,10)=-BHO30 A(6,11)=F6 A (7,6) = 1. D+0 A (7,7)=-OA31 A (7,8J =- (T3-T1) /2. D+0 A (7, 11) =F7 A (8,7)=-FK31 A (8, 8) -BETA88 A (8,9) =1.D + 0 A (8,11) =F8 A ( 9 , 6 ) = 1.0+0 A (9, 7) =-OB32 A (9,8)=- (T3-T2) /2- D+0 A (9,11) =F9 A (10,7)=-FL32 A(10,8)=BET108 A (10, 9)= 1. D+0 A (10, 11) =F10 BETOBN END  72  C C C  SUBROUTINE 0P EN ( A , Z ,XOLD,K1,K 2,K 3,1,1V Al,ISP,111,ND 1,TEMP3,RH03,A3,X3,VOID3) IMPLICIT REAL*8(A-H,0-Z) . i . I N THIS ROUTINE PT 3 IS UNKNOWN AND THE PTS 1, 2 KNOHN. EFFOBTS HAVE BEEN HADE TO REDUCE REDUNDANCE. DIMENSION A{ND,ND) ,Z (IVAR,ISP,ITI) , XOLD (ND) X=4.D*0 ...THE RESERVOIR PRESSURE IS TAKEN ATMOSPHERIC... PR=1.D5 DO 4 11=1,10 DO U JJ=1,11 a A ( I I , JJ) =0.D+0 IF (I.NE. 1) GO TO 5 Z1=Z(1,K2 -1,K1) T1=Z(2,K2-1,K1) U1=Z(3,K2-1,K1) ;  P1=Z (4,K2-?1,K1)  H1=Z(5,K2 -1,K1) TEMP1=Z (6, K2-1,K1) RH01=Z (7,K2-1,K1) A2=Z (8,K2-1,K1) X1=Z(9,K2-1,K1) VD1=Z (10,K2-1,K1) Z2=Z (T,K2,K1-1) T2=Z {2,K2;K1-1) U2=Z (3,K2,K1-1) P2=Z (4,K2,K1-T) H2=Z (5,K2, K1-1) DELZ21=?Z2-Z1 DELT21=T2-T1 DELU2t=U2-U1 DELP21=P2-P1 DELH21=H2-H1 ;  CALL SS(P1,TEMP1,X1,VD1,C,BPH,RHP,V) A1=C R1PH=RPH R1HP=RHP RH01=1-D0/V ALF1-ALFA (P1,TEMP1,X 1,U1) CALL SS(P2,TEHP2,X2,VB2,C,RPH,RHP,V) A2=C R2PH=RPH R2HP=RHP RH02=1. DO/V ALF2=ALFA(P2,TEMP2,X2,U2) FK1=FK(P1,TEMP 1,X1,U1,A1,R1HP,ALF1,RH01) FL2—FL(P2,TEMP2,X2,U2,A2,S2HP,ALF2, RH02) Z3=Z2 T3=T2+.002D+0 P3=P2 H3=H2 TEMP3=TEMP2 X3=X1 VOID3=VD1  73  C  C  TEHP0=TEMP1 X0=X1 U3=-(P3-P1)/{A1*BH01) P3=Z (5,K2,K1-T) H3=HG* {P3-P0)/RH01 TEMP3=TSAT(P3) .-.IS PT 3 IN THE CGHPBESSED LIQOID STATE? IF SG TEH=T{HF). IF (X3.NE.0. DO) GO TO 7 TEMP3=TCH3) GO TO 8 7 X3=QUTY{P3,TEHP3,H3) 8 CALL SS {P3,TEHP3,X3, v"B3,C,BPH, BHP,V) A3=C VGID3=7D3 0A31={(03+A3J + (01 + A1))/2.D+0 T3= (Z3-Z1)/UA31+T1 ..^.INITIALIZING XOLD ASSAY. XOLD{1) =Z0 XOLD (2) =T0 XOLD (3) =00 XOLD(4)=P0 XOLD(5)=H0 XOLD C6) =Z3 XOLD {7)=T3 XOLD (8) =03 XOLD (9)=P3 XOLD < 10) =H3 5 CONTINUE. A{1, 1) -1.D+0 A (1,4)=-DELZ21/DELP21 A (2, 1) =1.D + 0 A (2,3)=-DELZ21/DEL021 A(3,1)=1.B+0  A 0,5)=-DELZ21/DELH21  A f4, 2) =1.D + 0 A (4,4)•=-DELT21/DELP21 A (5, 1)=-1.D+0 IF (I.EQ. 1) GO TO 6 Z3=XOLD(6) T3=XOLD(7) U3=XOLD(8) • P3=XGLB &ff* H3=XOLB(10) TEHP3=TSAT (P3) IF {X3i HE. 0. BO) GO TO 17 TEHP3=T{H3) GO TO 18 17 X3=QUTY(P3,TEHP3,H3) 18 Z0=XOLD(1) T0=XOLD{2) 00=XOLB 13) P0=XOLO (4) H0=XOLD{5) TEHPO=TSAT (PO)  IF{XG.JE.O.BO) GO TO 27 TEMF0=T{H0) GO TO 6 27 X0=QOT¥ (PG,.TEMPO,HO) 6 ALFO=ALFA(P0,TEMPO,XO/UO) CALL SS (PO,TEMPO,XO,VD0 C,RPH,RHP,V) A0=C R0PH=BPH B0HP=RHP RHO0=1.DO/V FM0=FM(PO,TEMPO,XO,00,AO,ROHP,ALFO,RHOO) F1=Z0~Z1- (DELZ21/DELP21) * (P0-P1) F2=2G~Z 1-(DELZ21/DEL021) * (00-01) F3=Z0^Z1-(DELZ21/DELH21) *(H0-H1) F4=TG-T1-(BELT21/DELP21)*(PO-P1) A (1,11) =-Fl A (2, 11) =-F2 A(3,11)=-F3 A (ft, 11)=-F4 X3=QOTYfP3,TEMP3,H3) ALF3=ALFA(P3,T EM P3,X3,U 3) CALL SS (P3,TEMP3,X3,VD3,C,RPH,RHP,V) A3=C R3PH=EPH E3HP=EHP RH03=1.DO/V VOI03=VD3 FK3=FK(P3,TEMP3,X3,03,A3,E3HP,ALF3,EH03) FL3=FL (P3,TEMP3,X3,G3,A3,E3HP,ALF3,RHG3) FM3=FM (P3,TEMP3,X3,03,A3,13HP,ALF3,RH03) 030=(03+00)/2.D+0 EH030=(RHO3+RHO0)/2.D+0 FH30= (FM3+FM0)/2.D + 0 0A31= {{03+A3) + (01 +A 1) )/2.B+0 ARH031=(A3 *RHO3+A1*.RHO1)/2.D+0 FK31= (FK3+FK1)/2.D+0 IF (I.EQ. 1) FK31=FK1 0B32=( (03-A3) + (02-A2) )/2.D+0 IF (I.EQ- 1) 0B32=U2-A2 AEH03 2= (A3*RH03 + A2*RHQ2)/2. D + 0 IF (I- EQ. 1) A1H032=A 2*RH02 FL32= (FL3+FL2) /2- D + G IF (I. EQ. 1)iFL32=FL2 BETA63-(3. D+0/2.D+0)*ALF0*BHO0*00** 2*(T3-T0)/2.D+0 BETA68=(3.D+0/2.D+0)*ALF3*1H03*U3**2*(T3-T0) /2.D+0 BETA88=ARH031+({3.D+0/2.D+0)*A3**2*ALF3*03**2*R3HP 1 +SH03*A3*ALF3*03) * (T3-T 1) /2. D+0 BET108 =-ARH032+((3.D+0/2.D+0)*A3**2*ALF3*03**2*R3HP 1-RH03*A3*ALF3*03) * (T3-T2)/2.D+0 F5=Z3-Z0-O30* (T3-T0) F6=P3-PO-BH030* (H3r H0)r-FM30* (T3-T0) F7=Z3-Z1-UA31*(T3TT1) F8=P3-P1*ARH031* (03-01)-FK31* CT3-T1) F9=Z3-X r  75  ...AT THIS POINT CHECK IS HADE ABOUT CHOKING. IF 0<A, P=PB. F10=03-A3 IF (DABS (F10) .GT. 1. DO) GO TO 1 A (10,8) = 1. DO GO TO 2 F10=P3-PB ft (10,9) =1. DO A (5,2 = U30 A (5,3 = (T0-T3) /2.D+0 A (5,6 -1.D+0 A (5,7 =-U30 A{5,8 = (T0-T3) /2,D+0 A(5,1 )=-F5 A (6,2 =FH30 A (6, 3 =BETA63 A (6,4 =-1. D+0 A (6, 5 =BHO30 A (6,7 =-FH30 A (6, 8 =BETA68 A (6,9 = 1. D+0 A(6,10)=-RHO30 A (6,1 )=-F6 A (7, 6 = 1. D + 0 A (7,7 =-OA31 A (7, 8 = (T1-T3) /2.D + 0 A (7,1 )=-F7 A(8,7 =-FK31 A (8,8 =BETA88 A (8,9 =1.D+0 A (8,1 )=-F8 A (9,6 = 1-D+0 A (9,1 )=-F9 A (-10, 1)=-F10 RETURN END  76  C C C C C C  C  SUBROUTINE CLOSED (A,Z,XOLD,K1,K2,K3,1,IVAR,1SP,ITI,ND 1, T E M P 3, R HO 3 , A 3 , X3 , 70 ID 3) IMPLICIT REAL*8(A~H,0-Z) ...PTS 0,2,3 ARE USED IN THIS ROUTINE- PROPERTIES OF PT 3 ARE UNKNOWN. WHEN FIRST TIME USED IN AN INTE1ATI0N, THE PROPERTY VALUES OF PT 0 ARE ASSIGNED TO POINT 3. AFTER THAT THE ITERATED VALUES OF XOLD RESPECTIVELY ARE USE UNTIL I T CONVERGES, PTS 0 AND 1 ARE ONLY ASSIGNED AT FIRST ITERATIONDIMENSION A (ND , ND) ,Z (IVAR ,ISP,ITI) , XOLD (ND) DO 4 11=1,5 DO 4 JJ=1,6 4 A(II,JJ)=0-D*0 I F f l . N E . 1) GO TO 5 Z0=Z (1, 1,K1-1) T0=Z (2, 1,K1-1) U0=Z (3, 1,K1-1) P0=Z (4, 1,K1-1) H0=Z (5,1,K1-1) IF (K1.NE.2) GO TO 3 Z2=Z0 IF (K1.EQ.2) Z2=Z2-4.D0 Z0=0.DO T2=T0 U2=00 P2=P0 H2=H0 TEM P2=TEMP0 X2=X0 GO TO 2 3 Z2=Z (1,2,K1-1) T2=Z (2,2,K1-1) U2=Z (3,2,K1-1) P2=Zf 4, 2,K1-1) H2=Z (5, 2,K 1-1) TEMP2=Z (6, 2,K1-1) RH02=Z (7,2, K1-1) A2=Z(8,2,K1-1) X2=Z (9,2,K1-1) VD2=Z(10,2,K1-1) 2 CALL SS (PO,TEMPO,X0,VDO,C,RPH,RHP,V) AO=C ROPH=SPH ROHP=RHP RHO0=1. DO/V ALFO=ALFA(PO,TEMPO,X0,UO) CALL SS(P2,TEMP2,X2,VD2,C,RPH,RHP,V) A2=C R2PH=RPH R2HP=RHP RH02=1.DO/V ALF2=ALFA(P2,TEMP2,X2,02) F80=FM(PO,TEMPO,XO,U0,AO,ROHP,ALF0,RHO0) FL2=FL (P2,TEMP 2,X2,U2,A2,R2HP,ALF2,RH02)  TEMP3-TEMP0 X3=X0 ?OID3=?B0 XOL1>(T)=Z3 XOLD <2)=T3 XOLD (3) =03 XOLD (4) =P3 XOLD (-5) =H3 5 COBTINOE IF (i.EQ- 1) t50 .TO 6 Z3=X0LD(1) T3=XGLD(2) 03= XOLD (3) P3=X0LD{4) H3=XOLD{5) TEHP3=TSST (P3) IF{K1»EQ.2)SQ TO 6 IF { X 3 * H E » 0 w D 0) GO TO 7 TEMP3=f(H3) GO TO 6 7 X3=Q0TT{P3,TEMP3,H3) 6 &LF3=ALFA(P3,TEMP3,X3,U3) CALL SS {P3,TEMP3,X3,VB3,C,:IPH,BHP,?) A3=C R3PH=SPH R3HP=BHP BHO3=1.D0/7 VGIB3=VD3 FH3=F»(P3,T1MP3,X3>03,A3,R3HP,ALF3,BH03) FL3=FL (P3,T1MP3*X3>33^ A3,R3HPSiF3,1H03) 0B32= tf03-A3) • C02-A2) j/2.B+0 ABH032= CA3*RH03*A2*RHO2)/2.D+O FL32= (FL3+FL2) /2.D+0 BHO30= {BHO3+1HO0) /2. D+0 FM 30= (FH3*Ffl0} /2. D + 0 BET&23=-ARH032+- { (3. D+0/2.D+Q) *A3**2*ALF3*03**2*13HP1 RHG3*a3*ALF3*B3) * (T3-T2)/2.D0 BETA33= {3. D+0/2. D+0) *BH03*ALF3*03**2* (T3-T0) /2.D0 Ft=-OB32*{T3-T2) +Z3-Z2 F2=P3-P2-ARH032*(03^02)-FL32*(T3-T2) F3= P3-P0*1HO30* (H3- BO) -FM30* (T3-T0) F4=Z3 • F5=03 AC1,1)=1^D0 A(T;2)=-UB32 A f1,3| = (T2-T3) /2.D+0 A {1,6)=-Ft A(2,2)=-FL32 A(2,3)=BETA23 A (2, 4) =1.D + 0 A{2,6)=-F2 A(3,2)=-FM30 A (3,3) =BETA33 A (3, 4) =1.D+0  &(3,5)=-BHO30 A (3,6)=-F3 A (a,1)-1.D+0  A(4,6) =-F4 A (5,3) =1.D + 0 A (5,6)=-F5 BETTJBN BHD  79  FUNCTION SIHUL (N,A,X,EPS, INDIC, NIC) .....THIS ROUTINE COMPUTES THE SOLUTION OF THE N X N+1 MATRIX A. AUGMENTED VALUES ARE IN N+T COLUMN.  C C C  C  C  5  C  C  8 9  11 C  C C  13  C C C C  14  IMPLICIT RFAL*8 (A~H,0-Z) DIMENSION 1101 (50) , JCOI (50) , JGRB (50) , Y (50) , A (NBC, NEC) ,X (N MAX=N IF (INDIC. GE.O) MAX=N+1 ... IS B LARGER THAN 50. IF (N.LE.50)GO TO 5 WRITE (6,200) SIMUL =Q.D+0 RETURN ......BEGIN ELIMINATION PROCEDURE....... DETER =1.D+0 DO 18 K = 1,N KM1 = K - 1 .......SEARCH FOR THE PIVOT ELEMENT PIVOT= 0.D+0 DO 11 1=1, N DO 11 J=1,N ..SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCR IF ( K . EQ. 1) GO TO 9 DO 8 ISCAN =1,KM1 DO 8 JSCAN =1,KM1 IF (I.EQ. I ROW (I SCAN)) GO TO 11 IF (J. EQ. JCOL (JSCAN) ) GO TO 11 CONTINUE IF (DABS (A (I, J) ) . LE. DABS (PIVOT) ) GO TO 11 PIVOT= A (I ,3) IEOW(K)=I JCQL(K) =J CONTINUE .......INSURE THAT SELECTED PIVOT IS LARGER THAN EPS . IF(DABS(PIVOT).GT.EPS) GO TO 13 SIMUL=0.D+0 RETURN .UPDATE THE DETERMINANT IRG¥K=IRO¥ (K) JCOLK=JCOL(K) DETER=DETER*PIVOT  VALUE  ......NORMALIZE PIVOT ROW ELEMENTS...... DO 14 J=1,MAX A(IROWK,J) =A (IR01K, J) /PIVOT .....V.CARRY OUT ELIMINATION AND DEVELOP INVERSE...... A (IROWK, JCOLK) = 1./PIVOT DO 18 1=1,N AIJCK=A(I,JCOLK) IF(I.EQ.IROWK) GO TO 18 A (I,JCOLK)=-AIJCK/PIVOT DO 17 J=1,MAX  80  C C  17 18  IP (J. HE. JCOLK) A (I,J)=A (I> J)-AIJCK*A (IHOSK, J) CONTINUE  20  .......ORBER SOLUTION VALUES (IF ANY)AND CREATE JOED AREA DO 20 1=1,N IROfI=IROI (I) JCOL 1= JCOL (I) JORB(I10HI)=JCQLI IF(INBICiGE.O) X{JCOLI)=A(IROSI,HAX)  C C  22 C C  C C C  i  26  27 C  28  29 C C C C  ADJUST SIGN OF DETERMINANT. ... . . INTCH=0 NH1=N-1 DO 22 I=1,NM1 IPt=I*1 DO 22 J=IP1 ,N I F (JORB (J) .GE. JOED (I) ) GO TO 22 JTEMP=JORD (J) JOES (J) = JOBD(I) JOED f l ) =JTEMP INTCH=INTCH+1 CONTINUE IF (INTCH/2*2-NE. INTCH) DETEB=-DETER • • .•.'.•v-.'»- -JIF -INDIC IS POSITIVE aET01S',SlT& RESULTS.:.,.^*, , IF(INBICLE-O) GO TO 26 SIMUL= DETER RETURN  30  200  I F INDIC IS NEGATIVE OS ZERO, UNSCRAMBLE THE INVE FIRST BY B O f S i . . . . DO 28 J=1,N DO 27 1=1,N IROWI=I10W(I) JCOLI=JCOL(I) Y (JGOLI) =A (IROWI, J) DO 28: 1=1, N A(I,J)=Y(I) ......... ..THEN BY COLUMNS.-- > . - . . DO 30 1=1,N DO 29 J=1,N I1OWJ=IR0U{J) JCOLJ=JCOl(J) Y (IS01J) =A (ly JCOLJ) DO 30 J=1,N A(I,J)=Y(J) > i i . .:.EETUR8 > FOR INDIC 'NEGATIVE OE 2ERO. . . - . i . SIMUL=BETER RETURN . . v . . wviiv. ...FORMAT* FOR- OUTPUT STATEMENT.:. . , FORMAT(10HON TOO BIG) END  81  C  C C C C  FUNCTION RE (P,T, X,U) -...COMPOTES THE REYNOLD*S NUMBER. IMPLICIT REAL*8{&-ff,0-Z) I I I L * 8 MU,MUF,MOS&T 6=9^8094D+O D=.032D+0 -.l.E=0.D'+0 • IF.C>O.Lfe.02D+O-):RETIIRNM0=M0F (P,T) IF {X.LE.0.1D+0 )GO TO 1 fiflO='T. D+0 / V (P,T X)MU=HOSAT(P,T) GO TO 2 1 RHG=1.D + 0 /VFfT) RE=RHO*0*D/N0 RETURN END :  ;  #  C C  82  C  C C  FUNCTION TASTAR( P,T,X,U) ..-.TWO PHASE MULTIPLIER (HANCOX-NICOL) . 0.<TASTAR<1. IMPLICIT REAL*8(A-H,0-Z) ; REAL*8 MUF,MUFF,MUG,HOGG G=9.8094D+0 TASTAR=1.D+0 IF (X.LE.0. 1D+0 . OR-X.GE.0.9D + 0 )RETURN VMAS =VELHAS (P,T,X,U) VMAS=DABS(VHAS) VGG=VG(P,T) VFF=VF(T) MUFF=MUF(P,T) MUGG=MUG (P,T) GAMMA=VM AS /DABS ( (G*MOFF* f t . D+0 /VFF-1.D+0 /VGG) )/VFF) **. 33 D=MUGG/MUFF BETA= (VFF/VGG) * (DABS (D) ) **.2D+0 TASTAR=1-D + 0 +X* (BETA-1. D+0 ) * (1.D+0 +3-57D + 0 *DEXP (-.00994 1*GAMMA) * (1. D+0 -DEXP(-4.96D + 0 *(1.D+0 -X)) )) TASTAR=DABS(TASTAR) RETURN END  83  FUNCTION FL (P»T,X, 0,»ft,RH,Al. BO)f ... THE CONSTANT IN DP-RADO=LDT (ALONG THE 0-A DIRECTION.) IMPLICIT BEAL*8(A-H,0-2) D=.032D+0 DELY=0.D+0 DELZ=4.D+0 /10.D+0 6=1*00+0 AG=9.8D0 FL=0.D+0 IFCO-LE.0. 02D+0 ) BETOBN QU=0.D+0 BHO=BO QO=QO*G SHP=BH/G**2 ALF=AL ALF1=ALF/G A=AA DYDZ-DELY/DELZ DYDZ=-T.D+0 ALF=ALF/2.D+0 ALF 1= ALF 1/2. D+0 FL=-(A**2*(QU+ALF*DABS(0)**3)*RHP-RHO*A*(ALF1*0*DABS(D) 1+AG*DYDZ)) BETORN END #  C  C C  8H  C  C C  FUNCTION FH(P,T,X, U,AA,RH,AL,80) ...THE CONSTANT IN DP-R DH=M DT (ALONG THE 0 DIRECTION.) IMPLICIT REAL*8 (A-H,0-Z) DELY-O.D+0 DELZ=4.D+0 /10.D+0 G=1.OD+0 AG=9.8D0 FM=0.,D + 0 IF (U.LE.O.02D+0 )RETURN QU=O.D+0 RH0=RO ALF=AL/G**2 ALF=ALF/2.D+0 FM=—RHO* (QU+ALF*DABS (U) **3) RETURN END  85  FU NCTION F K (P , T ,X, U ,ftA ,.B H , AL , SO) -- - THE CONSTANT IN DP*HA DO=K DT (ALONG THE U+A DIBECTION.) IMPLICIT REAL*8 (ft-H,0-Z) DELY=0.D+0 DELZ=4.D+0 /10.D+0 G=1.0D+0 AG-9.8D0 FK=0.D+0 IF (O.LE.,Q2D+0 )SETOBN QD=0.D+0 BHO=10 CO=QD*G BHP=RH/G**2 SLF=AL ALF1=ALF/G A=AA -  DYDZ=DELY/DE1Z DYDZ=-1.D+0 ALF=ALF/2.D+0 ALF1=ALFl/2-D+O FK=- (A**2* (QO+ ALF* DABS (0) **3) *BHP+BHO* A* (ALF1 *0*DABS (0) 1+AG*DYDZ)) BETOBN END  C  90 2  10 2 100 101 103  10 4 105  106 C  FUNCTION T(E) ...FINDS THE TEMP,GIVEN THE ENTHALPY. IMPLICIT EEAL*8fA-H,0-S) ; COMMON ENTH ENTH=E IF (BiGT.Oi D+0) GO TO 2 WSITE(6,90) FOBMAT(2X, * NEGATIVE ENTHALPY*) STOP XLO=Q.D+0 XN=0.D+0 SOLN=0.D+0 EP5= 1.B-3 H=2-D+0 FL0=FUN (XLO) DO 102 1=1,1000 XB0=XLO*H*(I) FSO=F0N (XBO) IF(FLO *FRO)100,101,102 CONTINUE XN= fXLO + XBG)/2.D+0 FH=FUN (XN) IF(FLO*FN) 103,104,105 SOLN=XHO GO TO 106 X!G=XN FBO=FN IF(DABS(XBO—XLO); .GT. EPS) GO TO 100 SOLN=XSO GO TO 106 SOLN=XN GO TO 106 XLO=XN FLO=FN IF(DABS(XBO-XLO).GT.EPS)GO TO 100 SOLN=XLO T=SQLN BETUBN END ...FUNCTION APPBOXIHATING H LIQUID. FUNCTION FUN (T) IMPLICIT BEAL*8(A-H,0-Z) COMMON ENTH A=6.0653D+0 B=3.S872D+0 C=.21752D-02 D=-*103390-04 E=.27801D-07 FUN=A+B*T*C*T**2+D*T**3 +E*T**4 F0N=FUN*1. D+3-SNTH BETUBN END  FUNCTION TSAT(P) iV^SlTUHATION LINE.. IMPLICIT REAL*8 (A-H,0-Z) COMMON Y P=P*9.869D-6 Y=218.1670*0 /P XLO=0.OD+0 XN=0.D+0 XLO=XLO+273.16D*0 SOLN=0.D*0 EPS=1.D-3 H=2.0D+0 FEO=FUNC(XLO) DO 102 1=1,1000 XRO=XLO*H*fT) FRO=FUNC(XBO) IF(FLO*FRO)100,101,102 10 2 CONTINUE 10 0 XN= tXLO*XBO)/2. D+0 FN=FUNC(XN) IF (FLO*FN) 103, 104, 1 05 101 SOLN=XBO TSAT=SOLN-273.16D+0 GO TO 106 103 XBO=XN FRO=FN IF(01BS(XBO-XLO) iGT.EPS)GO TO 100 SOLN=XBO TS A T=SOL N-27 3.16 D+0 GO TO 106 10 4 SOLN=XN TSAT=SOLN-273.16D+0 GO TO 106 10 5 XLO=XN FLO=FN IF (DABS(XRO-XLO).GT.EPS)GO TO 100 SGLN=XLO TSAT=SOLN-273. 16D+0 106 P=P/9.869D-6 BETURN END C LOG FUNCTION OF SATURATION PRES. OF TEMPERATURE. FUNCTION FUNC(T) IMPLICIT SEAL*8 (A-H,0-Z) COMMON Y a=3.3463130D+0 B=4.T4113D-2 C=7.515484D-9 D=1.3794481D~2 E=6.56444D- 11 PC=218.167D+0 TC=647.27D+0 X=TC-T FUNC= (X/T) *((A + B*X+e*X**3 +E*X**4 )/( 1.D+0 +D*X) ) F4JRC=FUNC-DLOG 10 (Y) RETURN END C  FUNCTION PSfiT(T) IMPLICIT HE At* 8 (&-H O-Z) T=T + 273. 16D+Q A-3.346313OD+0 B=4.14113D-2 C=7.515484D-9 D=1.3794481D-2 E=6. 56444 D-11 PC=218.167D+0 TC=647.27D+0 X=TG-T r  Y= tX/T)*( ta*B*X+C*X**3 +E*X**4 )/(1.D+0 *D*X) )  Y=10-D+0 **Y PSaT=PC/I T=T-273.16D+0 PSaT=PSlT/9.869D-6 1ETUBN END  89  C C  C  SOB BOOTING SS (P,T,X,VD,A,RPH,RHP, V) ..4 SOUND SPEED=1/(BPH+RHP/R)**2. VOID FRACTION IS ALSO FOON IMPLICIT BEAL*8 (A-H,0-Z) .,.PAS., ABE FED. A BETOBNS IN M/S. DELP=1.D-2 DELT=1.D-3 PDELP=P+DELP TBELT=T+DELT G=9.8094D+0 VFTP= (VF (TDELT) -VF (T) ) /DELT HFTP= (HF (TDELT) -HP (T) )/DELT I F (X.LE.0.0B+0 )GO TO 1 VGPT=(VG(PDELP,T) - VG(P T))/DELP VGTP=(VG(P,TDELT) -VG(P,T))/DELT HGPT= (HG (PDELP, T) -HG (P,T) )/DELP HGTP= (HG (P,TDELT)-HG(P,T) ) /DELT VGF=VG (P T) VFT=VF(T) V=X*VGF+(1.D+0 -X)*VFT RHO=1.DO/V B1=X*VGPT B2=X*VGTP+(1.D+0 -X)*VFTP B3=X*HGTP+ (1-D + O -X) *HFTP B4=X*HGPT BPH=-(1.D+0 /V**2) * (B1-B2*B4/B3) RHP= - (1. D+G /V**2) * (B2/B3) VD=X/ (X+ ( 1-X) *VFT/VGF) GO TO 2 V=VF(T) VD=0iDO VFF=-(1-D+O /V**2) RPH=0.D+0 BHP=VFF*(VFTP/HFTP) Y=DABS (RPH+V*RHP) A=1.D*0 /DSQRT (Y) RETURN END r  f  C C C C C C  1  FUNCTION VF{T) ....SAT VALUE USED FOR LIQUID SPECIFIC VOLUME. IMPLICIT REAL*8 (ft-H,0-Z) A=t.0142D+0 B=-.46549D-03 C=.10871D-04 B=-.41385D-07 E=.82105D-10 VF=a+B*T+C*T**2+D*T**3+ E*T**4 VF=VF/1. D+3 RETURN END  FUNCTION HF{T) . . < . . : SAT VALUE USED FOB LIQUID SPECIFIC ENTHALPY.IMPLICIT BEAL*8(A-H,0-Z) A=6.0653D+0 B=3.9872D+0 C=.21752D-02 D=-.103399-04 E=.2780lD-07 HF=A+B*T+C*T**2+D*T**3 +E*T**4 HF=HF*1.B+3 BETUBN E.ND  92  C C  FUNCTION HG (P,T) ..... SUPERHEftT SPECIFIC ENTHALPY... IMPLICIT REAL*8 (A-H,0-Z) .. ..m FEED PAS. AND DEG C. HG RETURNS IN J/KG. Tl=T+273. 16 D+0 P1=P*9.869D-6 S=1/T1 A=8087O.D+O *S**2 A=10.D+0 **A A1=A*DLOG (10. D+0 ) A1 = A1 *S ** 3 * 1617 40. D +0 A2=2.D + 0 *S*A A3=A1+A2 A3=A3*(-264T.62D + 0 ) F0=1.89D+0 + A3 C1=-2641w62D+0 *S*A C=1.89D+0 +C1 A4=S*A A4=A4*A1 A4=A4*(-2641;.62D+Q ) D=A4/S E1 =82.546D+0 *S-1*6246D5*S**2 Gl=82.546D+0 -3.2492D5*S Fl=2.D+0 *S*C**2*ET+2.D+0 *S**2*C*E1*D+S**2*C**2*G1 E3=.21828D + 0 -1.2697D5*S**2 G3=-2. 5394D5*S F3=4.D+0 *S**3*C**4*E3+4.D+0 *S**4*C**3*D*E3*S**4*C**4*G3 E12=3.6 35 D-4-6,768D64*S**24 G12=-1.6243D66*S**23 F12=-(13* D+0 *S**12*C**T3*ET2+T3.D+0 *S**13.D+0 *C**12 *B*E12+S**13*C**13*G12) H=F0*Pl+F1*Pl**2/2.D+0 +F3*P1**4y4.D+0 +F12*P1**13/13.D+0 H=H*. 101325D0 CP=CPT (T)/1.D+3 H1= (H•CP * 2 50 2.36D•0)*1.0003 DO HG=H1*1.D+3 RETURN END  93  C C  FUNCTION VG (P,T) ..SUPERHEAT SPECIFIC VOLUME.-.. WE FEED PUS. aND DEG C. VG RETURNS IN M**3/KG. IMPLICIT REAL*8(a-H O-Z) P1=P*9.869D-6 a5=.2088543D+5 T1=T+273.16D+0 S=1/T1 a=8G870*S**2 f  a=io.D+0 * * a ai=a*DLOG(10.D+0 ) A1=A1*S**3*161740.D+0 a2=2.D+o * s * a A3=AT+A2 A3=A3* (-2641i62D+0 ) FO=1.89D+0 +A3 CT=-2641. 62D+0 *S*A C=1.89D+0 +C1 E1 = 82.546D+0 *S-1.6246D5*S* *2 E3=.21828D+0 - 1.26 97D 5 * S ** 2 E12=3.635D-4-6-768D64*S**24 B=C+C**2*E1*S*P1+C**4*E3*S**3*P1**3-C* VG=4.55504D+0 /(P1*S)+B VG=VG/1.D+3 RETURN END  FUNCTION 8UF(P;T) i . - .LIQUID VISCOSITY, . v . IHPLICIT RBAL*8{A-H,O-Z) ...PAS. AND DEG. C R1E FED. HUF EETUIN IN KG/M SEC. EEAL*8 HUF/HULSAT T1=T+273. 16D+0 A6=T1-140.D*0 A6=1.D+0 /A6 A6=247.80*0 *A6 A6=10.D+0 **A6 A6=24T.4D*Q *A6 HULSAT=A6*1.D-6 A7=T1-305.D+0 A7= 1.0467D+0 *A7 A7= A7* {P-PSAT fT) ) *1.D- 6* 1.01 325D+0 * 9 . 869D-6 A7=1.D+0 +A7 A7=A7*A6 HUF=A7*1. D-6 RUF=HUF/10.D+0 RETURN END :  95  C C  FUNCTION MUG (P,T) VAPOR VISCOSITY. IMPLICIT REAL*8 f A~H,0-Z) ...WE FEED PAS. AND DEG. C. MUG RETURNS IN KG/M SEC. REAL*8 MUG,MU1 MU1=.4G7D+0 *T+80.4D+0 VGG=VG (P,T) *1. D+3 HUG=MUl-{1.D+0 /VGG)* {1858. D+0 -5.9D+0 *T) MUG=MUG*1.D-6 MUG=MUG/10.D+0 RETURN END  96  C C  FUNCTION HUSAT {P,T) .v;.iSATURATION VISCOSITY... IMPLICIT HEAL*8(A-H,0-Z) . . . m S .  AND  DEG.  C  REAL*8 HUF,HUSAT T 1=T+ 273. 16D+ 0 A6=T1-140.D+0 A6=T.D+0 /A6 A6=247.8D+0 *A6 A6=10.D+0 **A6 A6=241.4D+0 *A6 SUSAT=A6*1.D-6 HUSAT=HUSAT/10.D+0 RETURN END  ARE  FED.  HO  SAT  1  SET  DEN  IN  KG/H  SEC.  97  FUNCTION' QUTY(P,T,H) ....EVALUATES THE QUALITY IMPLICIT REAL*8(A-H,0-Z) C WE FEEB PAS., DEC C, AND J/KG. HFF=HF(T) HGG=HG (P,T) Q0TY= (H-HFF) / (HGG- HFF) C WRITE(6,100)QUTY C 100 FORMAT (1X, *QUTY=' , G14. 7) QUTY=DABS (QUTY) RETURN END C  98  C  FUNCTION HC(P,T,X,U) IH.PLICIT REAL*8 (A-H,0-Z) .PAS-, DIG C, AND S/SE A l l FED IN. HC RETURN IN J/M**2 S C HEAI*8 HUF,BUFF D=.032D+0 G=9.8094D+a flUFF=MUP (P> T) CRK=CK <P,T) VHAS = VELHAS (P T,X,U) CP=CPT(T) HC=.023D+0 * (CKK/D) * (VM AS *B/H0FF) **- 8D + 0 HC=HC* (CP*HUFF/CKK) **. 33333D+0 RETURN END #  99  C C  FUNCTION CK(P,T) ....EVALUATES THE COEFFICIENT OF CONDUCTIVITY. IMPLICIT REAL*8 (A~H,0-2;) .^WE FEED PAS. AND DEC. C.THE FUNC GIVES J/SEC. M. DEG C. TQ=273iTSD+0 T1=T+T0 X=T1/T0 A0=-922.47D+O AT=2839*5OD+0 A2=-180O.7D+0 A3=525.77D+0 Att=-73.attD+0 B0=-.9473D+0 B 1=2. 5186D+0 B2=-2.0012D+0 B3=v51536D+0 C0=i.6563D-3 Cl=-3. 8929D-3 C2=2.9323D-3 C3=-7.1693D-4 CK=A0 + A1*X+A2*X**2 CR=CK+A3*X**3+AU*X**4 CK1=B0+B1*X+B2*X**2+B3*X**3 P1=P-PSAT(T1-T0) P1=P1*9.869D-6 P1=P1*1.01325D+0 CK1=CK1*P1 CK2=C0+C1*X+C2*X**2+C3*X**3 CK2=P1**2*CK2 CK=CK+CK1+CK2 CK=CK*1.D-3 RETURN END  FUNCTION VELMAS (?,T,X,W) COMPUTES THE MASS VELOCITY.IMPLICIT HEAL*8(A-B,O-Z) VV=V (P,T,X) VELMAS =U/VV HETUHN END  FUNCTION ALFA CP >T, X ,U) • »•.•..<••. FB©«-'- F=XL*A*0**2 FOB FBICTIOH. . •IBPfclCIT• BSSL*8 (A-H,0-2) 'ALFA--- *D0 GO TO 3 ALFfi-0.0+0 IF (Us LE ...02D+0 )RETURN 1=RE (P,TyX,U) TSTAR=TASTSR (P,T,X>0) ALFA=0.0D+0 IF (B.I.E. 1. D+0 ) RETURN IF{R.LE.,200*0.0+O ) GO TO 2 F=.046D+0 /R**0.2B+0 GO TO 1 2 F=64.D+0 /B 1 CONTINUE G= 9.8094D+0 D=.032D+0 AKAL=0.D+0 ALFA= (H.D+0 *F*TSTAR/D+AKAL) /2.D+0 3 RETURN END  102  C C  FO NCTIG K Q { H C , 7 , T) IMPLICIT REAL*8 (A-H,0-Z) ... J/H**2 S CV M**3/KG, AN0 C ABE FED. FROM THE WALL RETURNS IN J/KG S. TW=253.D+0 D=.032D+0 Q=0. D+0 IF(HC.EQ.0.D+0 )BETUBN Q= ( D + 0 *HC*V/D) * (TW-T) RETURN END  Q, HEAT FLUX  103  C  FUNCTION V (P,T,X) EVALUATES THE SPECIFIC VOLUME.„*... IMPLICIT REAL*8(A-H 0-Z) VFF=VF (T) IF{X.LE.. 1D+0 )RETURN V=X*?G <P,T) + (1. D+0 -X)*VFF RETURN END #  C  FOHCTIOH --CPT (T) IMPLICIT REAL*8 {A-H,0-2) ....GlfEH T m BEG. C, CPT BETBBHS AS 3/KG. Tl=T+27 3. 16D+0 A=1.4720 B=7,5566B-4 C=47w8365B+0 B=273.16B+0 CPT=A* (Tl^-B) + (B/2. D+0 }*<T1**2 -D**2 ) 1 +C*tBL0GfT1)-BLOGCD}} CPT=CPT*1. B+3 RETUBI EKB  105  Figure  1  A P i p e a t Q u i e s c e n t S t a t e w i t h Some I n i t i a l Void Oriented V e r t i c a l l y f o r Upflow Blow-down  106  0  J  ,  0  2  ,  I  4  I  Time, ms Figured  , I  6  , I  , T  8  10  .  P r e s s u r e h i s t o r y a t c l o s e d end s h o r t l y a f t e r p i p e burst  107  Figure 3  C h a r a c t e r i s t i c s Diagram i n Z-T  Plane  108  R-WAVE A  z  A P  3f 31  31 31 31 31 31 31  ( u + a )  + ( a R )  A t  A u  = K  A t  L-WAVE P-WAVE  fc 32 " 32 32 32" 32 32 32 32 z  =(u  a)  A P  A Z  A P  30 30 30 30" 30 30 30 30 = U  A t  R  A h  = M  ( a R )  A u  = L  A t  A t  A Z  A z  h  At  01 21  A t  A t  01 21  A u  A u  01 21  A P  A P  Ah 01 21 " Ah j Q 1  2  t  z  Figure 4  G e n e r a l Two-Phase a l g o r i t h m t o advance s o l u t i o n dz/dt=u±a c h a r a c t e r i s t i c s  along  109  L-WAVE At =z /(a-u) 3 2  P-WAVF A P  A P  2  32"  ( a R )  3 2  32 2 32 32l u  = L  A t  30" 30 30 30 30 R  A h  = M  A t  C l o s e d end a)  C l o s e d End  R-WAVE A  Z  A P  3f 31  + ( a R )  P  31 31 31 31 31 31  ( u + 3 )  A t  A U  = K  or  A t  U  a 3  i f  =  a  u  3 3 < a  3  P-WAVE A z  A P  A 2  01  A t  A t  01 21  A P  01  Au Au  Q 1  2 1  = U  A t  _ R  A h  = M  Open end  oi Ah Ah  2 1  b) Figure 5  30 30 30 30 30 30 30 30  Open End Wave i n t e r a c t i o n s a t t h e b o u n d a r i e s  A t  Figure 6  E x p a n s i o n f a n a t t h e end o f a d u c t s u d d e n l y opened t o atmosphere  ~1  Projectile  'Breech h a)  E x p a n s i o n Fans C e n t e r e d a t p o i n t o  b)  E x p a n s i o n Fans f o r Z e r o Mass B u l l e t  Figure 8  E x p a n s i o n Fans f o r H e l i u m  Gun  Q_  * i-  *  2.-I  1  0  20  Figure 9  • 1 40  1 60 Time,  Pressure at Closed ^  Upflow  •  Downflow  1 80  4 * 100  T 120  ( P 0 = 3 . 5 MPa  )  T  mS End  1— 140  0.0  Figure  0.5 Void F r a c t i o n , a  10  1.0  C h o k i n g e f f e c t a t Open End; on Each S u c c e e d i n g E x p a n s i o n F a n , V o i d F r a c t i o n I n c r e a s e s and L o c a l V e l o c i t y Draws C l o s e r t o L o c a l Sound Speed  1  —I  20  40  1  1  1  1  60  80  100  120  Time, F i g u r e 11  Void U p f 1 ow Downf1ow  1—  140  mS  Fraction at Closed  End  (po=-12 MPa  )  1.0. c o  +J  •r—  o «o  it  0.5.  "O •I—  o  O.Oi  —I  1—  20  60  40  1  80 Time,  F i g u r e 12  Void  •  Upflow  •  Downflow  1  1  100  120  1—  160  140  mS  Fraction at Closed  E n d (P0=.18 MPa  )  en  • P0=.18 MPa  0  o  o  o  o  o  A ~1~  40  20 F i g u r e 13  ~1  60  A  80 Time,  mS  Pressure at the Closed  •  Upflow  O  •  Downflow  A  10 End  O A  O A  P0=.12 MPa A  C/1  CD  01  +-> ca  o  T i m e , mS Mass Flow Rate a t t h e Open F i g u r e 14 P0= 1 8 M P a P0=.12MPa '• Upflow A Downflow —  End  MAIN PROGRAM  C a l l EVAL ( f o r CLOSED)  3 C a l l EVAL ( f o r 11 PHASE)  C a l l EVAL ( f o r OPEN)  Figure  15  Computer Flow  Chart  SUBROUTINE EVAL  C a l l OPEN  Deter=SIMUL "„  J = 1  »  1 0  X N ( J ) > XOLD(J)?  1  Store returned r e s u l t s i n Z-Arraj  Return  F i g u r e 15  Computer Flow C h a r t  SUBROUTINE IIPHAS  S e t A. .=0  Computef , a , x, Al, etc. for p o i n t s 0 and 3 Assign values to p o i n t s 1 and 2  I n i t i a l i z e points 0 and 3  E v a l u a t e Az 21 etc.  I n i t i a l i z e the solution , X(j)  < •  I n i t i a l i z e K, L, M f o r 1=1, o t h e r w i s e assign values Compute augmented matrix comprized o f -F's  S e t up A.. i = l,10 a n d j = l,ll 1  Return  A s s i g n new s o l u t i o n s to p o i n t s 0 and 3 •  End  F i g u r e 15  Computer Flow C h a r t  122  Table 1 Initial  P a r a m e t e r s f o r t h e Gas Dynamic Problems  Conditions: u=0.0 m/s p=7.0 MPa T=243.0 C  Boundary C o n d i t i o n s : a t z=0.0 m: u=0.0 m/s a t z=4.0 m: p=p  i f u<a  e  o t h e r w i s e u=a ( c h o k i n g  condition)  Parameters: a=804.2 mm dA/dz=0.0 m dZ/dz=0.0 D =32.0 mm -1 k/l=0.0 m q=0.0 W/m 2  e  1  C o n s t i t u a t i v e Laws: F=(4f D  +k/l)u =2fu /D 2 2  2  e  e _ where f = f r i c t i o n f a c t o r = 1 0  6  T h e o . ____ • . N u m e r . Sheo. burner. Velocity P 0.0 0.0 ~ 07298?"540'E- 0 2 : ~ 6 T 2 . 9 9 l " 7 7 6 E - 0 2 0 . C 0. 1 3 3 6 9 9 9 E + 0 4 0 . T 4 3 3 83 9F + 0 0 0 . 1 4 4 0 1 5 0 E + 0 0 0 . 3 0 9 5 3 4 0 E - 0 2 0 . 3 0 97 6 7 5 E - 0 2 0 . 6 8 9 1 5 26=+ 0 2 0 . 1 3 1 3 94 2E +C4 0 . 2 9 9 5 2 1 3 5 + 0 0 C. 3 0 0 2 8 1 5E + 0C 0 . 3 2 0 5 E 4 8 E - 0 2 0 . 3 2 0 8 8 2 6 E - 0 2 0 . 1 3 7 8 2 9 6 E + 0 3 0.1291088E+04 •C.4695316F+00 0.4705089E+00 0.3322399E-C2 0.3326067E-02 0.2067453E+C3 C.1268228E+04 Z  Z  S o u n d  S  e e d  7j77J54^4TWnni70T7J^^  0 . 8 5 - 6 2 1 2 3 E + 00 0. 8 5 7 9 2 9 9 E + 0 0 G . 3 5 7 5 4 1 2 E - 0 2 0 . 3 5 8 0 6 7 1 E - Q 2 0 . 3 4 4 5 7 9 5 5 + 03 0 . 1 2 2 2 4 9 2 E + 0 4 O..l'n75738E + 01 0 . 1 0 7 8 0 0 3 E + 0 1 0 .37 130 2 7 E - 0 2 0. 3 119C74E-0 2 0. 4 1 2 4 9 1 IE + C2 0. 1 1 9 9 6 1 6 E + 04 0 . ~ T 3 T 4 T 7 % T C 1 O . T 3 T 7 8 ' 2 8 F + a i C .3 85 8 7 C 0 E - 0 2 0 . 3 8 6 5 6 5 4 5 - 0 2 0 . 4 8 2 4 1 4 3 E + 0 3 C l 1 7 6 7 2 5E + 0 4 0 . 1 5 7 5 4 7 3 E + 01 0 • 1 5 7 9 2 6 2 E + 0 1 0 . 4 0 1 3 1 3 6 E - C 2 C. 4 0 2 1 0 7 0 E - 0 2 0 . 5 5 1 3 3 0 3 E + 03 0 . 1 1 5 3 8 5 C E + 0 4 0 . 1 8 5 9 5 6 5 E + 01 0 . 1 8 6 4 3 6 9 5 + 0 1 0 . 4 1 7 7 0 5 3 E - 0 2 0 . 4 1 8 6 G 4 4 E - 0 2 0 . 6 2 C 2 4 4 6 E + G3 G. 1 1 3 0 9 6 1 E + 0 4 0 T 2 T 6 9 4 ? 5£ + G 1 ~ 0 T 2 " T T 5 T 4 5 P + 0 1 0 . 4 3 5 1 2 ^ e T ^ C 2 ^ 4 " 3 T O 7 3 T 0 7 ^ W r " 5 " c T E T 0 T ~ C V r i O W c ^ + l ? 4 _  r  _  0 7 7 3 6 8 T 3 4 E - 2 3 0 . 0 0 . 3 2 0 7 3 35E-G 2 0 . 3 21124 2 E - 0 2 - 0 . 3 3 8 8 1 2 2E- 2 G 0 . 1291 1 82E+04 0. 1 5 6 7 6 8 2 E + 0 0 0 . 1 5 6 9 5 3 2 E + 0 0 0 . 3 3 2 6 0 0 0 E - 0 2 0 . 3 3 3 0 6 6 9 E - 0 2 0 . 6 8 9 1 6 7 8 E + G 2 0 . 1 2 6 8 3 2 8 E + 04 0.3281616E+C0 0.3286345E+0C C.3451372E-02 G.3456880E-02 0,1378351E+G3 0.1245472E+04 0T5lT32~7Trr+Wl^TT^ 0.72C3540E+C0 0.7217692E+00 0.3724398E-02 0.3731794E-02 0.2756753E+03 0.1199747E+04 C o 9 4 4 3 O 0 5 E + G 0 C . 9 4 6 4 C l 9 E + G0 0 . 3 8 7 3 2 5 2 E - G 2 0. 3 8 8 1 7 C 5 E - 0 2 0 . 3 4 4 5 9 6 2 E + 0 3 0 . 1 1 7 6 8 7 7 E + 0 4 0 . 1 1 8 9 2 0 4 E + 0 1 0 . T 1 9 2 1 6 1 E + 0 1 0 . 4 0 3 1 2 2 ? E - Q 2 0 . 4 0 4 0 8 2 4 E - 0 2 0 .4 1 25 1 T C * + C 3 C . 1 T 5 4 C C 2E + 0 4 0. 1 4 5 7 1 0 6 E + C1 0. 14 6111 6E+ 0.1 0 . 4 1 9 9 0 7 5 E - G2 0.4 2 0 9 9 2 0 E - 0 2 0 . 4 8 2 4 3 6 8 5 + 0 3 0 . 1 1 3 1 1 2 6 E + 04 Oc1750304E+01 0.1755579E+01 0.4377655E-C2 0 . 4 3 £ 9 8 4 6 E - 0 2 0.5513550E+03 0.1108244E+04 0. 207134-OE + O l 0 . 2 0 7 f f l 3 0 E + 0 l 0 . 4 5 6 7 W 5 T a 2 ~ 0 T 4 T 8 l 3 4 3 E - 0 2 0 . 6 2 0 2 7 C 6 e + C 2 G. 1 0 8 5 3 5 9 E + G 4 r  0. 8 0 7 2 2 4 3 E - 2 3 C . C 0 7 3 45 29 6 I E - 0 2 0 . 3 4 5 9 0 5 7 E - 0 2 - C . 3 3 8 8 1 2 2 E - 2 0 0 . 1 2 4 5 5 7 6 E + 0 4 0 . 1 7 2 1 3 0 3 E + 0 0 0 . 1 7 2 4 5 9 2 E + 0 0 0 . 3 5 8 7 9 6 3 E - 0 2 0 . 3 5 9 5 0 0 1 E - 0 2 0. 6 8 9 2 1 18E +C2 0 . 1 2 2 2 7 1 6 E + G 4 0 . 3 6 1 1 1 9 3 E + 0 0 0 . 3 6 1 9 0 7 2 E + Q0 0 . 3 7 3 1 0 1 8 E - 0 2 0 .3 7 39 1 0 4 E - 0 2 0 . 1 3 7 8 4 3 7 E + C3 0 . 1 1 9 9 8 5 6 E + Q 4 0. 5 6 8 6 0 1 9 E + 0 0 0 . 5 7 0 0 C C 6 E + 00 0 . 3 8 8 Z T T Z F ^ C 2 C. 3 892~0~325-02 0. 2 0 6 7 6 71 E • 03 0 . 1 1 7 6 9 9 4 E + 0 4 0. 7 9 6 4 C 1 2 E + 00 0 . 7 9 8 5 8 4 6 E + 0 0 0 . 4 0 4 4 0 1 6 E - 0 2 0 . 4 C 5 4 5 1 8 E - 0 2 0 . 2 7 5 6 9 1 2 E + C 3 0 . 1 1 5 4 1 29E + 0 4 0. 1 0 4 6 5 4 9 E + C1 0 . 1 Q 4 9 7 1 8 E + 01 0 . 4 2 1 5 4 9 0 E - 0 2 0 . 4 2 2 7 3 7 3 E - 0 2 0 . 3 4 4 6 1 5 1 E +03 0.1 1 3 1 2 6 1E + 0 4 0TI3213'2'3'F + 0 1 0 . 1 3 2 5 7 0 2 E + 0 1 0 V 4 3 9 8 1 1 1 E - G 2 0 . 4 4 1 1 4 9 8 E - 0 2 0 . 4 1 3 5 3 8 4 E + 0 3 0 . 1 1 0 8 3 8 9 E + 0 4 0 . 1 6 2 3 2 6 6 E + 0 1 0 . 1 6 2 9 1 1 9 E + 0 1 0 . 4 5 9 2 8 5 7 E - C 2 0 . 4 6 0 7 8 9 1 E - 0 2 0 « 4 8 2 4 6 C I E + C 3 0. 1 0 8 5 5 1 4 E + 0 4 0 . 1 9 5 5 2 3 7 E + 01. 0. 1 9 6 2 8 6 4 E + 0 1 G . 4 8 0 0 8 2 5 E - C 2 0 . 4 0 1 7 6 6 2 E - 0 2 0 . 5 5 1 3 7 9 4 E + 02 0. 1 0 6 2 6 2 6 E + 04 -—  .,  .  ,  ^  ro  CO  124  Table 3  Initial  P a r a m e t e r s f o r P r e s s u r i z e d Water P r o b l e m  Conditions: u=0.0 m/s P=3.5 MPa a=0.001 ( v o i d f r a c t i o n ) h=l .0459 MJ/kg T=243.0 C  Boundary C o n d i t i o n s : a t z=0.0 m: u=0.0 m/s a t z=4.0 m: P=P i f u<a o t h e r w i s e u=a ( c h o k i n g  condition)  Parameters: A=804.2 mm dA/dz=0.0 m dZ/dz=0.0 D =32.0 mm, k?l=0.0 m" q=0.0 W/m  2  1  Constituative Equations: F=(4fx*+k.)u e D  1  where f = f r i c t i o n f a c t o r = 0 . 0 4 6 R ~ * f o r R >2000 =64/R f o r R <2000  2  2  e  2  e  x*=l x ( 0 - l ) { l 3 . 5 7 +  +  e  -°-  (1-e" - ^-^)}  0 0 8 8 4 Y  4  9 6  1/3 where  y=G/{gp (p -p )y } f  f  g  f  0 2 and here  8= ( v / u ) ' P / p g  f  f  g  G i s t h e mass v e l o c i t y  (kg/m s )  e  e  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0094637/manifest

Comment

Related Items