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 Carolina State University, 1971 H. S., Oregon State University, 1974 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES in the Department of Mechanical Engineering He Accept This Thesis as Conforming to the Required Standard THE UNIVERSITY OF BRITISH COLUMBIA December, 1978 (c) Hooshang Mozaffari-Nejad, 1978 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 degree a t the U n i v e r s i t y o f B r i t i s h Columbia, I agree 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 st u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y purposes 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 u n d e r s t o o d t h a t c o p y i n g 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 w r i t t e n p e r m i s s i o n . 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 Vancouver, Canada V6T 1W5 Date A p r f T ?, 1979 ABSTRACT Present research was undertaken to study the effect of pipe orientation during discharge of saturated water through a sudden opening. A n a l y t i c a l investigations have been made for various ranges of i n i t i a l pressure, including high pressures previously studied by others. A U m long 32 mm in diameter pipe has been used e n t i r e l y for t h i s investigation to f a c i l i t a t e comparison with available 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 above atmospheric, body forces were obviously expected to play an important r o l e . The e f f e c t became espe c i a l l y pronounced when d i r a c t l y comparing upflow versus downflow. Even at higher pressures, i t has 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 horizontal pipe. When the pipa i s oriented v e r t i c a l l y (up or down), t h i s w i l l no longer be true since turbulent mixing occurs. The research was undertaken because of i t s direct relevance to the loss-of-coolant accident (L3CA). In order to obtain licenses f o r operating reactors, certain safety standaris oust be met. The reactor must be shown to be safe i n normal operation as well as i n the case of a hypothetical accident. Different models have been used to analyze tha LOCA, simplest of which i s the equal velocity-equal temperature (EVET) model, which considers both phases i n equilibrium at a l l times. the EVET model (used e n t i r e l y i n t h i s work) was thought to be a good model for v e r t i c a l flow. Several, computations were made using the method of c h a r a c t e r i s t i c s to study the e f f e c t of gravity on pressure and void f r a c t i o n h i s t o r i e s at the closed end. Void f r a c t i o n , the r a t i o of vapor area to pipe area at any cross section, was investigated because of i t s direct relevance to mixture density and speed of sound. Also, t o t a l discharge or mass flow rate, a combined measure of density and v e l o c i t y , was calculated at the open end. Runs were made for water at i n i t i a l pressures of 1.2, 1.8, and 35 atm (at the top of the pipe). A l l cases started just at saturation with a small void f r a c t i o n of .001 i n a 4 m loog .032 m diameter pipe (Figure 1). The hydrostatic head of water column exerting an approximate pressure of 40000 Pa was comparable to the former pressures, but quite d i f f e r e n t from the l a t t e r of 35 atm. For 1.2 atm., or .12 MPa, the r a t i o of gravity force to pressure force (re c i p r o c a l of Froude number times Euler number) defined as H3P'" was about 33%; i t was 25% for the .18 MPa case, and only 1.1 Jt for the case when i n i t i a l pressure was 3. 5 MPa. The "3P" r a t i o was shown to be an important measure of the gravity e f f e c t . When GP=1.1%, there was n e l i g i b l e difference between the pressures along the pipe for the upflow and downflow i v case- In both lower i n i t i a l pressure cases, however, there seemed to be a clear difference between the pressures for upflow and downflow at the closed end. V TABLE OF CONTENTS Page ABSTRACT -. i i LIST OF FIGURES v i i LIST OF TABLES v i i i LIST OF SYMBOLS ix ACKNOWLEDGEMENTS xi I. INTRODUCTION - - 1 1.1 General .................................... 1 1.2 Previous Work 5 1.2.1 General 5 1.2.2 Previous experiments .............. 6 1.2.3 EVET Model 7 1.2.4 EVUT Model 8 1.2.5 UVOT Model 9 1.3 Scope of Present Work ...................... 10 II. GOVERNING EQUATIONS: 12 2.1 Conservation Equations ..................... 12 2.2 Compatibility Relations .................... 13 II 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 for Air 4.3 Results for Helium 4.4 Results for 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 Page 1 A Pipe at Quiescent State with 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 ..---......105 2 Pressure History at Closed End Shortly 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 ....................107 4 General Two-phase Algorithm to Advance Solution along C h a r a c t e r i s t i c ........................... 108 5 Have Interactions at the Boundaries .....................109 6 Expansion Fan at the End of Duct Suddenly Opened to Atmosphere .110 7 Pressure History at Closed End of Duct f o r A i r Problem ....................................111 8 Expansion fans for Helium Gun 112 9 Pressure at Closed End (P0=3.5 MPa) 113 10 Choking E f f e c t at Open End ....114 11 Void Fraction at Closed End (P0=0. 12 MPa) 115 12 Void Fraction at Closed End (P0=0.18 MPa) ......116 13 Pressure History at Closed End (Low Pressure Case) 117 14 Mass Flow Rate at the Open End 118 15 Computer Program ........................................119 v i i i i LIST OF TABLES Table Page 1 Parameters for Gas Dynamic Problems 122 2 Results f o r Helium 123 3 Parameters for High-pressure Water Problem .............. 124 NOMENCLATURE A = duct cross section flow area (m2) a = sound speed (m/s) aO = i n i t i a l sound speed (m/s) Al = volumetric concentration of vapor at a section d = pipe diameter (m) I = f r i c t i o n a l force per unit mass of f l u i d (N/kg) g = g r a v i t a t i o n a l acceleration (m/s2) h = mixture enthalpy = xhg * (1-x)hf (J/kg) 1 = pipe length (m) P = pressure (Pa) q = rate of heat added per unit mass of f l u i d (H/kg) R = mixture density = AlRg • (1-Al)Rf (kg/m3) s = the known dire c t i o n S = entropy (J/kg K) t = time (s) T = temperature (K) u = mixture v e l o c i t y (m/s) x = guality y = c h a r a c t e r i s t i c direction Z = elevation above an a r b i t r a r y datum (m) z = space coordinate (m) Subscripts g = vapor phase f = l i q u i d phase i = i n t e r f a c e w = wall ACKNOWLEDGEMENTS Author wishes to express h i s deepest appreciation to his supervisor. Dr. E. G. Hauptmann without whose help and advice this work would not have been possible. Useful discussions with Dr. P. H i l l at the i n i t i a l part of the research proved to be e s s e n t i a l and i s acknowledged. Many thanks to Dr. B. McDonald, Dr. D. Kawa, and Dr. W. Mathers of AECL (Whiteshell, Manitoba) for t h e i r many useful explanations and understanding throughout the research. Thanks to Gary Ford f o r proof-reading part of the thesis and Alan Steeves f o r double-checking some of the mathematical derivations; and to many others for having given t h e i r help and support whose names are too numerous to mention here. Support for t h i s research from the National Research Council of Canada i s g r a t e f u l l y acknowledged. Computing f a c i l i t i e s were provided by the University Computing Center. .1 I. INTRODUCTION 1. 1 General Because of i t s direct r e l a t i o n s h i p to the so-called l o s s of coolant accident, (LOC&), the analysis of transient steam-water flow and associated heat transfer phenomena (flow-boiling) i s focus of ouch current research. Interest i s centered around the development of more general conservation laws from which the governing equations for flow-boiling are derived and improved numerical techniques f o r solution of these governing equations. Presence of time as an additional variable i n the system makes transient two-phase flow an even more challenging area to explore than that of steady-state flow. For the transient problem the flow usually can not reach thermodynamic equiliorium, i t i s not f u l l y developed and the two phases are not uniformly d i s t r i b u t e d . Thus, i n t e r f a c i a l transport processes have to be considered. Unsteady three-dimensional problems are most representative of the complete flow-boiling situations i n reactor safety studies, however these are very d i f f i c u l t to handle. By use of su i t a b l e c r o s s - s e c t i o n a l area averaging techniques, a range of 2 unsteady one-dimensional analyses of two-phase flow can be considered. Models of the actual f l u i d flow regime are then formulated and applied to the problem to-obtain solutions. The simplest model used to handle such problems i s the equal velocity-equal temperature (EVET) model, which assumes the flow medium as being a homogeneous, bubbly mixture (pseudo-gas) . As a r e s u l t of t h i s assumption the vapor and l i q u i d phases are assumed to have equal v e l o c i t i e s and temperatures. A more complicated model i s the equal velocity-unequal temperature (EVUT) model, which assumes the 1 l i q u i d and vapor phases have equal v e l o c i t i e s but d i f f e r e n t temperatures. A r e l a t i v e l y sophisticated model (which has met with l i t t l e success to-date) i s the unequal velocity-unequal 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 i s conversion of part of the l i q u i d to vapor due to depressurization, occurs when l o c a l pressure drops to or below the saturation value for any given temperature. The latent heat needed to vaporize the l i q u i d i s supplied to the two-phase interface from the bulk of the l i q u i d . The time taken for f u l l nucleation to occur from the i n i t i a l depressurization state depends on the model chosen. Models such as the homogeneous equilibrium (EVET) model do not allow any time for 3 departure from thermal equilibrium. Nonegulibrium (such as EVOT) models, however, allow a few milliseconds "relaxation time". The present work uses the method of c h a r a c t e r i s t i c s or wave-tracinq technique, which has been used extensively to analyze gas dynamic problems, to analyze one-dimensional two-phase flow problems. This method i s based on a mesh of c h a r a c t e r i s t i c s dz/dt=u*a, u-a. I t greatly reduces the numerical d i f f i c u l t i e s present i n other f i n i t e difference c a l c u l a t i o n s , making i t useful f o r benchmark solutions and model te s t i n g . The grid of points at which the solution i s c a r r i e d out i s not the usual space-time grid but points of i n t e r s e c t i o n of c h a r a c t e r i s t i c l i n e s . These are l i n e s along which small-amplitude pressure disturbances propagate; thus the idea of wave tracing. The main advantage of a technique i n which the c h a r a c t e r i s t i c paths are followed throughout i s i t s lack of propagation of errors. Since d i s c o n t i n u i t i e s i n the sol u t i o n must propagate along the c h a r a c t e r i s t i c s , methods which do not trace the c h a r a c t e r i s t i c s w i l l d iffuse any such d i s c o n t i n u i t i e s . Since speed of sound i n two-phase mixture i s a very s e n s i t i v e function of void f r a c t i o n , staying within the scope of the above algorithm i s not always possible. In addition a c h a r a c t e r i s t i c changes i t s slope (in the z-t plane) abruptly i n 4 moving from the l i q u i d to the two-phase region. The reason for such an abrupt change in slope i s the sudden drop i n the l o c a l sound speed from the l i q u i d to two-phase value. In regions where such rapid but continuous changes occur as a res u l t of change i n speed of sound, one needs to relocate the advanced point on the saturation l i n e ( i . e . enthalpy changes from the l i q u i d value to two-phase value) or to use a very fine mesh. Interpolations are performed along the c h a r a c t e r i s t i c s to su b s t a n t i a l l y reduce the spurious numerical d i f f i c u l t i e s . Hancox (1) has shown that a wave-tracing technique has much less numerical d i f f u s i o n than f i n i t e difference methods. However i t i s much slower and more d i f f i c u l t to program. In addition because of the nature of i t e r a t i v e technique (Newton-Baphson method) required to obtain the solution at each point, i t i s extermely sensitive to i r r e q u l a r i t i e s i n the f l u i d properties such as £, SR/SP, 3R/8h, h, or externally s p e c i f i e d quantities such as dA/dz, dZ/dz, d, and K / l . I f such terms are discontinuous, d i f f i c u l t i e s can be expected. I n i t i a l l y , the research set out to study the p o s s i b i l i t y of applying the method of hodograph transformations (used successfully to analyze i s e n t r o p i c gas-dynamic problems) to two-y phase problems (appendix C and D). However the p a r t i a l d i f f e r e n t i a l equations proved to have undetermined c o e f f i c i e n t s . 5 Hence, closed-form a n a l y t i c a l solutions could not be obtained. The present research has focussed on the e f f e c t of orientation of the body force (gravity) at a range of i n i t i a l pressures. 1.2 Previous work 1. 2. 1 General Solving a complete unsteady three-dimensional flow-boiling problem i s a formidable task. Thus averaging methods have been attempted by Panton (2), Delhaye(3), Vernier (4) , Banerjee(5), Mathers (6), Agee(7), Lyczkowski (8) to make the problem simpler. These averaging methods change the problem from an unsteady three-dimensional to an unsteady one-dimensional problem. A set of constituative relationships are then formed, which include transfer terms (mass, momentum, and energy) between the phases and between each phase and the boundaries. These are complementary sets of equations to the governing equations (continuity, momentum, and energy) , and are c a l l e d "external* 1 i n t e r f a c e c o n s t i t u t a t i v e r e l a t i o n s h i p s . To handle the flow problem, one needs to solve the derived one-dimensional system of equations numerically. Predicted r e s u l t s from numerical c a l c u l a t i o n s are then compared with the 6 experiments to check accuracy of the obtained s o l u t i o n . 1.2.2 Previous Experiments Frontier work i n flow-boiling experiments were done by Edwards and 0 , B r i e n ( 9 ) f who studied the rapid depressurization of subcooled l i q u i d from long, horizontal pipes. The pressure range for these experiments was from 500 to 2500 psig, and the temperatures varied from 467 F to 636 F. A l l pipes chosen for these experiments were 4 m long (with various diameters). For a l l runs, the pressure gage at the closed end of the pipe recorded a small dip 1 ms a f t e r rupture had occurred. Thermal noneguilibrium was explained to be the reason f o r the dip. Edwards and Hather(10) did another set of blow-down experiments with 4.0 meter long pipes (with various diameters). Most of th e i r tests started from a subcooled condition, 35 bars above saturation pressure. Another experiment was done by Uancox and Baner jee (11) , which was e s s e n t i a l l y comparable to that of Edwards and Mather for 32 mm diameter pipe. The experiments mentioned above showed a s i m i l a r depressurization trend at the closed end. The blow-down time i n a l l cases was l e s s than a second. Numerical models such as the equal velocity-equal temperature model, or the equal velocity-unequal temperature 7 model have predicted lower depressurization times than those measured experimentally. The reason for t h i s discrepancy has been attributed to the assumption of equal v e l o c i t i e s for both phases. 1.2.3 The Equal Velocity-Equal Temperature (EVET) Model: This i s the simplest model used to predict solution to reactor blow-down problems. B r i t t a i n (12) , Martini (13), Turner(14), Mathers(15), Karwat(16), and Moore(17) incorporated t h i s method i n their solutions. In t h i s model phases have equal v e l o c i t i e s and temperatures (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 transfer of momentum and heat from the wall to the f l u i d . Edwards and Mather (10) compared t h e i r "EVET" model with the i r own experiments. Generally speaking comparison between the two was not good. F i r s t the i n i t i a l pressure dip at the closed end (measured experimentally) was not predicted numerically because phases were assumed to have equal temperatures at a l l times. Second the predicted long-term depressurization occurred sooner that i n the actual s i t u a t i o n since phase v e l o c i t i e s were considered equal i n the model. In the actual blow-down experiments (oecause of flow 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 predicted the rig h t trend, but improvements were needed to make t h i s model more accurate. In p a r t i c u l a r , extensions to allow for thermal nonequilibrium and unequal v e l o c i t i e s of phases vere needed. I t was clear that improvements could not be made by simply changing of c o n stituative r e l a t i o n s . 1.2-4 The Egual Velocity-Unequal Temperature (EVUT) Model This model assumes that vapor and l i q u i d phases have the same v e l o c i t i e s but different temperatures. There are ce r t a i n situations where t h i s model seems necessary. Early depressurization of blow-down phenomenon may be mentioned as an example. During tiie blow-down, vapor and l i q u i d v e l o c i t i e s are i n i t i a l l y s i m i l a r but the l i q u i d may be extremely superheated. Edwards and O'Brien (9) showed that during rapid depressurization of high enthalpy water, pressure undershoots the saturation values. Banerjee(18) and Ferch(19), with a simple rela t i o n s h i p of i n t e r f a c i a l heat transfer, have shown close agreement between t h e i r numerical code and the experimentally measured undershoot (Figure 2) . They chose two values of time constants for heat transfer between phases, the smaller of which gave better agreement with the experiments. 9 So far i t has become clear that simpler representations of transient flow-boiling, (EVET and EVDT models) do predict the trend of the experiments guite well and could be s u f f i c i e n t f o r many cases. Another l e v e l of sophis t i c a t i o n i s reguired for more complex predictions. This l e v e l i s met with unequal the velocity-unequal temperature (UVUT) model, which produces other d i f f i c u l t i e s . 1.2.5 The Unequal Velocity-Unequal Temperature (UVUT) Model Sets of equations of t h i s type (separate continuity, momentum, and energy equations f o r each phase) produce imaginary c h a r a c t e r i s t i c s . I t was demonstrated by Bryce(20) that whatever the school of thought, the simplest UVUT model (using a form of " s l i p " ) did not give converging solutions as mesh size and time step were refined. Thus although numerically dispersed solutions are av 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 (UVUTUP) model that has real c h a r a c t e r i s t i c s and gave r e a l i s t i c propagation speeds for separate flows. Stuhmiller(21) also suggested a somewhat d i f f e r e n t approach with i n t e r f a c i a l pressure d i f f e r e n t from phases. Under some conditions t h i s model also led to a system of equations with r e a l c h a r a c t e r i s t i c s . 10 None of these unequal systems with r e a l c h a r a c t e r i s t i c s are e n t i r e l y convincing since there are conditions where d i f f i c u l t i e s show up. More work i s needed not only to advance the model, but also to derive a r e l i a b l e set of constituative r e l a t i o n s . 1.3 Scope of Present Work The present research inte r e s t has been to see the e f f e c t of gravity in the blow-down problem, using EVET model and assuming that a l l problems started off with a small (figure 1), say .001, void f r a c t i o n (void f r a c t i o n i s r a t i o of the area of void to t o t a l area at any pipe cross-section). I t would seem proper to use a more sophisticated and complicated model such as EVUT (which allows departure form thermal eguilibrium). However from the r e s u l t s obtained by Hancox, and Mathers (1), the gr a v i t y e f f e c t s prevailed around 200-300 ms from the time of rupture and were unimportant to b o i l i n g e f f e c t s at the closed end. At high pressure regions, horizontal flow i s s t r a t i f i e d . S t r a t i f i c a t i o n , however, i s not expected to be s i g n i f i c a n t i n v e r t i c a l l y upward or downward flows. Hence the gravity e f f e c t s were expected to p r e v a i l even at high pressures. The research, thus, encompassed solving a range of problems including the 11 reactor f a i l u r e problems. Tnermodynamic derivations have been made to obtain the sound speed as a function of pressure and temperature instead of pressure and enthalpy (since the sound speed i s a function of density and density a function of pressure and enthalpy). The steam-water properties used are from continuous functions i n Keenan and Keyes (22). 12 IX. GOVERNING EQUATIONS 2. 1 Conservation Equations: The homogeneous equilibrium (EVET) model has been adopted to i l l u s t r a t e the method of c h a r a c t e r i s t i c s solution technique, used extensively in gas-dynamic problems. One-dimensional transient f i e l d equations r e s u l t when flow variables are averaged over the cross-section, and hence vary only i n the d i r e c t i o n of flow. In t h i s model i t i s also assumed that the l i q u i d and vapor v e l o c i t i e s and temperatures are the same at any cross-section. The resultant conservation equations may be rearranged into the following:* (from Appendix A) Du/Dt+(aP/a»z)/R = -F - g(dZ/dz) = C l , (2.1) Dh/Dt+ a 2 (du/dz) = a 2 [ (q*uF) dR/dP| h - (u/A) (dA/dz) J = C2, (2.2) DP/Dt*a 2R(iu/8z) = a 2 t - (q+uF) 2>R/dh| p- (uR/A) (dA/dz) j=C3, (2.3) where D/Dt = 8/2>t • u(3/tfz) , a 2 =1./ £ SR/dP|h • (dB/dh | P)/R J, 1 dfi/ap|h = p a r t i a l d e r i v a t i v e of density with respect to pressure at constant enthalpy 13 Equation (2.1) i s the momentum, (2.2) the continuity, and (2.3) the energy equation. In deriving these equations density i s used i n the form B = fi (P,h). Equations (2.1) throuqh (2.3) are quasi-linear p a r t i a l d i f f e r e n t i a l equations of the hyperbolic type. From the theory of hyperbolic equations (for example see Appendix B) , the c h a r a c t e r i s t i c directions can be derived: dz/dt = u, u*a, u-a. These represent repectively the p a r t i c l e path, the pressure disturbance propagation in the d i r e c t i o n of flow, and the pressure disturbance propagation opposite to d i r e c t i o n of flow. 2.2 Compatibility Relations The ordinary d i f f e r e n t i a l equations, c a l l e d compatibility equations, which apply along these c h a r a c t e r i s t i c s are: (Appendix B) along dz/dt = u+a, dP + fia du = (C3 • B a d ) dt = K dt, along dz/dt = u-a, dP - Ba du = (C3-BaC1)dt = L dt. 14 along dz/dt = u: dP - Bdh = (C3 - BC2)dt = M dt where K, L , and H are constants. P h y s i c a l l y , the f i r s t two eguations represent a step change in pressure as related to ve l o c i t y change i n crossing a c h a r a c t e r i s t i c l i n e (a form of momentum equation). The f i n a l equation r e l a t e s a step change in pressure as a resu l t of enthalpy change in crossing a c h a r a c t e r i s t i c l i n e (a form of energy eguation). 15 III- NUMERICAL PfiOCEDUEES 3.1 General A f i n i t e difference algorithm was used to advance solutions along the c h a r a c t e r i s t i c s i n the z-t plane from known nodes (Figure 3). C h a r a c t e r i s t i c s l i n e s (along which small disturbances travel) chosen were the l e f t running, the r i g h t running, and the p a r t i c l e path (or trace of the particles) . Geometry chosen for a l l the computer runs was a pipe 4 a long and .032 m in diameter to compare with available experiments. I n i t i a l l y , the pipe contained high pressure water with a .001 void f r a c t i o n - At time t=0* one end of the pipe was suddenly ruptured and the subsequent blow-down occurred (Figure 1). Special f i n i t e difference algorithms were also used to advance solutions on the boundaries. The boundaries considered were the open end and the closed end. The task of advancing the solution was greatly reduced because of s t a r t i n g the blow-down with small i n i t i a l Al (void t r a c t i o n ) . As a r e s u l t of t h i s s t a r t i n g two-pnase condition, a s p e c i a l t r a n s i t i o n algorithm (from one-phase l i q u i d to two-phase mixture) was not needed. 3.2 Two-phase Algorithm 16 The algorithm used f o r solving the point of int e r a c t i o n of opposite-running waves inside the pipe was c a l l e d "IIPHAS". There were 10 unknowns z,t,u,P, and h related to points 3 and 0 (Figure 4) . Ten nonlinear equations were used to solve for these unknowns: six of these equations were the c h a r a c t e r i s t i c directions and t h e i r corresponding compatibility r e l a t i o n s ; the other 4 were obtained using interpolations between known nodes 1 and 2- The Newton-fiaphson method was used to solve these nonlinear equations. Only the solutions f o r f i v e properties at point 3 were kept; solutions f o r point 0 were discarded ( i . e . , the p a r t i c l e s were not mapped continuously i n z-t plane as the algorithm progressed). 3-3 Open-End Algorithm: For the open end boundary the routine "OPEN" was used, which was exactly the same as MIIPHAS M with the difference that a choking c r i t e r i o n replaced the left-running c h a r a c t e r i s t i c and i t s compatibility eguation. This required that P3=Pfi i f u3<a3, otherwise u3=a3 and P3 was found i t e r a t i v e l y (Figure 5-b). 3.4 Closed-End Algorithm For the closed end of the 4 m long test tube, the routine MCL0SED w set up the c o e f f i c i e n t matrix to solve 5 equations and 17 5 unknowns for point 3, using the Hewton-Baphson method of solution (Figure 5.a). Here only left-running and p a r t i c l e path c h a r a c t e r i s t i c s were used to determine the advanced s o l u t i o n (meaning the properties at point 3). 3.5 Expansion Fans For high i n i t i a l pressure cases, an i n i t i a l expansion fan (consisting of 8 waves for water and 11 waves f o r a i r and helium) emanating from the rupture end at time zero was used to analyze the problem and solutions were advanced at grids of these waves in t e r s e c t i n g one another (Figure 3). To set i n i t i a l conditions on the waves, pressure dropped i s e n t r o p i c a l l y from i t s o r i g i n a l value to atmospheric at equal increments on 8 waves. The vel o c i t y u was found from dP*Ba2du=0, using the pressure drop dp from one fan to the one immediately succeeding i t . Enthalpy was found from dP-Bdh=0 (Figure 6). The formulae were ite r a t e d u n t i l convergence on u=a occurred (within a tolerance of .5 m/s). The expansion fan method did not seem to work s a t i s f a c t o r i l y at low starting pressures (when the gravity e f f e c t was important), because hydrostatic head was not accounted for appropriately. As a re s u l t of t h i s shortcoming the c h a r a c t e r i s t i c s crossed one another and caused d i f f i c u l t i e s , producing meaningless r e s u l t s . 18 Since i t was d i f f i c u l t at these pressures to i n i t i a l i z e the problem properly, the idea of using an expansion fan was dropped and 8 wavelets (gradual pressure reduction) were chosen i n i t i a l l y along the t=0 axis. This method was successful because i t enabled posing the hydrostatic e f f e c t s c o r r e c t l y . 19 IV RESULTS AND DISCUSSIONS 4.1 General Tiie method of c h a r a c t e r i s t i c s algorithm was used to make computer runs for a i r , helium, and water. Hater was considered at an i n i t i a l l y saturated condition with 0.001 void f r a c t i o n . This eliminated the spurious d i f f i c u l t y in going from one to two-phases ( i . e . relocating advance point on saturation line) . Accuracy of the program was checked by making preliminary runs with a i r and helium. Eleven fans were chosen to analyze the a i r and helium problems, and eight waves to analyze the high pressure water problem. These expansion waves emanated from the open end at t=0«- due to a sudden rupture. The expansion fan method worked well f o r high pressures since gravity did not play an important r o l e (ratio of gravity to pressure ,GP, was l e s s than 1.5%). Using more than 8 waves for water did not seem to make much difference i n the accuracy of re s u l t s . For the water runs at lower pressures, however, the expansion fan method did not prove to be the correct way of analyzing the problem. The r a t i o GP was 33% fo r P0=1.2 atm, and 25* f o r P0=1.8 atm (meaning measurable r a t i o of gravity to the pressure f o r c e ) . In both cases the strategy was to i n i t i a l i z e the problem along the t=0 20 axis to account for hydrostatic head. The pressure at the open end was allowed to drop gradually, so that a centered expansion fan did not e x i s t . Several rates of depressurization were t r i e d . I t was found that for the .12 MPa case l e s s than 5 ms depressurization rate did not a f f e c t the numerical c a l c u l a t i o n s , when the expansion fan method was used for P0=1. 2 atm, subsonic conditions prevailed because the ve l o c i t y at the open end was much les s than the l o c a l speed of sound. The expansion fan did not seem organized. Some of the fans crossed over one another causing computational d i f f i c u l t i e s . I n i t i a l i z i n g the problem from t=0 a x i s , however, seemed to work guite well. DZ/dz=1 defined upflow, and dz/dz=-1 defined downflow. Each run was carried for 30 time advances of about 140 ms. 4.2 fiesults for Air An i s e n t r o p i c problem with a i r was chosen to check the accuracy of the algorithm. A 4 meter long pipe with 32 mm diameter was depressurized rapidly from an i n i t i a l pressure of 7. MPa to .738 of i t s o r i g i n a l value (Table 1). The flow was subsonic. Isentropic relations including state eguations were used to i n i t i a l i z e the problem along the fans at t=0 + . To 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 can not be completely i s e n t r o p i c . Hence a very small amount of f r i c t i o n 21 was introduced as an approximation of isen t r o p i c flow. The f r i c t i o n factor chosen was f=1.E-6 (making f r i c t i o n force F=(2f/d)u^) very small to model isentr o p i c s i t u a t i o n . The algorithm progressed i n z-t plane with interactions of re f l e c t e d waves from the closed end with oncoming waves emanating from the open end. Results compared quite well with those of Owczarek (23) for i s e n t r o p i c flow (Figure 7). 4.3 Results f o r Helium Isentropic problems may generally be solved a n a l y t i c a l l y by hodograph transformations (using Biemann Invariants). The solutions are hypergeometric functions (Appendix C). In the case of helium, with k=5/3, the analytic solution f o r a ruptured tube i s quite simple. Parkinson (24) showed that the space and time (z and t) can be expressed as ve l o c i t y and sound speed (u and a) as follows: z(u,a) = 2tu(9a2-3u2*27a02 J/(27a 3) t(u,a) = 6 (9a 2 -u2*9a02)/(27a3) where aO and a are the i n i t i a l and l o c a l speed of sound, respectively (Figure 8)-The same i n i t i a l condition as for the a i r problem was used to solve the helium problem numerically and results were compared with a n a l y t i c a l values. Accuracy was to within 0.1% 22 (Table 2). The f r i c t i o n f a c t o r was again chosen as 1.E-6 to 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 . 4.4 Results for Sater Three runs were made with water. For the high pressure condition (P0=3.5 HPa, and T0=243 C), eight expansion fans emanating from the open end were chosen. I n i t i a l conditions were set on leading edge of the fan (Table 3). On the t r a i l i n g edge of the fan ( v e r t i c a l l i n e i n z-t plane), choking occurred. Values f o r intermediate waves were it e r a t e d i s e n t r o p i c a l l y . Hence i n i t i a l conditions were set up on the expansion fans. Boundary conditions were considered at closed and open ends. At the closed end the velocity was zero while at the open end pressure was set equal to ambient i f the flow was subsonic. I f choking occurred (u=a), l o c a l pressure was solved f o r i t e r a t i v e l y . A two-phase m u l t i p l i e r between w 0 n and M 1 " was chosen to determine f r i c t i o n factor (25). There did not seem to be much difference between upflow and downflow i n the pressure hi s t o r y at closed end (Figure 9). Hence even though flows were fundamentally d i f f e r e n t from horizontal to v e r t i c a l orientations (because of horizontal flow being s t r a t i f i e d , but v e r t i c a l flow being turbulent mixing), very l i t t l e gravity effect was sensed i n the pressure h i s t o r y . 23 In t h i s case GP was small (1.1%). The hump in the pressure curve at about 90 ms was f e l t to be due to the i n t e r a c t i o n of compression waves in the medium. The speed of sound was a very sen s i t i v e function of void f r a c t i o n and depeneded l i t t l e on pressure. Its minimum value was about 110 m/s when void f r a c t i o n was .5 (Figure 10), Pressure decreased on each successive fan as did enthalpy, temperature, and density. Sound speed i n i t i a l l y decreased on each successive fan (as void f r a c t i o n increased to .5) and then increased to 177 m/s (when void f r a c t i o n increased from -5 to 1) . Quality throughout the run was less than 0.3. Void f r a c t i o n , though, increased quickly to 0.9 and higher values. For P0= 1.2 atm ( s l i g h t l y more than atmospheric pressure), eight wavelets set the i n i t i a l conditions from the t=0 axis. The hydrostatic head of water column exerted a pressure of approximately 40000 Pa, comparable to i n i t i a l pressures at the top of the pipe i n some cases. For .12 MPa, the r a t i o of gravity force to pressure force, GP, was 33%; i t was 25% f o r .18 MPa case. The «*GPM r a t i o showed to be an important measure of the gravity e f f e c t . In both lower pressure cases, i . e . , .12 and .18 MPa, the void f r a c t i o n increased gradually with time at the closed end for both up and down flow (Figure 11 and 12). 24 Pressure at the closed end f o r downflow increased s l i g h t l y to above i t s o r i g i n a l value at around 90 os for both .12 and .18 MPa pressure cases (Figure 13) . The hump was thought to be due to wave i n t e r a c t i o n of compression waves coming from the high pressure region (bottom of pipe). 25 V CONCLUSIONS It was found during the course of the present research that the method of c h a r a c t e r i s t i c s i s a very accurate numerical technigue for analyzing two-phase problems . This was shown by comparing numerical r e s u l t s obtained f o r helium with the available closed form a n a l y t i c a l solution (hypergeometric ser i e s solution). Tne method of c h a r a c t e r i s t i c s , however, i s a c o s t l y technigue since convergence f o r each advanced point i n the z-t plane i s required. Hence i t i s a good method f o r benchmark comparisons with other numerical techniques. Reynolds number f o r the blow-down water problems was of the order of 10 6. Hence as many other authors have done, the vi s c o s i t y e f f e c t s were considered n e g l i g i b l e i n the flow medium and were neglected from the analysis. Viscosity e f f e c t s prevailed, however, at the wall i n the form of f r i c t i o n . Comparison between experimental r e s u l t s of v e r t i c a l l y oriented pipes and the predicted numerical values should be made as a next step to see i f the EVET model gives f a i r prediction f o r such flow regimes (since v e r t i c a l flows involve turbulent mixing and are not s t r a t i f i e d as are horizontal flows). In general body forces were found to s i g n i f i c a n t l y a f f e c t the flow, especially at lower i n i t i a l pressures. There were pronounced differences i n void 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 the closed end of the pipe. Due to gravity the mass flow rates f o r upflow and downflow, at the open end of the pipe were quite different for the case of .12 HPa, whereas they were rather s i m i l a r for the case of .18 HPa (when pressure was higher). Accurate property eguations are needed otherwise the r e s u l t s f o r sound speed are meaningless (this was underlined in t h i s work while trying to calculate sonic speeds from available tabulated thermodynamic property values for Freon-114). Generally for high i n i t i a l pressure conditions, the expansion fan method was found to be adeguate. This i s because at high i n i t i a l pressures gravity e f f e c t was rather unimportant. Hence at any orientation the i n i t i a l pressure was assumed uniform along the pipes. Expansion fans could not be used, however, for low i n i t i a l pressure cases, since i n i t i a l pressure i n the pipe was not uniform because of hydrostatic head. Instead of an expansion fan, the pressure at the open end was allowed to gradually drop. This was found adeguate for simulating a rapid opening. I t was found that less than 5 ms for the opening did not make a difference i n f i n a l r e s u l t s i n the case of .12 HPa. For future work, further research of the unfinished work i n hodograph transformation f o r non-isentropic problems i s 27 suggested (Appendix D). This B u s t be done by mixed experimental and perturbation techniques giving r i s e to semi-closed form solutions. I t i s evident from the l i t e r a t u r e survey that so far in a l l previous work, the two-phase and s u p e r c r i t i c a l regions have been well studied. The n e a r - c r i t i c a l region could be an i n t e r e s t i n g area for future research. Similar to the work done by Hancox(1), in the present research, the f i n a l form of the governing equations are derived making use of density i n the form R=R (P,h). The reason for such a choice i s because in the compressed l i q u i d state, density i s a rather i n s e n s i t i v e function of pressure and enthalpy and hence the d i f f u s i o n errors resulted from forming the density derivatives are quite small. It may prove, with further studies, to be advantageous to use the density i n alternate 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 also be pursued i n the area of developing a more sophisticated technique for handling complicated pipe net-works with associated changes of area and with the shock-waves present. 28 VI- REFERENCES 1 Hancox, W. T., Bathers, W. G., and Kawa, D. , Analysis of Transient Flow-Boiling; Application of the Method of Ch a r a c t e r i s t i c s , AIChE Paper 42, 15th National Heat Transfer Conference, August 1975. 2 Panton, B., "Flow Properties f o r the Continuum Viewpoint of a Non-Equilibrium Gas-Particle Mixture, J. F l u i d Mech., 31, 273-303, 1968. 3 Delhaye, J.M. , and Achard, J. L. ,On the Averaging Operators Introduced i n Two-Phase Flow Modelling, Invited Paper presented at 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Toronto, August 1976. 4 Vernier, P. , and Delhaye, J. H. ,General Two-Phase Flow Equations Applied to the Thermo-Dynamics of Boiling Water Reactors, Energie Primaire 4, No- 1-2, 1968-5 Banerjee, S., Ferch, R., Mathers, W.G., McDonald, B.H., Dynamics of Two-Phase Flow i n a Duct, Proc. Sixth International Heat Transfer Conference, Toronto, August 1978. 6 Mathers, S.G., Ferch, R.L., Hancox, W.T., and McDonald, B-H., Equations for Transient Flow-Boiling i n a Duct, Invited Paper Presented at 2nd CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Paris, 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 S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Paris, June 1978. 8 Lyczkowski, R.H., Theoretical Bases of the Drift-Flux F i e l d Eguations and Vapor D r i f t - V e l o c i t y , Proc. Of Sixth 29 International Heat Transfer Conference, Toronto, August 1978-9 Edwards, A. B. , and O'Brien, T. P., Studies of Phenomena Connected with Depressurization of Mater Beactors, Journal BNES Vol. 9, pp. 125-135, A p r i l 1970. 10 Edwards, A. B., and Mather, D. J . , Some UK Studies Belated to Loss of Coolant Accident, Topical Meeting on Mater Beactor Safety, Salt Lake City OSAEC Report CONF-730304 (1973.) 11 Banerjee, S., and Hancox, M.T., On the Development of Methods f o r Analyzing Transient Two-Phase Flow, Int. J. Multiphase Flow, i n press, 1978-12 B r i t t a i n , I., and Fayers, F.J., A Beview of Becent Mork i n the Modelling and Solution of the Hydrodynamic Eguations i n the LOCA Code, EELAP-OK, Invited Paper presented at the 1st CSNI S p e c i a l i s t s 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 One-dimensional Transient Two-Phase Flow Model and Its I m p l i c i t F i n i t e Difference Solution, Invited Paper Presented at the 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Toronto, August 1976. 14 Turner, H.J., and Trimble, G.D., Calculation of Transient Two-Phase Flow, Invited Paper presented at the 1st CSNI S p e c i a l i s t s Meeting on Transient 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 Difference Solutions to the Transient Flow-. Boi l i n g Eguations, Invited Paper presented at the 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Toronto, August 1976. 30 16 Karwat, H., BBOCH-S, A Code to Investigate Blowdown of Boiling Hater Beactor Systems, Hucl. Eng. Design 9, 363-381, 1968. 17 Moore, K.V., and Bettig, W.H. , RELAP4-A Computer Program for Transient Thermal-Hydraulic Analysis, Aerojet Nuclear Company Beport ANCB-1127, 197 3. 18 Banerjee, S., Ferch, B . L., and Mathers, I . G. , the Dynamics of Two-phase Flow i n a Duct, Paper Submitted to 6th International Heat Transfer Conference, August 1978. 19 Ferch, B.L., Method of C h a r a c t e r i s t i c s Solutions for Noneguilibrium Transient Flow B o i l i n g , Submitted to Int. J . Multiphase Flow, 1978-20 Bryce, W.H., Determination of the Hyperbolic Begime of the Hydrodynamic Equations Modelled i n the LOCA Code EELAP-OK, Invited Paper Presented at the 1st CSNI S p e c i a l i s t s Meeting On Transient Two-Phase Flow, Toronto, August 1976. 21 Stuhmiller, J.H., The Influence of I n t e r f a c i a l Pressure Forces on the Character of Two-Phase Flow Model Eguations, Int. J. Multiphase Flow, 3, 551-560, 1977. 22 Keenan, J.H., and Keyes, F.G., Thermodynamic Properties of Steam, Wiley, N.T., 1958. 23 Owczarek,J.A., Fundamentals of Gas Dynamics, International Textbook Co., p 330, 1964. 24 Parkinson,G.V., The Short-Chambered Helium Gas Gun, Defence Besearch Board of Canada, Technical Memo. TMAB-60, 1960. 25 Hancox, 8. T., and N i c o l l , W. B. , Unpublished Westinghouse Canada Ltd., Beport, 1972. 31 26 Carnahan, L.W., Applied Numerical Method, Wiley, p 321, 1969. 27 F r i t z , John, "A Note on •Improper' Problems in P a r t i a l D i f f e r e n t i a l Equations," Communication on Pure and Applied Math., 7, 591-594, 1955. 28 fiichtmeyer, fi. D. And Morton, K. W. , "Difference Methods for I n i t i a l Value Problems," Interscience, New York (1967). 29 Payne, L. E., Bounds in the Cauchy Problem f o r the Laplace Equation , Arch, f o r Bat. Mech. and Anal., 5, 35-45 (1960). 30 Schaefer, P. 8 . , On Uniqueness, S t a b i l i t y and Pointwise Estimates in the Cauchy Problem for Coupled E l l i p t i c Equations, Quart. Of Appl. Math., 321-328, 1973. 31 Lax, P. D., D i f f e r e n t i a l Equations, Difference Equations and Matrix Theory, Comm. Pure and Appl. Math., 11, p 175-194, 1958. 32 J a r v i s , Stephen J r . , "On the Formulation and Numerical Evaluation of a Set of Two-Phase Flow Equations Modelling the Cooldown Process," BBS Tech. Note No. 301, January 4, 1965. 33 Siegmann, E. B., Theoretical Study of Transient Sodium B o i l i n g i n Beactor Coolant Channels U t i l i z i n g a Compressible Flow Models" Argonne Nat. Lab., ANL-7842, July 1971. 34 Boure, J. A., "GEVATBAN: A computer Program to 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. , Giot, H. And Reocreux, M. , Choking Flows and Propagation of Small Disturbances, 32 European Two Phase Flow Group Meeting, Brussels, Paper FL, 1973. 36 Gidaspow, D i m i t r i , Lyczkowski, B. W. , Solbrig, 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 Society Transactions, 17, p 249-250, November 1973. 37 Boure, J. A., Dynamics des Ecoulements Diphasiques -: Propagation de Petites Perturbations, Grenoble CEA-R-4456, October 1973. 38 In h i s writeup Professor Gidaspow shows t h i s a n a l y t i c a l l y for the incompressilble l i m i t . Numerical r e s u l t s show that the same i s true i n the compressible case; see Gidaspow et. A l . , Trans. ANS Winter Meeting (San Francisco, November 1973) , p 249. 39 Hancox, W.T., and Banerjee, S-, Numerical Standards i n Flow B o i l i n g Analysis, Nuclear S c i . and Eng. 64, p 106-124, 1977. 40 Wallis, G. B., One-dimensional Two-phase Flow, Ch. 9, McGraw B i l l , *969. 41 Shapiro, A. H., The Dynamic And Thermodynamics Of Compressible F l u i d Flow, Vol. 2, Ch. 24, Bonald Press Co. 42 Necmi, S. , and Hancox, W. T. , An Experimental and T h e o r i t i c a l Investigation of Blow-Down from a Horizontal Pipe, 6th Int. Heat Transfer Conference, Toronto, 1978. 43 Bousseau, J . C , Blowdown Experiments and Interpretation, European Two-phase Flow Group Meeting, Haifa, I s r a e l , 1975. 44 Mathers, W. G., Zuzak, 8 . W., McDonald, B. G., and Hancox, W. T., on F i n i t e Difference Solutions to the 33 Transient Flow-Boiling Equations, Proceedings of the OECD NEA S p e c i a l i t s Meeting on Transient Two-phase Flow Toronto, August. 1976 45 Baker, 0., the O i l and Gas Journal, Vol. 53, p. 175-195, 1954. 46 Banerjee, S., Hancox, H.T., J e f f r i e s , fi.B., and Sulatisky, M-I.,. Transient Two-Phase Flow during Blowdown with Heat Addition, AIChE Symposium Series, Heat Transfer Vol. 174, 1978. 47 Premoli, A., and Hancox, V.T., An Experimental Investigation of Subcooled Blow-down with Heat Addition, Invited Paper presented at the 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow. 48 Arrison,N.L., Hancox,H.T-, Sulatisky,M.T., and Banerjee,S. , Blowdown of a Recirculating Loop with Heat Addition, Heat and F l u i d FlOw i n Water Beactor Safety, I n s t i t u t i o n of Mech. Engs. London, p 77-83, 1977. 49 Hadamard,J., Lectures on Cauchy's Problem i n Linear P a r t i a l D i f f e r e n t i a l Eguations, Yale Dniversity Press, New Haven, Conn. , 19 23. 50 Boure,J.A.,and Latrobe, A., On well-Posedness of Two-Phase Flow Problems, Paper No. 76-CSME/ CSCHE-9 presented at the 16th National Heat Transfer Conference, St. Louis, Missouri, August 1978. 51 Bound Table Discussion on Two-Phase Flow, F i f t h International Heat Transfer Conference, Tokyo, 1974. 52 Soo, S.L. , Multiphase Mechnics of Single Component Two-Phase Flow, Physics F l u i d s , 20, p 568-570, 1977. 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 Equations of Components i n a Mixture, Int. J . Multiphase Flow, i n press. 54 Boure, J.A., Fritte,A.A., Giot,M.H., and Beocreux,M.L., Highlights of Two-phase C r i t i c a l Flow, Int. J . Multiphase Flow, 3, p 1-22, 1976. 55 Boure,J.A., On a Unified Presentation of the Non-Equilibrium Two-Phase Flow Models, Non-Equilibrium Two-Phase Flows ASME, N.Y., p 1-10, 1975. 56 Truesdell,C., and Toupin,B., The C l a s s i c a l F i e l d Theories, Handbuch der Physik, Vol. 3/1, Springer Verlag, B e r l i n , 1960. 57 Standard,G., The mass, momentum, and energy equations f o r heterogeneous flow systems, Chem. Engng. S c i . , 19, p 227, 1964. 58 Solbrig,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 Describe Steam-water Behaviour i n Nuclear Safety Studies, Paper Presented at 15th National Heat Transfer conference, San Francisco, 1975. 59 Banerjee,S., A Surface Benewal Model for I n t e r f a c i a l Heat and Mass Transfer f o r Two-Phase Flow, Int. J . Multiphase Flow, i n press. 60 Yanenka, N.M., The Method of F r a c t i o n a l Time Steps: The Solution of Problems of Mathematical Physics i n Several Variables, ed. M- Holt, Springer-Verlag, N. Y., 1971. 61 Marchuk,G.I., A Theoretical Heather Forecasting Model, Dok. Akad, S c i . USSB, 155, p 10-12, 1965. 35 62 Douglas, J . , Alternating Direction Methods for Three-Space Variables, Numerische Hathematik, 4, p 41-63, 1962. 63 Lax,P.D., Beak Solutions of Nonlinear Hyperbolic Eguations and Their Numerical Computation, Comm. Pure Appl. Math., 7, p 159-193, 1954. 64 Lax,P.D. , and Hendroff,B., Difference Schemes with High Order of Accuracy for Solving Hyperbolic Equations, Comm. Pure Appl- Math., 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 for Hydrodynamic Calculations, LASL Beport No. LA-2139, N. Mexico, 1957. 66 Gentry,B.A., Martin,B.E., and Daly, B.J., An Eulerian Differencing Method f o r Unsteady Compressible Flow Problems, J . Coop. Phsics, 1 , 87-118, 1966. 67 HcGormack,B.W., Numerical Solution of the Interaction of a Shock . Have with a Laminar Boundary Layer, i n Proc. 2nd Int. Conf. on Num. Methods i n Fluid Dynamics, ed. By M. Holt, Vol. 8, Springer-Verlag, N.Y., 1971. 68 Porching, T.A., Murphy,J.H., and Bedfield, J.A., Stable Numerical Integration of Conservation Eguations for Hydraulic 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 Calculations, submitted to AIAA Journal, 1978. 70 Beichenback,V.H-, and Dr e i z l e r , H., Ober den Querschmitt Sanderungen and Gitter i n Kanalen anf Stosswellen, Z e i t . Agewandte Physik, 12, p 62-7 1, 1960. 71 Von Neumann, J . , and Bichtmeyer, B.D., A Mehtod for the Numerical Calculation of Hydrodynamic Shocks, J . Appl. 36 Physics, 21, p 232-257, 1950. 72 McDonald, B.H., Hancox,8. T. , and Mathers, H.G., Numerical Solutions of the Transient Flow-Boiling Equations, Invited Paper Presented at 2nd CSNI S p e c i a l i s t s Meeting o Transient Two-phase Flow. Paris, June 1978. 73 Courant,B., Isaacson, E., and Bees, H., On the Solution of Nonlinear Hyperbolic Equations by F i n i t e Differences, Comm. Pure Appl. Math., 5, P 243-255, 1952. 74 Cu r t i s , A.C., and Beid, J.S., FOBTRAN Subroutines f o r the Solution of Sparse Sets of Linear Equations, AERE-B-6844, U.K. Atomic Energy Research Establishment, Harwell, 1971. 75 Harlow,F.H., and Amsden, A.A., Numerical Calculation of Almost Incompressible Flow, J. Comp. Physics, 3, p 80-93, 1968. 76 Hirt,C.W., and Oliphant,T.A., S0LA-PLO0P: A Non-Equilibrium D r i f t Flux Code f o r Two-Phase Flow i n Networks, Invited Paper Presented at the 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Toronto, August 1976. 77 Graf, 0., Bomstedt, P., and Merner,H., Application of the ASHE Method to Two-Phase Flow Problems, Invited Paper Presented at 2nd CSNI S p e c i a l i s t s Meeting on Transient Two-Phase 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 for the Solution of F l u i d Eguations, 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 Nucleaires 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., Depressurization Studies, Phase 2: Besults of Tests 115 and 130, UKAEA Beport SBD-B-115, 1978. 81 Ibid, Depressurization Studies, Phase 3: Besults of Tests 144 and 145, OKAEA Beport SBD-B-77, 1977. 82 Ibid, Depressurization Studies, Phase 3: Besults of Tests 142 and 143, OKAEA report SBD-B-76, 1977. 83 White,E.P., and Duffey, B.B., A Study of the Unsteady Flow and Beat Transfer in the Beflooding of Water Beactor Cores, 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., Contribution a 1'etude des debits c r i t i q u e s en ecoulement diphasique eauvapeur, Thesis, Univ. Sc 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 at the 1st CSNI S p e c i a l i s t s Meeting on Transient Two-Phase Flow, Toronto, August 1976. APPENDIX A DERIVATIONS 39 DERIVATION OF EQUATIONS USEFUL FOB THE CHABACTEBISTICS: I t i s desired to transcribe the following u n i - d i r e c t i o n a l governing eguations in t o a new set of eguations useful for application of 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 : dB/dt+(B/A) (4Au/3z) = 0, (Continuity) du/dt • (3P/tfz)/R = —F- g(dZ/dz), (Momentum) dh/dt - d(P/R)/dt • (P/BA) (3Au)/3z) = g • uF, (Energy) where q i s J/Kg/Sec, and F i s N/Kg or M/Sec*. Using the r e l a t i o n B = B(P , h ) , (because density i s f a i r l y independent of P,h in the compressed l i g u i d s t a t e ) , and using an experession for the speed of sound as 1 a* = 1./ ( 3E/3P|h + (1. /B) 2>B/dh | P ), the above set may be developed into a d i f f e r e n t set of eguations. Application of the chain rule gives: dB/dt = 3B/ap|h (dP/dt) • 3B/ah|P (dh/dt), (A.1) where from the speed of sound r e l a t i o n , (dB/iP)Jh = l./az - (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./a2 - (1./B) (»B/Sh) |P) (dP/dt) • (3R/»h)|P . dh/dt. Upon substitution for dh/dt from the energy eguation, 1 3E/aP|h = p a r t i a l of density with respect to pressure at 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 c o l l e c t i n g the c o e f f i c i e n t terms f o r dB/dt, and dP/dt, i t follows: (dB/dt) ( H (P/B2) (SB/ah) |P) = (1./a2) (dP/dt) • (dB/3h) | P [q+uF-(P/EA) 3Au/3z) 3-Now substitution for dB/dt from the continuity eguation r e s u l t s : -(B/A) (3Au/az) (1.*(P/B 2) ( a f i / a h)|P) = (1/a*) (dP/dt) «• (d£/3h)|P rLg*uF-(P/BA)) (3Au/az) . Opon cancelling the egual -(P/BA) (3Au/2z) (dB/3h)P terms from both sides, we a r r i v e at the following: — (B/A) (aAu/Sz) = (1-/a?) (dP/dt) + (3B/3h) |P[q + uF], - (Bu/A) ( d A / a z)-B3u / a z = ( 1 . / a z ) ( d P / d t ) • (3B/ah) | P Qq*uF], or f i n a l l y : or 1 41 Equation (A.4) i s the new form of the energy equation. Hext to obtain the new form of the continuity equation, we substitute for dP/dt from the energy equation: d B / d t = (SR/aP) |h dP/dt + 3fl/ah|P dh/dt = (3B/ a P)Ih 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) (3Au / a z ) ] | +B(1-/a 2--dB/aP|h) dh/dt. Col l e c t i n g c o e f f i c i e n t s of dfi/dt and dh/dt terms gives: d E / d t J.1. + (P/B) ( a f i/aP|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 substituting from continuity eguation for dfi/dt gives: -(B/A) (3Au/az) (1. • (P / E 2 ) a £ / a h | P ) = B / a 2 . dh/dt -R (3R / a P|h)[q + uF-(P/RA) (3Au/az)3# or (1./a 2) dh/dt • 1/A (cJAu/az) = afi/3P|h [q + uFJ or dh/dt • a. a u / a z = a 2 (3R/aP|h) [q+uF] -u/A aA/Sz"^. (A. 5) (A.5) i s the new form of the continuity equation. Momentum equation remains the same as before. Thus continuity, momentum, and energy equations are respectively : 42 dh/dt+a2 du/dz = a 2 dE/dP|h [g+uF]-u/A «A/dzJ = C2 (A.5) du/dt*1./B aP/Sz = -F-gdZ/dt= C l (A.6) dP/dt * Ra A du/az = -a4I>B/ah| p[g+up3+ (Bu/A) av azj = C3 (A.4) C1, C2, and C3 are c a l l e d the constants of integration. The l e f t hand side of the above hyperbolic eguations are the p a r t i a l derivatives of u, p, and h with respect to space and time ( z , t ) . These equations may be readily translated into c h a r a c t e r i s t i c d i r e c t i o n s . APPENDIX B DERIVATION OF THE OBDINABY DIFFERENTIAL EQUATIONS 4 U DERIVATION OF THE ORDINARY DIFFERENTIAL EQUATIONS: In Appendix A the following set of hyperbolic equations were derived: dh/dt • a 23u/3z = C2 (Continuity) du/dt+ (3P/az)/fl = C1 (Momentum) dP/dt • Ra* 3u/3z = C3 (Energy) From here one may pursue to f i n d the ordinary d i f f e r e n t i a l eguations along the c h a r a c t e r i s t i c s . The above can be written as: h',t • uh«,z • a*u»,z = C2 u»,t • uu«,z + P',z/R = C1 P«,t + uP«,z tRa2- 3u/3z = C3. Now transformation of above to another set of orthogonal coordinates s (along which the properties and t h e i r derivatives 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 2u*,y y',z • h*,s s»,t * uh*,s s',z • 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 f #y V,z/z = c i P'#y y*»t • uP«,y y»,z • Ba*u«,y y*,z • P« ,s s»,t • uP»,s s«,2 • fia^.uffs s*,z = C3. i h» , t = dh /a t , etc. 45 These equations can be written i n matrix form as following: !"aV#z f*,t*uy«,z y , t * u y , z (y',z)/fi o Ba^y •, z y' #t*uy',z 0 u' ,y P',y h«,y = Known. Since u>,y , P',y , h ,,y (the disturbances across the c h a r a c t e r i s t i c s ) are indeterminate, then i n order that the above i s true, the determinant of matrix must vanish. Expanding the determinant gives : (y',t + uy,z) [(ySt+uy^z)*-a 2-(y«,z) 2] = 0. From dy = 3y/3z dz+dy/dt dt =0, i t can be written: y'#t/y«,z = -dz/dt = -u, (3.1) or y'#t/y« #z = -dz/dt = -(u±a), (B.2) Eguations (B. 1) and (B-2) are the set eigenvalues or c h a r a c t e r i s t i c d i r e c t i o n s . Left eigenvectors can be obtained using these eigenvalues. Substitution from (B. 1) i n the matrix of c o e f f i c i e n t s r e s u l t s the respective l e f t eigenvector i n the following way: 46 £61 62 63] a* 1-/B = 0, Ba 2 0_ or expanding, 61 a« • Ba« 63 = 0 62 (1./H) - 0 0 = 0 , and the eigenvector can he chosen as [-B 0 1], S i a i l a r l y for dz/dt = - (y» , t) / (y 1 ,z) = u+a: [6 1 62 63] -a -a 1-/B = 0, Ba* -a or expanding the above gives: 61 a* - a6 2 • Ba2- 63 = 0, 62 /f i - a 63 = 0, and -a 51 = 0. The respective l e f t eigenvector could be Co Ba l"3« (B.3) Similarly f o r dz/dt = - (y* , t ) / ( y ' ,z) = u-a, the eigenvector i s i.0 -Ha 1], (B.4) I t i s desired now to transform the governing eguations into these dir e c t i o n s . For the dz/dt = u*a, multiply the continuity equation by *0*, the momentum equation by •Ba•, and the energy equation by * 1•, we get: Ba Iau/3t *udu/dz +(SP/dz)/S J • a P / a z *u(dP / a z ) • Ba^au/az = C3*BaC1, or ap/at+(u*a) a p / a z * aa £ a u / a t + (u+a>au/az ] = C3*aad or dF/dt+Ba du/dt = C3*BaCl. (B. 5) Si m i l a r l y 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 , we s i m i l a r l y w i l l have: -B Qh/dt*a 23u/az3I • dP/dt*Ba 2du /az = -BC2+C3, 48 -B ah/at-Buah/az-Ba^au/az+ap/at-map/az •Ba 2 a u / a z = -BC2+C3, or -B ( a h / a t + u a h / a z ) * a p / a t + u a p / a z = -BC2*C3, or f i n a l l y dP/dt-Bdh/dt = C3-BC2. (B-7) The equations (B.5)-(B.7) are set of ordinary d i f f e r e n t i a l equations (called the compatibility r e l a t i o n s ) . These equations can be written i n incremental form useful for f i n i t e difference numerical computations. APPENDIX C SOLUTION TO ISENTBOPIC PROBLEM 50 ISENTROPIC SOLUTION: Similar but simpler analysis than tor a nonisentropic problem (Appendix £) can be done to arrive at compatibility r e l a t i o n s for an isentropic problem. For the case of an isentropic problem, the governing equations are: (a£ / a t)/fl • aR/az = 0 (Continuity) ( C 1) du/dt • (dP/fiz)/R = 0 (Momentum). (C.2) The equation of state i s : ^ dP/aRJS=a2' (where a i s speed of sound and S, the entropy) . Thus ap/az = (ap/afl) (aR/az) = a 2 aR/az. Then, the equations (C .1) and (C.2) may be written as:? R» #t • u R«,z • R R',z = 0, R (u»#t • u u»,z) • a 2 fi«,z = 0. It follows from s i m i l a r analysis as previously (Appendix fi) that 1 aP/aR|S = p a r t i a l of pressure with respect to density at constant entropy 2 R»,t = d f i/at , e t c . 51 (a/B)£R' #t * (u±a) R\z2 i{,a*rt «• (u±a) a*rz] = 0, or (a/B) (dB/dt) du/dt = 0 along dz/dt = uia. (C.3) Hence (a dB)/B • u = 2r along dz/dt = u+a, and (a dfi)/B -u = 21 along dz/dt = u-a. Por polytropic gases a dB/B = 2a/(k-1), where k i s the s p e c i f i c heat r a t i o . The terms r and 1 are c a l l e d the Biemann Invariants and since they are constant along each c h a r a c t e r i s t i c , one may consider them to be the new coordinates and i n v e r t the space and time coordinates to express as z = z ( r , l ) and t = t ( r , l ) . Then from dz = (u*a) dt, dz = 3z/3l dl= (u+a) 3t/31 d l or 3z/3l = (u+a) St/31. (C.4) Si m i l a r l y 3z/3r = (u-a) 3t/3r (C.5) Then assuming that both z, and t are continuous functions i n r and 1, the eguations (C.4)-(C.5) may be combined as follow: z«,lr = (u#a)«,l t»,l • (u*a) t ' , l r z»,rl = (u-a)«,l t»,r • (u-a) t«,rl. 52 Setting the above equal, we get: 2a t»,rl • (u+a)«,l t«,l - (u-a)*,l t»,r = 0, or since Alfa • u = 2r, and Alfa-u=21, then u=r-l, and Alfa = r * l where Alfa = 2a/(k-1) ; 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 • t«,r)/(r+l) = 0, (C.6) where L * (k + 1)/2(k-1) and k = (2L+1)/(2L-1) . For a i r , L=3, and for He, L=2. The solution to the equation (C . 6 ) i s hyper-geometric functions. The hyper-geometric solution f o r the helium-gas gun i s used to check the numerical routine used f o r the present work. APPENDIX D ANALYTIC SOLUTION TO NONISENTHOPIC PBOBL 54 HONISENTBOPIC PROBLEM: I n i t i a l l y , i t was desired to Obtain exact results f o r the nonisentropic problem along the same l i n e s as isentropic problem. To do t h i s , the compatibility r e l a t i o n s were written as before (as derived i n Appendix B) : dP*Hadu = (C3+BaC1)ds along dz/dt=u+a (D-1) dP-Badu = (C3-BaC1)ds along dz/dt=u-a (B.2) dP-Bdh = (C3-BC2)ds along dz/dt=dt=u. (D.3) Integration was then performed f o r both sides of the above equations. Renaming A l f a = ^ Radu, Beta =^Bdh, 2n = P+Alfa, 21 = P-Alfa,and m = P-Beta, we f i n d that P - n+1, A l f a = n - l # and Beta - n«-l-m. From there solving f o r n, n, and 1, we get: 2n = 55 Now transforming z and t to a,n, and 1 di r e c t i o n s , gives:* z',1 = (u+a) t«,l or z* ,ln = (u+a)' ,n t»,l • (u*a) t« ,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 Z',DD = U',n t ' a • u t',mn Next taking the derivative of the above eguations with respect to a,m,l, respectively, we get: z*,lnm = (u*a) ' f n i B t» ,1* (u*a) •,n t',lBt(uta) t • , l n * (u*a) t',lnm z' ,nlm = (u-a) • ,1m 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 1 , " ! • u',1 t',»a •• u t«,anl Setting once the r i g h t hand side of the f i r s t and second of the above eguations egual and next the f i r s t and t h i r d , we get: 2a t«,anl • 2 a»," t»,ln + (u*a) »,n t ' , l a - (u-a) »,1 t*,mn + (u*a)',an t«,1 - (u-a) *.,•! t*,n = 0, i z ' ,ln = a*z/aian , etc. 56 or c o l l e c t i n g terms, i t follows: a t',mnl • a»,n t f , l m • (uia)',niii t* ,1 -u*,nl t»,» + (u*a)«,m t ' , l n - u«,l t,,mn = 0. Form there cancelling the t',onl term gives: (u-a)»,n t* ,1m • 2 u*,m t ' . l n * (u+a)«,l t»,mn- (u+a)*,an t',1 -(u-aj'^ml t # , n •-2 u« ,nl t*,m = 0. There i s no hope cf solving 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 useful to define a=a (P,T) instead of a=a(P,h), because, then, the property tables can readi l y be applied. Using chain rule r e l a t i o n s , the followings can be written: da=3h/aP|T dP*dh/dT|P dT (E.1) dv=3v/dP|T dP+av/aTJP dT- (E.2) Then i r o n (E.2) , av/apjh=av/ap|T*av/aT| p a i/apjh. (E.3) But from (E.1) 3T/aP|h=-(dh/aPlT)/(3h/aT|P) , (E.4) and also from (E.2) av/3h|P = (2v/aT 1 P) / (3h/3TJ P) . (E.5) Terms such as dB/dP)h and dB/dhjP need to be derived in a more meaningful way. He may write: 3R/3PJh = 3(1/v)/aP|h = -(av/aP|h)/v 2, and (E.6 3R/3hJP = d(1/v)/ah|P = - (Sv/dh | P)/v2: , (E-7) and since v=x vg • (1-x) vf, then combining equations (E.3), (E.4), (E.5), (E.6), and (E-7), we get: aR/ap|h = -{ • x avg/ap|T+(i-x) avf/ap|i j - i ax avg/aup • (1-x) a v f / d T J P ] [x Shg / a p j h * (1-x) a h f/3Pih] / i x a h g / a T | P+(1-x) 3hf/3T|P]j/v 2 and 3R/3h|P = - { I x a v g / S T I P+(1-x) 3vf/3T | P]/[x 3hg/dTJ P x) ahf /ar j P3]/V 2. where the qu a l i t y i n terms or void f r a c t i o n i s : x = A l / [ A 1 * {1-A1) vf/vf] , and the void f r a c t i o n i n terms of quality i s : Al = x/£ x* (1-x) vf/vg J . \ • APPENDIX F THE COMPUTER PROG BA M 61 The program IIPHAS for handling two-phase flow water proolem was written i n Fortran Language and run on Amdahl 470/V6-II computer. Double precision was used for better accuracy. A t y p i c a l run f o r th i s program took 400 seconds for complete depressurization of 120 deg. C high enthalpy water f i l l e d , 4 meter long (.032 m i n diameter), tube at saturation pressure of about 1.8 atm. I n i t i a l l y , the program was mainly debugged with help of Interactive Program (IF). System units MKS ( i . e . meter, kilogram, and second) were used throughout the analysis. That means, length i n m, time i n s, velocity i n m/s, pressure i n pa, enthalpy in J/kg, temperature i n C, and density i n kg/m3. Thermophysical properties were also used i n metric units ( i . e . , kg m/s f o r absolute v i s c o s i t y , J/kg B f o r c o e f f i c i e n t of conductivity , etc.) A three dimensional array Z(I,K1,K2) stored a l l calculated nodal values throughout the run and array values were saved i n a f i l e f o r l a t e r use. 1=1, ... ,10 indicated variables z, t, u, P, h, T, B, a , x, and void f r a c t i o n . K1 was time index f o r r e f l e c t i v e waves from closed end of the tube. K2 was the space index representing the encounter of ref l e c t e d wave from closed end with on-coming subsequent waves. 62 DESCRIPTION OF BOOTINES: HAIN Contained body of the program. Routine "FAN" contained i n t h i s routine i n i t i a l i z e d the expansion waves emanating from open end at time zero. Other routines contained i n MAIN were START and EVAL. START Read i n i t i a l values generated by Routine FAN from a f i l e . EVAL N-R i t e r a t i v e routine used to evaluate properties using SIHUL, the Nonlinear Equation Solver. IIPHAS Set up c o e f f i c i e n t matrix f o r intermittant waves. OPEN Set up c o e f f i c i e n t matrix for open end condition. This d i f f e r e d from the IIPHAS one only i n containing the choking condition f o r the open end. CLOSED Set up c o e f f i c i e n t matrix for closed end condition. SIMOL Newton-fiaphson Nonlinear Equation Solver. Having received the c o e f f i c i e n t matrix A, i t returned the solut i o n matrix X. RE Evaluated Reynolds number. TASTAR Evaluated two-phase m u l t i p l i e r for the f r i c t i o n a l force. FL Evaluated the right hand side constant of the equation of left-running c h a r a c t e r i s t i c . Evaluated the r i g h t hand side constant for the eguation of p a r t i c l e path. Evaluated the right hand side constant for the equation of riqht-running c h a r a c t e r i s t i c . Evaluated temperature from the enthalpy i n compressed l i q u i d region. 64 IMPLICIT HEAL*8(A-H,0-Z) DIMENSION A (30, 30} ,Z(10, 8,30) ,XOLD (30) ,XINC (30) C THIS PBOGBAM TO ADVANCE THE SOLUTION OF C TRAVELLING* PBESSURE MV1S F.BOH THE RUPTURE END OF HIGHLY C PRESSURIZED 4 H LONG TUBE BY THE METHOD OF CHARACTERISTICS. C THE N-RAPSON NONLINBAB EQUATION SOLVES IS USED TO OBTAIN XO C THE SOLUTIONS. THE ROUTINE EVAL HANDLES THIS. K5=1, CALL C 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) C K1 INCREMENT IS IN TIME AND K2 ALONG THE PIPE. C 50 0 CONTINUE DO 50 K1=2,ITIME C ...K3 IS A SWITCH. , K3=0 MEANS EVEN TIME ADVANCE AND K3=1 O C ... 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 8 Z (I, 1,K1) = XOLD(I) 4 WRITE {6, 101) (Z (I1,1,K1) ,11=1, 1X>) yKt,K 2, K5 101 FORM AT (6G1 4 . 7 / 4 G1 4 . 7 , 313) C .... USING IIPHAS 10UTINE TO CALCULATE^ THE INTERMITTANCE POI 9 DO 100 K2=M,N C ....CALLING THE IIPHAS-*. CALL EVAL(A,Z,XOLD,XINC,K1,K2,K3,1,1VAB,ISP,ITI,ND) DO 6 1=1,5 6 Z (I,K2,K1) =XOLD fI+5) K5= 1 WRITE(6, 101) (Z (I1,K2,Kl)iIT=1, 10) ,K1,K2,K5 100 CONTINUE K2=NN C -...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 7 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 50 CONTINUE STOP END 65 SUBROUTINE START (Z,XOLD,IVAR,ISP,ITI,ND,X,NN) IMPLICIT '.REAL*8 (A-H,0-Z) C ......INITIALIZES THE XOLD AND Z (I,K2,K1) ARRAYS. K1=1. C 1=1 IS Z, 1=2 TIME, 1=3 U, 4 P, AND 5 H,6 TEMP, 7 DENSITY,8 C SOUND SPEED,: 9 THE QUALITY, AND 10 THE VOID FRACTION. C 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 SUBROUTINE EV1L (1,2,XOLD,.XING,KT, K 2 ,K3, K5,1?HE,ISP,ITI, ND) IMPLICIT RE1L*8(A-H,0-Z) C .....K5=1 CALLS THE IIPHAS ROUTINE, =2 THE IIPHAS, =3 THE C -.CLOSED, AND =4 THE OPEN ROUTINE. C IT AUTOMATICALLY STORES THE NEW VALUES OF TEMP AND VOID C IN Z(6,K2,K1) AND Z (10, K2 ,K1) RESPECTIVELY. DIMENSION A (ND, ND) ,Z (IVAR ,ISP,ITI) , XOLD (ND) , XIHC (80) EPS 1=1-00-10 EPS 2= 1.0 DO ITMAX=6 C - . ... .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) 5 N=5 14 CONTINUE C — ..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 C . .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 C —-.ASSIGNING COMPUTED VALUES TO THE ARRAY Z. 2 Z(6,K2,K1)=TEMP3 Z (7,K2,K1) =RH03 Z (8,K2,K 1)=A3 2 (9,K2,K1) =X3 Z (10rK2,KT)=VOID3 IF(K5„NE.3) GO TO 1 GO TO 6 10 CONTINUE 6 BETUBN END 68 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) C ..THIS SOLUTION IS THE APEX Of Ai: TBIINGLE, PT 3, WHOSE BASE C POINTS 1 AND 2 ARE KNOWNf WE KNOl Z,TyU,P,H FOR THEM.) THE C AFTER THE FIRST TIME CALCULATION IN AN ITERATION ALL C CALCULATIONS CONSISTED OF FIXED VALUES FOB 1, AND 2 ABE C BYPASSED. DIMENSION A (ND,ND) >Z {IVAE,I S P 11 TI) , X OL D (ND) PR=T.D5 C ...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 C .... AFTER THE FIRST ITERATION MANY STEPS ARE BYPASSED. IF (I. HE. 1) GO TO 5 KK2=K2M C .....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= a L F A C P 2 , T E M P 2,X2 , U 2 ) FK1=FK (Pl,TEMPT,X1,U1,At,R1HP,aLF1,RH01) FL2=FLfP2,T E M P 2,X2 , U 2 , A 2 i B 2 H P , A L F 2,SH02) Z3=Z2 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 IF (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 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) C . i . I N THIS ROUTINE PT 3 IS UNKNOWN AND THE PTS 1, 2 KNOHN. C EFFOBTS HAVE BEEN HADE TO REDUCE REDUNDANCE. DIMENSION A{ND,ND) ,Z (IVAR,ISP,ITI) , XOLD (ND) X=4.D*0 C ...THE RESERVOIR PRESSURE IS TAKEN ATMOSPHERIC... PR=1.D5 DO 4 11=1,10 DO U JJ=1,11 a A (II, 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 TEHP0=TEMP1 X0=X1 U3=-(P3-P1)/{A1*BH01) P3=Z (5,K2,K1-T) H3=HG* {P3-P0)/RH01 TEMP3=TSAT(P3) C .-.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 C ..^.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,VD0rC,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 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 A (5,3 A (5,6 A (5,7 A{5,8 A(5,1 A (6,2 A (6, 3 A (6,4 A (6, 5 A (6,7 A (6, 8 A (6,9 A (6,1 A (7, 6 A (7,7 A (7, 8 A (7,1 A(8,7 A (8,8 A (8,9 A (8,1 A (9,6 A (9,1 A (-10, RETURN END = U30 = (T0-T3) /2. -1.D+0 =-U30 = (T0-T3) /2, )=-F5 =FH30 =BETA63 =-1. D+0 =BHO30 =-FH30 =BETA68 = 1. D+0 D+0 D+0 A(6,10)=-RHO30 )=-F6 = 1. D + 0 =-OA31 = (T1-T3) /2. )=-F7 =-FK31 =BETA88 =1.D+0 )=-F8 = 1-D+0 )=-F9 1)=-F10 D + 0 76 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) C ...PTS 0,2,3 ARE USED IN THIS ROUTINE- PROPERTIES OF PT C 3 ARE UNKNOWN. WHEN FIRST TIME USED IN AN INTE1ATI0N, C THE PROPERTY VALUES OF PT 0 ARE ASSIGNED TO POINT 3. C AFTER THAT THE ITERATED VALUES OF XOLD RESPECTIVELY ARE USE C UNTIL IT CONVERGES, PTS 0 AND 1 ARE ONLY ASSIGNED AT C FIRST ITERATION-DIMENSION 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 IFfl.NE. 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 C 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*13HP-1 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) C .....THIS ROUTINE COMPUTES THE SOLUTION OF THE N X N+1 C MATRIX A. AUGMENTED VALUES ARE IN N+T COLUMN. C 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 C ... IS B LARGER THAN 50. IF (N.LE.50)GO TO 5 WRITE (6,200) SIMUL =Q.D+0 RETURN C ......BEGIN ELIMINATION PROCEDURE....... 5 DETER =1.D+0 DO 18 K = 1,N KM1 = K - 1 C .......SEARCH FOR THE PIVOT ELEMENT PIVOT= 0.D+0 DO 11 1=1, N DO 11 J=1,N C ..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 8 CONTINUE 9 IF (DABS (A (I, J) ) . LE. DABS (PIVOT) ) GO TO 11 PIVOT= A (I ,3) IEOW(K)=I JCQL(K) =J 11 CONTINUE C .......INSURE THAT SELECTED PIVOT IS LARGER THAN EPS . IF(DABS(PIVOT).GT.EPS) GO TO 13 SIMUL=0.D+0 RETURN C C .UPDATE THE DETERMINANT VALUE 13 IRG¥K=IRO¥ (K) JCOLK=JCOL(K) DETER=DETER*PIVOT C C ......NORMALIZE PIVOT ROW ELEMENTS...... DO 14 J=1,MAX 14 A(IROWK,J) =A (IR01K, J) /PIVOT C C .....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 17 IP (J. HE. JCOLK) A (I,J)=A (I> J)-AIJCK*A (IHOSK, J) 18 CONTINUE C C .......ORBER SOLUTION VALUES (IF ANY)AND CREATE JOED AREA DO 20 1=1,N IROfI=IROI (I) JCOL 1= JCOL (I) JORB(I10HI)=JCQLI 20 IF(INBICiGE.O) X{JCOLI)=A(IROSI,HAX) C C ADJUST SIGN OF DETERMINANT. ... . . INTCH=0 NH1=N-1 DO 22 I=1,NM1 IPt=I*1 DO 22 J=IP1 ,N IF (JORB (J) .GE. JOED (I) ) GO TO 22 JTEMP=JORD (J) JOES (J) = JOBD(I) JOED fl) =JTEMP INTCH=INTCH+1 22 CONTINUE IF (INTCH/2*2-NE. INTCH) DETEB=-DETER C • • C .•.'.•v-.'»-i-JIF -INDIC IS POSITIVE aET01S',SlT& RESULTS.:.,.^ *, , IF(INBICLE-O) GO TO 26 SIMUL= DETER RETURN C C IF INDIC IS NEGATIVE OS ZERO, UNSCRAMBLE THE INVE C FIRST BY BOfSi.... 26 DO 28 J=1,N DO 27 1=1,N IROWI=I10W(I) JCOLI=JCOL(I) 27 Y (JGOLI) =A (IROWI, J) DO 28: 1=1, N 28 A(I,J)=Y(I) C ......... ..THEN BY COLUMNS.-- > . - . . DO 30 1=1,N DO 29 J=1,N I1OWJ=IR0U{J) JCOLJ=JCOl(J) 29 Y (IS01J) =A (ly JCOLJ) DO 30 J=1,N 30 A(I,J)=Y(J) C C > i i . .:.EETUR8 > FOR INDIC 'NEGATIVE OE 2ERO. . . - . i . SIMUL=BETER RETURN C . . v . . wviiv. ...FORMAT* FOR- OUTPUT STATEMENT.:. . , 200 FORMAT(10HON TOO BIG) C END 81 FUNCTION RE (P,T, X,U) C -...COMPOTES THE REYNOLD*S NUMBER. IMPLICIT REAL*8{&-ff,0-Z) IIIL*8 MU,MUF,MOS&T 6=9^8094D+O D=.032D+0 C -.l.E=0.:D'+0 C • IF.C>O.Lfe.;02D+O-):RETIIRN-C M0=M0F (P,T) C IF {X.LE.0.1D+0 )GO TO 1 fiflO='T. D+0 / V (P,T#X)-MU=HOSAT(P,T) C GO TO 2 C 1 RHG=1.D + 0 /VFfT) RE=RHO*0*D/N0 RETURN END 82 FUNCTION TASTAR( P,T,X,U) C ..-.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 C TASTAR=1.D+0 C 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 C ... 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 C FL=0.D+0 C 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 8H FUNCTION FH(P,T,X, U,AA,RH,AL,80) C ...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 C FM=0.,D + 0 C 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 , ft A ,.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 FUNCTION T(E) C ...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) 90 FOBMAT(2X, * NEGATIVE ENTHALPY*) STOP 2 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 10 2 CONTINUE 100 XN= fXLO + XBG)/2.D+0 FH=FUN (XN) IF(FLO*FN) 103,104,105 101 SOLN=XHO GO TO 106 103 X!G=XN FBO=FN IF(DABS(XBO—XLO); .GT. EPS) GO TO 100 SOLN=XSO GO TO 106 10 4 SOLN=XN GO TO 106 105 XLO=XN FLO=FN IF(DABS(XBO-XLO).GT.EPS)GO TO 100 SOLN=XLO 106 T=SQLN BETUBN END C ...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) C 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 FUNCTION PSfiT(T) IMPLICIT HE At* 8 (&-HrO-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 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 SOB BOOTING SS (P,T,X,VD,A,RPH,RHP, V) C ..4 SOUND SPEED=1/(BPH+RHP/R)**2. VOID FRACTION IS ALSO FOON IMPLICIT BEAL*8 (A-H,0-Z) C .,.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 C IF (X.LE.0.0B+0 )GO TO 1 VGPT=(VG(PDELP,T) - VG(P rT))/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 f 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) C GO TO 2 C 1 V=VF(T) C VD=0iDO C VFF=-(1-D+O /V**2) C RPH=0.D+0 C BHP=VFF*(VFTP/HFTP) Y=DABS (RPH+V*RHP) A=1.D*0 /DSQRT (Y) RETURN END 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 FUNCTION HG (P,T) C ..... SUPERHEftT SPECIFIC ENTHALPY... IMPLICIT REAL*8 (A-H,0-Z) C .. ..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 FUNCTION VG (P,T) C ..SUPERHEAT SPECIFIC VOLUME.-.. C WE FEED PUS. aND DEG C. VG RETURNS IN M**3/KG. IMPLICIT REAL*8(a-HfO-Z) P1=P*9.869D-6 a5=.2088543D+5 T1=T+273.16D+0 S=1/T1 a=8G870*S**2 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 FUNCTION MUG (P,T) C VAPOR VISCOSITY. IMPLICIT REAL*8 f A~H,0-Z) C ...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 FUNCTION HUSAT {P,T) C .v;.iSATURATION VISCOSITY... IMPLICIT HEAL*8(A-H,0-Z) C . . . m S . AND DEG. C ARE FED. HO SAT1 SET DEN IN KG/H SEC. 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 97 FUNCTION' QUTY(P,T,H) C ....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 98 FUNCTION HC(P,T,X,U) IH.PLICIT REAL*8 (A-H,0-Z) C .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 9 9 FUNCTION CK(P,T) C ....EVALUATES THE COEFFICIENT IMPLICIT REAL*8 (A~H,0-2;) C .^ WE FEED PAS. AND DEC. C.THE 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 OF CONDUCTIVITY. FUNC GIVES J/SEC. M. DEG C. 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 FO NCTIG K Q { H C , 7 , T) IMPLICIT REAL*8 (A-H,0-Z) C ... J/H**2 S CV M**3/KG, AN0 C ABE FED. Q, HEAT FLUX C 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 103 FUNCTION V (P,T,X) C EVALUATES THE SPECIFIC VOLUME.„*... IMPLICIT REAL*8(A-H#0-Z) VFF=VF (T) C IF{X.LE.. 1D+0 )RETURN V=X*?G <P,T) + (1. D+0 -X)*VFF RETURN END 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 F i g u r e 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 V o i d O r i e n t e d V e r t i c a l l y f o r Upflow Blow-down 106 0 J , , , , , I I I I T 0 2 4 6 8 10 Time, ms . F i g u r e d 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 b u r s t 107 F i g u r e 3 C h a r a c t e r i s t i c s Diagram i n Z-T P l a n e 108 R-WAVE  A z 3 f ( u + a ) 3 1 A t 3 1 A P 3 1 + ( a R ) 3 1 A u 3 1 = K 3 1 A t 3 1 P-WAVE A Z 3 0 = U 3 0 A t 3 0 A P 3 0 " R 3 0 A h 3 0 = M 3 0 A t 3 0 L-WAVE fcz32=(u"a)32At32 A P 3 2 " ( a R ) 3 2 A u 3 2 = L 3 2 A t 3 2 h t A Z 0 1 A t 0 1 A u 0 1 A P 0 1 A h Q 1 A z 2 1 A t 2 1 A u 2 1 A P 2 1 " A h 2 j z F i g u r e 4 General Two-Phase a l g o r i t h m to advance s o l u t i o n a l o n g dz/dt=u±a c h a r a c t e r i s t i c s 109 P-WAVF A P 3 0 " R 3 0 A h 3 0 = M 3 0 A t 3 0 C l o s e d end a) C l o s e d End L-WAVE A t 3 2 = z 2 / ( a - u ) 3 2  A P 3 2 " ( a R ) 3 2 u 2 = L 3 2 A t 3 2 l R-WAVE A Z 3 f ( u + 3 ) 3 1 A t 3 1 A P 3 1 + ( a R ) 3 1 A U 3 1 = K 3 1 A t 3 1 P a i f u 3 < a 3 o r U 3 = a 3 P-WAVE A z 3 0 = U 3 0 A t 3 0 A P 3 0 _ R 3 0 A h 3 0 = M 3 0 A t 3 0 A 2 0 1 A t 0 1 A P 0 1 A u Q 1 A h o i A t 2 1 A u 2 1 A h 2 1 Open end b) Open End F i g u r e 5 Wave i n t e r a c t i o n s a t the b o u n d a r i e s F i g u r e 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 P r o j e c t i l e '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 Zero Mass B u l l e t F i g u r e 8 E x p a n s i o n Fans f o r Helium Gun Q_ * * i- 2 . - I • 4 * 1 1 1 1 T T 1— 0 20 40 60 80 100 120 140 T i m e , mS F i g u r e 9 P r e s s u r e a t C l o s e d E n d ( P 0 = 3 . 5 MPa ) ^ U p f l o w • D o w n f l o w 0.0 0.5 1.0 V o i d F r a c t i o n , a F i g u r e 10 Choking 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 Fan, 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 — I 1 1 1 1 1 1— 20 40 60 80 100 120 140 T i m e , mS F i g u r e 11 V o i d F r a c t i o n a t C l o s e d E n d (po=-12 MPa ) U p f 1 ow Downf1ow 1.0. c o •r— +J o «o i t 0.5. "O •I— o O.Oi 1 — 20 —I 1 1 1 1— 60 80 100 120 140 T i m e , mS 40 160 F i g u r e 12 V o i d F r a c t i o n a t C l o s e d E n d ( P 0 = . 1 8 MPa ) • U p f l o w • D o w n f l o w en • 0 o o o o A o A ~1~ 20 40 ~1 6 0 80 1 0 T i m e , mS F i g u r e 13 P r e s s u r e a t t h e C l o s e d End • U p f l o w O • D o w n f l o w A P0=.18 MPa O A O A P0=.12 MPa A C/1 CD 01 +-> ca o F i g u r e 14 P0=.12MPa ' • U p f l o w A D o w n f l o w T i m e , mS M a s s F l o w R a t e a t t h e O p e n E n d P0= 18MPa — 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) F i g u r e 15 Computer Flow C h a r t SUBROUTINE EVAL C a l l OPEN Deter=SIMUL "„ J = 1 » 1 0 XN(J)> XOLD(J)? 1 S t o r e r e t u r n e d r e s u l t s i n Z- A r r a j Return F i g u r e 15 Computer Flow C h a r t SUBROUTINE IIPHAS Se t A. .=0 A s s i g n v a l u e s t o p o i n t s 1 and 2 I n i t i a l i z e p o i n t s 0 and 3 E v a l u a t e Az e t c . 21 I n i t i a l i z e the s o l u t i o n , X ( j ) 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 Computef , a, x, A l , e t c . f o r p o i n t s 0 and 3 < • 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 a s s i g n v a l u e s Compute augmented m a t r i x c o m p r i z e d o f -F's Set up A.. i = l,10 a n d 1 j = l , l l Return • End F i g u r e 15 Computer Flow C h a r t 122 T a b l e 1 Parameters f o r the Gas Dynamic Problems I n i t i a l C o n d i t i o n s : 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 e i f u<a o t h e r w i s e u=a ( c h o k i n g c o n d i t i o n ) P a r a m e t e r s: a=804.2 mm2 dA/dz=0.0 m dZ/dz=0.0 D =32.0 mm e -1 k/l=0.0 m 1 q=0.0 W/m C o n s t i t u a t i v e Laws: F=(4f + k / l ) u 2 = 2 f u 2 / D D 2 e e _6 where f = f r i c t i o n f a c t o r = 1 0 Z T h e o . ____ • Z.Numer. Sheo. b u r n e r . V e l o c i t y S o u n d S P e e d 0.0 0.0 ~ 07298?"540'E- 0 2:~6T2.99l"776E-020. 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 67 5 E-0 2 0 . 6 8 9 1 5 26=+ 0 2 0 . 1 313 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 . 1 2 9 1 0 8 8 E + 0 4 • C . 4 6 9 5 3 1 6 F + 0 0 0 . 4 7 0 5 0 8 9 E + 0 0 0 . 3 3 2 2 3 9 9 E - C 2 0 . 3 3 2 6 0 6 7 E - 0 2 0 . 2 0 6 7 4 5 3 E + C 3 C . 1 2 6 8 2 2 8 E + 0 4 7j77J54^4TWnni70T7J^^ 0.85-62123E + 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 + 04 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 94? 5£ + G 1~0T2"TT5T4 5 P + 0 1 0 . 4 3 5 1 2 ^ e T ^ C 2 _ ^ 4 " 3 T O 7 3 T r 0 7 _ ^ W r " 5 " c T E T 0 T ~ C V r i O W c ^ + l ? 4 -0 7 7 3 6 8 T 3 4 E - 2 3 0 . 0 0 . 3 2 0 7 3 3 5 E - G 2 0 . 3 2 1 1 2 4 2 E - 0 2 - 0 . 3 3 8 8 1 2 2 E - 2 G 0 . 12 9 1 1 82E+04 0. 1 5 6 7 6 8 2 E + 00 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 1678E+G2 0 . 1 2 6 8 3 2 8 E + 04 0 . 3 2 8 1 6 1 6 E + C 0 0 . 3 2 8 6 3 4 5 E + 0 C C . 3 4 5 1 3 7 2 E - 0 2 G . 3 4 5 6 8 8 0 E - 0 2 0 , 1 3 7 8 3 5 1 E + G 3 0 . 1 2 4 5 4 7 2 E + 0 4 0 T 5 l T 3 2 ~ 7 T r r + W l ^ T T ^ 0 . 7 2 C 3 5 4 0 E + C 0 0 . 7 2 1 7 6 9 2 E + 0 0 0 . 3 7 2 4 3 9 8 E - 0 2 0 . 3 7 3 1 7 9 4 E - 0 2 0 . 2 7 5 6 7 5 3 E + 0 3 0 . 1 1 9 9 7 4 7 E + 0 4 C o9443O05E+G0 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 + 03 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 TC* + C 3 C . 1 T 5 4 C C 2E + 04 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-02 0 . 4 8 2 4 3 6 8 5 + 03 0.1 1311 26E + 04 O c 1 7 5 0 3 0 4 E + 0 1 0 . 1 7 5 5 5 7 9 E + 0 1 0 . 4 3 7 7 6 5 5 E - C 2 0 . 4 3 £ 9 8 4 6 E - 0 2 0 . 5 5 1 3 5 5 0 E + 0 3 0 . 1 1 0 8 2 4 4 E + 0 4 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 r 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 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 + 04 0 . 1 7 2 1 3 0 3 E + 00 0. 1 7 2 4 5 9 2 E + 00 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 + 00 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 + 00 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 + 04 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 51E +03 0.1 1 3 1 2 6 1E + 04 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 84E + 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 + 04 0 . 1 9 5 5 2 3 7 E + 01. 0. 1 9 6 2 8 6 4 E + 01 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 T a b l e 3 Parameters f o r P r e s s u r i z e d Water Problem 124 I n i t i a l C o n d i t i o n s : 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 c o n d i t i o n ) P arameters: A=804.2 mm2 dA/dz=0.0 m dZ/dz=0.0 D =32.0 mm, k?l=0.0 m"1 q=0.0 W/m C o n s t i t u a t i v e E q u a t i o n s : F=(4fx*+k.)u 2 where f = f r i c t i o n f actor=0.046 R e~* 2 f o r R e>2000 D e 1 2 =64/R e f o r R e<2000 x*=l + x ( 0 - l ) { l + 3 . 5 7 e - ° - 0 0 8 8 4 Y ( 1 - e " 4 - 9 6 ^ - ^ ) } 1/3 where y = G / { g p f ( p f - p g ) y f } 0 2 and 8= ( v g / u f ) ' P f / p g here G i s the mass v e l o c i t y (kg/m s) 

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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0094637/manifest

Comment

Related Items