TH E U N I V E R S I T Y OF M I C H I G A N COLLEGE OF ENGINEERING Department of Aerospace Engineering Technical Report A NUMERICAL STUDY OF THE EFFECTS OF INTENSE TRANSIENT HEATING IN A LAVAL NOZZLE James M. Noyak DRDA Project 011116 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GI-36248X WASHINGTON, D.C. administered through: DIVISION OF RESEARCH DEVELOPMENT AND ADMINISTRATION ANN ARBOR February 1974

ACKNOWLEDGMENTS The author takes this opportunity to express sincere appreciation to Professor Richard L. Phillips for his kind patience, understanding, and able direction during the last three years. His friendship and support have made this dissertation research particularly memorable. I would like also to thank the other members of my committee, Professors Thomas C. Adamson, Arthur F. Messiter, Martin Sichel, and James O. Wilkes for their interest and support. The support of all the people associated with the Gasdynamics Laboratory and Aerospace Engineering Department of The University of Michigan is hereby gratefully acknowledged. I must single out Mr. John Van Roekel for the countless times he helped me in understanding the idiosyncrasies of the MTS computing system. Special thanks go to Mrs. Eunice Sherbrook and Mrs. Christine Long for their capable typing of this manuscript. Mrs. Susan Tack also assisted in the typing of earlier drafts. The financial assistance provided me through an NDEA Title IV Fellowship and National Science Foundation research grant GK-14757 is sincerely appreciated. And finally, I acknowledge the immeasurable support of my iii

wife, Carmen, who during the last six years has handled the toughest assignment of all (and passed with honors)-being mother to our three children. iv

TABLE OF CONTENTS Page DEDICATION A CKNOWLEDGMENTS LIST OF ILLUSTRATIONS NOMENCLATURE 1. INTRODUCTION 1.1 The Gas Blast Circuit Breaker 1.2 Numerical Solution of the Conservation Equations of Gasdynamics 1. 3 Overview of Present Research 2. QUASI-ONE-DIMENSIONAL MODEL 2. 1 Introductory Remarks 2.2 The Quasi-One-Dimensional Governing Equations 2. 3 Numerical Technique 2. 4 Quasi-One-Dimensional Results and Discussion ii iii vi xiii 1 2 9 20 22 22 25 29 58 3. AXISYMMETRIC MODEL 3. 1 Introductory Remarks 3.2 The Axisymmetric Governing Equations 3. 3 Numerical Technique 3. 4 Axisymmetric Results and Discussion 178 178 182 183 198 4. CONCLUSIONS APPENDIX REFERENCES 222 225 243 v

LIST OF ILLUSTRATIONS Figure Page 1-1 Arc-Nozzle Configuration of Interest in Study 4 1-2 Electrical Conductivity of SF6 Versus Temperature 5 2-1 Illustration of CFL Stability Criterion 51 2-2 Alternate Interpretation of CFL Criterion 54 2-3 Transient Temperature Profiles for Isentropic Nozzle Flow 62 2-4 Transient Massflow Profiles for Isentropic Nozzle Flow 65 2-5 Temporal Density Variation for Isentropic Nozzle Flow 67 2-6 Temporal Velocity Variations for Isentropic Nozzle Flow 68 2-7 Temporal Temperature Variations for Isentropic Nozzle Flow 69 2-8 Spatial and Temporal Temperature Variations for Isentropic Nozzle Flow 71 2-9 Spatial and Temporal Massflow Variation for Isentropic Nozzle Flow 72 2-10 Spatial and Temporal Density Variation for Isentropic Nozzle Flow 73 2-11 Spatial and Temporal Velocity Variation for Isentropic Nozzle Flow 74 2-12 Spatial and Temporal Pressure Variation for Isentropic Nozzle Flow 75 vi

Figure Page 2-13 Duct Configuration for Rayleigh Flow 77 2-14 Sketch of Rayleigh Line on T-s Diagram 81 2-15 Static Temperature Distribution for Increasing Values of Q 84 2-16 Mach Number Distribution for Increasing Values of Q 85 2-17 Transient Mach Number Profiles for Q = 0. 6 88 2-18 Temporal Variations of Duct Exit Temperature 89 2-19 Density Distribution for Q = 0. 5535 90 2-20 Velocity Distribution for Q = 0. 5535 91 2-21 Static Temperature Distribution for Q = 0. 5535 92 2-22 Mach Number Distribution for Q = 0. 5535 93 2-23 Assumed Electrical Conductivity-Static Temperature Dependence 99 2-24 Initial Transient Massflow Profiles for I = 100 amps (DC) 101 max 2-25 Initial Transient Pressure Profiles for I = 100 amps (DC) 103 max 2-26 Overall Massflow Variations in I = 100 amps (DC) max 105 2-27 Overall Pressure Variation for I = 100 amps (DC) max 106 2-28 Temperature vs. Time at Several Axial Stations for I = 100 amps (DC) 108 max 2-29 Static Temperature for I = 500 amps (DC) 109 2-30 Steady-State Heat Addition Rate Distribution for I = 100 amps 111 max vii

Figure Page 2-31 Steady-State Hot (Ima = 100 amps) and Cold Flow Static Temperature Distributions 112 2-32 Steady-State Total Temperature Distribution for I = 100 amps 114 max 2-33 Steady-State Hot (Imax = 100 amps) and Cold Flow Mach Number Distributions 116 2-34 Steady-State Hot (Imax = 100 amps) and Cold Flow Pressure Distributions 117 2-35 Steady-State Hot (Ima = 100 amps) and Cold Flow Density Distributions 118 2-36 Steady-State Potential Difference (in volts) for I = 100 amps 119 max 2-37 Transient Temperature Response at x = 0. 2 for Different Values of J and I = 100 amps (DC) max 121 2-38 Transient Temperature Response at x = 0. 8 for Different Values of J and I = 100 amps max (DC) max 122 2-39 Three Nozzle Geometries for Comparison Study 128 2-40 Massflow Versus Current for Three Nozzle Geometries 131 2-41 Steady-State Heat Addition Rate Profiles for Increasing Current-Conical Nozzle 133 2-42 Steady-State Static Temperature Profiles for Increasing Current-Conical Nozzle 135 2-43 Steady-State Total Temperature Profiles for Increasing Current-Conical Nozzle 136 2-44 Steady-State Potential Difference Profiles for Increasing Current-Conical Nozzle 137 viii.. lrlll

Figure Page 2-45 Steady-State Mach Number Profiles for Increasing Current-Conical Nozzle 138 2-46 Steady-State Mach Number Profiles for Increasing Current-Hyperbolic Nozzle 139 2-47 Heat Addition Rate Profiles for Three Nozzle Geometries (I = 500 amps) 141 max 2-48 Potential Difference Profiles for Three Nozzle Geometries (I = 500 amps) 143 vmax 2-49 Static Temperature Profiles for Three Nozzle Geometries (I = 500 amps) 144 max 2-50 Total Temperature Profiles for Three Nozzle Geometries (I = 500 amps) 145 max.1 2-51 Mach Number Profiles for Three Nozzle Geometries (I = 500 amps) 146 max 2-52 Static Temperature Surface for Initial Moments of Free Decay in the Conical Nozzle 148 2-53 Static Temperature Surface for Initial Moments of Free Decay in the Hyperbolic Nozzle 149 2-54 Static Temperature Surface for Initial Moments of the Decay in the Logarithmic Nozzle 150 2-55 Initial Static Temperature Decay at Four Axial Locations in the Conical Nozzle 152 2-56 Initial Static Temperature Decay at Four Axial Locations in the Logarithmic Nozzle 153 2-57 Relative Rate of Decrease (Dimensionless) of Static Temperature for Three Nozzle Geometries 155 ix

Figure Page 2-58 Full Cycle Density Response in the Conical Nozzle (j = 6) 159 2-59 Full Cycle Velocity Response in the Conical Nozzle (j = 6) 160 2-60 Full Cycle Static Temperature Response in the Conical Nozzle (j = 6) 161 2-61 Full Cycle Rate of Energy Addition Response in the Conical Nozzle (j = 6) 162 2-62 Full Cycle Density Response in the Conical Nozzle (j = 13) 163 2-63 Full Cycle Velocity Response in the Conical Nozzle (j = 13) 164 2-64 Full Cycle Static Temperature Response in the Conical Nozzle (j = 13) 165 2-65 Full Cycle Rate of Energy Addition Response in the Conical Nozzle (j = 13) 166 2-66 Full Cycle Density Response in the Conical Nozzle (j = 31) 167 2-67 Full Cycle Velocity Response in the Conical Nozzle (j = 31) 168 2-68 Full Cycle Static Temperature Response in the Conical Nozzle (j = 31) 169 2-69 Full Cycle Rate of Energy Addition Response in the Conical Nozzle (j = 31) 170 2-70 Full Cycle Massflow Response in the Conical Nozzle (Throat) 174 2-71 Full Cycle Massflow Response in the Hyperbolic Nozzle (Throat) 175 x

Figure Page 2-72 Full Cycle Massflow Response in the Logarithmic Nozzle (Throat) 176 3-1 Illustration of Coordinate Transformation 187 3-2 Mach Number Isoline Evolution-Effect of Wall Boundary Condition 201 3-3 Mach Number Isolines in a 45~-15 Conical Nozzle with Circular-Arc Throat 204 3-4 Axial Mach Number Distribution in a 45~-15 Conical Nozzle 206 3-5 Axial Static Pressure Distributions in a 45~15~ Conical Nozzle 207 3-6 Axial Density Distributions in a 45~ 15~ Conical Nozzle 208 3-7 Axial Static Temperature Distributions in a 45~-15~ Conical Nozzle 209 3-8 Temporal Static Temperature Response in Mid-Subsonic Nozzle Region 212 3-9 Temporal Static Temperature Response in Transonic Nozzle Region 213 3-10 Temporal Static Temperature Response in Mid-Supersonic Nozzle Region 214 3-11 Mean Total Temperature Axial Distributions for Non-Adiabatic Flow Case 216 3-12 Mean Density Axial Distributions for NonAdiabatic Flow Case 218 3-13 Mean Static Temperature Axial Distributions for Non-Adiabatic Flow Case 219 3-14 Mach Number Isolines for Both Adiabatic and Non-Adiabatic Flow 220 xi

NOMENCLATURE a speed of sound, m/sec 2 A nozzle cross-sectional area, m; also, matrix of coefficients in conservation equations 2 A arc cross-sectional area, m arc A percent axial rate of change of A b coefficient in Eq. (1. 11); also, semi-conjugate axis c specific heat at constant volume v c specific heat at constant pressure C coefficient matrix CSF safety factor for time step control e specific internal energy E electric field strength, volts/ m f, g,h general dependent variables; also, column vectors representing flow variables F vector representing finite-difference solution values F defined by Eq. (1.8) j current density, amps/m2 i, j spatial finite difference indices I current, amps; also, total number of radial difference mesh segments J total number of axial difference mesh segments L nozzle length, m rn massflow, kg/sec xii

M n p p q Q Q r r c R R s S t T T cut -off u V V V Mach number temporal finite difference index operator representing general overall difference scheme 2 static pressure, newtons/mn or atm defined by Eq. (3. 2) order of truncation error net rate of heat addition per unit volume, watts/m heat addition per unit mass between two locations, Joules/kg radial coordinate radius of curvature, m gas constant; also nozzle wall radius, m percent axial rate of change of R specific entropy column vector of "source-like" terms time, sec static temperature, K static temperature at which gas electrical conductivity first becomes significant, OK axial velocity component, m/ sec radial velocity component, m/ sec velocity vector, m/sec voltage defined by Eq. (3. 16) xiii

w velocity component in z direction, m/sec x axial coordinate x percent axial distance from nozzle throat y Cartesian coordinate direction z Cartesian coordinate direction Z stretching function a conical nozzle half-angle; also, inlet flow angle f defined by Eq. (3. 14) E truncation error y ratio of specific heats Amd mass flow decrement dec Ar difference mesh distance in r direction At difference mesh distance in t direction'Atlag time lag of flow response to current input Ax difference mesh distance in x direction X At/ Ax (fb At/ Ar v frequency of AC current input p density, kg/m3 a electrical conductivity, mho/m T time constant of static temperature decay xiv

Subscripts crit E hf isentr I max. s. ref t tan thr r,t,x,y "critical" conditions (M = 1) nozzle or duct exit "heat front" location isentropic conditions nozzle or duct inlet maximum "imaginary" node; also, initial conditions steady state value reference quantity "total" conditions tangency locations in conical nozzle wall contour nozzle throat location denotes partial differentiation with respect to respective variable Superscripts * () A ( ) denotes non-dimensional quantity "predicted values"; also, mean values denotes state at which M = 1 in Rayleigh process XV

1. INTRODUCTION The present study has as its primary goal the numerical solution of the non-linear partial differential conservation equations that govern the unsteady, inviscid, mixed subsonic-supersonic flow within a Laval nozzle. The problem of interest involves an intense, time varying and spatially varying, bulk heat addition within the nozzle, which acts as a transient "driving force" on the internal gasflow. In what follows we consider both a quasi-onedimensional model and an axisymmetric model based on Euler's equations. Each of these initial, boundary-value problems is discretized by means of a two-step finite-difference technique and subsequently solved on a digital computer. The ensuing analysis places particular emphasis on the dynamic convective heat transfer processes within the nozzle and on the role that nozzle geometry plays in altering these processes. This work was initially motivated by a desire to better understand some of the physical phenomena underlying the operation of the so-called gas blast circuit breaker, which finds wide use in the power industry as an effective interrupting device for high power AC electric systems. It is believed that convective heat transfer plays an important role in arc extinction, and nozzle design criteria for this type of breaker have long been sought. 1

2 Since the "art" of numerical analysis of non-linear partial differential equations is in its infancy (relative, say, to the progress that has been made in the mathematical solution of linear equations), a concurrent goal of this study is to evaluate the feasibility of using modern numerical techniques to effectively model the circuit breaker problem. In particular, we draw some conclusions regarding more realistic multi-dimensional models, including those that account for real gas effects, such as thermal and electrical conductivities and thermal radiation, in the working fluid. Before describing the numerical models and associated analyses, we briefly discuss the physical problem of interest and provide some background in the area of numerical analysis of partial differential equations. 1. 1 The Gas Blast Circuit Breaker 1 2 Lee and Rieder both review the state of circuit breaker technology, discuss the various types of breakers in use, and point out areas in which research would be beneficial to a better understanding of the physical phenomena involved in their performance. Gas blast breakers find wide use in ultra-high voltage AC systems with voltage magnitudes up to 750 kV and current values as high as 50, 000 amperes. As their name implies, they depend on the "blast" of some gas from a high pressure reservoir to rapidly cool and then

3 extinguish an electric arc that completes the external circuit to be interrupted. Although there exist several mechanical configurations for such breakers, the arc-nozzle combination shown in Fig. 1-1 illustrates the principles involved in axial-flow breakers and forms the basis of the ensuing analyses. After the separation of the two electrodes and the subsequent formation of the arc along the nozzle centerline, a gas enters the nozzle from an upstream reservoir and exhausts to a low pressure receiver. The pressure drop across the nozzle is so great that the back (receiver) pressure offers no resistance to the interior gas flow, which freely expands throughout the converging and diverging nozzle regions. Gas blast breakers use either air or sulfur hexafluoride (SF6) as their working gas. The latter substance has certain special 2 3 properties that are conducive to arc extinction. Figure 1-2 shows the dependence of the electrical conductivity, a, of SF6 on static temperature and pressure, according to the recent analytical 4 results of Liebermann. The abrupt increase of a with temperature at about 4000 K, a behavior found in most other gases, allows for a clear definition of an electric arc "boundary", as illustrated in Fig. 1-1. As the AC current input to the arc undergoes a complete temporal cycle, this boundary defines an arc crosssection that, on the one hand, may fill the containing nozzle during peak current transition or, conversely, may shrink to a narrow

NOZZLE CONTOUR f THROAT GAS FLOW N RESERVOIR ELECTRODE / t INLET PLANE RECEIVER (P 0) r L ELECTRODE EXIT PLANE Figure 1-1. Arc-Nozzle Configuration of Interest in Study

5 100. 80.'E I 0 E 60. p 16 atm (/4 20.x 103 b 40. 20. 0 0 4. 8. Ia Temperature (~K) 16. Figure 1-2. Electrical Conductivity of SF6 Versus Temperature.

6 filament along the nozzle centerline as the current magnitude becomes zero. This "current zero" period, during which the power input to the arc instantaneously vanishes, provides the optimum conditions for arc extinction. To convert the electrically conducting plasma into a non-sconducting gas that will withstand rapidly increasing circuit recovery voltages without thermal breakdown (or "re-ignition"), an enormous amount of energy must be removed from the arc region within a few microseconds. The energy transfer mechanisms available to do this include laminar and turbulent thermal conduction, convection, and radiation. An understanding of these mechanisms is, of course, essential for circuit breaker design. And yet, their nonlinear and highly coupled nature render the problem almost entirely insoluble. The possibility of significant electrode and gaseous nonequilibrium effects further complicate the task of investigators. For these reasons, previous analytical and experimental studies have attempted to isolate only one of the important transfer processes. The majority of previous analytical work in this area, beginning 5 with that of Cassie, deals with the idealized, cylindrically symmetric, "asymptotic arc column" in which radial conduction, and possibly radiation, are the dominant heat transfer mechanisms. One imagines that the arc is confined within a long, well-cooled, constant-area tube so that far away from the electrodes all axial gradients become negligible. Through various mathematical and physical simplifications,

7 a generally comprehensive understanding of dynamic electric arc behavior has been obtained and the results have given some insight into the role of thermal conduction near current zero passage in gas 6 blast breakers. Anderson presents an excellent summary of this work. In the few cases for which radial convection effects were included in the analyses, either an assumed velocity behavior (e.g. a linear dependence on the radial coordinate) or numerical solution techniques 7 were necessary. Phillips recently demonstrated that the momentum conservation equation must be solved to properly account for radial flow effects for the cylindrical arc. He also showed that under certain conditions induced radial velocity and the associated compressibility effects can dominate thermal conduction effects in a SF6 circuit breaker arc. The addition of axial convection effects (which are obviously important to the gas blast circuit breaker) to the cylindrical arc problem necessitates the use of numerical solution methods, unless far-reaching and generally unrealistic idealizations are incorporated into the analysis. Stine and Watson, for example, added the axial convective transfer term to the steady-state Elenbaas-Heller equation, which was used in the conduction dominated cylindrical arc studies discussed above. They decoupled this energy balance relation from the momentum and continuity equations by assuming a pseudo-onedimensional, uniform pressure, "slug" flow, and thus obtained a

8 closed form solution for the enthalpy distribution. In a later in9 vestigation, Watson and Pegot expanded this work and performed what is perhaps the most comprehensive numerical study of constricted arc-heater flow. Their analysis included both radial and axial convection and real gas effects, but was limited to subsonic flow velocities and steady state conditions (i. e. a DC current input to the arc). Standard finite-difference techniques were used for both symmetric and asymmetric conditions. Numerical analyses that include constrictor (nozzle) geometry effects generally assume the flow is one-dimensional or piece-wise 10, 11 12 one-dimensional. Swanson and Roidt, however, were able to include turbulence and radiation effects in a solution of the conservation equations for an arc core lying on the centerline of a breaker nozzle. The arc geometry was determined through an approximate integral formulation of the thermal boundary layer equation applied at the arc boundary. The nozzle geometry only determined the axial pressure gradient in the arc core and conditions were again limited to those of the steady state and subsonic flow velocities. In the present work we propose to go beyond the above investigations and concentrate directly on the gasdynamic effects involved in gas blast AC circuit breakers. Such an analysis requires the use of sophisticated numerical techniques for solving the unsteady conservation equations. Since the field of numerical gasdynamics is a relatively

9 new one, a review of its terminology, a brief historical development, and consideration of the present "state of the art" will expedite the discussions of following sections. 1. 2 Numerical Solution of the Conservation Equations of Gasdynamics Although the field of research dealing with the computational solution of systems of non-linear partial differential equations has been a viable one for only a quarter century, an extensive literature on the subject has emerged. The general area of fluid dynamics has probably 13 14 been the greatest impetus inthis direction. Sichel and Emmons present overviews of previous work in this area. Owing to a lack of a rigorous mathematical foundation, however, numerical fluid dynamics remains more of an "art" than a science. It is hardly an exaggeration to state that there exist nearly as many numerical methods as there are physical problems. In a manner similar to that of mathematical solution procedures for linear equations, these methods (and physical problems) may be categorized according to the type of equations involved-elliptic (e.g., Laplace's equation), parabolic (the heat equation), or hyperbolic (the wave equation). In what follows, we will be concerned with the area of gasdynamics, in which the working fluid is compressible. The governing conservation equations then form a system of hyperbolic partial differential equations, or a mixed hyperbolic-parabolic system for the case of the compressible

10 Navier-Stokes equations, which contain second order viscous and 15 thermal conduction terms. Cheng gives an excellent discussion of the many numerical considerations involved in solving these equations. 1. 2. 1. Definitions and Terminology To further limit our considerations of this extensive field, we also restrict the discussion to finite-difference methods for the Eulerian form of the conservation equations. Among other things, this rules out our use of other possible numerical techniques, such as the method of characteristics (for purely hyperbolic systems), the method of integral relations, and finite-element theory. Only the first of these offers a practical alternative in the present study, but it has the disadvantage of requiring a much more complicated computer program logic (particularly in the axisymmetric case with 16 three independent variables ). In addition, the inclusion of viscous and thermal conduction effects is not possible for future extensions of the models. The Eulerian viewpoint is preferred here, because it represents a more straightforward approach to the problem of unsteady nozzle flow. Lagrangian coordinates find application to problems that include free-surface boundaries, two or more fluids with different state equations, or certain processes of interest that are associated with fluid elements. Two classical references on finite-difference techniques are 17 18 those of Richtmyer and Morton and Forsythe and Wasow. The

11 former has particular application to the present problem. In general, the technique of finite-difference modeling involves replacing the physical domain of interest by a finite set of nodal points at which the solution to the continuous equations are to be approximated. Wherever the exact solution is analytic, the derivatives are replaced by finite differences which follow from a Taylor series expansion of the exact solution, denoted as f, about any node at (x,t) (say, for the unsteady, one-dimensional case). For example, from 2 af 1 2a f f(x ax,t) =f(x,t)~ ax +.A. (Ax)1.1) ax 2' x,t ax xt we obtain the forward difference approximation, 2 =A f-A"x +o[(Ax)2 (1.2) ax x 2 2 ax x x)t where A f f(x + Ax,t) - f(x,t) (1.3) x Ax and the backward difference approximation, af a 1 2 axf = f+ Ax a 2f +. [(X) (1. 4) ax x,t x,t where V f = f(x,t) - f(x - Ax,t) (1 5) x Ax

12 The addition of (1. 2) and (1. 4) results in the central difference operator, 6 f x f(x+ Axt) f - -x,t) (16) x 2Ax Similar expressions hold for differences for the variable t. There exist many other approximations for first and higher order derivatives, some of which are mentioned below. They all are based on the Taylor expandability of the exact solution of the continuous differential equations to be solved, and thereby assume that the sum of the Taylor series is approximated to a high degree by its leading terms. The remaining terms then constitute the truncation error associated with a given difference scheme. We can define more precisely the truncation error of a differencing scheme and simultaneously explain the notions of consistency, of explicit versus implicit techniques, and of the "order of accuracy" 19 of any scheme. Following the ideas of Moretti, we denote the approximate finite difference solution at x = j Ax and t = n t as F.. j Consistent with the previous notation, f then denotes the exact solution at the same node. The overall differencing procedure is thus n one that arithmetically operates on F. at all nodes in the domain of interest to produce approximate values of Fn+1 (at t = (n+l)At). This procedure may be specified by a global operator, (', such that

13 F" =a!(FJ. (1.7) Now, if F1 is calculated from exact values of the solution, that is, if lF - F(+1 ^fn) (I (1.8) then the truncation error at each node (j,n) is given by: n+l E = F - f. (1.9) Note that if one begins with the exact values, f, which are continuous and satisfy the differential equations, then the ratio F - f J At should approach the same limit as fn+l f3 j At when At approaches zero. That is, F - fn+l lim - ) (litm0 (1.10) At -O0 At - 0 " Difference techniques that satisfy Eq. (1. 10) when applied to a differential equation are said to be consistent. Evaluation of the truncation error E for a given technique (e. g., a given operator, C ) requires the replacement of each f'n in Eq. (1. 8) with its Taylor series j

14 equivalent about the node (j, n). The partial time derivatives that appear in these series can be expressed as spatial derivatives by means of the associated partial differential equations. The truncation error may then be written in general form as = b(Ax)q + [(Ax)q+] (1.11) where At is given in terms of Ax, based on stability considerations (to be discussed below). The order of accuracy of the given differ19 encing scheme is then defined as (q-1). However, Moretti notes that for this definition to be precise: (1) all the Taylor series involved must converge rapidly and monotonically and (2) the factor b, which depends on local values of the dependent variables and their derivatives (and is unique to each difference operator, ), must be of the order of one. In the case of non-linear governing equations, Moretti demonstrates that these criteria are not always met. Schemes that may be represented by Eq. (1. 7) are known as explicit techniques, because the calculation of F+ depends only on J known values of the dependent variables at previous times. As such, these techniques are particularly convenient to use, since they involve a simple "marching forward" procedure from one time line to the next. They have the disadvantage, however, of imposing a limit on the time step of integration to maintain a stable solution procedure. Implicit schemes, on the other hand, which may be represented as

15 F^^l -1(:, Fj), (1.12) ] J J for example, usually need no limit on At to ensure numerical stability. Their major shortcoming stems from the need to efficiently solve a large set of simultaneous algebraic equations. In general, explicit techniques find wider use in gasdynamics. There is reason to believe, in fact, that even with the possibility of larger time steps, implicit schemes are ill-suited to the solution of hyperbolic equations. We discuss this topic in more detail in Sec. 2. 3. 3. A great deal of work has recently been devoted to the question of stability of finite-difference schemes. We will not delve into this extensive topic in detail here; most references on finite-difference methods, such as Ref. 17, devote a great deal of discussion to this concept. In general, any differencing method is said to be stable if any errors (truncation, round-off, boundary errors, etc.) introduced into the computation procedure do not grow in an unbounded manner with increasing time and eventually "swamp" the solution. 20 Unfortunately, stability analyses, such as that due to Von Neumann apply only to linear equations. In the non-linear case the usual procedure is to linearize the governing equations and assume the resulting stability criteria apply locally over the domain of interest. However, for all hyperbolic partial differential equations and systems of equations, such as those studied here, it is known a priori that the

16 Courant-Friedrichs-Lewy stability condition must be met no matter what form of differencing scheme is used. We defer a more thorough discussion of this topic to Sec. 2. 3. 3 on time-step control. 1.2.2 Concepts of Numerical Gasdynamics Computer solution of the conservation equations of fluid dynamics had its start around 1950 with Von Neumann and Richtmyer's work on the ENIAC computer. In their attempts to solve for inviscid flows containing shock waves, they proposed the use of an "artificial viscosity", explicitly incorporated into the governing equations, to numerically account for the discontinuous change in fluid properties 21 across the shock. In a manner similar to the true physical characteristics of a real gas, the flow then undergoes a continuous transition through a "numerical shock wave" of finite thickness. Hence, the use of the Rankine-Hugoniot relations is avoided and finite-difference methods applicable to continuous flows suffice for the entire computational domain. Since that time, numerous investigators have used this device not only for flows containing shock waves, but also as a means for damping numerical instabilities that occur near boundaries, steep gradients, etc. 22 In 1954, Lax22 extended the above ideas by proposing the use of the so-called "conservation form", or "divergence form", of the equations of gasdynamics, as opposed to the better known Euler equations. In this way, the unsteady one-dimensional conservation equations,

17 for example, have the form ag _ ah Da~g a=~ ah a(1. 13) at ax where g represents a vector whose elements are the three dependent variables of interest and h is a vector function of g. Euler's equations, on the other hand, have the form af af f - A Of (1. 14) at ax in which f is the vector of dependent variables and A is a matrix whose elements depend, in general, on f. The innovative feature of form (1. 13), as put forward by Lax, lies in the fact that in the steadystate the Rankine-Hugoniot equations for conditions across a shock 22 wave are automatically satisfied. Moretti points out that this advantage does not carry over to unsteady flow calculations. A further advantage attributed to the conservation form, as opposed to Euler's equations, is that the fluid mass, momentum, and energy are conserved during a finite-difference calculation. As a consequence, this technique has attained great popularity among researchers in computational physics, who have applied its principles to a wide range of fluid dynamic and gasdynamic problems. Nevertheless, there exists some skepticism and convincing evidence that it offers few, if any benefits19 23 if any benefits

18 The use of an artificial viscosity and the conservation form of equations to implicitly "fit" shock waves into inviscid flows entails two major disadvantages. One is the fine mesh size (and corresponding large computational time) needed to maintain a reasonably narrow (but still unrealistic) shock transition region. A large viscosity relieves some of these difficulties, but causes significant errors throughout the overall flowfield. The second disadvantage is the non-physical oscillations that appear on the high pressure side of the shock and become more intense with decreasing viscosity. (See Ref. 24, for example. ) In an effort to relieve some of the difficulty, Lax 25 and Wendroff proposed a new difference technique in 1960, which is second order accurate in both the time and space variables. The scheme involves approximating the dependent variables, f, at each new time line by a second-order Taylor series expansion; hence, 2 af i 2at2f f(t + at) = f(t) + At + - (At) 2 (1. 15) ~a~~t 22 The governing equations in conservation form provide expressions for the above time derivatives in terms of spatial derivatives, which are than replaced with centered differences. No artificial viscosity was used in the original work. Although the technique still produced oscillatory behavior of the solution in the vicinity of the shock, the Lax-Wendroff technique has become quite popular because of its high

19 26,27 overall accuracy. Some impressive results have been obtained The technique also led to the development of several multi-step Lax-Wendroff schemes that avoid the cumbersome form of the full 17 scheme in Eq. (1.15). One of the first of these is due to Richtmyer 28 29 30 31 More recently, Gourlay and Morris28' have developed and studied several variations of these schemes. The final topic relating to computational fluid dynamics that we mention here concerns the time dependent approach, which has found increasing usage in recent years. In general, it simply entails the numerical solution of the unsteady form of the conservation equations. 32 However, in 1965, Crocco advocated this approach as a means of solving fluid flow problems in the steady state. The flow would thus evolve from some specified initial condition and approach the desired equilibrium state asymptotically in time. Since the steady form of the conservation equations are elliptic, this procedure avoids the difficulties often encountered with specifying conditions at boundaries that lie at infinity. Instead one solves an initial-value problem, which is hyperbolic or parabolic in form, with greater freedom in the 33, 34 choice and specification of boundary conditions33. Time-dependent techniques have been used in numerous fluid dynamics problems and all significant references cannot be mentioned here. (See Refs. 33 and 34 for brief surveys.)

20 1. 3 Overview of Present Research As we mentioned above, the primary intention of this work is to numerically calculate unsteady flows in Laval nozzles and investigate the effects of an intense time-varying bulk heat addition on the internal gasdynamics. To accomplish this, we use a modern two-step finitedifference technique which has second-order accuracy in both time and 15 space. Cheng offers convincing arguments that second-order schemes are optimum for time-dependent problems. Our overall technique may be termed time-dependent in the general sense described above, as well as in the more restrictive definition, since we use it to calculate steady-state solutions as well. Neither the artificial viscosity nor conservation form concepts are used. Flow throughout a Laval nozzle experiences large gradients in fluid properties that would cause ar tificial dissipative mechanisms to produce significant errors in the assumed inviscid flow. And since we are not concerned with shock wave phenomena in this analysis, the natural stability characteristics of the chosen difference scheme suffice. In addition, our governing equations cannot be written in a true conservation form, because they contain "source-like" terms due to heat addition, nozzle area variation, and terms due to the axisymmetric form of the conservation equations. The numerical study, as a whole, is significant in that it considers transient internal flows whose velocity ranges from near zero to highly supersonic values in the domain of interest.

21 In regard to the gas blast circuit breaker problem, the analysis focuses on convective energy transfer processes and how they are influenced by nozzle geometry and heat input. As in previous studies, we omit non-equilibrium, turbulence, magneto-hydrodynamic, and electrode effects in order to obtain a tractable problem formulation. The implications of such simplifications are thoroughly discussed 1 6 8 9 72 elsewhere'''' In addition, we consider the working gas to be inviscid, so that numerical difficulties associated with nozzle wall boundary layer effects are avoided. As such, the ensuing analysis implicitly assumes that the internal gasflow follows the specified wall geometry variations; that is, boundary layer separation does not occur. This assumption is expected to be valid for the large inlet-toreceiver pressure ratios, high Reynold's numbers, and nozzle geometries studied here. Furthermore, the results of Watson and Pegot indicate that viscous forces in a plasma-generator nozzle and constrictor have a negligible effect on flow properties when internal passage dimensions are comparable to those used in the present study.

2. QUASI-ONE -DIMENSIONAL MODEL 2. 1 Introductory Remarks The quasi-one-dimensional conservation equations of gasdynamics, which are widely used in the area of internal gasflows, form the basis for the initial approach to numerically modeling the present problem. Although the assumption of one-dimensionality is restrictive in its basic implications, these equations have been found, in many cases, to yield reliable quantitative, as well as qualitative, results; flow through ducts, heat exchangers, turbo-jet engines, and wind tunnel and rocket engine nozzles represent typical examples. The possibility of accounting for (in a manner consistent with the one-dimensional assumption) the effects of heat addition, a variable duct geometry, wall friction, body forces, variable gas composition, and mass injec34 tion, further enhances the applicability of this approach. Crocco presents a detailed discussion of the steady form of the quasi-one35 dimensional equations, while Shapiro's text is perhaps the classical 36 reference in this area. Sentman's work gives a complete and mathematically rigorous derivation of the general, unsteady equations and was the author's chief reference in formulating the governing equations of this section. With regard to the present investigation of flow in a Laval nozzle with a general time-varying bulk heat addition, a one-dimensional 22

23 model offers several benefits. First of all, it provides the opportunity to evaluate a given numerical technique before it is applied to the more comprehensive axisymmetric physical model. We use the term ^numerical technique" here and throughout the remainder of this thesis to signify the overall method used to discretize the physical model and produce an approximate computer solution for a given set of problem parameters. Such a technique would therefore include: (1) a finite difference scheme that accurately integrates the governing partial differential equations within the interior of the domain of interest, (2) a suitable numerical scheme (which is compatible with that of the domain interior) to impose the boundary conditions that define the problem at hand and determine a unique solution, (3) a method to dynamically control the size of the time step used in integrating the equations and thereby ensure a stable solution procedure, and (4) initial values of the dependent variables that either represent the true physical initial conditions or, in the case for which only a steady-state condition is of interest, prevent numerical instability. Other more subtle ingredients of the numerical technique are the forms in which the differential equations are written (e.g. conservation versus non-conservation form), the size and shape of the finite-difference grid used, and methods for reducing numerical data and displaying the final results. Although the additional spatial (radial) coordinate of the axisymmetric governing equations introduces numerical considerations not present

24 in the one-dimensional model, virtually all procedures developed in the latter case carry over to the higher dimension model. For example, the manner in which the inlet boundary condition is incorporated into the numerical model proves to be more critical in determining a stable numerical solution than conditions at the nozzle exit, axis of symmetry, or wall. In addition, the numerous exact analytical solutions available for quasi-one-dimensional flows allows us to demonstrate the validity of the technique that is ultimately adopted. In short, the one-dimensional model is a logical "numerical stepping-stone" to the axisymmetric problem. Another benefit to be derived from the one-dimensional approach concerns the physical implications that may be drawn from the resulting solutions pertinent to the circuit-breaking potential of Laval nozzles. The unsteady quasi-one-dimensional equations yield significant qualitative information pertaining to the gross effects of an "arc-like" bulk heat addition on the mean flow properties within the nozzle, and to the roles played by axial convection and nozzle geometry. One parameter of interest to the circuit-breaker designer, for example, is the thermal time constant associated with the decay of the static temperature (and hence, with the decay of electrical conductivity) both near current zero passage and during the "free decay" from a high temperature steady state. Also, the transient phenomenon of

25 reverse flow near the nozzle inlet during the first moments of energy addition (which is one undesirable form of circuit-breaker`clogging"), the validity of the quasi-steady assumption commonly used in analytical treatments of the "near peak current" period of an AC arc current input, the effects of the time lag between heat input and gasflow response, and the importance of supersonic flow velocities represent additional phenomena of importance that may be examined through a one-dimensional analysis. And finally, we should mention that the numerical solution of the quasi-one-dimensional conservation equations governing mixed subsonic-supersonic flow in a Laval nozzle is a non-trivial problem. Some of the more recent investigations along these lines were mentioned in the Introduction. Most entail considerations of the steady state only; and, to this author's knowledge, the case including an intense time-varying bulk heat addition has not been previously investigated. In the remainder of the present section, we detail the form of the governing equations used in the one-dimensional model, the overall numerical technique involved in their solution, and the subsequent results obtained from this approach. 2. 2 The Quasi-One-Dimensional Governing Equations The quasi-one-dimensional equations of gasdynamics express the laws of conservation of mass, momentum, and energy of a fluid flowing

26 within ducts of varying cross-sectional area. The prefix "quasi-" generally reflects the inherent inconsistency between the assumption of one-dimensionality (with the coordinate, x, directed along the duct axis) and the inclusion of area variation effects. One further assumes that the effects of any mechanism, such as heat addition, are "felt" by the fluid at each x instantaneously and uniformly over the local duct cross-section. Consequently, the fluid properties that appear in the following equations represent mean values at each axial station. The equations can be derived either by considering a control volume fixed 35 36 in space3, or by following the motion of a fixed-mass fluid element36 They do not, in general, follow from the three-dimensional conservation equations as do the true (exact) one-dimensional equations, which hold along a streamline in an inviscid three-dimensional flowfield. Proceeding with the fomulation of the working equations, we write the continuity equation in general form as (pA) + (puA), (2.1) at ax assuming there are no mass sources or sinks present. The fluid density (p) and velocity (u) are functions of both time (t) and the single spatial coordinate. Since the duct cross-section area, A, is a function only of x, we have by the above equation that ap ap a+ pu+ + = (2. 2) at ax ax with A- d( n A)/dx.

27 If we now assume that the fluid is inviscid and non-heat-conducting, there are no body forces acting on the fluid, and a state of thermodynamic and chemical equilibrium exist at each x and t (so that, in the usual manner, the normal stress within the fluid is equated with the equilibrium thermodynamic pressure, p), then the momentum equation takes the familiar Euler form in which at au + u a (2. 3) at ax pax and u represents the fluid velocity. It is interesting to note that the duct area, A(x), does not appear in Eq. (2. 3); hence, the equation holds regardless of the spatial (or temporal) variation of A. Finally, we write the energy equation in a form that is analogous to the First Law of Thermodynamics for a fluid element; De Q D(/p) (2 Dt p Dt where e denotes the specific internal energy of the fluid and Q represents the net rate of heat addition per unit volume (watts/m3) at any point along the duct. In general, Q is a function of both x and t and accounts for duct wall heat transfer effects, radiation, and all other forms of energy addition or removal. A discussion of the specific form of Q used in the present study will be deferred to the section on the one dimensional results.

28 If the working fluid is now assumed to behave as a perfect gas with p = pRT (2. 5) and e c T (2.6) where cV is the specific heat at constant volume, then the foregoing conservation equations may be written with p, u, and T as the dependent variables. For the remainder of this thesis, we adopt the standard subscript notation to denote partial differentiation. Equations (2. 2), (2. 3), and (2.4) then become for a perfect gas: Pt = UPx - P (U + u) at -uux R (Tx + x. (2.7) ut x x p Tt =-uT -(y- 1) T (u +) +(Q/pc) y represents the ratio of specific heats cp/Cv. Equations (2. 7) may now be normalized with respect to a characteristic length, L, a characteristic velocity, V, and some reference ref' values of the fluid properties. We therefore introduce the following dimensionless variables: P*-,u u-, T*- T P ref V T ref Vref Tref *p P A, Q P* p, ref - ref

29 where e ref. ref v -1 Qref — L-' c ~* - - = (y -1) ref L v R V t* ef t and x* L L The reference velocity is, for convenience, given by the relation V,=(p /p )1/2 ref (Pref/Pref/ The working equations now take the following non-dimensional form: = - up - p (u + uX) Ut= -uu - T -(T/p) Px (2.8) Tt =- uT - ( - 1) T(u +uA) + (- )(Q/p) in which we choose to drop the asterisk superscript to simplify the notation and will re -introduce it only when it is necessary to reference dimensional values in relation to their dimensionless counterparts. In the numerical results to follow, L is the nozzle length and the fluid upstream reservoir conditions are used as the reference state. 2. 3 Numerical Technique The task of choosing a numerical scheme to integrate a nonlinear, coupled set of partial differential equations, such as (2. 8) above, is complicated by the fact that the literature abounds with numerical techniques applicable to fluid dynamic problems, almost all of which purport

30 to offer certain advantages that make them preferable to others. We described some of these in the Introduction. The Lax-Wendroff (L-W) technique perhaps comes closest to being a "universal" scheme in fluid dynamics. It offers second order accuracy in both time and space (in the sense previously defined), numerical stability (with the aid of an additional artificial damping term in certain problems involving stagnation points, sonic lines, and/or shock waves near an axis of sym39 metry ), and a certain proven versatility in that it has been successfully applied to a wide range of inviscid flows. However, the inclusion of diffusive second order terms in the working equations (most notably, those terms involving viscous and heat-conduction effects) renders the L-W differencing scheme virtually useless, because of the great number and complexity of the resulting terms in the differenced equations. Even in the case of multi-dimensional inviscid flow in which an auxiliary mapping of the physical domain onto a rectangular computational domain is necessary (such as the mapping used in the axisymmetric model considered in the present study), one must question the practicality of the method (see, for example, Ref. 27). In these respects, then, the Lax-Wendroff scheme lacks efficiency and flexibility. In the present investigation, in which we perform a partial parametric study of the effects of nozzle geometry and type of heat input, an efficient numerical technique is essential. In addition, the technique should be flexible enough to handle both single and

31 multi-dimensional flow cases, as well as diffusive fluid mechanisms, such as thermal conduction, for future extensions of the axisymmetric model. In light of these requirements, the two-step L-W integration schemes previously discussed offer a workable alternative to the full Lax-Wendroff method while maintaining second order accuracy. 37 MacCormack recently developed such a scheme to solve the unsteady, axisymmetric, Navier-Stokes equations describing the fluid-like cratering of hypervelocity cylindrical aluminum pellets on semi-infinite aluminum targets. Consequently, the equations contained a coefficient 19 of viscosity of large magnitude. Moretti conducted a very rigorous comparison of this scheme with several others prominent in the literature. He considered the initial compression period of the one-dimensional "shock-tube" problem in which a piston is accelerated linearly from rest into a quiescent gas. By terminating the solution before the compression waves coalesce into a shock, he avoided the problem of consistency of the finite difference schemes when a discontinuity is present in the flow. After comparison of the results with the known exact solution, MacCormack's scheme proved to be superior from the standpoints of accuracy and ease of programming. This finite difference method has since been applied to a variety 38 of problems. Anderson used it in the time-dependent solution of Laval nozzle flow with chemical and vibrational relaxation effects

32 40 present-as applied to gasdynamic laser devices. Moretti studied both the steady state and transient asymptotic approach to the steady state for inviscid compressible flow past a circular cylinder. Moretti 41 42 43 44 and Pandolfi', Moretti et al., Kutler and Lomax, and Thomas 45 et al., all applied MacCormack's scheme to the steady, inviscid, three-dimensional conservation equations describing supersonic flow past various body and wing-body configurations at angle of attack. (These works offer interesting contrasts as to the manner in which discontinuities and severe flow property gradients are handled numerically, and represent excellent examples of the differing philosophies present in the modern computational fluid dynamics literature.) And finally, MacCormack has applied his original scheme, as well as drastically altered forms of it, to the time-dependent computation of the flowfields involved in the shock wave-laminar boundary layer inter46 47 action problem and supersonic flow over a double wedge. We now apply the original form of this scheme to the problem at hand with the ultimate test of its applicability here being the quantitative and qualitative character of the subsequent results. 2. 3. 1 Interior Points In order to discuss in a systematic manner, the overall numerical technique used in this study, we distinguish between those points, or nodes, in the computational domain that constitute either the domain

33 boundary or its interior. In the one-dimensional model, the equations constitute an initial-, two-point boundary-value problem with boundaries located at x = 0 and 1. Correspondingly, the interior comprises all nodes of the finite-difference mesh, except two. These interior nodes are discussed first. To facilitate the illustration of MacCormack's integration scheme, we write Eqs. (2. 8) in an equivalent vector form as: f - Cf +S. (2.9) t x The column vectors, f and S, represent the dependent variables and "source-like" terms in the governing equations, respectively; Pl f= u (2.10) and - puA S 0. (2. 11) ( - 1)gQ/p) - Tu] The coefficient matrix, C, is a function of f; u p 0 C = T/p u 1.(2.12) 0 (y - 1)T u

34 We now divide the domain of interest uniformly in the axial direction and seek approximate values of f at a discrete set of J + 1 nodes spaced Ax apart at each value of time t = nA t(n = 0,1, 2,..), where n 0 < x < 1 and Ax = 1/J. In the usual manner, denote the finite difference value of f at any axial station along the nozzle, where x = (j - 1) Ax (j = 1,2,... J, J + 1), by f. Then, given initial values of f. at 3 e time t, the first, or predictor, step of MacCormack's scheme estabn n+l lishes temporary values of f., denoted fJ, at time t + A t by means n of a simple first order Euler method with forward differences in both x and t. Thus, fn+1 n n n f. f.1 f" -f. n j+1 n T] — Cn j+l fJ + Sn (2. 13) At Ax Defining X = A t/Ax and re-arranging (2.13), we write the predicted n values explicitly as: fn = f. - XC. (f1 f) + tS. (2.14) j+1 3 n 3 n+1 The final corrected values at t + A t, f., follow from what amounts a t a o t a t b nning o t t to a temporal average of those at the beginning of the time step, fj, and values at t = (n + 2) Ant, denoted here by f., which are based on the predicted values given by (2. 14) and first order backward differences in x. Hence,

35 ~n+l 1 n n+2) n+f (f + f (2.15) 3 j 33 where n+2 n+l n+l n+1 n+in f +2 f. (f1 f -1 + +A t S.. (2. 16) 3J 3 J 3 ~3-1 n ] The manner in which the time step, Ant, varies with n will be discussed below. In practice, Eqs. (2.15) and (2.16) are combined so that, for example, the first component of (2.9), the continuity equation, appears in the corrector step as: Pn+l 1 n n+ n+1 n+1 n+1 n+n+ n n+ 2ivit'( -u l P = P + P - AL (j - Pj-1 ) + P (u - uj-1 - A t p u. A.. (2. 17) n 3 3' Before proceeding further, we introduce a more convenient notation than the obviously cumbersome one displayed in Eq. (2.17) above. Since each finite difference equation refers to computations about the general node (j, n) in the (x, t) plane, only the indices different from this central node need be specified. Therefore, f. may simply be written as f, f j as fl. and so forth. There should be no confusion between j+1 j+1' the exact and approximate (differenced) values of f, because the context in which they appear will make the distinction obvious. Furthermore, let the predicted values of f be denoted as f, and A t as simply n

36 At. Then, for example, the general,vector form of the corrector step, Eqs. (2.15) and (2.16), becomes: 2 l + f - XC (f- fj ) +At S (2.18) This notation will be particularly convenient for the discussion of the axisymmetric model. We now write out the full form of the finite-difference equations corresponding to the domain interior (2 < j < J). For the predictor step, Eqs. (2. 10) through (2. 12) and (2.14) yield: P = p - X[u(pj+ - P) + P(uj - U) - At pu u= u - X[u(u+ - u) + (+T. - T) + T( - p)/p] (2.19) T = T X[u(T -T) + (y - 1) T(u - u)] + At (y - l)[(Qd/p) - TuA] For the corrector step, Eqs. (2.18) yields: n+1 I - P = -{P +P - xP - -1 + p( - u1) At pu n+1^ I {u+u-T + T( 1) +(-T + u = 2-u + u - x[u(u - ujz) +(T < - T + - +j 1)lP] (2. 20) Tn+l= (T + T (- [u T. ) + (Y - 1) T(i - u )] + At ( ) / - Y)[(Q/) - Tu]}

37 2. 3. 2 Boundary Points Boundary conditions, of course, determine the uniqueness of any mathematical solution. In the present, one-dimensional, case we must consider the two permeable boundaries of the computational domain corresponding to those at which the gas flow enters (inlet; j = 1, x = 0) and leaves (exit; j = J + 1, x = 1). It is interesting that texts on the fundamentals of numerical analysis of partial differential equations, such as Ref. 17, offer little guidance in the imposition of boundary conditions; and yet these are the very conditions that distinguish a given physical model. The bulk of the literature puts little emphasis on this aspect of numerical techniques and there seem to be as many methods for numerically simulating conditions at the boundaries of the computational domain as there are solutions. One must turn to the more recent literature for logical guidance 48 19 in this area. Moretti is perhaps the most vehement critic of many of the recent trends in computational fluid dynamics. He stresses the importance of boundary conditions to a numerical technique and suggests the manner in which the various types of boundaries (rigid wall, permeable boundaries, symmetry axes or planes, and flow discontinuities) should be handled. His ideas are based on an adherence to the physical principles underlying the fluid flow in a given problem, as well as to the corresponding mathematical

38 characteristics of the associated equations. As an example, consider the case, which is pertinent to the axisymmetric model of the present study, of a fixed rigid boundary in an inviscid flow. Mathematically, it is well known that the proper boundary condition to be applied is the vanishing of the normal component of the velocity. Hence, to determine the density, temperature, and the tangential velocity component at the wall (which are needed in each step of the difference solution), Moretti notes that one must use information only from the domain interior plus this one boundary condition; otherwise, the problem will be over-specified and therefore will not be consistent with the physical problem considered. And yet, it is most common in the literature to use the so-called "reflection rule" at wall boundaries. (See, for 27 26 example, Serra and Burstein.) In this technique, an imaginary row of finite-difference nodes is added "inside" the wall and the density, temperature, and tangential velocity values there are related to those within the interior such that their gradient normal to the wall vanishes. (That is, they are locally even functions with respect to the wall boundary.) In this manner, the same finite-difference relations conveniently apply at both the boundary and interior points. However, these boundary conditions do not correspond to the true physical situation, except under very special circumstances. Nevertheless, investigators have obtained seemingly valid results through their use.

39 44 Limited studies of the reflection rule have both confirmed and 48 denied its value. This apparent contradiction, and others associated with the numerical solution of non-linear partial differential equations (for example, the use or non-use of the conservative form of the governing equations), is possibly best explained by the relatively large degree of error allowed for in finite difference approximations. Any error in the boundary values of the same or higher order as the truncation error of the difference scheme will generally be inconsequential. Even in the case of lower order errors, the pseudo-diffusivity implicit in the differencing scheme may sufficiently damp these errors at points away from the boundary before they affect the overall quality of the results. For reasons such as these it is often difficult to draw an analogy between the differential and differenced problems in order to determine whether a computational model is well-posed; it will generally appear to be over-specified. The case of first-order partial differential equations differenced by means of a second-order scheme, such as in the present study, represents a good example. Cheng5, who discusses some of the above ideas, notes that by using a second-order difference scheme one assumes the first and second derivatives of the Taylor series expansion of the difference quotient are non-zero. The subsequent

40 evaluation at each point of the second-order derivative implicit in the scheme therefore requires values of the dependent variables at two neighboring points and provisions for two boundary conditions, as opposed to the one condition normally expected for a well-posed differential problem. The point to be made by the above discussion is that errors due to the numerical domain boundaries are unavoidable. The proper imposition of extra boundary conditions requires a knowledge of the exact solution; and this is, of course, unknown. These errors need not be disastrous as long as their magnitudes are small and they can be damped by the natural dissipation associated with the differencing 15 scheme. In regard to the latter point, it is well-known that the high wave number components of an error spectrum are damped at a greater rate; therefore, the use of extrapolated boundary values, whose related errors change sign every few integration steps, is preferable to specifying fixed values of the dependent variables at the boundaries (whose related errors are not damped at all). In fact, the reflection rule previously mentioned probably owes its "success" to these same principles, since the associated error is not a fixed value and likely changes signs frequently during the course of the integration. With the above ideas in mind, we proceed to specify the boundary conditions for the one-dimensional model.

41 Consider now the downstream boundary as being located at the nozzle exit. As mentioned in the introduction, we consider only the case of negligible back pressure so that the flow is unimpeded and may freely seek supersonic velocities in the diverging region of the nozzle. As a consequence, the dependent variables at the exit plane are independent of downstream conditions and may be determined by means 48 of a simple extrapolation of upstream (interior) values4. (Note that the interior finite-difference equations (2.19) do not hold here, because 27, 38 they require a forward spatial difference.) In similar studies the authors have used a simple linear extrapolation on the grounds that any perturbation caused by the associated truncation error cannot travel upstream in a supersonic flow. While this is true in the continuous differential problem, a consideration of finite-difference equations, such as (2. 19) and (2. 20), shows that any error will indeed propagate "upstream" at a rate of one mesh node per time step. Initial numerical results of the present model based on a linear extrapolation at the exit confirmed this by displaying oscillatory spatial distributions in the dependent variables near the exit boundary at each time step. We therefore adopt a second-order extrapolation formula, which is consistent with the accuracy of the difference scheme used and has been found to yield good results. Hence, the dependent variables at the exit (j = J + 1) are calculated in both steps from values at the three adjacent upstream nodes;

42 J+1 J-2 -3 (Tl - J) (2. 21) in the predictor step, and f n+t fn+ 3 n+- fj+ (2. 22) in the corrector step. The use of Eqs. (2. 21) and (2. 22) implies that the third partial derivative of each dependent variable with respect to x is equal to zero. Finally, we note here that this boundary condition has always resulted in supersonic flow velocities at the nozzle exit. The inlet proves to be the most interesting and, correspondingly, the most critical boundary in the determination of stable and physically reasonable results in the present study. Several possibilities exist in the treatment of this boundary, more than one of which were actually tested here by numerical experiment. Each warrants mentioning, because it helps to illustrate the importance of this aspect of the numerical model and the considerations involved in developing a technique that is applicable to both the one and two-dimensional problem. The simplest method for specifying the upstream state entails essentially no restriction on the flow. One assumes the inlet boundary to be located in the converging portion of the nozzle. After determining the flow properties within the interior domain, the model evaluates the dependent variables at the upstream boundary by means of an extrapolation of adjacent interior values. A consideration of the nondimensional difference equations shows that with this method the

43 interior flow obtains no information about an upstream reservoir, or lack thereof. Hence, if the initial conditions correspond to true isentropic flow and no heat is added (Q 0) for t > 0, then the numerical computation results in the flow finally coming to rest in a portion of the converging nozzle region and reversing its direction there,while maintaining a positive axial velocity in the diverging region. The computer run eventually terminates when the nozzle is vacated and the density values become negative. One concludes from this somewhat naive first approach that some sort of inlet flow condition must be specified in order to simulate the existence of an upstream reservoir. And yet, to be accurate, any such specification of flow properties necessitates a knowledge of the true time-varying solution at the inlet. 19 48 Moretti suggests a solution to this dilemma, which is both mathematically and physically reasonable. For the present case of an internal flow, one assumes the nozzle is fed by an infinitely long constant area duct. Since the initial perturbation formed by the transient flow adjustments within the nozzle can never reach the inlet (where x = - oc), the flow properties there do not change with time and may, therefore, be fixed at their initial values. To allow for the semiinfinite physical space (- c < x <, t > 0) computationally, one maps it onto a rectangle by means of an appropriate stretching function, Z(x).

44 The governing equations are then modified only by the inclusion of a coefficient, proportional to dZ/dx, in each of the spatial derivative terms. Without presenting the details, we mention that the above scheme was briefly tested during the initial development of the one-dimensional model. Several disadvantages associated with its use became apparent. The most significant involves the large number of computational nodes needed for a stable solution-two to four times the number necessary if only the nozzle interior is considered. There are two reasons for this. First of all, the node immediately adjacent to the inlet boundary must be sufficiently displaced from the nozzle so that the first transient disturbance does not reach it in the time to which the computation 48 extends. Also, good resolution must be maintained in transition 50 regions, such as the nozzle inlet at which the area begins decreasing and the heat front at which Q first becomes non-zero. As discussed below, the time step involved in any explicit differencing scheme is limited by stability considerations and, for the present problem, becomes smaller as the static temperature increases with an increase in the amount of heat addition. Hence, the need to "carry along" many computational nodes, at which the information is essentially irrelevant, soon becomes intolerable. An extension of these ideas to the non-adiabatic axisymmetric model is entirely unrealistic if

45 computer processing times are to remain within practical bounds. The notion of mapping the inlet plane to - co was therefore rejected. The nozzle inlet plane represents the most logical and convenient choice for the upstream boundary, and is the one that we use to complete our model. In his quasi-one-dimensional representation of 51 the flow in a gasdynamic laser nozzle, Anderson assumed that the nozzle inlet was essentially in the reservoir (that is, the inlet area ratio was greater than six) and that the fluid density and temperature therefore remained fixed with time at their stagnation values. The inlet velocity, however, was allowed to vary. Its magnitude was computed from a simple linear extrapolation of the values at the two adjacent interior nodes. We note that, based on our initial considerations of numerical boundary conditions, these restrictions should be quite inaccurate, because an error of constant magnitude that cannot be damped by the difference scheme is introduced into the interior at each time step. Also, a linear extrapolation is inconsistent with the use of a second-order numerical technique. In what follows,therefore, we extend the above ideas in a manner that allows the inlet static properties of the fluid, as well as the velocity, to vary with time. To simulate an upstream reservoir, the total, or stagnation, pressure and temperature are assumed fixed at reservoir values. This specification of inlet conditions is precise for steady-state nozzle flow. During transient flow periods it is accurate in the limit of small flow (pressure) disturbances

46 at the inlet plane. However, whatever errors are generated at the inlet have a random character and can be damped away from the inlet by the natural dissipative mechanisms associated with the differencing scheme. In fact, results to follow indicate that variations in inlet flow conditions have only a local effect on the interior nozzle flow. We consider a physical model of the circuit-breaker arc-nozzle configuration equivalent to that represented by Fig. 1-1 of the Introduction. The upstream electrode is located just downstream of the inlet plane, such that any heat addition due to Joulean dissipation within the arc begins affecting the gasflow at x = 0. In this manner, since the fluid is nonheat-conducting, the inlet flow remains adiabatic. Now, visualizing isentropic flow upstream of and including the inlet, we assume that the stagnation, or "total", values of the fluid properties at the boundary (specifically, the total density, temperature, and pressure) remain constant. Theoretically, the total conditions at any point in a flow are those that would exist if the flow were brought from its given state to one of zero velocity by means of a steady isentropic process. By definition then, the total pressure and temperature for a perfect gas are given, in dimensionless terms, by: Pt =p +-1 27/1 (2.23) and

47 Tt T(I + 21 M2) (2.24) in which M2 =u/(TT). (2.25) Since the reference state corresponds to that of the upstream reservoir, as mentioned previously, pt and Tt assume fixed values of unity and the above equations reduce to a pair of relations for the inlet static pressure and temperature. We have, therefore, that T= ( - )u2 (2.26) and p [+ ( u2/T 17 (2.27) at the inlet. Given a value for the velocity, u, Eqs. (2. 26) and (2. 27) and the perfect gas law, p=pT, (2.28) determine the flow at the upstream boundary. The inlet flow velocity needed in the above equations is established by means of the differenced momentum equation for j = 1. The predicted value, ui, follows directly from the second of Eqs. (2. 19): u = ul - [u (u2 - u) + (T2 - T1) + T (P2 - P1)/P1] (2. 29)

48 In order to take the backward spatial difference required in the corrector step, Eqs. (2. 20), we imagine an additional "virtual" node upstream of the inlet, where x = - Ax, and compute values of the dependent variables there with a second-order extrapolation. If we denote these values by f, it follows that T 3(f -f2) + 3 (2. 30) o 12 3 This extrapolation is consistent with the overall accuracy of the numerical technique, and is further justified by the generally smooth spatial distributions of the dependent variables near the nozzle inlet. n+I The corrected value of the inlet velocity, u1, is now given by the equation: uI1 1={u1 + u - [ - o + (T1 - To) + T (P1 - o/l]} (2. 31) n+l Given values of ul, or u, Eqs. (2. 26) through (2. 28) determine both the predicted and corrected values of the inlet static temperature and density needed to continue the numerical integration. One final consideration will complete the overall formulation of the inlet boundary condition. The preceeding relations are sufficient to determine the flow conditions at the nozzle inlet plane as long as the massflow there remains positive for the timespan considered. In this manner, only the upstream "cold" flow traverses the boundary. However, in the case of a large amount of heat addition to the

49 downstream flow, the fluid velocity in the nozzle inlet region may momentarily become negative; that is, the nozzle may "clog". When this occurs the total energy, and hence, the total temperature, at the inlet no longer remain constant, but increase due to the resulting convection of the heated gas from the nozzle interior. To allow for these circumstances in the numerical model, we let the dimensionless inlet total temperature increase above a value of one whenever Tu, in the predictor step, or u1, in the corrector step, becomes negative (as calculated in either Eq. (2.29) or (2. 31)). The static temperature value necessary to evaluate Tt1 by Eq. (2. 24) is obtained from the differenced energy equation similar to (2. 19) and (2. 20). Hence, with Q1 = 1 = 0 we have that T1 T1 - [u1 (T2 T1) +(- 1) T1 (u2 ul)l -(y -1) At T1 u1 (2.32) and n+l I T1 = 2{T1 + T1 - U I (T1 - - 1) T (u 1 )] -(y-1) AL /tT1uT II} (2.33) whenever ul, u1 < 0 and/or Ttl, Ttl > 1.0. The inlet total pressure, Ptl, remains fixed at the reference value as before to simulate the presence of an upstream reservoir.

50 The above boundary conditions as a whole have been found to yield good results, both qualitatively and quantitatively, and they are consistent with a stable numerical technique as long as the incremental time step, At, is properly controlled. 2. 3. 3 Time Step Control As mentioned in the Introduction, any explicit finite-differencing scheme, such as MacCormack's, has associated with it some limitation on the time step used. For the case of a purely hyperbolic system of differential equations, a necessary condition for convergence of the finite-difference representation to the exact solution is the so-called CFL, or domain of dependence, condition first put forward by Courant, 52 Friedricks, and Lewy in 1928. In mathematical terms, this criterion requires that the domain of dependence of the partial differential equations be contained within that of the corresponding difference equations. For the present one-dimensional model with fixed Ax, this translates into the relation: < Ax At<(u ). (2.34) (u + a) max We illustrate the significance of the CFL criterion by considering a general interior mesh node located at (j, n + 1) in the x-t plane, as shown in Fig. 2-1. The finite-difference solution at this node (point A in the figure) given by Eqs. (2. 14) through (2. 16), depends on values

Slope =At/Ax Chorocteristic: Slope= (u+ + n + I l a A(j,n+l) a t t = (n I )Att = nAt — / At /r, ~-,or ~ m -,,u.m.. \\I x O BI D D E'C Domain of Depe of Point A mdence - ---- 6Ax Finite- Difference Mesh i-i i j+l Figure 2-1. Illustration of CFL Stability Criterion

52 of the dependent variables at only the three additional nodes on the previous time line t - nAt (shown as darkened circles). These four nodes constitute the "computational molecule" of MacCormack's method. The numerical domain of dependence of the point (j, n + 1), therefore, corresponds to the triangle ABC. On the other hand, a consideration of the characteristic curves converging to point A, which locally have slopes of (u + a) in the x-t plane, shows that its domain of dependence in the differential formulation is the triangular region ADE. (Actually, the characteristic curves are not straight lines as shown, but may be represented as such to a good approximation in the limit of small At.) It is clear that Eq. (2. 34) requires ADE to lie within ABC, given that the sum (u + a) is a maximum at point A for the time line t = (n + 1) At. In practice, the characteristic slopes must be evaluated at t=n At, because the values of u and a at the succeeding time line are a priori unknown. To allow for curvature of the characteristics, a "safety factor, CSF is incorporated into the criterion, which now takes the altered form: Ax nC Ax At = C uAx... SF (u[ +a,)' SF u +(7T1/2; 3 3 j 0 < CSF < 1 (2.35)

53 in which the speed of sound has been non-dimensionalized with respect to Vref The CFL criterion may now be given an alternative interpretation that is perhaps more intuitively appealing. Figure 2-2 shows the characteristics that issue from the nodes at B and C. It follows from the theory of hyperbolic partial differential equations that data along BC determines the solution at points within the region BCD alone, the solution at other points in the x-t plane being undetermined by data only on this interval of the t = nAt time-line. Hence, Eq. (2. 35) assures the consistency of our numerical technique by restricting the location of point A to the region BCD. To conclude this discussion of time step control, a few additional notes are in order. First of all, a further consideration of Fig. 2-1 indicates that for CS < 1 the numerical scheme makes use of more SF information than it needs in order to calculate the solution at point A, since the segments BD and EC of the t = nAt time-line are outside its 19 interval of dependence. Moretti mentions that this inaccuracy is not detrimental to the numerical scheme, because values along the true interval of dependence, DE, are interpolated from those at the three nodes on BC. One would expect, however, that a greater accuracy results from values of CSF closer to (but not greater than) one. This was, in fact, verified by the numerical results of the present study both quantitatively, for the steady-state isentropic and

Slope n n = (u? + a. j j -I ) Slope CSF( en 4b. t = nt -- j j+l j+2 Figure 2-2. Alternate Interpretation of CFL Criterion.

55 "Rayleigh process" cases for which an exact solution is known, and qualitatively, for all cases which as a whole indicated better stability characteristics for larger values of CSF. Other studies seem to 53 verify these ideas as well. For example, Rubin and Burstein show a definite improvement in the qualitative aspects of their results for a standing shock wave as the value of CSF is increased up to a value of 0. 98. For C = 1. 0, however, they found complete instability SF of their techniques. With the above discussion in mind, we interject here some related comments regarding the relative merits of implicit and explicit numerical techniques for solving hyperbolic partial differential equations. Generally, implicit differencing schemes result in a set of finitedifference equations for which the domain of dependence of a mesh point, such as point A, is the entire preceding time line. (See, for example, Refs. 54 (pp. 505-512) and 55 (pp. 49-54).) As a result, the domain of dependence of the differential equations is always contained within that of the difference equations; that is, the CFL criterion is satisfied for all values of At. Although the ability to choose a large time increment would be very desirable, we note that the use of an implicit scheme with hyperbolic equations is similar to having a very small value of CSF in an explicit method. Because most of the information used to calculate the solution at each point of the succeeding

56 time-line is extraneous, this should result in a very inaccurate solution. In this sense then, explicit differencing methods are more consistent with the mathematical nature of hyperbolic equations. On the other hand, implicit schemes correspond more closely to the nature of parabolic and elliptic equations. It is not unusual, therefore, that explicit techniques are predominantly used in the area of numerical gasdynamics. Finally, note that the CFL criterion is a necessary condition for stability and convergence of the numerical scheme. Equation (2. 34) represents an upper bound on the value of At used during the course of the integration. Several explicit differencing schemes, such as the full Lax-Wendroff method, have stability properties that are more restrictive than the CFL condition. Although the numerical technique used here, based on MacCormack's scheme, was found to be stable for values of CSF as high as 0.95 in certain cases of a "well-behaved" solution, values from 0. 75 to 0. 85 were used for the general "hot flow" runs in which large transient variations of the dependent variables existed and the characteristics,such as those in Fig. 2-2, had large curvatures. Furthermore, since the sum of the fluid velocity and speed of sound is always greatest in the diverging portion of the nozzle, we evaluate Eq. (2. 35) at only one or two nodes in this region of the domain. In fact, it has been found that in the majority of cases the minimum time step is determined at the nozzle exit.

57 2. 3.4 Initial Conditions A specification of the initial values of the dependent variables, f. (1 < j < J + 1), necessary to start the integration, completes the formulation of the quasi-one-dimensional numerical model. If one desires only a steady-state solution and the preceding transient flow processes are irrelevant, virtually any initial distributions of the flow variables that do not cause numerical instabilities will suffice as initial conditions. Experience with this numerical model has shown that instabilities do not appear if the f. are consistent with the physical assumptions inherent in its development. For example, the initial exit Mach number should be greater than one as assumed by the extrapolation boundary condition there. For the majority of cases considered in this study, however, the initial flow transients are of interest. Consequently, a typical run is "started" with the exact steady-state isentropic flow values as initial conditions and, with Q = 0, the integration allowed to proceed until the corresponding finite-difference isentropic solution is established. (The criteria used to define the achievement of a "steady-state" in the numerical results are outlined in the results section.) This latter state then acts as the initial condition for cases in which an initial isentropic flow is assumed to exist. Although the difference in the values of the dependent variables between the exact and finite-difference solutions is less than one per cent,

58 as demonstrated by the results to follow, this procedure eliminates any non-physical initial transients caused by the truncation error. For cases in which some "hot flow" steady-state, corresponding to a given heat input in the nozzle, is required as the initial condition, the dependent variable distributions established by the difference equations themselves are used. 2.4 Quasi-One-Dimensional Results and Discussion 2.4. 1 Validation of Technique It is not uncommon in the literature for the author of a numerical study of some physical or engineering problem involving the solution of non-linear partial differential equations to simply assume his technique produces results of a specified order of accuracy. In fact, the investigators in some experimental studies often refer to numerical solutions, to which they compare their results, as "exact". However, one can rigorously demonstrate only the consistency of a finite-difference scheme when the equations to be solved are non-linear. As discussed in the Introduction, information concerning the order of accuracy, stability, and convergence properties of the scheme follow only after the working equations have been linearized and, as such, this information is not always reliable when applied to the actual numerical solu19 tion. Moretti mentions several examples of this and shows that

59 even the simple first-order Euler differencing scheme, which is unconditionally unstable in the linear case, produces relatively good results compared to other first-order techniques when applied to the one-dimensional conservation equations of gasdynamics. 37 46 MacCormack' performed the linearized analysis of the present scheme and we do not repeat it here. To demonstrate the value of the results produced by the numerical technique set forth in the preceding sub-sections we recall the previously mentioned studies that involved the use of MacCormack's method (for both steady and non-steady problems), particularly those of Moretti, and carry out our own comparison with the exact analytical solutions of two wellknown steady-state problems. First consider the case of isentropic flow in a Laval nozzle; that is, with Q = 0. 0 for all t and given some initial distributions of the dependent variables we allow the solution to approach a steady state condition asymptotically in time. Recalling that the present model assumes the back pressure to be negligible, we expect this steady state to correspond to the well-known supersonic branch of the isentropic area ratio-Mach number relation56: 2 1 2' - 1 2 ^A =+- L+1 -r-+ M2)] (2. 36) A =2 7+ 1 + 2M

60 in which, for convenience, the nozzle area has been non-dimensionalized with respect to the throat area, Athr. Besides demonstrating the credibility of our numerical technique, this simple first case illustrates some of the transient response characteristics of Laval nozzle flow. It also serves as a reference case against which we may later compare cases of non-zero heat addition. The nozzle geometry considered corresponds to a wall radius that varies hyperbolically with axial distance and is symmetric about the nozzle throat. The area ratio then follows from the equation: A(x) = + (x -thr)/b2 (2. 37) in which the semi-conjugate axis, b, of the hyperbola is given in terms of the inlet area ratio, AI = A(x=0), and the distance to the throat, xthr; thus, 2 2 b =xthr /(A -1). (2.38) Figures 2-3 through 2-12 present the results for the case of linearly varying initial distributions of p, u, and T. We have arbitrarily -3 chosen a range of values from 10 to 1.0 with density and temperature decreasing from nozzle inlet to exit and velocity increasing. (Note that for y = 1.4 this corresponds to an initial supersonic value for the exit Mach number; ME = 8.44.) With xthr = 0.5, the inlet and exit area ratios are identical; they are chosen in this example to have a value of 5. 0.

61 We should mention here that themajority of graphs presented in this thesis have been produced by computer on a remote Computek 401 storage-tube graphic display terminal. Spatial distributions, such as those in Figs. 2-3 and 2-4, are generated as a sequence of straight lines connecting the actual computed values at each node and intermediate points that are added to enhance the smoothness of each curve. The intermediate values result from a third order interpolation scheme 57 due to Akima, which has been found to yield very satisfactory curve-fits to our numerical data even in regions of large spatial variations. In the above-mentioned figures, for example, there is a total of 81 points-21 nodal points and 60 interpolated points. On the other hand, all temporal distribution curves, such as those in Figs. 2-5 through 2-7, comprise a sequence of straight lines connecting only the actual calculated data points. For most quasi-one-dimensional and axisymmetric runs, the dependent variable values at every node are stored at intervals of At = 0.01. In a few cases At = 0.02. And finally, the curves shown in the "hidden-line" surface plots, such as Figs. 2-8 to 2-12, do not in general pass through the computed data (nodal) points upon which they are based, but result from an interpolation implicit in the plotting subroutine. Figure 2-3 shows the evolution of the static temperature distribution from its initial linear profile to the final "steady-state" profile at which t = 3. 0. With a nozzle length, L, of 10 centimeters,

RUN NO. 1 21 NODES, NOZ GEOM COMB: 1 - 1 NOZ LEN(M) = 0.1000, XM(THR) = 0.5000 R (INLET) 5.0000, AR(EXITI = 5.0000 1.0 o\ " k 0.6 - OU L.u CL- 0*4 p^ on TIMEf(1) = T IME (2) = TIME (3) = TIME (4) = TIME (5) = 0.0 0.2000 0.4000 0.6000 3.0000 AXIAL DISTANCE PLOT u 26 Figure 2-3. Transient Temperature Profiles for Isentropic Nozzle Flow

63 Pref= 10 atm, T = 300 K, and air as the working fluid, a unit of Pref ref time corresponds to 0. 34 msec. It should be noted that for the present example a specification of L and the reservoir conditions is only needed, as above, to refer values of the dimensionless variables to corresponding dimensional quantities. As evidenced by the working equations (2.8), the isentropic nozzle flow problem scales in such a manner that the solution depends only on the specific heat ratio, y, of the gas and the spatial derivative of the logarithm of the nozzle area, X(x). The abscissa in the figure (and all figures in this section with axial distance as the abscissa) represents a percent distance from the nozzle throat, x; hence, = (x - thr)/Xthr (2. 39) and x = 0.0 therefore corresponds to the throat of the nozzle. The numbers shown on each curve in this and other graphs serve only as identification with respect to the parameter values (t, in this case) listed in the legend above the graph and have no relation to the number or location of nodal points. Note that the first four curves correspond to the early moments of transient flow adjustment. Since little change occurs in the temperature distribution after a time of 0. 6, only the final curve is added for comparison. The figure as a whole demonstrates a gradual, generally uniform relaxation of the static temperature;to its equilibrium state. The accuracy of the numerical scheme

64 is demonstrated in the resulting values of the dependent variables which, as a whole, differ from the exact values by less than 0. 1%. The exception to this occurs in the density values in the exit region of the nozzle, which become relatively small in magnitude (e. g. p 1 = 0. 0633) and are in error by as much as 4% at the exit itself. J+I Increasing the number of axial nodes, J + 1, from 21 to 31, however, reduces this error to less than 1%. In contrast to the static temperature behavior, the adjustment of the massflow distribution to its uniform steady state value throughout the nozzle is significantly more dynamic. Figure 2-4 presents the massflow profiles for the same five values of time as in the previous figure. The peculiar shape of the initial profile is artificial and results from the initial linear variations assumed for p and u and from the symmetric nozzle area-ratio distribution. Because of the reference quantities defined previously, the dimensionless massflow plotted in the figure is non-dimensionalized with respect to the quantity m, = p, V. A,. (2. 40) mref Pref ref thr (2.40) The figure clearly demonstrates the existence of a lag between the transient massflow response in the supersonic and subsonic regions of the nozzle. This increased resistance to changes in the state of the fluid for higher Mach numbers is obviously attributable to the greater

RUN NO. 1 21 NODES, NOZ GEOM COMB: I - I NOZ LEN(M) = 0.1000. Xw(THR) = 0.5000 AR (INLET) = 5.0000, R (EXIT) = 5.0000 TIME (1) TIME (2) TIME (3) TIME (W TIME (5) =0.0 = 0.2000 = 0.4000 =0.6000 = 3.0000 1.4 1.2 a (c Cr) A: 1.0 0.8 0.6 0.4 CD 0.2 0.0 -1oi AXIAL DISTANCE PLOT # 27 Figure 2-4. Transient Massflow Profiles for Isentropic Nozzle Flow

66 58-60 momentum (or "inertia", as Zuckler refers to it) of the fluid. As one might intuitively expect, this effect exists whether there is heat addition to the flow or not. It can be detrimental to the interruption of AC circuits, because it prevents the flow from returning to its original cold-flow state near the moment of current reversal, which is the critical time for arc interruption. The fluid, therefore, retains a certain "residual temperature" that makes it less conducive to dielectric recovery. The effect of nozzle geometry on reducing this effect is discussed below. Before terminating the discussion of Fig. 2-4 we note that the massflow attains a uniform steady-state value within 1% of the exact value of 0. 684, except near the nozzle exit. Again, this is caused by the more inaccurate density values there and is alleviated if ten additional nodes are added along the x axis. Figures 2-5 through 2-7 show the temporal variation of each of the three dependent variables at five representative axial locations within the nozzle- the nozzle inlet (j = 1), mid-subsonic region (j = 6), throat (j = 11), mid-supersonic (j = 16), and exit (j = 21) regions. The figures illustrate the attainment of steady-state conditions throughout the nozzle. Although it is not entirely obvious from these graphs, the density of the fluid in the subsonic region converges toward an equilibrium condition less rapidly than the other variables and acts,

RUN NO. I 21 NODES. NOZ GEOM COMB: I - I NOZ LEN(M) = 0.1000, Xw(TTHR) = 0.5000 R (INLET) = 5.0000. Rt(EXIT) = 5.0000 JPOSN(1) = 1= JPOSN 2) = 6 JPOSN 3) = 11 JPOSN( 4) = 16 JPOSN( 5) = 21 1.2 1.0 >(f) Lu 0.8 0.6 0.4 0.2 0.0 O.i 0.5 1.0 1.5 2.0 2.5 3.0 TINE PLOT u 22 Figure 2-5. Temporal Density Variation for Isentropic Nozzle Flow

RUN NO. 1 21 NODS. NOZ GEOM COMB: 1 - I NOZ LEN(M) = 0.1000. X (THR) = 0.5000 Rw (INLET) = 5.0000, lg (EXIT) = 5.0000 JPOSN( I) JPOSN( 2) JPOSN 3) JPOSN( 4) JPOSN( 5) = 1 = 6 - 11 = 16 = 21 2.4 2.0 I —. ( —C: ---- 1.6 1.2 0.8 0. 4 0.O 0.< 0.5 1.0 1.5 2.0 2.5 3.0 TIME PLOT u 23 Figure 2-60 Temporal Velocity Variations for Isentropic Nozzle Flow

RUN NO. 1 21 NODES, NOZ GEOM COMB: 1 - 1 NOZ LEN () = 0. 1000, X (THR) = 0.5000 RN(INLET) = 5.0000. RN(EXIT) = 5.0000 1.0~... LL o.e K — 0.6 L / s C /. n y JPOSN( 1) 1 JPOSN( 2) = 6 JPOSN( 3) = 11 JPOSN( 4) = 16 JPOSN 5) = 21 VW r 0.5 1.0 1.5 2.0 2.5 3.0 TIME PLOT F 24 Figure 2-7. Temporal Temperature Variations for Isentropic Nozzle Flow

70 therefore, as a "gauge" of the overall progress toward the steady state. During runs of the computer program, values of the dependent variables at two specified nodes are monitored from the remote computer terminal. A steady state is assumed to exist whenever the -4 subsonic density value converges to within 10. This condition is admittedly arbitrary and the results of both hot (Q > 0) and cold flow runs are examined graphically to verify the achievement of equilibrium. (The program contains provisions for restarting previous runs if convergence appears to be incomplete.) Although graphs, such as those in Figs. 2-3 and 2-4, are useful for demonstrating some aspects of the temporal variations in the dependent variable distributions, Figs. 2-8 through 2-12 display essentially all of the results produced by the numerical technique and, as such, represent a very effective means for visualizing the flow dynamics in the nozzle. Figure 2-8 shows the rather undramatic evolution of the static temperature distribution from its initial linear profile to that at t = 1.0 where it is very near its numerical steadystate (as verified by Fig. 2-7). Figure 2-9, on the other hand, is significantly more interesting. It readily demonstrates the more dynamic massflow behavior in the inlet region of the nozzle ando the approach to a uniform axial distribution for greater values of time. The remaining "surface plots" in Figs. 2-10 through 2-12 show the

RUN 1 I 1.000 1.2000- 0.800^^s 0 LIJ [~ 0.800 0.600 *(a~~~~~~~~~~0.400 LLj 0.400 Zl' o.. ot Z 021- y0.200 LU 0.0 0e be be.0 X-XC THRE PLOT a 2 Figure 2-8. Spatial and Temporal Temperature Variation for Isentropic Nozzle Flow

RUN # 1 PLOT # 28 1.500 ~'\ "t O. 00 1 000- 0 0 Cr) 0.500 0.0 0.0/ 0,~,~. O,. ~ X-XC THR) Figure 2 -9. Spatial and Temporal Massflow Variation for Isentropic Nozzle Flow

RUN 1 PLOT # 30. 200 1. 000 1.000 () 0.800 1Z- ~~~~~~~~~0. 400. 0.400 ~'o ToOD > <^ oO ~-~ ~ 0.200 0.200 0.0, 0.0,,. I. C. P. /Cd' O' q~~~~~%1 X-XC THR3 Figure 2-10. Spatial and Temporal Density Variation for Isentropic Nozzle Flow

RUN a I 1.000 2,400 0- - o. 0.00 1.900 0.400 -.. 0.800 illT~~~~~~~ ~o o.200 0. 0. X-XC THRJ PLOT * 4 Figure 2-11. Spatial and Temporal Velocity Variation for Isentropic Nozzle Flow

LRUN 1 PLOT # 31 1.000 lO Or (n Q 1.2001.000 - o. Boo - 0. 8000. 6000.400 0.2000.C.O mo V VV 4; X-XC THR) Figure 2-12. Spatial and Temporal Pressure Variation for Isentropic Nozzle Flow

76 corresponding behavior of the other dependent variables, velocity and density, and also the static pressure. In the last figure, note the tendency of the pressure distribution to "over-shoot" its steadystate profile and develop a maximum upstream of the nozzle throat. The introduction of a large amount of energy addition in the nozzle makes this transient effect much more pronounced and the resulting adverse pressure gradient can cause the inlet flow to momentarily reverse itself. This form of nozzle clogging is investigated further below. As a final note on the isentropic case, we mention that the integration from t = 0. 0 to 3.0 for J = 20 involved less than twenty seconds of computer processing (CPU) time. Before considering the general case of heat addition to the nozzle flow, we compare the solution produced by the present numerical technique to another steady-state problem for which an exact analytical solution is available-the so-called "Rayleigh process" in which heat is added to a frictionless constant area duct. For this purpose, we simply attach a duct to the hyperbolic Laval nozzle of the previous example, as illustrated in Fig. 2-13. A prescribed amount of heat is then added to the duct alone, so that the nozzle acts only as a device for supplying supersonic flow to the duct inlet (x = 0. 5). The throat of the overall system is now located where x = 0. 25 and, as before, the back pressure is assumed to offer no impedance to the internal flow.

B 40 IB b L O 0.25 (= Xthr) 0 0.5 0.5 I1O 1.0 (:XE) x L -1.0 1.0 1.0 I 3.0 A J —,OF'q Isentropic Nozzle (0=o) "Royleigh Duct" (Q>0) Figure 2-13. Duct Configuration for Rayleigh Flow

78 This example problem proves to be particularly interesting, because it offers insight into several aspects of the one-dimensional numerical model. First of all, it represents another exact test with which the accuracy and overall validity of the model can be demonstrated. Besides providing a case in which heat is added to the flow, it also illustrates the importance of the boundary conditions in determining a stable solution. In addition, this problem demonstrates the ability of the finite-difference scheme to handle various types of discontinuities such as: (1) the discontinuity in dA/dx which occurs between the nozzle and constant area duct, (2) the "heat front" at which Q first assumes a non-zero value, and (3) a gasdynamic shock wave. To properly present the numerical results of this problem, it is necessary to first detail some background information pertaining to 36 Rayleigh processes. Sentman thoroughly discusses the particular case for which a Rayleigh duct is fed by a converging-diverging nozzle. Ideally, the length of the Rayleigh duct is immaterial in terms of the resulting flow pattern. Instead, the steady-state fluid properties at any value of x are determined by the flow conditions at the duct inlet (which also corresponds to the nozzle exit) and the amount of heat added to the flow between the inlet and x per unit mass of fluid, denoted here by Q. In the steady state, Q is equal to the difference in total

79 enthalpy between these two points and is related to the rate of heat addition per unit volume at each x, Q, by the following dimensionless equation: Q(x) = Q / (pu) (2.41) hf in which we take Q -P /p v (2242) ref Pref/Pref ref (2.42) xh in Eq. (2. 41) refers to the axial position of the heat front defined hf previously. We can also express the total temperature at each x in terms of Q; Tt(x)=l + (- )/y] Q(x) (2.43) with the specific heat at constant pressure, cp, referenced to the perfect gas constant. In what follows we consider the case of an initially isentropic flow throughout the entire nozzle-duct system to which increasing amounts of energy are subsequently added in the duct alone. The flow is allowed to equilibrate before the value of Q (and hence, of Q) increases to a new incremental value. In order to interpret the numerical results, it is useful to introduce the concept of the "Rayleigh line", which represents the locus of all possible steady states of a Rayleigh

80 process on a temperature entropy (T-s) diagram for given duct inlet conditions. The curve corresponding to the present example is sketched in Fig. 2-14. The Rayleigh line is determined by the steadystate continuity, momentum, and state equations and is given in dimen36 sionless terms (with s c = R) by the equation36 ref p ref y-1 s5-. ^ lIn (T/) I 2 -(/ (2.44) 7 + 1) ~ (+)2 -4y (T/T in which the symbol (A) represents that particular flow state, reached by means of a Rayleigh process, at which the Mach number is exactly unity. The energy equation determines the relative location of two points along the Rayleigh line. On the T-s diagram the initial steady state isentropic process assumed for the nozzle-duct system corresponds to the vertical line, AB, shown in Fig. 2-14. The arrows indicate the direction of flow from nozzle inlet (point A) to duct exit (point B). The static temperature is noted to decrease past its "critical" value at the nozzle throat, denoted as Tit to its final value at point B. Since no heat is added in the duct and the area is constant there, the flow conditions remain unchanged. Hence, point B initially corresponds to conditions at both the nozzle exit and all points within the Rayleigh duct. Now, if Q is

81 C) /I T 33 Royleigh Line T r - -E - -A ~ r<z Thf -- A I S Figure 2-14. Sketch of Rayleigh Line on T-s Diagram

82 increased throughout the duct to some small value, the steady-state temperature in the duct increases monotonically from point B until it reaches a certain value at the duct exit. The entire nozzle-duct flow process would now correspond to path ABC in the figure. The corresponding duct Mach number, on the other hand, decreases monotonically to an exit value, ME, which is greater than 1.0 at point C. Further increases in the amount of heat addition to the duct flow result in similar changes in the distributions of the dependent variables with the steady-state duct exit conditions advancing along the Rayleigh line until point D in Fig. 2-14 is reached and the exit Mach number becomes exactly equal to unity. Up to this point, flow conditions throughout the nozzle remain unchanged. The value of Q corresponding to sonic duct exit conditions, Qmax depends only on the Mach number at the duct inlet and follows from Eq. (2.43). Increasing Q above this critical value causes a shock wave to form at the duct exit. However, a shock is unstable in a Rayleigh duct and will tend to "jump" across the duct into the diverging portion of the nozzle. The reason a shock wave does not assume a steady-state position in a Rayleigh duct may be explained, physically, by the fact that it represents an adiabatic process and therefore does not provide a mechanism by which the flow can adjust itself to different amounts of heating. Now, with a shock wave in the nozzle and the back pressure

83 still negligibly small, the Mach number would be subsonic at the duct inlet and increase monotonically with x up to M - 1. This process would be represented by a different Rayleigh curve than the one shown in Fig. 2-14; however, we need not discuss this possibility further since the numerical model under consideration is not expected to handle the presence of shock waves in the flow. Explicit procedures to do so (derived from the Rankine-Hugoniot relations, for example) have not been incorporated into the analysis and the working equations are not written in the conservation form, which other investigators claim is necessary for the implicit handling of shocks. With the above information regarding the expected qualitative flow behavior of the nozzle-Rayleigh duct configuration, we now discuss the results of the associated numerical experiment, which are presented in Figs. 2-15 through 2-22. The starting initial conditions correspond to the actual numerical results for isentropic flow in the nozzle (obtained in the previous example) and uniform values of the dependent variables throughout the duct equal to those of the nozzle exit. Curve number one of Figs. 2-15 and 2-16 show the static temperature and Mach number distributions associated with this initial state. The second curve on each figure represents the numerical isentropic condition for the nozzle-duct system as a whole; that is, the equilibrium condition reached after the integration is

FUN NO. 29 41 NODES. NOZ GEOM COMB: 1 - 1 NOZ LEN M) = 0.200. XN(TIB) = 0.2500 R (IINLET) = 5.0000. R(EXIT) = 5.0000 TIME () = TIME (2) = TIME (3) = TIME q) = 0.0 1.0000 3.0000 5.0000 1.2 Lu Lu A? I 1.0 0.8 0.6 0.4 0.2 00 4:^ 0.0 L -1.0 AXIAL DISTANCE Figure 2-15. Static Temperature Distribution for Increasing Values of Q

RUN NO. 29 41 NODES, NOZ GEOM COMB: I - 1 NOZ LEN(M) = 0.200. XN(THR) = 0.2500 FR (INLET) = 5.0000, fI (EXI T) = 5.0000 TIME (1) TIME (2) T IME (3) TIME (q) 0.0 = 1.0000 = 3.0000 5.0000 3,5 3.0 2.5 2.0 (2: qClj 1.5 1.0 0.5 06 0 -L -1.0 1.0 1.0 15I I 2 I 1.5 2.0 2.5 3.0 -0.5 0.0 0.5 AXIAL DISTANCE PLOT m 16 Figure 2-16. Mach Number Distribution for Increasing Values of 0

86 allowed to proceed from the above starting conditions to t = 1.0 with Q equal to zero everywhere. An adjustment of the duct flow to the numerical error generated at its inlet by a discontinuity in the value of dA/dx causes curves one and two to differ. This error could be reduced by increasing the number of axial nodes used, which in turn would increase the resolution of the numerical scheme in the vicinity of the duct inlet. However, the results to be discussed below for Q > 0.0 in the duct confirm that the present value of J = 40 is sufficiently accurate for the present purposes. For simplicity, we now assume that the value of Q increases impulsively in time to the same positive value throughout the Rayleigh duct; that is, Q = constant for 1.0 < x < 3. 0. Since the term (pu) is constant in the duct (for steady-state conditions), the critical value of Q necessary to bring M to 1.0 follows directly from Eqs. (2.41) A and (2.43). Hence, with M(x = 1.0) = 3.175, we have that Tt= 1.56 and Q = 0. 564. In the steady state, this amount of heat addition max corresponds to path ABD in Fig. 2-14. In the unsteady case, however, this limit must be approached gradually if M is to remain greater than unity during the transient flow adjustments. The remaining results illustrate this point. Curves 3 and 4 in Figs. 2-15 and 2-16 correspond to the numerical steady state associated with values of Q = 0. 3 and 0. 5. The exact

87 analytical values at each node are not shown, because their difference from the numerical results would be indistinguishable. The error in the dependent variables is again generally less than 1% in both the nozzle and duct, except for the three nodes in the transition region at x = 1.0. The error here is due to both the discontinuity in wall slope and the inaccuracy in density values mentioned previously. Figure 2-17 illustrates the temporal behavior of the Mach number profile soon after the value of Q is again incremented such that Q = 0.6 throughout the duct. Since this amount of energy addition is more than that necessary to thermally choke the duct, a shock wave should form at the duct exit as discussed above and attempt to quickly traverse the duct. The numerical technique does indeed account for this phenomenon as evidenced by the unstable behavior of the Mach number profiles shown in the figure near the exit. The numerical integration terminates as the density becomes negative three nodes upstream of the exit for t = 5. 838. Figure 2-18 shows the associated variation of the exit static temperature (j = 41) with time for the entire time range from t = 0.0 until the shock occurs and causes instability. Note the achievement of steady-state conditions for the lesser values of Q, as well as the final unstable behavior when t > 5. 0. This verifies the conjecture that the present model cannot handle the presence of shock waves in the flow without producing unstable results.

RUN NO. 29 41 NODES, NOZ LEN (H ) Rw (INLET) - 6. 5s C] 4. 2. 1. 40Z GEOH COBt 1I - I - 0.2000 a XM (THR ) 0.2500 5.0000 * RN(EXIT) = 5.0000 TIME (1) - TIME(2) = TIME (3) = TIME i4) =. 5.0000 5.7000 5.8000 5.8200 4 00 CO 1 3 O..T-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 AXIAL DISTANCE PLOT a 17 Figure 2-17. Transient Mach Number Profiles for Q = 0. 6

JPOSN( 1) = L4 RUN NO. 29 41 NOOES, NOZ GEOM COMB: I - I NOZ LEN(M) = 0.2000. XM(THR) = 0.2500 Rm (NLET) = 5.0000. N(EXIT) = 5.0000 2.5 Lu Lii 2.0 1.5 1.0 0.5 0.0, 0. A _0.3, 1.0 0.5, 3.0 0.6, < t < 1.0 <t< 3.0 -t < 5.0 t 5.0 0.0 0. 1. 2. 3. 4. 5. 6. TINE PLOT a 19 Figure 2-18. Temporal Variation of Duct Exit Temperature

RUN NO. 29 41 NODES, NOZ GEOM NOZ LEN(M 0.2000 Ri (INLET) - 5.0000 1.0 0.8 0.6 c(n 0. C] TIME(1) = 10.0000 COMB: 1 - 1 X XMTHR) = 0.2500. Rm(EXIT) = 5.0000 CO 0o 0.2 1.0 1.5 2.0 2.5 3.0 AXIAL DISTANCE PLOT u 7 Figure 2-19. Density Distribution for Q = 0. 5535

RUN NO. 29 41 NODES, NOZ GEOM COMB: I - 1 NOZ LEN(M) = 0.2000, XN(THR) = 0.2500 RM (INLET) = 5.0000. R n(EXIT) = 5.0000 TIME(l) = 10.0000 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.Q 0.2 CD 0.0 L -1.0 -0.5 0.0 0.5 1.0 1.5 3.0 AXIAL DISTANCE Figure 2-20. Velocity Distribution for Q = 0. 5535 PLOT * 8

RUN NO. 29 41 NODES, NOZ GEOM COMB: I - 1 NOZ LEN(M):= O2000. X~(THR) = 0.2500 R- (INLET) = 5.0000, R(EXIT) = 5.0000 TIHE (1) = 10.0000 1.Li LJ C, lu Lu 1.2 1.0 0.8 0.6 0.q 0.2 0.0 L -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 AXIAL DISTANCE Figure 2-21. Static Temperature Distribution for Q = 0. 5535 PLOT a 9

RUN NO. 29 41 NODES. N02 GEOM COMB: 1 - 1 NOZ LEN(M) = 0.2000. X (THIB) = 0.2500 R* tINLET) = 5.0000.,RI (EXIT) = 5.0000 TIME () = 10.0000 3.5 3.0 2.5 2.0 1.5 1.0 0.5 1.0 1.5 2.0 3.0 AXIAL DISTANCE Figure 2-22. Mach Number Distribution for Q = 0. 5535 PLOT a 11

94 It was found that the value of Q could be increased gradually up to a value of 0. 5535 without causing numerical instabilities in the solution. The resulting exit Mach number is 1.0006. Thus, the model predicts the proper amount of heat addition to choke the flow to within 2%. If desired, a larger value of J would improve the accuracy. Figures 2-19 through 2-22 present the final distributions of the density, velocity, static temperature, and Mach number for the case of thermally choked flow. The slight inaccuracy at the nozzle exit (x = 1. 0), caused by the discontinuous wall slope, is more perceptible in these figures. 2.4. 2 Nozzle Flow with Intense Heat Input In order to apply the quasi-one-dimensional numerical model to the circuit-breaker problem, we now specify the form of the heat addition rate per unit volume, Q(x, t), which appears as a source term in the energy conservation equation and acts as the "driving function" for the hot flow cases to follow. In reality, there are several energy transfer mechanisms of importance during the operation of a gas-blast circuit-breaker, as mentioned in the Introduction, which could conceivably be incorporated within the expression for Q. In what follows, however, we assume that all the energy liberated within the arc plasma through Joulean dissipation remains in the gasflow. In light of the qualitative nature of the one-dimensional model and the large

95 magnitude of the Joulean heating term, wall energy transfer effects are considered negligible. We also ignore radiation losses, although they could easily be included in the model at a future date, if the assumption of an optically thin gas is made. Assuming that there are no externally applied magnetic fields and that those induced by the presence of the arc are negligible, we may write that Q(x,t) = Qref j2(xt)/U[T(xt)] (2.45) in which j is the current density (amps/m), a is the electrical conductivity of the gas (mho/m; and shown in the equation as a function of static temperature), and Q is dimensionless as before. With the further assumption that the arc is sustained by some external "current source", represented by I(t), Q takes its final form; 2 Q(x, t) = Q f 2IM-t)] (2.46) a[T(x,t)] A (x) arc with A (x) as the local arc cross-sectional area. arc In what follows we consider two forms of current input for the above "arc-like" heat addition model: 1) the "step-modulated" (DC) case with associated "free decay", in which

96 Ia (amps);0.0 < t <t max -- max I(t) = (2.47) 0.0; t> t max and 2) a simple sinusoidal (AC) current variation where I(t) Im (amps) sin (27r v t) (2.48) and a frequency, v, of 60 Hz is assumed as standard. For the DC case, therefore, the current is a step-function in time and is applied with the nozzle flow conditions initially at a steady-state corresponding to either isentropic flow (Q = 0.0) or to a previously applied lower value of I. Similarly, the AC current waveform is applied when initially max isentropic flow conditions exist in the nozzle. The assumption of an ideal current source to replace the external electrical circuit, which the arc completes, creates difficulties in the numerical technique that have no physical counterpart. The problem occurs as a result of the extremely small values of the electrical conductivity associated with relatively cold gas conditions; that is, for T less than about 3000 K. Such conditions prevail both initially, when the flow is isentropic and the current is suddenly increased, and near current-zero passage for an AC current input. Consider, for example, SF6 at a temperature of 1000 K and a pressure of 1 atm, which has

97 -28 an electrical conductivity of the order of 10- mho/m (Ref. 61). The corresponding value of Q, g iven by Eq. (2.46), would be astronomical when t = 0 in the DC case, as well as approaching current-zero in the AC case (with a small but non-zero value of current), for any reasonable value of the arc cross-sectional area. The rapid changes in the dependent variables caused by this unrealistically large heat input would, as one might expect, cause instabilities in the numerical solution procedure. Of course, in an actual gasblast circuit-breaker system the arc represents an important component of the overall circuit and its dynamic behavior couples non-linearly with that of the remainder of the circuit. Hence, the current flow depends strongly on the arc's presence, particularly near current zero conditions when the arc impedance is greatest. There exist several possibilities for circumventing these difficulties. For instance, one could somehow model the dynamic behavior of the arc geometry such that A = A (x, t) -a function of time as arc arc well as axial distance. However, such a model would of necessity be almost entirely artificial, since the physical mechanisms underlying arc dynamics in the present problem are little understood. In addition, these mechanics are clearly of a multi-dimensional nature and would probably not be warranted in the present one-dimensional analysis. Another possibility is to simply assume a constant value for a that is

98 representative of the gas as a whole. But we expect the temperature of the gas to vary over a wide range of values from a few hundred (if the reservoir temperature is taken as 300 K) to several thousand degrees Kelvin for currents up to 500 amps. Because the electrical conductivity depends so strongly on the temperature, it would not only be difficult to choose a representative fixed value, but the model would lack an important physical aspect of the gas behavior. In any case, the numerical results would still not reflect the true behavior of circuitbreaker nozzle flow under the special low temperature conditions noted above. We choose, therefore, to define an ersatz gas whose electrical conductivity varies linearly with temperature for 0 < T < 5 x 10 K, 3o such that a(T = 5 x 10 ~K) = 35. 56 mho/m. The latter value corresponds to that recently calculated by Liebermann4 for SF6 at p = 4 atm. 3o For temperature values above 5 x 10 K, the electrical conductivity is based on a curve-fit of the same data of the above reference. Since a varies with pressure to a much lesser degree than with temperature, particularly for temperatures below 2 x 104 K (see Fig. 1-2 ), a pressure dependence is not included. For the reservoir pressure of 10 atm used in the following runs, the static pressure value of 4 atm used here corresponds approximately to conditions in the majorportion of the nozzle about the throat region. Figure 2-23 shows the resulting variation of electrical conductivity with static temperature.

3000. 2500. C) 2000. t looo. 1000., 500..0... I 0. 2000.4000. 6000. 8000. 10000. TEMPERA TURE (CDEG K)J PLOT e 11 Figure 2-23. Assumed Electrical Conductivity-Static Temperature Dependence

100 For the remainder of the quasi-one-dimensional results the gas constant corresponding to SF6 and a value of y = 1. 2 are used in order to be consistent with our model for the electrical conductivity and to lend as much relevance as possible to the results. Also, unless specified otherwise, the nozzle is assumed to be 10 centimeters in length with a throat radius of one centimeter. These dimensions are generally representative of circuit-breaker nozzles currently used in industry. As an initial demonstration of the above arc-like heat addition model and the effects of an intense energy input on nozzle flow conditions, consider the same symmetric hyperbolic Laval nozzle of the previous examples to which a DC current of 100 amperes is suddenly added when t = 1.0. As usual the flow is initially isentropic and reaches a new higher temperature steady state, in this case when t = 6.0. Figures 2-24 through 2-36 demonstrate some of the pertinent aspects of the resulting flow behavior. Figure 2-24 shows the transient behavior of the massflow profile during the initial stages of flow readjustment to the intense heating. While the gas velocity initially increases in the supersonic flow region of the nozzle, the lower velocity fluid near the inlet rapidly decelerates and reverses direction just prior to t = 1.3. This "clogged" flow condition is noted to be more severe (in terms of the minimum massflow value) near a time of 1. 5. The profile actually begins to rise when

RUN NO. 31 31 NODES. NOZ LEN M) FW (INLET) l.C NOZ GEOM COMB: 1 - 1 = 0.1000, XxTIHR) = 0.5000 5.0000. F.(EXIT) = 5.0000 TIME () = TIME (2) = TIME (3) = TIME (4) = TIME (5) = 1.0000 1.1000 1.2000 1.3000 1.4000 TIME (6) = TIME (7) = 1.5000. 6000 ) () -(J LL 0.5 0.0 -0.5 i. -1.0 L -1.0 0.0 0.5 t.0 AXIAL DISTANCE PLOT * 21 Figure 2-24. Initial Transient Massflow Profiles for I = 100 amps (DC) max

102 t = 1.48 and becomes positive throughout the nozzle again as t = 1. 76. Since a unit of dimensionless time, t, now corresponds to 0. 765 milliseconds, the period of reverse flow persists for about 0. 214 milliseconds. For times greater than 1. 76 the massflow profile continues to oscillate about its new steady-state value for several cycles, but in the present case remains positive. When t = 6.0 the massflow has a value of 0. 313, less than half the numerical isentropic value of 0. 649 for this nozzle geometry. A decrease in massflow under these conditions is a familiar phenomenon in the gasdynamics of internal flows and is associated with the existence of any irreversible physical mechanism (such as heat addition, wall friction, etc.) that acts to increase the entropy of the working fluid. Therefore, the presence of the electric arc will always lead to massflow values less than those when the arc is absent, except for brief periods during the early transients of arc ignition, such as in the diverging portion of the nozzle as shown in Fig. 2-24. We investigate the possibility of minimizing this effect in a circuit-breaker nozzle below. The reason for this clogging phenomenon is best explained by considering the initial transient behavior of the axial pressure distribution. Figure 2-25 presents the pressure profiles corresponding to the first four curves in the previous figure. In response to the sudden heat input, a pressure maximum quickly develops just upstream

RUN NO. 31 31 NODES, NOZ GEOM COMB: 1 - 1 NOZ LEN[M) = 0.1000, XN(THR) = 0.5000 RM(INLET) = 5.0000, FN(EXIT) = 5.0000 1.4 l 1.2 _- - TIME (1) = TIME (2) = TIME (3) = TIME (t) = TIME (5 = 1.0000 1.1000 1.2000 1.3000 1.4000 Cu (n (or) CL 1.0 0.8 0.6 O.4 CO co 0.0 -1 AXIAL DISTANCE PLOT * 38 Figure 2-25. Initial Transient Pressure Profiles for I = 100 amps (DC) max

104 of the nozzle throat. The resultant adverse gradient in the nozzle inlet region causes the fluid velocity to decrease and eventually come to rest and reverse its direction. On the other hand, the fluid downstream of the pressure maximum now experiences a more favorable pressure gradient which resuits in a temporary increase in the massflow there. Again, note the lag that exists between the pressure response and that of the massflow. The pressure peak begins to subside when t = 1. 3, while the inlet massflow continues to decrease until t = 1, 5. For values of t greater than those considered in Fig. 2-25, the pressure distribution gradually relaxes to one of a more smooth axial variation in which the gradients are negative throughout the nozzle. We again utilize the advantages of surface plots to demonstrate the overall spatial and temporal characteristics of the nozzle flow in the present example. Figures 2-26 and 2-27 show the massflow and static pressure surfaces for a range of t from 1.0, when the current is pulsed from 0 to 100 amperes, to t = 2. 6, when the flow transients have essentially ceased. The inlet clogging effect is clearly evident, as well as the associated pressure peak upstream of the nozzle throat. The wavy pattern near the nozzle inlet in both figures is caused by the heat front located at the second node, which in turn results from the discrete nature of the numerical model. These oscillations in the axial distributions of the dependent variables occur in the vicinity of steep gradients or discontinuities whenever real gas effects are not present. In the

RUN # 31 PLOT # 19 1.000 zQ 0.500-.j L o. 0 - (f) % -0.500 oo -1. 000, I* 2.600 0 Cl un 4 I it Figure 2 <~~ t h ab~!C 1. uuu g O <;, -s, I X-XCTHR) -26. Overall Massflow Variation in I max - 100 amps (DC)

RUN, 31 PLOT # 21 2.600 1.500- \ 2.200 Lu 1.000 - 1. oo Lu 0 500 1.400 O. o.0~ o " i~1.i000 X-XC THR) Figure 2-27. Overall Pressure Variation for I = 100 amps (DC) max

107 present case, they are somewhat exaggerated by the numerical interpolation procedure of the plotting program. Note that they dissipate for larger values of time. The gasflow in the nozzle rapidly adjusts to the intense heat input by means of a series of pressure waves which are initiated by the first pressure pulse upstream of the throat and which subsequently traverse the nozzle in a downstream direction. Figure 2-28 shows the temporal response of the static temperature at five axial locations within the nozzle; x =0.0, 0. 2, 0. 5 (= xthr), 0.8, and 1.0. The relative positions of the most prominent peak of each curve indicate the downstream motion of the initial stronger wave. For greater values of Ima, this max' initial pressure pulse is more intense, as are the subsequent transient responses of the flow variables. Figure 2-29, which is inserted here somewhat out of context in order to emphasize this point, shows a portion of the static temperature surface for the case of I = 500 max amperes (DC) in a different hyperbolic nozzle geometry for which A(x = 0) = 10.0 and Xth = 0. 6. Note that time is now the horizontal axis in contrast to previous surface plots. Only the early transient period encompassing the motion of the initial pressure wave is shown. The figure clearly demonstrates the downstream motion of the wave and the increase in its velocity as it approaches the higher Mach number flow just upstream of the throat. We note further that the

RUN NO. 31 31 NODES, NOZ LEN CM) fR (INLET) = Li4.C NOZ GEOM COMB: 1 1 = 0.1000, Xx (THR) = 0.5000 5.00000. R(EtXIT) = 5.0000 JPOSN( 1) JPOSN ( 2) JPOSN C 3) JPOSN ( 4) JPOSN ( 5) 1 7 = 16 = 25 = 31 ) 3.5 L4J LUj Lu QL 3.0 2.5 2.0 1.5 1.0 0.5- - 1.0 TIME PLOT 8 Figure 2-28. Temperature vs. Time at Several Axial Stations for I = 100 amps (DC) max

16.000 — ^t oT" l^^^^^/^~ / ~ ~0.400 ) 12.000 t — c 8.oo -000o.00ooo Lu Q, -o. 2oo t 4. oo _02000, --- --— ~-\ — Lj..j. /0.400.o -0.600 TI/ME PLOT * 15 Figure 2-29. Static Temperature for I = 500 amps (DC) nmax

110 temperature "peak" is more pronounced than in the 100 ampere case and, of course, the overall temperature values and the differences between those of the transient "peaks and valleys' are greater. The irregularity of these peaks in the diverging region of the nozzle is due to poor resolution in the plotting scheme and does not reflect the exact behavior of the numerical results. Figures 2-30 through 2-36 present the steady-state "hot flow" axial distributions of each flow variable, as well as the initial isentropic distributions for comparison, for the symmetric Laval nozzle in which I = 100 amperes. Qualitatively, these curves are representative of the manner in which the present heat addition model alters the conditions throughout Laval nozzles as a whole. Figure 2-30 shows the final distribution of Q, which is seen to decrease away from the nozzle throat in both the subsonic and supersonic regions. (The heat addition rate reaches its maximum slightly upstream of the throat.) Two conditions are responsible for this behavior. Perhaps the more important effect is the arc cross-sectional area, which in this case corresponds to that of the nozzle itself and therefore is a minimum at the throat. While this alone would account for the axial variation of Q shown in the figure, the static temperature of the gas, and therefore its electrical conductivity, vary significantly through the nozzle. Figure 2-31 shows the steady-state temperature profiles for both the cold and hot flow cases.

RUN NO. 31 31 NODES, NOZ GEOM COMB: 1 - 1 NOZ LEN(M) = 0.1000, X (THR) = 0.5000 R (INLET) = 5.0000, FA(EXIT) = 5.0000 TIME(1) = 6.0000 12. 10. 8. LI lu 6. 4. I-a h"-A pr 2. -1'-I. -0.5 0.0 0.5 1.0 AXIAL DISTANCE PLOT # 16 = 100 amps. Figure 2-30. Steady-State Heat Addition Rate Distribution for I max

PN NO. 31 31 NODES, NOZ LEN (M) AR (INLET) = 4.0 3,5 TIME (1) = TIME (2) = 1.0000 6.0000 40 GEOM COMB: 1 - 1 0.1000, XN(THR) = 0.5000 5.0000, AR(EXIT) = 5.0000 s CQ 3.0 2.5 2.0 1.5 1.0 P-I ISo 0.5 o.o L -1.0.-4.0 I I I I -0.5 0.0 0.5 1.0 AXIAL DISTANCE PLOT 422 Figure 2-31. Steady-State Hot (I = 100 amps) and Cold Flow Static Temperaturmeafistributions

113 The rapid increase of the static temperature in the converging region of the nozzle tends to counteract the effect of A in increasing Q toward arc the throat, while the more uniform temperature in the diverging nozzle region causes the electrical conductivity to have a much less significant effect in this respect. The maximum heating rate of 11. 27 at the nozzle 3 3 throat corresponds to 1.49 x 10 watts/cm for the present problem parameters. Figure 2-30 also shows the steep axial gradient in Q at the nozzle inlet, which results from the heat front discussed above. The static temperature shown in Fig. 2-31 reaches a maximum value of 3. 533 (1060~K for T = 300~K in the present case) downref stream of the nozzle throat where x = 0. 36 and gradually decreases toward the exit. The steady-state total temperature, shown in Fig. 2-32, on the other hand, increases monotonically as it must for Q > 0. 0 throughout the nozzle. We note that a decrease in the static temperature, which can be beneficial in a circuit-breaker nozzle, since it entails a lower electrical conductivity, can occur in the steady state only if the Mach number of the flow increases with x. This fact is easily verified by a consideration of the static/total temperature relation for a perfect gas, Eq. (2. 24), Tt/T= 1 +[(y-1)/2]M2, t

RUN NO. 31 31 NODES. r NOZ LEN (M) Rm (INLET) = 5.5 5.0 CLL o.. 3.5 — J 3.0 I, 2.5 - 2.0 1.5 TIME (1 = 6.0000 NOZ GEOM COMB: 1 - 1 =0.1000,l X(tTHl) = 0.5000 5.00000. Am(EXIT) = 5.0000 1.0 j -1.1 -0.5 0.0 0.5 1.0 AXIAL DISTANCE PLOT * 19 Figure 2-32. Steady-State Total Temperature Distribution for I = 100 amps max

115 which implies that with Tt increasing and T decreasing with x, M must increase. Figure 2-33 compares the hot and cold flow Mach number distributions. In addition to an overall decrease in M, note from the figure that the steady-state sonic point shifts downstream of the throat when heat is added to the flow. It may be shown from a consideration of the onedimensional governing equations that this displacement of the sonic point increases with increasing Q and is also dependent on the nozzle geometry and type of gas. (See Refs. 62 and 63, for example.) For the case shown X(M = 1.0) 0. 107 and the exit Mach number decreases from 2. 785 in cold (isentropic) flow to 2. 311 when I = 100 amperes. max Figures 2-34 and 2-35 further illustrate the overall effects of heat addition on the nozzle flow. The static pressure increases above its cold flow value throughout the nozzle, while the density is noted to decrease rapidly in the subsonic region. Given values for the current, I, and electrical conductivity of the gas, a, at each x and t, we may express the one-dimensional electric field strength, E (volts/m), by means of Ohm's Law as: E(x*,t*) - j(x*,t*) (t*)(249) - a(x*, t*) - ((x*, t*) A (x*) arc The associated potential difference, V(volts), at any point in the nozzle is then

RUN NO. 31 31 NOOES NO N GEOM COMB: 1 - 1 NE LEN(M) = 0.1000, Xw(THR) = R,(INLET) = 5.0000, ft(EXIT) = TIME(1) = 0.0 TIME(2) = 6.0000 0.5000 5.0000 3.0 2.5 2.0 1.5 1.0 0.5 __5 LI Q.5 LjO -0.5 AXIAL DISTANCE PLOT * 44 Figure 2-33. Steady-State Hot (I = 100 amps) and Cold Flow Mach Number I5iaributions

RUN NO. 31 31 NODES. NOZ GEOM COMB: 1 - 1 NOZ LEN(M) = 0.1000. XN[(THR) = 0.5000 RA(INLET) = 5.0000, Rln(EXIT) = 5.0000 TIME(1) = 1.0000 TIME (2) 6.0000 1.0 0.8 lu LLJ LU Q_ 0.6 0.4 0.2 0.0 I 1. AXIAL DISTANCE PLOT a q43 Steady-State Hot (I = 100 amps) and Cold Flow Pressure Disributions Figure 2-34.

RUN NO. 31 31 NODES, NOZ GEM COMB: 1 - 1 NOZ LEN(M) = 0.1000, X (THR) = R (INLET) = 5.0000, R (EXIT) = 1.0;_ 0.8 \), TIME (1) = TIME (2) = 1.0000 6.0000 0.5000 5.0000 (r 0 0.6 0.4 I 0.2 AXIAL DISTANCE PLOT a 41 Figure 2-35. Steady-State Hot (I = 100 amps) and Cold Flow Density DisEtlSitions

TIME: (l = 6.000 RUN NO. 31 31 NOOES. NOZ GEOM COMBi 1 - I NOZ-LtN-M) 1 0.1000. X*(THR) = 0.5000 FW (INLET) = 5.0000. R (EXIT) = 5.0000 3600. 3000. Lu (.0 -...J "400. 1800. 1200. 600. 0. i -. -0.5 0.0 0.5 1.0 AXIAL DISTANCE Figure 2-36. Steady-State Potential Difference (in volts) for I max PLOT a 32 - 100 amps

120 X* V(x*,t*) = L E(x*, t*) dx*. (2.50) 0 The asterisks are included in the above two equations only to denote the dimensionless nature of the independent variables; E and V are more meaningful in their dimensional form. For the present case of a 100 ampere input, the potential difference across the nozzle, shown in Fig. 2-36, is 3.079 x 10 volts. Although it is not possible to rigorously prove the accuracy of the above and following numerical results for unsteady Laval nozzle flow with heat addition, we can demonstrate the self-consistency of the technique and thus infer its validity. First of all, we note that the results of the hot flow runs as a whole are in qualitative agreement with the very limited experimental results relevant to this problem, such as those of Kogelschatz and Schade6, who compared the axial Mach number and pressure distributions in a Laval nozzle with and without an arc present. Secondly, we compare the transient and steady-state solutions of the previous example in which different values of J + 1, the number of axial finite-difference nodes, are used. The static temperature variations at nozzle locations where the flow is subsonic and supersonic, x = 0. 2 and 0. 8, presented in Figs. 2-37 and 2-38, are representative of the overall results. It is noted from the first

J = 50 1.9 - 30 LU 1.7 1.6 - C, 1.4 - / Lu L 1.3 1 / 1.0 L 0. 1.0 2.0 3.0 4.0 5.0 6.0 PLOT u 7 Figure 2-37. Transient Temperature Response at x = 0. 2 for Different Values of J and I.. = 100 amps (DC) IImAx

QU l — i, 0; 4a.0 3.5 3.0 2.5 2.0 1.5 1.0 - 0.5 I__L 1.0 Figure 2-38. 2.0 3.0 4.O 5.0 6.0 TIME PLOT * 14 Transient Temperature Response at x = 0. 8 for Different Values of J and I = 100 amps (DC) max

123 of these figures that the difference between steady-state values in the subsonic region for J = 20 and 50 is less than 2%. Also, the greatest discrepancy between solutions occurs at the transient "peak" near t = 3. 0, but when J = 30 it is still within 2% of the higher resolution curve. From the figure, we conclude that solutions with larger numbers of nodes are better able to resolve transient oscillations at each node; however, 30 axial nodes seems to be quite adequate. In fact, the transient solutions in regions of the nozzle other than near the inlet are barely discernible from each other for this range of J. Figure 2-38 shows the three resultant temperature variations in the diverging nozzle region. Again, the lower curve corresponds to J = 20 and the upper curve to J = 50. Similar curves for the nozzle throat are indistinguishable in graphical form. The apparent decrease in accuracy of the solution in the subsonic flow region (if we assume the case with J = 50 is the more accurate) could result from the relatively poor matching of numerical and physical domains of dependence. Recall that the time step used, At, is determined in the supersonic portion of the nozzle where the term (u + a) takes its greatest value. Hence, our previous discussion of Fig. 2-1 would indicate that the domains of dependence are more closely matched near the nozzle exit. If this factor does indeed significantly affect the accuracy of the numerical technique, then a decrease in the

124 value of CSF should change the solution throughout the nozzle. To test this conjecture, we varied CSF from 0. 1 to 0. 9 in the present example with J = 30 and compared the transient response at several axial stations within the nozzle. The results were virtually identical for the major portion of the nozzle. The only variation occurred near the nozzle inlet where a better resolution of the first transient peak was obtained for greater values of CSF; the supersonic flow region was not affected. We conclude that the inaccuracy in the numerical solution that exists at lower Mach numbers is a matter of spatial, rather than temporal, resolution. We also note that compared to the isentropic cold flow case with J = 20, which required less than 200 integration steps and a final -2 At = 1. 58 x 10, the corresponding solution with the same number of axial nodes and I = 100 amperes required 850 time steps and max -2 At = 0. 613 x 1Q on the final step. For larger numbers of nodes At -2 decreases proportionately as 1/J; for example, At = 0. 245 x 10 for J = 50. It also decreases, in general, as I and, hence, the gas max static temperature increase. To keep the CPU time within reasonable limits, therefore, it is necessary to keep J as low as possible while maintaining sufficient accuracy. A value of J = 30, based on the above results and the author's experience with this numerical technique in general, was considered sufficient and used in the remaining numerical experiments discussed below.

125 We proceed now to consider the effects of nozzle geometry on the transient and steady-state flow properties in a Laval nozzle with arclike heat addition and to draw some conjectures regarding circuitbreaker design. While the overall hot flow results of the present one-dimensional model indicate that the nozzle inlet geometry has little effect on flow conditions in the diverging section of the nozzle, its influence on the intensity and duration of the nozzle clogging phenomenon is significant. Results show that clogging is minimized if the inlet area ratio is large and if the area variation throughout the inlet region is strongly convergent; that is, if lA(x) I is large. For comparison with the previous example, consider the piecewise-hyperbolic nozzle geometry, used below in a further study of geometry effects, in which the inlet area ratio is doubled [A(x=0) = 10.0] and xthr is decreased 10% to 0.4. With A = At in the diverging portion of the nozzle and I arc thr max = 100 amperes, the total energy added to the flow in this case is about 50% greater than that for the symmetric hyperbolic nozzle (based on the resulting steady-state exit total temperatures). If, for convenience, we define the duration of clogging as that period from the time at which current is initially applied to the first moment at which the massflow becomes positive again throughout the inlet, then its value is 16% less

126 for the nozzle with the shorter, more rapidly converging inlet. In addition, the massflow decreases to a minimum value at x = 0.0 of - 0. 72, compared to - 0. 92 for the symmetric nozzle. We now investigate the influence of the diverging portion of a Laval nozzle on flow characteristics by appending three basically different divergent geometries onto the more efficient hyperbolic inlet in which AI = 10.0 and thr = 0. 4: thr 1. The familiar hyperbolic nozzle in which A(x) = 1 + (x - Xthr)/b and the semi-conjugate axis is now given in terms of the exit area ratio, AE= A(x=l), such that b =(1 - th/(AE -) (2. 51) 2. A conical nozzle with a circular throat wall variation of radius of curvature r, in which A(x) = Ra + m L(x - t) / Rthr2, > xtan (2.52) where x tan represents the axial location of the point of tangency between the circular throat region (xthr < x < xtan) and the conical divergent region of half-angle x; R is the wall radius at any x, so that Rtan = R(xtan) and Rthr = R(thr); m = tan co (2.53) and L is the nozzle length.

127 3. A diverging geometry that has an elongated throat region, which may be represented analytically by the expression: A(Xr) = 1+ X-Xthr A(x) I + + (2.54) Xthr in which the exponent N may be used to determine the nozzle exit area. It should be noted that all three of the above geometries result in a continuous wall slope, dA/dx, throughout the nozzle. However, 2 2 d A/dx is discontinuous at the nozzle throat for all cases, and at xtan in the conical nozzle. tan For purposes of comparison, each nozzle has an exit area ratio, AE of 4.0. This corresponds to a = 9. 5941~ with rc = thr (= 10 m) for the conical nozzle and to N = 1. 3215 in Eq. (2. 54) for the logarithmic nozzle geometry. Each nozzle is 10 centimeters in length. Figure 2-39 shows the resulting wall radius variation for each. The conical geometry is more easily fabricated and, therefore, more commonly used in circuit-breakers and related experimental studies. There is some evidence, based on an examination of the one-dimensional steady-state equations, that suggests possible difficulties in achieving a continuous subsonic-to-supersonic flow transition in the 62 presence of an arc-like heat input for a conical nozzle. The following results tend to refute this, however. The elongated throat geometry

3.5 2: — 0:: 3.0 2.5 2.0 1.5 1.0 0.5 0.0 -1.0 )olic 0O 00 0.0 1.0 1.5 DI~STANCE FRON THROAT PLOT s 7 Figure 2-39. Three Nozzle Geometrics for Comparison Study

129 9 is more representative of that used in arc-heaters where the electric arc is constricted over much of its length so that the entire gasflow must pass through the high-temperature plasma. The object is to produce a high enthalpy gas stream for testing of materials, to simulate hypersonic reentry conditions, to produce thrust, etc. The hyperbolic diverging geometry, on the other hand, represents a compromise between these two extremes. To study both the transient and steady-state flow characteristics of the three geometries, we step the current up to 500 amperes in increments of 100 amperes, allowing the flow variables to reach a steady state at each step. The initial isentropic massflow value for each nozzle is the same, since their throat areas are identical. To enhance the validity of this comparison, we assume the arc crosssectional area is the same for all three cases. Hence, for x Xt x is equal to the nozzle throat area, 7Rthr. For x < thr, the arc fills the nozzle inlet; that is, A (x) = A(x). T he latter stipulaarc fills the nozzle inlet; that is, A (x) = A(x). The latter stipulaarc tion on the inlet arc area variation also diminishes the unrealistic effects of a severe heat front by allowing Q to increase more gradually downstream of the nozzle inlet. (See Eq. (2.46).) We first obtain an indication of the effect that the diverging portion of a Laval nozzle has on the clogging phenomenon. Using the previous definitions of the intensity and duration of clogging we find

130 that, in response to the 100 ampere step input in current, the conical geometry results in a minimum inlet massflow of - 0. 69 compared to - 0. 72 for the hyperbolic and a- 0. 77 for the elongated throat geometry. More significantly, the respective durations of clogged flow operation correspond to t = 0.49, 0. 64, and 0. 68. One may therefore conjecture that a conical divergent nozzle geometry is most beneficial in minimizing reverse flow, or clogging, effects in circuit-breaker nozzles. As the current is stepped to increasing values, the overall massflow throughout each of the three nozzles tends to decrease. Unlike the initial transient response following heat input to the isentropic flow with I = 100 amperes, the massflow profile in the succeeding cases max with larger current inputs tends to oscillate in such a manner that the inlet flow may momentarily reverse direction several times before equilibrating to a new, generally smaller, uniform value. Figure 2-40 shows the steady-state massflow for each nozzle geometry, referenced to the isentropic value, as a function of current. For relatively small amounts of heating, there is a severe drop in the massflow that each nozzle can pass. Again, the conical geometry proves to be the most efficient. However, for cases of large values of current input, which perhaps better represent the conditions in an actual circuit-breaker, the effects of nozzle geometry become less significant. Between 400 and 500 amperes, the massflow for the case of the logarithmic nozzle

131 1.0 0.8 -E \\ Log \\\-Hyperi 0.4 tcConico 0.4 - 0.2 1 0 100 200 300 I (amps) -. 500 Figure 2-40. Massflow Versus Current for Three Nozzle Geometries

132 geometry variation actually increases to a small extent. It is not clear to the author whether this behavior reflects that of the true physical problem or is associated only with inaccuracies in the numerical technique. Figures 2-41 through 2-45 present the steady-state axial distributions of several flow variables for the conical nozzle. The five curves, such as those for the energy addition rate, Q, shown in Fig. 2-41, correspond to the five values of current. Wherever applicable, the initial isentropic distribution is also included for comparison. We choose the conical geometry for this demonstration of the influence of an increasing electrical current, because the results generally indicate it is the most efficient from a circuit-breaking standpoint. The rate of heat addition to the nozzle flow has its maximum value just upstream of the nozzle throat for all values of the current. The sudden drop of Q in the supersonic flow region for I = 300 and 400 max amperes (curves 3 and 4 in Fig. 2-41) indicates the achievement of the effective "cut-off temperature" at 5000 K (T = 16. 67 for Tref = 300 K) at which the electrical conductivity increases more rapidly with T. For I = 500 amperes we note that the entire diverging portion of max the nozzle has values of static temperature that relate to values of electrical conductivity consistent with those of SF6 at P = 4 atm, as assumed in our a - T model of Fig. 2-23. Figure 2-42 shows the

RUN NO. 22 31 NODES. NOZ GEOM COMB: 1 - 2 NOZ LEN(M) = 0.1000. XM(THR) = R (INLET) = 10.0000, RK(EXIT) = 0.000 4.0000 TIME (1) TIME (2) TIME (3) TIME (4) TIME (5) = 5.9000 = 1O.7000 = 15.6000 = 20.5000 = 25.3000 (100 amps) (200 amps (300 amps) (400 amps) (500 amps) 55. 50. s5. 110. 35. 30. Q: Ljji 14J 25. 20. 15. 10. 5. -0.5 0.0 0.5 1.0 1.5 AXIAL DISTANCE Figure 2-41. Steady-State Heat Addition Rate Profiles fo'r Increasing Current-Conical Nozzle

134 corresponding static temperature distributions with the cut-off value of 16. 67 indicated by the horizontal dashed line. Also note from this figure that the temperature increases monotonically with x, unlike that in Fig. 2-31for the symmetric hyperbolic nozzle, which had a shorter, more rapidly diverging, supersonic region than the nozzles considered here. We should note that the flow behavior displayed in Figs. 2-41 and 2-42 is indicative of all three nozzles; each attains temperatures greater than T o throughout the diverging section for I 500 cut -off max amperes. The total temperature distributions for the conical nozzle are shown in Fig. 2-43 and are similar to those of preceding examples. For I = 100 and 200 amperes Tt increases in a nearly linear fashion. max t It reaches an exit value of 29. 75 for I = 500 amperes. max While the potential difference along the nozzle increases somewhat in the inlet region with increasing current, the effects of an increase in electrical conductivity cause the electrical field to decrease downstream of the throat and, hence, the overall potential difference decreases. Figure 2-44 presents the resulting voltages. The substantial decrease in Mach number for I = 100 amperes max (shown as curve number 2 in Fig. 2-45) relative to the initial cold flow distribution is further evidence of the substantial alteration of the nozzle flow characteristics caused by relatively small amounts of

RUN NO. 22 31 NODES, NOZ GEOM COMBS 1 - 2 NOZ LEN(M) = 0.1000, XN(THR) = O.41000 R t(INLET) = 10.0000. RM (EXIT) = 4.0000 24. r TIME (1) TIME (2) TIME (3) TIME (q) TIME (5) = 1.9000 - 5.9000 = 10.7000 - 15.6000 = 20.5000 TIME(6) = 25.3000 lu Lev k-. lu 20. 16. 12. 8. CL-l Cen 40 AX I AL DISTANCE PLOT a 27 Figure 2-42. Steady-State Static Temperature Profiles for Increasing Current-Conical Nozzle

RUR NO. 22 31 NOES. N GEOM COMB1 I - 2 NOZ LEN(M) = 0.1000, X-(THR) = O.AOO RK (INLET) = 10 0000. R, (EXIT) = q. 0000 TIME1) = 1.9000 TIME(2) = 5.9000 TIME(3) = 10.7000 TIME () = 15.6000 TIME(S) = 20.5000 TIME(6) = 25.3000 30. 25. Lu 20. 15. LO. 5. CAO 0. L-1.0 AXIAL DISTANCE PLOT - 31 Figure 2-43. Steady-State Total Temperature Profiles for Increasing Current-Conical Nozzle

RU NO. 22 31 NODES, NOZ GEOM COMBS I - 2 NOZ LEN(M) = 01000 X (THR) RKr (NLET) 10.0000, Rt(EXIT) 3600. 3000. Lu 2400. - F —eioo. V. 1800. >. 1200. 600. - -1} -0.5 0. 00 qL 0000 TIME(1) = 5.9000 TIME(2) = 10.7000 TIME 3) = 15.6000 TIME () = 20.5000 TIME(5) = 25.3000 -I 0.5 AXIAL DISTANCE PLOT # 30 Figure 2-44. Steady State Potential Difference Profiles for Increasing Current-Conical Nozzle

RUN NO. 22 31 NODS, NOZ GEO COMS: I - 2 NOZ LEN(H) = 0.1000, XN(TfR) = RF (IINLET) = 10.0000, AR (EXIT) = o.0000 Iq. 0000 TIME (1) TIME (2) TIME (3) TIME (4) TIME (5) - 1.9000 = 5.9000 = 10.7000 = 15.6000 = 20.5000 TIME() = 25.3000 3.0 2.5 2.0 1.5 CA3 CO 00 1.0 0.5 -0.5 0.0 0.5 1.0 1.5 AXIAL DISTANCE Figure 2-45. PLOT * 28 Steady-State Mach Number Profiles Current-Conical Nozzle for Increasing

RUN NO. 21 31 NOMES NOZ GEON COM8. I I NoZ LEN (N = 0. 1000 X THRf) = R (INLET) = 10.0000 RF EXfXIT) = 3.0 2.5 J 2.0 1.5 1.0 0 5 TIME(l) = 0.0 TIHME2) =.0000 0. q000 TIME (3) = 8000 q.0000 TIME y) = 12"0000 TIME tS) = 1S.O5000 TIME6) [= 21.0000 i. Co cO) -aoI0 Figure 2-46. -0 5 0.0 0.5 1.0 AXIAL DISTANCE.5 Steady-State Mach Number Profiles for Increasin Current-Hyperbolic Nozzle PLOT i 7

140 heat addition. Further increases in current have less effect on the flow until the cut-off temperature is reached and the Mach number increases correspondingly. The conical nozzle displays a unique flow behavior in the supersonic region for lower values of the electrical conductivity (curves 2 and 3). In this case the influence of large heat addition rates balances the expansion effects of an increasing area ratio. This results in a nearly constant value of Mach number throughout the diverging conical nozzle region. This phenomenon is not evident in the hyperbolic or logarithmic nozzle geometries. For the purpose of comparison, Fig. 2-46 presents the corresponding Mach number profiles for the hyperbolic nozzle. Figures 2-47 through 2-51 compare the steady-state distributions of the flow variables throughout the three nozzle geometries for I = 500 amperes. The existence of real gas electrical conducmax tivity values downstream of the nozzle throat lends particular significance to this high current case. The energy addition rate profiles in Fig. 2-47 display several interesting features. The influence of the heat front at the nozzle inlet is obvious. It is noted that the distribution of Q throughout most of the converging region does not change significantly for different diverging geometries. This holds true for the other flow variables as well, except for the total temperature profile of the logarithmic nozzle. The heating rate profiles of

60. I = 500 amps max 50. - -2.0 oL -1..0 log -hyperbolic - conical I,^ I' -0.5 0.00 0.5 1.0 1.5 AXIAL DISTANCE PLOT # 3 Figure 2-47. Heat Addition R Geometrics (Im vmax {ate Profiles for Three Nozzle = 500 amps)

142 the hyperbolic and conical nozzles are quite similar, while that of the logarithmic nozzle reaches a maximum further upstream of the throat and has lower values throughout its constrictor region. Downstream of the constrictor, however, the nozzle diverges rapidly, producing lower static temperatures and, thus, higher values of Q. The potential difference across the entire nozzle holds particular significance for gasblast circuit-breakers. In this respect, Fig. 2-48 shows that the conical nozzle geometry with 2068 volts is the most effective in providing a large impedance to the external circuit. The corresponding values for the hyperbolic and logarithmic nozzles are 2041 volts and 1980 volts, respectively. These results for the voltage change across the nozzle may be partially explained through a consideration of the associated static temperature profiles in Fig. 2-49. The lower temperature values for the conical nozzle throughout its inlet region and downstream of the throat result in smaller values of the electrical conductivity and correspondingly greater voltages. The relatively colder flow in the exit regions of the logarithmic and hyperbolic nozzles is not sufficient, in relation to the conical geometry, to offset the trend established by the upstream gasflow. It is interesting that the total temperature distribution throughout the logarithmic nozzle, shown in Fig. 2-50, lies significantly below

I = 500 amps max Lu rC -J (0 C^ 2000. 1600. 1200, 800.O0.O -1,0 -0.5 0.0 0.5 1.0 1.5 AXIAL DIS TA NCE PLOT (k Figure 2 -48. Potential Difference Profiles for Three Nozzle Geometries (I = 500 amps) max

I = 500 amps max 24. ( — cl.. Lu 20. 16. 12. 8. 44 O0 -0.5 0.0 0.5 AXIAL DIS TA NCE PLOT, I Figure 2-49. Static Temperature Profiles for Three Nozzle Geometries (Im = 500 amps) max

30. I = 500 amps max 25. 20. 15. 10. S. 0. - -1.0 log -hyperbolic -conical -0.5 0.0 0.5 1.0 AXIAL DISTANCE PLOT S 5 Figure 2-50. Total Temperature Geometries (I max Profiles for Three Nozzle = 500 amps)

I = 500 amps max 1.5 s 1.0 0.5 - 00 " -1w0 log -0.O 0.0 0.5 1.5 AXIAL DISTANCE PLOT a 2 Figure 2-51. Mach Number Geometries (I ma Profiles for Three Nozzle = 500 amps) 1X

147 the other two curves, while Fig. 2-49 shows a higher static temperature through its inlet and constrictor regions. This indicates that although the electric arc liberates a lesser amount of energy within the flow, the elongated throat region of the nozzle causes a greater percentage of this energy to be converted into thermal energy of the gas. On the other hand, the more rapidly diverging conical and hyperbolic nozzles allow a greater proportion to be converted into directed kinetic energy of the flow. These facts are reflected in the Mach number profiles of Fig. 2-51. Note that the conical nozzle has the lowest exit Mach number. Its advantages with regard to circuit-breaker design result from the rapid flow expansion that exists immediately downstream of the throat. As a natural continuation of the above cases of a step-modulated current input, we now consider the associated free decay from the high temperature steady-state flow, corresponding to I = 500 amperes, max to isentropic flow conditions. Figures 2-52 through 2-54 present the static temperature surfaces for the initial moments of decay in the conical, hyperbolic, and logarithmic nozzles, respectively. For each case the high current steady-state initial conditions are allowed to persist from t = 0.0 to t = 0. 5, at which time the current is instantaneously set at zero. The overall decay characteristics are noted to be most rapid in the conical nozzle, while the logarithmic geometry leads to a more gradual response.

NIU 3 37 PLOT. 17 24.000-?/ 1.500 llJ 20.000 o k.' 16.00 0.500 Lii tC 08.000 0.0 _ Lu 4.000. 0.0 -1 000,,. TIME Figure 2-52. Static Temperature Surface for Initial Moments of Free Decay in the Conical Nozzle.

RUN # 39 PLOT. 19 24.000 1.500 LUi LU Li K — CD 20.000 16.000 12.000. 000 4. 0.0 -1. 000 O L TIME Figure 2-53. Static Temperature Surface for Initial Moments of Free Decay in the Hyperbolic Nozzle.

RUNR 38 PLOT 18 24.000 -.1500 20. o000 16.000 1 2.000 Q" 8.000'u 4. o oo- ~ ~ -o.500 0.0 I' -1.000'bo`i ~ o r4 5 4 TIME Figure 2-54. Static Temperature Surface for Initial Moments of the Decay in the Logarithmic Nozzle.

151 Of particular interest with regard to the circuit-breaking potential of each nozzle geometry is the initial rate of decrease of the static temperature of the gas. Greater decay rates imply a more rapid dielectric recovery of the gas, which is of critical importance during the current zero period in an AC breaker when only a few microseconds are available for arc extinction and prevention of subsequent reignition. Figure 2-55 compares the initial temperature decay at four equally spaced locations in the diverging portion of the conical nozzle; x = 0.4 (throat), 0.6, 0.8, and 1.0 (j = 13, 19, 25, and 31, respectively, in the figure legend). The decay rate appears to be nearly uniform throughout this nozzle region. The hyperbolic nozzle displays similar decay characteristics. However, the corresponding curves for the logarithmic nozzle, shown in Fig. 2-56, indicate a less rapid decrease in static temperature within the elongated region downstream of the throat; the throat itself, however, and the rapidly expanding exit region are noted to promote a rapid flow decay. These results can be made more quantitative by introducing the concept of a conductance "time constant", or characteristic time of response to temporal changes, which is frequently used to compare 65 the transient characteristics of circuit breakers. In view of the one-dimensional nature of the present model and the assumed temperature-only dependence of the gas electrical conductivity, we may define a dimensionless time constant, r*, as:

RUN NO. 37 31 NOES, NOZ GEOM COMB: I - 2 NOZ LEN(M) = 0.1000. X,(T7ft) = RK (INLET) = 10. 00 R (ElXIT) = JPOSN( 1) = 13 JPOSN( 2) = 19 O 4OOo JPOSN( 3) = 25 040000 JIPOSN( q) = 31 q.0000N Lj - Lu Lu k 20. 16. 12. 8. 4. C-, Cll t[' 0. L 0.0 0.2 O0. 0.6 0.8 TIME PLOT # 10 Figure 2-55. Initial Static Temperature Decay at FourAxial Locations in the Conical Nozzle.

RUN NO. 38 31 NODES, I NOZ LEN tM) Rw (INLET) = 22,. 20. NOZ GEOM COMB: I - 3 = 0.1000, Xn(THR) = 0, 4000 10.0000. RK(EXIT) = 4.0002 JPOSN( 1) JPOSN( 2) JPOSN( 3) JPOSN (4) = 13 = 19 = 25 = 31 141 cL LI 18. 16. 14. 12, 10. 8. 6.'4* crA Co4 2. 0. L 0.0 O.4 0.6 0.8 1.0 TIME Figure 2-56. Initial Static Temperature Decay at Four Axial Locations in the Logarithmic Nozzle.

154 -I aT (T*) -1 (2.55) T at* t * The right hand side of Eq. (2. 55) represents a relative, or percentage, rate of decrease of the gas static temperature. It was numerically calculated for the initial moment of flow decay (t0* = 0. 5) with the use of a simple first order approximation of the derivative, aT/at*, over an interval At* = 0.01. Figure 2-57 shows the resulting values (the reciprocal of the dimensionless time constant) as a function of axial distance for each of the three nozzle geometries. The maximum percentage rate of decrease, that is, the smallest time constant, occurs in the conical nozzle downstream of the throat at f = 0. 12. Since -4 t f 7.65 x 10 sec, this maximum value, (1/' *) = 11.9, corresref max ponds to a dimensional time constant of 64. 3 gjsec. Although no experimental data directly relatable to this study are available, the above value is of the proper order of magnitude. We might also note that the present model accounts for only convective energy transfer mechanisms. Radial conduction and radiation both play additional important roles in promoting dielectric recovery in a circuit breaker arc and would further reduce the above time constant value. Note further from Fig. 2-57 that the maximum percentage rates of decay for the hyperbolic and logarithmic nozzles occur, respectively, at the throat and just upstream of the throat. (The latter case

12. r0. - / \ /\// —conlcai 8. 2. 0. -1.0 -0.5 0.0 0.5 1. 1. X-XCTHRJ PLOT I 21 Figure 2-57. Relative Rate of Decrease (Dimensionless) of Static Temperature for Three Nozzle Geometries 5

156 represents only a local maximum,) Siddons and Heron66 have in fact reported experimental evidence using streak photographs that show arc interruption in a nozzle first occurs in the throat region. Our results therefore indicate that axial convection and nozzle geometry must play dominant roles in the interruption process. The effect of the exit region of the logarithmic nozzle is again evident in the time constant behavior. However, the relatively large values of r between the throat and exit would obviously be detrimental to the effectiveness of this geometry in circuit breaking. On the other hand, the time constant remains generally low in value throughout the supersonic flow region of the other two nozzles, especially the conical one. An important conjecture that follows from these results then is that supersonic flow velocities (and hence, converging-diverging nozzles) can be beneficial in gas blast circuit breakers. In the conical nozzle, for example, the arc would tend to first become non-conducting just downstream of the throat; however, the non-conducting zone should rapidly extend downstream through the entire diverging portion of the nozzle, thereby enhancing the ability of the breaker to withstand high recovery voltages in the external circuit that tend to cause breakdown of the gas and a subsequent arc reignition. We turn now to a consideration of the sinusoidal current source, given by Eq. (2.48), as a model of an AC current input to a circuitbreaker arc. For a more meaningful correlation with previous

157 results, we use the same three nozzle geometries as above and consider a full cycle of the current-time variation for which I max = 500 amperes. With a frequency, v, of 60 Hz, current zero occurs for (dimensionless) t = 0.0, 10. 9, and 21. 8. And I = 500 amperes max occurs for t = 5.45 and 16. 35. Initial conditions are again those of isentropic flow. Before presenting specific results of these numerical experiments, we can mention some general comments regarding clogging effects. The temporal increase in current is sufficiently rapid in these cases that the clogging phenomenon is very similar to that occurring for the step-modulated (DC) current input. The mechanisms leading to clogging are the same in both cases and they occur, as expected, immediately following current-zero (t = 0. 0 and 10. 9, in the present cases) as the magnitude of the current, and hence the heat addition to the gas flow, begins increasing. The conical nozzle again proves to be most efficient in minimizing this effect, while the logarithmic geometry leads to more intense clogging. Since the conical nozzle geometry displays the greatest overall potential for a gas blast circuit breaker of the three geometries considered in this study, we again use it as the prime example in presenting details of the results for an AC current input. Results corresponding to the hyperbolic and logarithmic nozzles are in general quite

158 similar. The more significant differences that exist in these two cases compared to the conical nozzle, are mentioned below. Figures 2-58 through 2-69 show the full-cycle temporal behavior of the density, velocity, static temperature, and heat addition rate at an inlet location (j = 6), the throat (j = 13), and the exit (j = 31) of the conical nozzle. The horizontal dashed line, included in each of these figures, designates the steady-state value of each variable, based on the previous DC current cases for I = 500 amperes. max As a whole the flow conditions within the nozzle follow the cyclic variation of the electric current input, which results in repeatable temporal patterns of the flow variables once the initial transients have decayed. The initial response differs from the overall temporal behavior because the flow is at its steady-state isentropic condition before the current increases; whereas succeeding current-zeros produce a different, higher-temperature flow state. As noted from the figures, there exists a time lag, Atlag, between the response of each dependent variable and the electric current input. If we measure this lag relative to the current-zero at t = 10. 9, the following approximate values are obtained:

JP'OSN( 1) - 6 RUN NO. 41 31 NOOES. NOZ GEOM COM1Brt -2 HNO LEN(4O = 0.1000. XM(TIHR) = 0.L000 t:INLET) = 10.0000, Fw(EXIT) =.0000 1..0 08 Cf) 0,7 0t6 0.5 04q cn CD 0,3 0.2 = 0. 1841 0. U, 8. 12,,6 2~ TIME Figure 2-58. OT Cycle Density Response in the ConicalNozz (j 6) Figure 2-58. Full Cycle Density Response in the Conical Nozzle. (j = 6)

JJOSN( 1) = 6 RUN NO. 41 31 NODES. NOZ GOM COMB 1I - 2 NOZ LEN(M) = 0.1 000, X~ (THR) = 0. 000 RA (INLET) =. 0000. R (EXIT) = 4. 0000 1 — C —~ L.U. 0.25 0.20 0.15 0.10 0.05 u =- 0. 2431 _ r ss (3 0.0 -0.05 0. TINE PLOT q45 Figure 2-59. Full-Cycle Velocity Response in the Conical Nozzle (j = 6).

JPOSN( 1) = 6 RUN NO. 41 31 NODES. NOZ GEOM COMB 1 - 2 NOZ LEN () 0. 1000, X (THR) = 0.4000 Rw (INLET) = 10.0000. Fl (EXIT) = 4.0000 -T = 6.459 7. b Lli 5. 4. 3. 2. I.r ~ —A 1. 0. 0. 8. 12. 16, 20. 24 TINE PLOT # 46 Figure 2-60. Full Cycle Static Temperature Response in the Conical Nozzle (j = 6)

JPOSN( 1) = 6 RUN NO. 41 31 NODES, NOZ GEOM COMBs 1 - 2 NOZ LEN(M) = 0 1000. Xm THR) = 0.4000 R* (INLET) = 10.000 R (EXIT) = 4.0000 9. 8. 7, 6. -sQ = 8.431 / ss - _ _ _ >LJ Lu 5. 4. 3. 2. 1. O. 0, TINE PLOT a 48 Figure 2-61. Full Cycle Rate of Energy Addition Response in the Conical Nozzle (j = 6)

RUN NO. t1 31 NOOES, N02 GEOM COMB: 1 - 2 NOZ LEN () = 0.1000, Xm (THR) = 0. 4000 f: (INLET) = 10.0000, I (EXIT) = 4.0000.JOSN( 1) = 13 o steady-state values corresponding to I = 100, 200, 300, and 400 amperes max 0.7 0,6 c) 0.5 0.4 0.3 0.2 COA = 0.0472 0.1 0.0 0. 4. 8. 16. TIME PLOT m 49 Full Cycle Density Response in the Conical Nozzl; (j = 13) Figure 2-62.

JPOSN( 1) = 13 RUN NO. 41 31 NODES, NOZ GEOM COMt I - 2 NOZ LEN(M) - O1000, X* (THR) = 04000 A (INLET) = 10.000, RK (EXIT) = 4.0000 4,0O -.0 o-,_ 2.5 -25 o steady-state values corresponding to I = 100 200, 300, and 400 amperes max c0 1.5 1,0 1 0. TINE PLOT a 50 Figure 2-63. Full Cycle Velocity Response in the Conical Nozzle (j = 13)

UN NO. 41 31 NODES, NOZ GEOM COMB: 1 - 2 NOZ LEN () = 0.1000. X(TIHR) = O. 000 Fx (INLET) = 10.0000. R (EXIT) = 4.0000 JPOSN( 1) = 13 o steady-state values corresponding to I = 100, 200, 300, and 400 amperes max 18. 16. -T = 16. 234 ss Lu Lu Lu k —.. 14. 12. 10. 8. 6. 4. 2. 0. I,0 TIME PLOT # 51 Figure 2-64. Full Cycle Static Temperature Response in the Conical Nozzle (j = 13)

JFOSN( 1) - 13 RUN NO. 41 31 NOtES, N02 GEOM COSBt I - 2 NOZ LEN(M) 0.1000, XN "(THR) 0,.4000 Rf (INLET) - 10.0000. F (EXIT) = 4.0000 60. 50. - - \ __ ~o. -Q =55. 36 ss >C3 Lu la 30, 20. 10. 0, 0. TIME PLOT * 53 Addition Response in the Conical Nozzle (j = 13) Figure 2-65. Full Cycle Rate of Energy

JOSN ( 1) = 31 RJN NO. 41 31 NOES, NOZ GEOMt COM8 1I - 2 NOZ LEN ) = 0.1000. Xw (THR) = 0.4000 R t(INLET) = 10.0000. R (EXIT) = 4.0000 0.06 0.07 0.06 k-. (f) Lu 0.05 0.l04 0.03 0.02 0,01 PSS = 0. 00472 0.0 0. 4. 8. 12. 16. 20. 24. TINE PLOT 59 Figure 2-66. Full Cycle Density Response in the Conical Nozzle'j = 31)

.JPOSN( 1) = 31 BUN NO. 41 31 NODES. NOZ GEOM COMB 1I - 2 NOZ LEN () 0.1 000. Xm(THR) = 0.4000 FR (INLET) - 10.0000. P (EXIlT) = 4.0000 = 9. 444 109 (2 t —-- 8. 7. 6. 5..4. 3. 2. O. 4. 12. 16. 20O 24. TIME PLOT # 60 Figure 2-67. Full Cycle Velocity Response in the Conical Nozzle (j = 31)

JOSN( 1) = 31 RUN NO. 41 31 NODES NOZ GEWO COMBt 1 - 2 NOZ LENQI) = 0.1000. X(tThIR) = 0.4000 tFINNLET3 -10.0000 AR(EXIT) = 4.0000 bi Lu 20. 16. 12. 8. 4. I' CD 0. 8. 12. 16. 20. 24. TIME PLOT # 61 Figure 2-68. Full Cycle Static Temperature Response in the CoiAical Nozzle (j = 31)

JPOSN( 1) = 31 BUN NO. 41 31 NODES. * NOZ LEN M) - RFi (INLET) = 12. 10. ~> _ 6. o 6 VOZ GEOM CGOMS I - 2 - 0.1000. Xi(THR) = 0.4000 10.0000. f (EXIT) =-.0000 Q =3. 835'(QSS 2. 0. 0. 4. 8. 12. 16. 20. TINE PLOT s 63 Figure 2-69. Full Cycle Rate of Energy Addition Response in the Conical Nozzle (j = 31)

171 Table I. Time Lag of Response of Dependent Variables j-6 j-13 j-31 p 0.56 0.365 0.65 u 1.10 0.68 0.55 T 0.36 0.24 0.0 In general, Atlag decreases with increasing Mach number, except in the case of density near the nozzle exit. This apparent discrepancy probably results from the numerical inaccuracy associated with the relatively low magnitude density values, as mentioned previously. We note, also, that the gas static temperature, and thermodynamic properties in general, respond most rapidly to the intense heat addition; whereas changes in the directed velocity of the gas at any nozzle location "wait" for the necessary pressure gradients to be established. Figures 2-62 through 2-64 contain eight additional data points that help to illustrate the "near quasi-steady" behavior of the nozzle gasflow when it is subjected to a 60 Hz sinusoidal current input. The points represent steady-state values of the dependent variables at the nozzle throat, based on the previous numerical results for the DC case. They correspond to values of I = 100, 200, 300, and 400 am s ad ae max amperes and are shown plotted at times corresponding to like values of the time-varying AC current input. The two intervals following

172 current zeros at t = 0. 0 and 10. 9 are shown. The static temperature values for the AC case fall closest to those of the steady state, while the fluid velocity displays the greatest discrepancy in this respect. A similar comparison of density values shows good agreement except for the case of I = 100 amperes (t = 0.7 and 11.59); that is, except max near current zero when the time rate of change of the current is greatest. As a whole, then, the results lend credence to the quasi-steady assumption that is often used in analytical studies of the near peakcurrent characteristics of AC breakers. We should mention, however, that as the frequency, v, increases above the present value of 60 Hz the temporal behaviors of the fluid properties are expected to lose the "rectified sine wave" appearance, displayed in Figs. 2-58 through 2-69, as dI/dt becomes larger overall and the gasflow is less able 72 to follow the driving function, I(t) The unusual energy addition rate variation at the nozzle exit, shown in Fig. 2-69, is typical of the diverging region of all three geometries. The sudden decrease in Q as current increases results again from the effective cut-off temperature present in the electrical conductivity/temperature model used in this study. This phenomenon may indeed be relevant to actual circuit breaker flow behavior with gases that display similar changes in a at given temperature values. 58-60 Zuckler discusses the significance of flow-inertia effects on circuit breaker performance. For very large AC current inputs

173 the massflow through the breaker nozzle can beconme very small. Then, as current-zero is approached, it rapidly increases but is unable to regain its full isentropic flow value due to the inertia of the high velocity gas stream. As a consequence, the electric arc cools less efficiently and dielectric recovery is impeded. One measure, therefore, of a given nozzle geometry performance is the massflow decrement, Amde, which may be defined as: Mdec misent - rm (t = 10.9 + Atl ) (2.56) dec isentr lag where rf (t = 10.9 + Atg) is the minimum massflow value in the vicinity of the half-cycle current zero. Of course, the massflow does not remain uniform throughout the nozzle durirg the transient flow variations, especially those near current-zero; hence, its temporal behavior will be different at different nozzle locations. However, it is logical to focus our attention on the nozzle throat region, since we have found that arc extinction is most likely to first occur there and it is generally representative of the overall nozzle flow response. Figures 2-70 through 2-72 show the temporal massflow variation at the throat of the conical, hyperbolic, and logarithmic nozzles, respectively. The massflow values are now referenced to that of the initial isentropic flow, mi., so that Amh follows immediately as shown isentrc' o edec in each figure. The conical nozzle geometry would again seem to be

Am /r. = 0. 145 dec isentr F —11 III 1 1 1 1 TIME LU 0.6 St___at__ tat POT2 0.2 ou __ 0.'. 8. 12.. 260. -TIME PLUOT.r 23 Figure 2-70. Full Cycle Massflow Response in the Conical No te~(Throat)

1.0 QC -) 0.8 0.6 0.4 0O2 Am d /m =0.152 dec isentr cn 01 0o0 0. Ii. 8. 12. 16. TINE PLOT a 22 Figure 2-71. Full Cycle Massflow Response in the Hyperbolic Nozzle (Throat)

1.0 etc 01 O.W8 O.6 01ff 0d(2 Arf /Mi. 0. 161 dec isentr Ac~ - -. __ r~~~~~~~~ ~. 4^^ ====^ _p^^:= =rr_~~~~~~~~~.. 0o0 0. 8. 12. 16. 20. TIME PLOT s 2N Figure 2-72. Full Cycle Massflow Response in the Logarithmic Nozzle (Throat)

177 most beneficial from a circuit breaking standpoint. Its massflow decrement has a value of 0. 145 compared to 0. 152 and 0. 161 for the hyperbolic and logarithmic geometries. That is, the conical nozzle geometry is most effective in overcoming the flow inertia and returning the fluid conditions nearer to those of true isentropic flow during current-zero passage. One of the most prominent features of the massflow response in Figs. 2.71 and 2-72 for the hyperbolic and logarithmic nozzles is the intense transient wave motions caused by the initial moments of heat input at t = 0.0 and 10.9 and, as mentioned previously, directly associated with nozzle clogging. This feature is noted to be much less pronounced in Fig. 2-70 for the conical nozzle. As a final note to the one-dimensional results, we mention that each of the above cases involving a full cycle AC current input required less than six minutes of CPU time on the IBM 360/67. As such, this represents the "worst case" of all the hot flow runs discussed.

3. AXISYMMETRIC MODEL 3.1 Introductory Remarks In the present section we consider the numerical solution of a multi-dimensional model of the circuit breaker nozzle problem. Although the previous quasi-one dimensional approach offers a great deal of insight into the gasdynamic characteristics involved, the basic physical problem at hand remains inherently three-dimensional. At the outset of this study the primary goal of the author was to develop a computer model, based on the assumption of symmetry about the nozzle centerline,that could predict the transient response of a coaxial "arc", or high-temperature gas core, to an electric current input. The inclusion in the model of real gas effects, such as thermal and electrical conductivity, was also considered desirable. The obvious advantages over the one-dimensional case of the additional radial coordinate would include the ability to define and monitor an arc geom - etry and to investigate the significance of radial variations in fluid properties on breaker performance. However, as this study progressed it became apparent that the practical and economic limitations imposed by the nature of the problem (e. g., the high static temperature values and associated decreased time-step of integration) were substantial. As a consequence, the 178

179 axisymmetric model presented here serves more as a study of the feasibility of developing a detailed, multi-dimensional, numerical model than as a tool to predict gas blast circuit breaker performance. Nevertheless, the model is significant because it can accurately predict the unsteady (or steady, in the limit of large t), rotational, inviscid, flow within a general Laval nozzle geometry that is subjected to a bulk heat addition. Furthermore, the numerical technique represents a substantial increase in efficiency (e. g., less CPU time and simpler program logic) over previous schemes presented in the literature, while maintaining second-order accuracy in time and space. The additional numerical implications associated with two spatial dimensions and time as independent variables turn out to be the overriding consideration in the development and extent of the axisymmetric model. One of the more important aspects of the problem concerns the substantially increased CPU time. Apart from other factors to be mentioned below, the computer time depends approximately on the number of spatial dimensions, m, according to: (CPU) (No. nodes per dimension)m Thus, with say 20 finite-difference nodes for each dimension, the total number of nodes at which the governing equations must be solved increases at least by a factor of twenty over the one-dimensional case. Other factors that increase the CPU time include the need for

180 solving the additional radial momentum equation and the extra terms that account for radial variations in each of the other conservation equations (i. e. radial convection, pressure gradient, etc.). Also, the need to apply conditions at two additional boundaries, the nozzle axis of symmetry and wall, requires the equations to be transformed and rewritten in new forms to ensure numerical stability. This increases the number of terms in each equation and, hence, the computer processing time. Finally, it is found that the criterion for controlling the time step of integration must be slightly more restrictive than that of the one-dimensional model. Few examples of the numerical solution of the unsteady, compressible, conservation equations for multi-dimensional internal fluid 68 flow exist in the literature. Kurzrock performed a very comprehensive numerical study of the compressible Navier-Stokes equations (including effects of viscosity and thermal-conduction) when he treated the two-dimensional shock tube problem. He reported computation times of 6. 6 to 24. 5 hours on an IBM 7044 digital computer. Numerical investigations concerning unsteady inviscid flow in a Laval nozzle have generally been in the form of time-dependent approaches to the 69 steady-state transonic flow problem. The works of Saunders Migdal t 50 27 fal70 Migdal et al., and Serra fall into this category, while Armitage

181 considered steady-state axisymmetric swirl flow. All of these studies used a variation of the full L-W technique with only Migdal solving the non-conservation form of the governing equations. Serra's study is the most recent, and perhaps the most complete in the area of internal gasflows. And yet, for an axisymmetric conical nozzle and adiabatic flow he reports a CPU time of 80 minutes on a UNIVAC 1108, which is comparable in speed to the IBM 360/67 used in the present study. (Results of the present technique for the same nozzle prove to be comparable to Serra's with a substantial reduction in computer time.) It is, therefore, apparent that any numerical study of the effects of unsteady heat addition, which involves an increase in the static temperature of the working fluid and an attendant decrease in the permissible time step of the finite-difference scheme (recall Eq. (2. 35)), must be limited in nature. Thus, in what follows, we exclude the influence of diffusive fluid mechanisms and concentrate instead on the effects of forced convection and nozzle geometry. The effects of viscosity and thermal conductivity are ignored because their inclusion would necessitate smaller time steps and finer computational grids in order to resolve the boundary layer effects at the nozzle wall and the arc boundary.

182 3. 2 The Axisymmetric Governing Equations Introducing the radial velocity component, v, we write the unsteady axisymmetric Euler equations in the following dimensional form: Continuity Pt - - VPr - UPX - p (div V) (3. la) Radial Momentum Vt = - v -uv -pr/P Vt r x r" (3. lb) Axial Momentum ut =-vur uu pX/P (3. 1c) Energy T =-vT - uT - P (divV) + Q/(pc ) r x (3. Id) in which P a [T(y - 1)/(ycv)] div V = u + (rv) /r x r (3.2) and (3.3) For a perfect gas, P in Eq. (3. 2) takes the form: P = (y - 1) T. (3. 4)

183 If the initial flow conditions assumed to exist throughout the nozzle are symmetric about the nozzle centerline (the x axis), then all subsequent transient flow disturbances will preserve axial symmetry. Hence, at r = 0 we have P =u T =p =v=0. (3. 5) r r r r After noting from l'Hopital's rule that at r = 0 div V = u +2v x r we write the equations governing flow along the axis as Pt = - UP - p (div V ) t= - uu -p/p (3. 6) T= - uT - P(div ) +Q/ (pcv) t X v 3. 3 Numerical Technique All of the ideas incorporated into the quasi-one-dimensional numerical model apply to the axisymmetric case, with some extensions at the nozzle inlet that are necessary to account for the fourth dependent variable, v. The numerical simulation of the additional boundaries at the axis and wall proves to be formidable and constitutes an integral component of the overall technique.

184 3. 3. 1 Nozzle Interior and Wall Points The application of any finite-difference scheme to the spatial domain defined by the nozzle inlet and exit planes, axis of symmetry, and wall results in non-uniform spacing of the differencing mesh at the nozzle wall, R(x). This requires the use of a complicated interpolation scheme to compute dependent variable values at wall nodes and, in turn, leads to smaller time steps (since At is generally determined by the smallest spatial increment) and, possibly, to numerical instability. Furthermore, any such interpolation is unique to a given nozzle geometry. To circumvent these difficulties, we map the above physical domain onto a rectangular computational domain in which the mesh spacing can be defined to be uniform in each spatial dimension. This transformation is effected through the definition of a dimensionless radial coordinate, r* = r/R(x). (3. 7) Noting that for any dependent variable, g, g ag/R(x) (3. 8) ar ar* and g=[g,r* R ag1/L (3.9) ax [x* a r*J /

185 where /dR d R - dx)R - d* (In ), (3 and using the same reference quantities as in the one-dimensional model, we re-write Eqs. (3. 1) in the following dimensionless vector form (again omitting the asterisk superscripts for convenience): ft = f -uf -. (3.11) t r x Again f represents the vector of dependent variables, P f v; (3.12) u T S is given by p (div V) S = ~ (3. 13) (px -rRpr) P (div V) - (y - 1) Q/p where:(x) L/R(x) (3.14) and now div V u - rRu + 3(rv) /r. (3.15) x r r/r

186 And finally, V represents a scalar variable; V - fv-rRu. (3.16) The underlined terms in Eqs. (3.13), (3.15), and (3.16), which result from the above mapping, represent an increase of about 35% in the total number of terms to be evaluated at each time-step of the numerical integration. The price to be paid for this convenient transformation, therefore, is a corresponding increase of equal proportion in the CPU time. We now divide the resulting computational domain, 0 < r < 1 and 0 < x < 1, into uniform radial increments, Ar= 1/I, (3.-17) and uniform axial increments, Ax =/J, (3. 18) with I and J as the total number of divisions in each direction, and seek values of f.n f(r.,x. t) (3. 19) ij i n that approximate the solution of Eqs. (3. 1). For the above we have r. = (i - 1) Ar (i = 1,2,..., I, I + 1), x = (j - 1) Ax (j = 1,2,... J, J + 1), and t n A t (n = 0, 1, 2,... ). Figure 3-1 shows the spatial portion of the finite-difference mesh in the rectangular spatial portion of the finite-difference mesh in the rectangular

187 // Nozzle Wall ( r = R(x)) 1.0 Inlet Plane (X=O) S, 1 4 4 _ _ -_ Exit Plane (x=L) / 0 0 ~Axis of Symmetry \A (r=O) 1.0 = r/R(x) r R(o) -! \;~~~ \** R (x) I -Lr_ * Constont _ _ e 0 0 L Figure 3-1. Illustration of Coordinate Transformation

188 computational domain and the corresponding mesh lines in physical space. Note that the transformation given by Eq. (3. 7) results in greater radial resolution of the numerical scheme in the nozzle throat region. It also has the obvious, yet notable property of mapping the nozzle centerline into itself, thus preserving the symmetry conditions that hold there. We can now formulate the difference equations for the nozzle interior, 0 < r < 1 and 0 < x < 1, by applying MacCormack's scheme to Eqs. (3. 11) in a manner directly analogous to that of the one-dimensional case. The predictor step is based on forward differences in all independent variables and follows from ij f ij -( i i+l, j i j-( int j+1- i() - A t(S')n (3. 20) n ij with A = Ant/Ar,(3. 21) and X= At/Ax. n S' denotes the vector S after the partial derivative terms of its components have been differenced according to MacCormack's scheme. Similarly, the corrector step follows with backward differences for the spatial derivatives:

189 ~n+1. nI. fn+-1 n~ n+~ f n+l _ (1/2) f.n + fiJ fiJ - f] (I i U 1] J ij - i-1,j n+f n+I n+1 n+l - (Xu -)(f f.i) - A t (S')ij (3. 22) - i] " l] 1,3-1 n''i]j With the simplified notation defined in the previous section f.. - f, fn. - f, etc.), the difference equations for the domain 13 ij interior, 2 < i < I and 2 < j < J, may now be formulated. For the predictor step, Eq. (3. 20) yields: P= P ( - (U((Pi - P) - ( - p) - p At (div V)' V = v-(( - (( - v) - (Xu)(vj+ - v) - (Pi+- P)/P u= u- ()(u -) (Xu) (u - u) (3. 23) [X(pj+l - p) - (i - ) AtR(pi+I - )] /p T = T- (<)(Ti++ - T) - (Xu)(Tj+ - T) - P At (div V)' + At (y - 1) Q/p in which At (div V)' - X(uj+l - u) - (i - 1) At R(ui - u) + i[() V. - (3. 24)

190 For the corrector step, Eq. (3. 22) yields: n+1 pn _ (1/2) p + p- ( T)(P- P ) - (Xu)( - ) - p t (div V)f] v " = (1/2) [v + v - ()R(v- -vi1) - (Xu)(v - j 1) - P Pi-)/p] n+1 un = (1/2) u + - - (. )(- i - (Xu)( -1) (3.25) 1-. J-[fx - p i) - (i - 1) At R(p - pi)L)] /P } - P At (div V)' + At (y - 1) Q/p in which At (div V)'- x (u -.1_1) - (i - 1) At R(u - il) J+ Av i- -)il] * v(3.26) The above sequence of forward-forward and backward-backward (ff/bb) spatial differences in the predictor and corrector steps, respectively, corresponds to the original and most common pattern of MacCormack's scheme. However, the other three possible sequences (fb/bf, bf/fb, and bb/ff all retain second-order accuracy in time and space. In fact, all four sequences can be used cyclically

191 and may produce different stability characteristics than the scheme used here. The interested reader should consult Refs. 37, 46, and 47 for a more detailed discussion of these ideas. We have chosen the most tested numerical scheme for the interior node calculation. Equations (3. 23) and (3. 25) do not, in general, apply to nodes along the nozzle wall where i = I + 1, mainly because application of the forward difference in the radial direction requires that the dependent variable values "inside" the non-permeable wall boundary be known. One way of obtaining such values is the reflection rule mentioned in the general discussion of numerical boundary conditions in section 2. 3. 2. We refer the reader to that section, rather than repeating here the reasons for our rejecting this physically and mathematically unappealing numerical scheme for handling solid boundaries. Another reason that the interior difference equations do not apply at the wall is because the velocity components there are coupled by the standard tangency condition for inviscid flow, i. e. v* = u* (dR/dx*)/L. (3. 27) Any mathematically proper boundary condition must be ba sed on Eq. (3. 27) and the difference equations must use information obtained only from the interior domain to establish flow conditions at the wall. Rather than employing the notion of characteristics, which represents one possible but more complicated approach, we solve the continuity,

192 axial momentum, and energy equations at the nodes where i = I - 1 using second-order backward differences for the radial derivatives in both the predictor and corrector steps, thus preserving the overall second-order accuracy of MacCormack's scheme. For any variable g (including "predicted" values, g), we have then that gr - (3gI+,j - 4gI,j +gI-11 )/(2Ar) (3.28) The difference equations that hold along the nozzle wall (i = I + 1, 2 < j < J) then follow directly from Eqs. (3-23) through (3. 26), except that Eq. (3. 27) replaces the radial momentum equation and relation (3. 28) is used to evaluate all radial derivatives. 3. 3. 2 Points on the Axis of Symmetry The finite-difference equations for nodes that lie along the nozzle centerline (i = 1; 2 < j < J) follow from Eqs. (3. 6) in a manner directly analogous to those of the nozzle interior and are therefore not detailed here. Note that the radial coupling of the centerline and interior difference equations occurs only through the term 23 v in the divergence of the velocity vector. One obtains the r associated backward radial difference needed in the corrector step by imagining a row of "virtual" nodes at r = - Ar and applying symmetry conditions about r = 0. Hence, the radial velocity component at a virtual node at any j is simply equal to - v2. 2?

193 3. 3. 3 Nozzle Inlet and Exit Points Calculation of the dependent variables at the nozzle inlet plane (j= 1, 1 < i < I + 1) proceeds in the same manner as for the onedimensional model. The inlet fluid properties, p, T, and p again follow from the assumption that any heat addition to the nozzle flow begins just downstream of the inlet plane. As a consequence, the incoming flow is adiabatic and the total conditions remain fixed at their dimensionless reservoir values of unity (since reference values corresto reservoir conditions). The inlet static temperature and density are then determined from the relations: T=1 ( 1) (2 + v2 (3.29) and P + (u2 + v)/ (3.30) once values for u and v have been established, as detailed below. The static pressure is obtained from the perfect gas law. We determine u from the axial momentum equation as in the onedimensional case. In the predictor step, which involves forward spatial differences, the resulting difference equation is identical to the axial momentum equations for the downstream nodes (j > 2). That is, one uses the second of Eqs. (3.6 ) for the nozzle centerline-inlet node

194 (i= 1, j = 1) and the third of Eqs. (3.23) for j = 1 and 2 < i < I. The nozzle wall-inlet node (i = I + 1, j = 1) calculation also follows from (3. 23) with second-order backward differences in place of the two radial derivatives. To carry out the backward spatial differences in the corrector axial momentum equation, we define a column of "virtual" nodes located upstream of the inlet plane at x = - Ax. The dependent variables (in particular, u and p) at these nodes are determined by a second-order extrapolation from the three adjacent downstream nodes, (i, 1), (i, 2), and (i, 3); hence, if these variables are denoted as f. 1, O it follows that fO 3(fi fi2) + f3 (3. 31) for 1 < i < I + 1. Then the second of Eqs. (3. 6 ) becomes (in abbreviated notation) [u + U ( -) - ( -o /] (3. 32) for i = 1, and the third of (3. 25) becomes u = u+u -()(U -i) - (x)(U - - [(p - o (i 1) t R(p - p i-)] / (3.33) for 2 < i < I. Again, for i = I + 1, one changes the backward radial differences in Eq. (3.33) from first to second order.

195 With values of u obtained from the above scheme, the specification of the inlet boundary conditions becomes complete through a determination of the radial velocity component, v. We accomplish this, in a 27 manner similar to that of Serra, by simply fixing the inlet flow angle, o., for 2 < i < I; then v. = u..tan c. (3.34) vti, 1 1 Since we know that a = 0.0 and OtI+1= tan -( dR 1 I+1 L dx* it is most convenient to assume that a.i varies linearly from r* = 0 to r* = 1.0; that is, Ui = (~I+1 (i - 1) Ar* (3.35) for 1 < i < I + 1. Although this assumption is somewhat arbitrary, numerical experiments indicate that a linear inlet flow angle variation closely approximates the true flow pattern in the nozzle geometries of interest here. If, for example, one assumes a parabolic radial variation of a. and computes the steady-state numerical solution for adiabatic flow in a given nozzle, then the flow angle variation for the remaining portion of the inlet region (beginning two nodes downstream of the inlet thplane) proves to be nearly linear. This result would further suggest that the

196 flow quickly adjusts itself to a pattern consistent with any nozzle inlet geometry and that ca does not significantly affect flow conditions downstream of the inlet plane. We therefore conclude that the assumption of a linear inlet flow angle variation is quite satisfactory. Finally, the exit plane conditions are again determined by means of a second-order extrapolation of the three adjacent upstream nodes (i, J-2), (i, J-1), and (i, J); thus, i, J+l i, J-2 il - fiJn+ for 1 < i < I + 1. As in the one-dimensional case, this implies that the flow is supersonic at the nozzle exit (PB 0). 3. 3. 4 Time-Step Control Since the governing axisymmetric equations(3. 11), are hyperbolic, the domain of dependence principle discussed in section 2. 3. 3 still implies a need to limit the size of the time step used in the numerical integration procedure. Experience with the one-dimensional model has indicated that this (CFL) criterion suffices in an essentially unmodified form when used in conjunction with MacCormack's scheme. That is, no additional, more restrictive, condition affects At. One would therefore expect this to hold true in the multidimensional case.

197 Although a more precise relationship can be derived to match numerical and physical domains of dependence in the axisymmetric case, it adds an unnecessary complication to the computational procedure. Instead, we determine the time step of integration with the following dimensionless expression, which approximates the CFL criterion and is similar to that used with the full L-W schem e = C Ar Ax1 (337) At= CSF min(rvl + a)' lul+a i. e. a safety factor times the minimum of the radial and axial "onedimensional CFL criteria". Numerical experiments suggest that CSF can be as high as 0. 9 (almost thrice the value of 1/v8 used with the full L-W method) without numerical instability. For the results described below, CSF ranged from 0.7 to 0.85 as an extra SF measure of safety and Eq. (3. 44) was evaluated every m (specified) integration step at odd-numbered nodes (i = 1, 3, 5,.., I + 1; j = 3, 5, 7,..., J + 1) to determine at which node At had a minimum. For intermediate steps, current values of u, v, and a at this node were then used to evaluate At.

198 3.4 Axisymmetric Results and Discussion 3.4.1 Validation of Technique Validation of any numerical procedure for solving axisymmetric nozzle flow must generally be limited to a comparison of its results with those of experimental studies and/or with approximate analytical solutions that, for example, consider the transonic flow regime in the vicinity of the nozzle throat. Unlike the one-dimensional case, there exist no exact analytical solutions for axisymmetric, rotational, flow throughout a Laval nozzle. Furthermore, previous applications of MacCormack's scheme to multi-dimensional fluid dynamics (see Sec. 2. 3) do not involve internal flows, so that it is difficult to infer from them the validity of the scheme's use in the present problem. We therefore simply cite these previous applications again, to lend credence to the technique generally, and offer the following example results to demonstrate its applicability in the present model. The first example illustrates both the overall characteristics of the results produced by this numerical technique and the importance of the wall boundary conditions in an inviscid, non-heat-conducting flow. Consider the symmetric hyperbolic nozzle wall variation of previous examples given by Eqs. (2. 37) and (2. 38) in which AI = A = 5 and xthr = 0.5. To simulate conditions in an experiI E tEr mental "cascade" nozzle, such as the one that was recently built at the Gas Dynamics Laboratories of The University of Michigan, in

199 which the nozzle consists of several water-cooled annular disks, one computer run was attempted with Ta (= TI for j = 1, 2,... wall T~1 J + 1) held fixed at the reservoir value of 300 K. Initial conditions for p, u, and T corresponded to one-dimensional isentropic flow in the given nozzle geometry with uniform radial distributions of each. The initial radial velocity component at each x was assumed to vary linearly with r from a value of 0.0 at i = 1 to vwall(x) = u1 D(x) dR(x)/dx at i = I + 1. The uppermost nozzle in Fig. 3-3 shows the resulting steady state Mach number isolines for adiabatic flow (Q = 0.0). The lines have the expected qualitative appearance except at the nozzle wall where the static temperature and, hence, the speed of sound of the gas remains at a value that is greater than that of the adjacent, expanding, interior flow. As a result, the Mach number isolines "jog" downstream in an unrealistic fashion to satisfy the wall boundary condition. This phenomenon represents a clear example of the difficulties one may encounter by imposing more conditions at a solid boundary than are required for a well-posed problem. The absence of any diffusive mechanisms in the working equations prevents the numerical solution from displaying a smoother adjustment of the fluid properties to the higher temperature wall. In short, a fixed wall temperature is inconsistent with the inviscid, non-heat-conducting, Euler equations.

200 The three nozzle sequence shown in Fig. 3-2 illustrates the early temporal evolution of the Mach number isolines when only the velocity tangency condition, Eq. (3. 27), is imposed at the wall and flow properties associated with the above "hot wall" steady state now serve as initial conditions. There are several points to note in regard to the behavior of the numerical results. First of all, the isolines rapidly adjust to the new "correct" boundary conditions at the nozzle wall and assume the expected shapes for an inviscid flow. Although the solution was carried out until t = 3. 0, a near steady state is reached at about t = 0. 5. Even when t = 0.1 the Mach number values along the wall have nearly completed their re-adjustment. Within the nozzle interior and along the axis of symmetry the curves do not vary from their initial positions. Since the initial conditions also corresponded to a numerical steady state, these results would tend to indicate that: (1) "extra" nozzle wall boundary conditions, such as a fixed temperature, have only a very local effect on the overall numerical solution and (2) the solution attained by the present technique represents the true steadystate solution of the governing equations, because it remains stable over a relatively long time interval. The problem illustrated in Fig. 3-2 (for t = 0.0 to 3. 0) required approximately 9 minutes of CPU time for 421 nodes (I = J = 20) and nozzle dimensions: L = 0. 1 meters and Rthr = 0. 01 meters. The thr

201 LLJ z I(n) -J 3 cr 1.5 1.0 0.5 0.0 1.5 1.0 0.5 0.0 1.5 1.0 0.5 0.0 10___ I I II i I I I I II I 1 - 0.75 0.0 0.75 AXIAL DISTANCE FROM THROAT Figure 3-2. Mach Number Isoline EvolutionEffect of Wall Boundary Condition

202 minimum time step in this case is established by the radial CFL condition two nodes upstream of the throat at the nozzle wall (i= 21, j-9). We now consider a more severe test of the axisymmetric model and compare results with those of experiment. Figures 3-3 through 3-10 present the numerical results for adiabatic flow in a conical nozzle with a circular arc throat in which the convergent and divergent half-angles are 45 and 15 degrees, respectively, the throat radius of curvature, r, equals 0. 5 in. and throat radius, Rthr, equals 0. 8 in. 73-75 Cuffel, Back, and Massier have performed extensive experimental studies of this particular nozzle. Because they are the only apparent investigators to measure static pressures within the nozzle interior, as well as along the wall and centerline, we make use of their 27 results to help verify the accuracy of our numerical technique. Serra makes a similar comparison to his full L-W numerical scheme, while Hopkins and Hill compare their analytical technique for the nozzle's transonic throat region. Cuffel et al. (in the first reference above) 77 also evaluated the analytical methods of Kliegel and Levine77 and 78 Shelton, which include method of characteristics solutions in the supersonic flow region, and the finite difference technique of Saunders69 (also based on the full L-W scheme).

203 For values of the parameter r /Rth less than unity (i. e. rapidly varying throat wall contours), as in the present nozzle where this ratio equals 0. 625, the steady-state inviscid flow field is characterized by strong radial gradients in the flow properties and interior shock waves formed by the coalescence of Mach lines that begin near the wall tangency point downstream of the nozzle throat. These Mach lines are generated by an "overturning" of the flow near the sharply curved throat wall and the subsequent compressive turning at a downstream wall region, which is necessary for the flow to satisfy the inviscid 79 velocity tangency condition there. Migdal and Kosson7 discuss such nozzle phenomena in more detail. Evidence of this onset in shock formation appears in Fig. 3-3, which shows several Mach number isolines computed by the present technique (solid lines) and the corresponding experimental data points. The large radial flow variations are apparent from the highly curved isolines, particularly those downstream of the throat. For example, M varies from a value of 0. 8 to 1. 4 in the nozzle throat plane. The line at which M = 1. 75 is sketched in the figure, based on the experimentally observed wall values. It is in t this region that the wall Mach number decreases (static pressure rises), due to the impingement of the overturned flow from the throat. Theory and experiment show excellent agreement, except within that

w L) z (1) cn 2.0 1.5 1.0 0.5 - 0 - -1.0 Contour M 1.75 CD -0.5 0,0 0.5 1.0 1.15 2.0 AXIAL DISTANCE FROM THROAT Figure 3-3. Mach Number Isolines in a 450-150 Conical Nozzle With Circular Arc Throat

205 portion of the nozzle interior at which the radial velocity component is large (i. e. the two data points adjacent to the nozzle wall for both M = 1. 6 and 1. 8). Cuffel et al., note that these points have an unknown error due to flow inclination to their axially extending pressure probe and subsequent lower-than-actual static pressure readirgs. They included the points in their data to show trends in the flow characteristics. In this regard, the odd shape of the theoretical isoline for M = 1. 8, which is a result of the large wall curvature at the throat, bears a striking resemblance to those found experimentally in this portion of the nozzle (see Ref. 73). Figures 3-4 to 3-7 show the resulting distributions of Mach number, static pressure, density, and temperature along the nozzle centerline (i = 1) and wall (j = 31), as well as the "quasi-streamlines" (i. e. lines of constant r*) at i = 16 and 24. The abscissa in each graph represents a percentage distance from the nozzle throat plane, given as x in Eq. (2. 39). The first of these figures also contains the experimentally determined values of the Mach number at the wall and axis. Agreement is generally very good, particularly on the axis. According to Ref. 73 the results of Saunders finite-difference scheme also agree well with experiment at the nozzle centerline. However, at the wall boundary, where he used a simple extrapolation of adjacent interior values to compute fluid properties, his solution varies markedly

RUN NO. 20 NODES M I X RXI L) = 31 X 31 (INJ/l CT R) 1.866.X(THR) = NOZZLE LENGTH = 1.500 THROAT RADIUS = 0.800 PSN(1 = I TIE - 1.500.5 F 2.0 _ 1.5 -- 1- 0 - 0.0 -1.0 -0.5 PSN ( PSN [ PSN 2) = (16) 3) = (2D) ) = (311 0.500 A Experimental Data (Wall) o Experimental Data (Centerline) 0 m^ 0.0 0.5 1.0 1.5 2.0 AXIAL DISTANCE PLOT 1 Figure 3-4. Axial Mach Number Distribution in a 45~-15~ Conical Nozzle

Nt ) = t16) fpSNt t3}= KMNoUH. 20 X RXSL) = 31 X 32 0OS~j 1 3^.Aml 0.500:1 (Ito/I:: tF = 1.868, I'0 NOZZLE LENGTH 1.500 THoFT RBROWIS = 0. 500 pS4t) = I ~ TIIE " 1.0 r LI:3 LLI CrL.OA 0.6 0.2 -1.0 -0.5 rT TA A N\F A II AXIAL UL."".-" PLOT. 3 i p u 45D-15~ Conical Nozzle Figure3-5. Axial Static pressure o

RUN NO. 20 NODES (RDi.RL X RXIAL) = 31 X 31 R Id /R (Trf) = 1.866,X THR) = 0.500 NOZZLELEENGTH = 1.500 THRORT RRUIUS = 0.80 PSN(l) = 1, TIME = 1.500 1.0 r PSN( 2) = (16) PSN[ 31 = (21) PSNt 1) = (31) 1jJ (n c 0.8 0.6 0.4 0 00 (00 0.2 0.0 L -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 AXIAL DISTANCE PLOT * 5 Figure 3-6. Axial Density Distributions in a 45~-15~ Conical Nozzle

RUN NO. 20 NODES (ROIRL X R (IN)/Rf (TR) = NOZZLE LENGTH THRORT RAlIUS PSNI1) 1, 1.0 r ( FXIA L) = 31 X 31 1.866.X (THt) = 0.500 = 1.500 O 0.800 TIME a 1.500 PSN( 2) = (161 PSN( 3) = (24) PSN( 4)= (31) LU cU L — - 4: ce LL LU 0.9 0.8 0.7 0.6 0,3 C)r 0.5 0.4 L -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 AXIAL DISTANCE PLOT T 6 Axial Static Temperature Distributions in a 45 -15 Conical Nozzle Figure 3-7.

210 from that of the present scheme and from experimental values. Hence, the present method, in which the conservation equations themselves are used to compute values at the wall, appears to be more accurate. We note further that our results for the supersonic nozzle regime 73 agree with those of Kliegel and Shelton who used method of characteristics techniques with transonic flow solutions as starting conditions. Such solutions are expected to be quite accurate. The present numerical results for the conical nozzle and those of Serra appear to be nearly identical; both show good agreement with the above experimental data at boundaries as well as in the nozzle interior. Serra's investigation considered essentially the entire nozzle length of approximately seven inches, whereas we chose to truncate the nozzle on each side of the throat, such that L = 1. 5 in. and x * thr = 1/3. In this manner, a more economical computer run was possible and the incipient shock wave that forms near the nozzle centerline about three inches downstream of the throat plane was avoided. With I = J = 30 and t* = 1. 5, our CPU time was just over five minutes. max A comparison with the 80 minute run time of Serra, who used 3000 mesh points, indicates that the present technique offers an increase in efficiency of about 75%. This figure may be conservative, since it does not account for the larger grid spacing used by Serra. (Recall that his computational domain was 7/1. 5 = 4.7 times larger.) A

211 coarser mesh should allow for an increased At, although the smaller value of CSF necessary with a full L-W scheme offsets this gain. While MacCormack's scheme, combined with the non-conservation form of Euler's equations, would probably not implicitly handle the presence of a shock wave, it could be supplemented with (preferably) 19,22 explicit procedures to do so if desired Since the quasi-one-dimensional results indicated that a conical nozzle geometry may be desirable for gas blast circuit breaker applications, it is significant that the numerical technique developed here can accurately predict the flow characteristics for the above case in an economical manner, relative to techniques previously described in the literature. A small radius of curvature ("tight") throat may indeed enhance the flow expansion effects, just downstream of the throat location, that were noted to lead to more efficient convective cooling and, hence, to smaller time constants of static temperature decay. To conclude the discussion of this example problem and demonstrate the transient characteristics of the numerical solution, we present the temporal static temperature response in the mid-subsonic (j = 6), transonic (j - 11), and mid-supersonic (j = 21) regions of the conical nozzle in Figs. 3-8, 3-9, and 3-10, respectively. For each of these axial locations, curves for i = 1, 16, and 31 are shown. The

Lj Q Lq L RUN NO. 20 NODES R X AXIL) 3 X 31 R (I N)/R (T = 1.866,X (THR) = 0.500 NOZZLE LEGTH 1.500 THRORT RPI = 0.800 0.98 - - 0.94 \ LI /<0.92 U \ 0.90 0.0 0.3 ( PSNt 1) = ( 1. 6) PSN( 2l = (16. 6) PSN( 3) = (31, 6) Ib3 IN3 ).6 0.9 1.2 1.5 TIME PLOT * 11 Temporal Static Temperature Response in Mid-Subsonic Nozzle Region Figure 3-8.

Q Lt Q Q S%% 4;^ RUN NO. 20 NODES (IRDIRL X F A (IN)/R THTF = NOZZLE LENGTH = THRORT RRDIUS = 0.95 j 0.090 "- 0.85 0.80 I- 0.75 0.70 - 0.0 iXIRL) = 1.866.: 1.500 0.800 31 X 31 X THR) = 0.500 PSN( 11 = ( 1.11) PSN(b 2) = (16.11) PSN( 3) = (31,11) 0.3 0.6 0.9 1.2 1.5 TINME PLOT # 12 Temporal Static Temperature Response in Transonic Nozzle Region Figure 3-9.

Li Q Q - RUN NO. 20 NODES DROIRLX X RXIRL) = 31 X 31 A(IN)/R(THR) 1.866,X(THR) = 0.500 NOZZLE LENGTH = 1.500 THRORT RRODUS 0.800 0.85 0.80 ~ - o.75 0.70 J - 0.65 -.J. 0.60 0.55 0.0 0.3 C PSN( 1) = ( 1.21) PSN( 2) = 16.21) PSN 3) = (31,21) TIME PLOT * 13 Figure 3-10. Temporal Static Temperature Response in Mid-Supersonic Nozzle Region

215 curves evidence the achievement of steady-state conditions. However, unlike the one-dimensional case, the supersonic region appears to require a greater amount of time to equilibrate. Each of the figures illustrate the more rapid flow expansion along the nozzle wall that results in the substantial radial gradients noted earlier. 3.4. 2 Flow with Heat Addition To show some of the effects that heat addition has on axisymmetric Laval nozzle flow and to demonstrate the performance of the numerical model in this respect, the final example involves a small amount of energy input to the standard symmetric hyperbolic nozzle. Initial conditions are those of the adiabatic (numerical) steady-state shown in the first example of this section and heat is added impulsively when t = 0.0, such that the heat addition rate per unit mass, l/p, has the same value throughout the nozzle. (In practice, the entire dimensionless "source-term", At (y - 1) Q/p, in the differenced energy -3 equation was simply set equal to 10 for 1 < i < I + 1 and 2 < j < J + 1.) Since p varies significantly from inlet to exit, this has the effect of a spatially varying Q, but is otherwise an arbitrary specification. One perhaps best appreciates the magnitude of this heat input by noting that in the ensuing steady state, the mean total energy of the working fluid increases almost 50% across the nozzle, as illustrated in Fig. 3-11. For axisymmetric flow, the "mean"

RUN NO. 10 NODESF(RAi L X RXIAL)] 21 X 21 R (IN)/Al(TH = 5. 000, X (T1R) = NOZZLE LENGTH - 0.100 TIHORT RADIUS = 0.010 PSN(1) = 0, TIME = 0.0 1.6 r TIME( 2) = 3.000 0.500 1.5 1.4 1.3 1.2 1.1 1.0 IN, 0.9 L-1.0 _L -1 -- _I -0.5 0.0 0.5 1.0 AXIAL DISTANCE PLOT a 4 Figure 3-11. Mean Total Temperature Axial Distributions for Non-Adiabatic Flow Case

217 value of any dimensionless variable g(x) follows in the usual manner from an integration across the local nozzle cross-section; thus, 1 g(x) = 2 gr dr (3.38) 0 (and there should be no confusion between g and the abbreviated notation for predicted quantities in the difference equations). Also note from the figure that the initial mean total temperature distribution, which was calculated numerically, has a uniform value of unity, as it must for adiabatic flow. Figures 3-12 and 3-13 show similar distributions for the fluid density and static temperature. The higher values of Q in the inlet region cause significant changes in p and T compared to the adiabatic case, even though the energy input is much less than the amount that would result from the presence of an electric arc. Figure 3-14 presents the Mach number isolines for both the adiabatic (dashed lines) and the hot flow steady states. The substantial effect of such a small amount of heating, illustrated by the downstream shift of the isolines in this figure, is consistent with the behavior predicted by the one-dimensional model. Figure 2-40, for example, indicated a large rate of decrease in the massflow through the nozzle

RUN NO. 10 NODES(RROIRL X RXIFL) = 21 X 21 R lIN)/R (TH S.,X0 (TH) = 0.500 NOZZLE LENGTH= 0.100 THRORT RAFIUS 0.010 PSN(I) = 0, TIME 0.0 1.0 0.8 \ V 0.6 (nr )-0. LU 0.4 TIME[ 2) = 3.000 00 AXIAL DISTANCE PLOT * 1 Figure 3-12. Mean Density Axial Distributions for Non-Adiabatic Flow Case

TIME (23 = 3.000 L ( I ( I NOZZLE LENGTH = 0.100 THOROT RROIUS = 0.010 PSN(1) 0, TIME = 0.0 1.4r 1.2 1 1. LJJ 0.6. —4 0.8 _ JI 0*0.2 0.2 I I I J-_ — ---— I —-A.-. - - I.A It 0.A0 L-1.0 BVOU1elrleb -0.5 0.0 0.5 XoU AXIAL DISTANCE PLOT * 2 Figure 3-13. Mean Static Temperature Axial Distributions for Non-Adi aatic Flow Case

w () z I-r cn 0 1.5 1.0 0.5 0-1.0 Wall Contour k3 0 -0.5 0.0 0.5 1.0 AXIAL DISTANCE FROM THROAT Figure 3-14, Mach Number Isolines for Both Adiabatic and Non-Adiabatic Flow

221 (and corresponding decreases in the Mach number distribution) for small increases in electric current input. Note also that the isolines shift downstream in a nearly uniform manner, thus retaining their original shapes. Although we carried out other hot flow runs with the axisymmetric model, the magnitudes of the energy input and the subsequent results were comparable to that of the preceding case. As noted earlier, increases in Q entail substantial decreases in the allowable time-step of integration. The above example required 12 minutes of CPU time for only 421 mesh nodes. However, a similar run, in which Q was fixed at a value of 0. 5 for all nodes, indicated that finer mesh sizes are needed when larger amounts of energy are added in the supersonic nozzle region. Inaccuracies caused by very low density values near the nozzle exit are again the reason for this. As a result, one obtains processing times of over 30 minutes (for a 31 by 31 grid) even for these small amounts of heating. Since the one-dimensional model shows that Q reaches values over 11 near the nozzle throat even for the low current case of 50 amperes (and values near 60 for currents of 500 amperes), CPU times of several hours appear to be necessary for a multi-dimensional analysis of the AC circuit breaker problem.

4. CONCLUSIONS In the preceding analyses we have numerically integrated the unsteady quasi-one-dimensional and axisymmetric Euler equations of gasdynamics for non-adiabatic flow through a Laval nozzle with negligible back pressure. Through the use of a second-order, explicit, two-step, finite-difference scheme and appropriate boundary conditions at the nozzle inlet and exit planes, we considered the entire nozzle interior domain for which the spatial variations in Mach number ranged from nearly zero to high supersonic values. In the axisymmetric case, a consistent set of boundary conditions along the nozzle centerline and wall allowed for a stable integration procedure. The overall numerical technique proved to be significantly more efficient than previous methods used for internal flow calculations, while retaining at least the equivalent accuracy. Artificial numerical dissipative mechanisms were not included in the difference procedure. The introduction of an intense, time-varying, bulk heat addition term into the quasi-one-dimensional energy conservation equation, which simulated conditions in a gas blast circuit breaker, led to useful qualitative information regarding such devices. For both stepfunction (DC) and sinusoidal (AC) current inputs, the numerical model predicted the transient reverse flow phenomenon which occurs in the converging nozzle region during the first moments of energy input. 222

223 Results showed that this is caused by the momentary development of a spatial maximum in static pressure just upstream of the nozzle throat. A comparison of three dissimilar nozzle geometries indicated that a conical nozzle may be most effective in minimizing this effect. The same type of geometry (conical) also proved to be the most efficient in dissipating the flow energy. That is, it had the smallest time constant associated with the decay of static temperature (and electrical conductivity) within the nozzle. In the numerical experiments, the spatial minimum of this time constant value occurred downstream of and near the nozzle throat. After this study was completed, the author learned that in recent experiments conducted at the Brown Boveri Company in Switzerland, with an arc-nozzle apparatus similar to that studied here, the arc first becomes nonconducting at a position 80 along its axis that lies just downstream of the nozzle throat. This agreement between experiment and the present model provides strong evidence of the importance of axial convection as an energy transfer mechanism in gas blast circuit breaking. The axisymmetric numerical model presented in the previous section offers a significant improvement in computational efficiency over past techniques for both steady and non-steady Laval nozzle flow. Because of the decreasing time step requirement as the static temperature of the gas becomes larger, this study indicates that a

224 physically reasonable, and numerically efficient, axisymmetric model of the circuit breaker problem cannot be achieved with present day numerical techniques, if cost is an important factor. The inclusion of thermal conductivity and/ or viscosity in the working equations would result in further reductions in At to maintain numerical stability, because of the need for a finer differencing grid to resolve boundary layer effects at the arc boundary and nozzle wall. From a numerical standpoint, the quasi-one-dimensional model holds the greatest potential for further investigations of gas blast circuit breaker performance, because it allows for a more realistic simulation of "arc-like" heating and requires less computer processing time than higher dimensional models. The energy addition term can easily be altered to include optically thin gaseous radiation and wall heat transfer effects (although the latter is expected to have little significance, because of the large magnitude of the Joule heating term). In addition, it would be of interest to study the effects of an external electrical circuit, perhaps with some idealized circuit model whose governing differential equation implicitly includes the arc impedance. In this manner, one might avoid the problem, experienced in the study, of an unusually large Joule heating value near current zero for real gas conductivity values. The above improvements would probably require a computation time on the order of one-half hour per run.

APPENDIX QUASI-ONE-DIMENSIONAL COMPUTER PROGRAM The quasi-one-dimensional model was programmed in the FORTRAN IV language for the University of Michigan's IBM 360/67 dual processor computer system. The general algorithm used to carry out the numerical solution is as follows: 1. Input all the necessary problem and program control parameters (e.g. tax CF' Pref T refetc.) max' SF ref ref' 2. Input the number of nodes, J + 1, y of the gas, the initial dependent variable distributions, pj, u., and T., and j' j j nozzle geometry, Aj; 1 < j < J + 1. 3. Initialize all necessary variables, including p, Q1 (= 0. 0 for all times), Ttl = 1. 0 and t = 0. 0. 4. Compute the pressure and heat addition distributions for 2 < j < J + 1. Also compute the speed of sound at three specified nodes to be used for the time-step calculation. 5. Depending on the current value of time, t, output the values of the dependent variables, pressure, and heat addition for 1 < j J+ 1; ift > t terminate run. - - - -- max 6. Calculate At at three specified nodes, choose the minimum value and compute X = At/Ax. 225

226 7. Compute the predicted values of the dependent variables, fY, 2 < j < J, using Eqs. (2. 19). 3 8. Calculate the predicted inlet velocity, ul, from Eq. (2. 30). a) If u > 0 and T = 1. 0, calculate T from Eq. (2. 27) and go to step 9 below. b) If u _ 0 or Ttl > 1. 0 (i. e. the nozzle is "clogged" or in the process of recovering from a clogged state), compute T1 from Eq. (2. 33) and then Ttl' from Eq. (2. 26). c) If Tt > 1. 0, set Tt = T' and go to step 9; otherwise, set Tt = 1. 0 and re-compute T1 by Eq. (2. 27). 9. Compute P1 and P1, the predicted inlet pressure and density, from Eqs. (2. 28) and (2. 29). 10. Compute the predicted heat addition values, Qj, for 2 j < J + 1. 11. Calculate the corrected values of the dependent variables, f., 2 < j < J, from Eqs. (2.20) n+1 12. Calculate the corrected inlet velocity, u1, using Eqs. (2. 31) and (2. 32), and determine the remaining inlet n+1 n+1 n+1 properties, p, T, and p1 by a procedure analogous to that in step 8. 13. Increment the current value of time; tn+ = t + A t. 4n+1 n n 14. Go to step 4 and continue the integration procedure.

227 The program is loaded and run from a remote, time-sharing teletype. terminal. In this manner, one may monitor the progression of the numerical technique by having values of the dependent variables printed-out to the terminal at specified time intervals. The alteration of program parameters, such as the safety factor, C, in the SF stability condition, sometimes mitigates instabilities before they completely distort the results and cause premature termination of the program. The algorithm allows for output of the results at specified intervals to a 9-track magnetic tape, and provisions for "restarting" previous runs are also available. The following pages contain listings of the major subroutines and variables used in the program along with descriptions of the function of each. Unless units are specifically noted, the variables are dimensionless and correspond to those in the equations of Sect. 2. 0. Note that most of the important variables have been made "common" to all subprograms. Finally, a listing of the computer program is presented.

228 SUBPROGRA MS BDRYS calculates inlet and exit values of the dependent variables; initial entry initializes the inlet static pressure and sets values of constants. BDRY1 entry to BDRYS that calculates the predicted values of density, velocity, temperature, and pressure at inlet and exit. BDRY2 entry to BDRYS that calculates the corresponding corrected values. IO handles all input and output to teletype, line-printer, magnetic tape, and disk files; initial entry accepts program starting and control parameters, reads initial conditions, and outputs them to line-printer. I01 subsequent entries to IO;outputs solution values to all devices as desired and terminates or restarts program when t > t max PROPS calculates fluid properties (p, E, and a) and heat addition Q at all nodes; initial entry fixes constants and initial values. PROPS1 subsequent entries to PROPS to calculate predicted or corrected fluid properties and heat addition. SIGMA calculates electrical conductivity as a function of static temperature. SPECIAL SUBROUTINES ENTPTS performs curve-fit of tabular data (see Ref. 57); initial entry inputs known (x,y) pairs. GETPTS subsequent entry to ENTPTS that provides vector of fitted values for specified values of x.

229 SETDSN system routine that assigns an I/O device to a given unit number. SKIP system routine that allows manipulation of magnetic tapes (e.g. positioning to given data file and/or record). TRACER system routines used to "trap" program interrupts PGMTP (caused by logical or execution errors) without terminating program execution. INPUT PARAMETERS CFL "safety factor", CSF, in CFL stability condition, Eq. (2.35). CRR gas constant; atm-m3/ (kg-OK). DTDF time increment at which solution values are to be written to magnetic tape. IMAX maximum or peak current, Imax (amperes for IQOP / 1). IPO number of time increments, DTDF above, to elapse between line-printer output of solution values (e.g., IPO = 10 and DTDF =. 1 implies printout every At = 1.0). IQOP indicates option for form of heat input; IQOP = 1, Ima for JHF < j < J + 1, if Ima < 0 Q = (j - JHF) Imax/ (J + 1 - JHF), if I > 0 = 2, DC current input, Eqs. (2. 46) and (2. 47) = 3, AC current input, Eqs. (2. 46) and (2. 48). IRES restart switch to allow restarting in middle or end of tape data file.

230 IRUN JCFL JHF JPO1, JPO2 PREF TMAX TREF run number. vector of 3 values of j (nodes) at which CFL criterion is to be tested. value of j at which "heat front" occurs (i. e. Q first becomes non-zero). two values of j at which dependent variable values will be written to the terminal in order to monitor program progress. reference pressure, P atmospheres. ref' maximum time, t ma, to which integration proceeds. reference temperature, T K. THER VARIABLES OTHER VARIABLES AARC ACFL AS AW CON DELX DT DTIM E FREQ GAM vector specifying axial variation of arc area, Aar; m. vector specifying speed of sound at JCFL. vector, Aj*. vector, A.. J vector, a.. distance between nodes, Ax. vector of 3 At values at JCFL. current time step of integration, At. electric field vector, E.; volts/m. frequency of AC current input, v; sec ratio of specific heats, y.

231 GEOM JJ1 JT Ll NREC NSTEP p Q R,U,T RB, UB, TB TT1 XTIM vector specifying nozzle geometry: GEOM(1)= Rth; m GEOM(2) =L; m GEOM(3) = xh GEOM(4) = A * GEOM(5) = A* E GEOM(6) and GEOM(7) specify type of converging and diverging geometry, respectively; 1 = hyperbolic, 2 = conical, 3 = logarithmic. GEOM(8) and above contain other parameters to specify the converging and diverging geometries depending on which of the three are used in each case. total number of nodes, J + 1. thr' A = At/Ax. number of records written to tape. current number of integration steps. static pressure vector, P.. J rate of heat addition vector, Q.. vectors of "corrected" dependent variables. vectors of "predicted" dependent variables. inlet total temperature, Tt1. current value of time, t.

232 PROGRAM LISTING REAL LI C0OM0N /REG/ SAVREG (18) COM0MN /PVARS/ RB(51),UBC51),TB(51) C0MM0N /CVARS/ R(51 ),U(51 )tT(51) C0MMON /PR0P/ P(51),Q(51),E(51),ACFL(3) C0MMON /SCRACH/ RS,US,TS,T 1,T2,3 T4 C0MMON /I02ALL/ JJJJl,DELXAW(5 1),AARC( 51 ),LI,GAM,GMI, & DDTGITTI,XTIM,DTIM,CFL,DT(3),JCFL(3) CALL TRACER(- 1) CALL PGWTP(SAVREG,&5) C C CALL ALL SU8PROGRAMS FOR INITIALIZATION C 1 CALL 10 CALL PR0PS CALL BDRYS DX IN=I./DELX C BEGIN INTEGRATION LOOP C 2 IS:I C C CALC FLUID PROPERTIES AT EACH NODE C CALL PR0PS (R,T,I S, &3) C C OUTPUT CURREtT VALUES OF DEP VARIABLES; C TERMINATE PR'0G IF TIME *GE. TMAX C 4 CALL O11(&1) CDT =CFL*DELX C C CALC MINIMUM TIME STEP BASED BN 3 SPECIFIED MIDES C DO lO I1t 3l 10 DT (I ) -CDT /U (JCFL (I) )+ACFL(I) ) DTISMIN. 1 (DT ( ) DT(2)) DTIA MIN 1 (DT(3),DTIM) LI lDTIM*DXIN DTG = TI MGMI C C CALC PREDICTED INTERI R VALUES C DO 20 J:,JJ JlJ+1 RS =R(J ) US3 U(J)

233 TS T(EJ) T I =R.(J I ).RS T2 sU 1)-US T3 T (Jl) - TS T4 US*AW (J) RBJ3)RS. *(US* T1S*T2 ) -DT IM*RS* T4 U (J)US -L *(,US*T2+T3+TS* TI /S ) 20 TB(J ).TS.L l,*(U S*T3+1GMI* TS*T2)+DTG l*((Q J)/RS ) -TS*T4) C CCALC CLC PEDICTED BDRY VALUES C CALL BDRY (&) 1S2 C CALC PREDICTED FUID PRBPERTIES & Q CALL PR@PS1 (RSBTBTIS&$) C C CALC CORRECTED INTERISR VALUES C D3 50 J2JJ RS sB(J) T4SUS*AW(J) T23 US.-UB (J! ) R)4.=S R(JAR )*(US*T S*2)-DT*RS* R (J)=.5, (.R (-JI)+.RS -L* ~ (& ST l +RS 1'2:) -DN l ~l, RS, 1;T4 ) U(J):,*(U(J)+US -L *(US*T2+T3+TS*T1/RS ) ) 30 T(J ) s5*(T(J )+'S- Il*(US*T3+GIStl*TS*T2)+DTGI*( (Q(J)/RS) & *TS*T4)) C C CALC CeRRECTED BDRY VALUES & INCREMENT TIME C CALL BDRY2(&3) XTI MXTI N+DTII G~ TO 2 9 P INT 10 SAVREG( ),SAVREG(2) XTI >1 100 FRHMATC0OPRIG INTRPT, PSW -'Z8,IXZ8tX,'TIPE =* & F9.4) CALL PGTP (SAVE~G, &5) 3 XTIM=1000. GE TD 4 END

234 SUBRBUTI E It REAL IMAX,Ll INTEGER WHERE INTEGER*2 LEN4 //8000/,LEN2 CBIMMN /REG/ SAVREG,(1B) COIMMN /CVARS/ R(51),U(51)tT(51) CIMMN. /PROP / P(51 ),Q(5 E ) E( 51 ),ACFL(3) C9M0aN /12P/ IMAX,I QBP,GEf((20),PREF,TREF,CR,IDF,TIMDF CMPIM1N /I 2ALL/ JJ,JJ DI,ELX,AW ( 51 ), ^ARC( 51),L l,GAM,GM & T DTG lTTI XTI'XTI-,MCFLDT (3),IJCFL(3) EQUIVALENCE (XT4GE~M(3)) DIMENSION TIM(2o00),AS(5 ) DATA YES/'Y~/ CALL PGMTP(SAVREG,&6) C C INPUT PRIG CBNTROL PARAMS FRiM TERMINAL C CALL SETDSN(17, **PRINT*',FDB7) I PRINT 100 100 FORMAT (OENTER TIM1E CSNTROL PARAMETERS:'/ & 5X'rTAX CFL,DEL-T(DF),IPt,& JCFL(I-3)*) REA6 10 lOlAXCPL,ITF,IP,.(JCFL(I"),!:1,5) I01 IFRMAT(3F1 5.6,51 5) PRINVT 102 102 FIRMAT('OENTER HEAT PARAMStIMAXCA'4P$) & IQ0P(lsQ CSNSTOI'/ & 5X,' 2:DC,3AC)') READ 103,IMAX,I QP 103 FRMAT (f 15.6,5I 5) PRINT 104 104'FfRMAT (OEP!TER REMAINING PARAMETERSt PREF(ATM),TREF(K),'/ & 5X "GAS COMST(A:TM-5M3/KG-K ),IRES,tiUN,JPr I,JPs 2*).REA l' 10,PREF,TREF,CRR,IRESIRUNJPIl,JP52 CRR:1./CRR IF(IRESEQ,.) GS Tb 2 C C IF NO RESTART{, READ INITIAL VALUES SF DEP VARIABLES C AND GEM VARIABLES FRIM FILE CALL SETSSN(1, -INPUT',FD ) READ(l) JJ1,ISETGAMTTI PAR,(GEWt..() I:Il,PAR), & (AW(J),J=1,JJI),(AARC(J),J:1,JJI),((ASJ)J =1,JJI), & (R(J),J =,l JJ I),(U(J ),J,JJ. ),.(T:.(J>,J l,JJ 1) REWIND 1 CALL SETBSN (, *-fINKx**,FDB l)

235 NREC:O G6 T $3 C C IF RESTART, PtSITItN TAPE TO LAST RECORD AND READ C INITIAL VALUES 2 CALL SETDSN(8,*T*'FDB88) CALL. SKIP( - -2,D 8 &204 &208,&21 2) CALL REgA(TI,(CI) EL2! I oLN2 9FD8 ) LREC =TIM(200) NREC LR EC REA^ (). JJ4 IRUN GA,TTI NPAR(GE@M(N), 1," PAR),(AW(J) & J =12 JJi ) (AARC.(J) J I,JJI )... (AS J),J Il,JJI ) CALL SKIP(0,93,tBS,&304,&308,&312) PRINT 115 115 FtR1AT(*ORESTART? AT LAST DISTR. RECRD?') READ l4,O9P 114 PmR.ATCA4) IF(0PoEQYES!) GI TO 8 PRINT 1t6 116 F RNTAT(-'OETER? VALUE *F NREC FIR STARTING DISTRS') READ 17tNREC 117 FPRM'AT(lO2 5) CALL SKIP(0,NREC LREC FDB8,&204,&208,&212) 8 READ8() (R(J) J 1,JJ ),(U(J),J =l JJ I ),(T(J)',J 1,JJ I ), & (P(J),J l JJ l) (Q(J),J =1JJl), (E(J),J =l,JJ ) CALL SET9SN"(8I *INK*,FDBB8) C C *UTPUT ALL INITIAL INF Te LINE PRINTER C 3 WRITE(7,118) IRUNIRES,ItMAXIBQPTMAX,CFL, & (JCFL (I),I I 3=,) voTDFIP0,ISET,JJ,GAM, & (GE~M(I),I I=lNPAR) 1 18 F.RMAT ("1,42X,'QUASI -NE-DIMENSIeNAL ARC-NMZZLE PROBLEM' & //57X,'RUN #',I4/6RESTART:', l2,5X,'MAX CURRENT(AMPS)=' &,F7'1!5X ~QD6T OPTION #*,12/ & OMAX TIME(INITIAL):',F62,25X,CFL COEFF -,F5.2,5X, & JCFLC(1-3): 3165X, DTDF',F7,45X,' P0,I 3/ &'-'/0',57X,'GE0METRY'//55X,*INPUT SET #',I3// &'-t AXIAL "N0DES 3 l 3,5X-,'GAMMA IF GAS.',F5,2/ &'OGEMETERY PARAMET~ERS FIR NOZZLE:'//(5EI 5,5)) WRITE(7,119) (AS(J),J:lJJl) 119 FlR AT ( (ONZ AREA VAR ATIIN(A(X)/ATHRBAT): //(5X,8E1 5.6)) WRITE( 7 20) (AW(J),J =9JJ ) 120 F0R AT OLOG AREA DERIVATIVE((DA/DX*)/A)t'//(5X,8E 5.6))

236 WRITE(7,121) (AARC(J),J 1,JJI) 121 F9RMAT(*OA??C AREA VARIATION.(METERS**2):'//(5X,8E 5,6)) C C INITIALIZE COEFFICIENTS,ETC. C JJ =JJ1-1 GMI:GAM-1, DELX JJ DELX =1./DELX IF(IRES.EQ.1*) GO TI 15 XTIM=0.O TIMPIOFOO. GO T0 16 15 XTIM-TIPM(REC) TIMDF:XTI M+DTDF 16 NSTEP:-t JT=.OOO+XT/DELX ICNT:IPt XIMAX I MAX I Q*PX =I QIP CALL SI MA(TREF,JJI) RETURN C C INTERNED INTRIES: *UTPUT BEP VARIABLES" T TAPE AND-0R C LINE PRINTER AS DESIRED C ENTRY II1(*) NSTEP:NSTEP+ I IF( IDF EQ..OA ND XTIILT TMAX). 6 TO S IF O-0 NREC.NREC+1 TIN(NREC) = XTI CALL SETDSN(t8,*T*,F88) VRITE(8) R) ((,Jl 1,JJl ) (U(J),J:1,JJI), (T(J) J l1,JJI), & (P(J),J 1,JJ. ) (Q (J),J l,JJJI),(.E(J ),J.l,JJ 1 ) CALL SETBSN(8,'*MSINK*, FDB8) TI!MDF =XTIM+DTDF CNT:I CNT+1 IF(ICNT LT.IP*.AND XTIM.LT o,T-AX) Of T8 5 ICNT:0 PRINT 107,NREC,NSTEP,XT Itl,TT,DTIM, (DT(I),I" l.,3 ),R(JPI ), & U(JPIl ), T(JPI ),R(JPS2),U(JP02).,T (JPe2) 107 FRMtAT('ONREC'*,I 52X,'NSTEP *=i 5,2X, TIME *,F10.4,2X, &'TTI =',-F6.3/*DTIM:',EI.4,*, DTt 3El 1.4/(SX,E1I 5.4)) WRITE(7,108) XTIM,NSTEPSTIM

237 108 F8RMAT(C"TIME',FO 4,5 X,'STEP #',IS5XDTIM =',E12.5/ 25 ICONXTINUE 110 FBfRM AT CE I55) I11 F/RMATr(o+* 14X *l' 14X,*lg4X,*`g,14X,** 14Xg$ ) 5 IF(XnTM*LY TMAX) XRETURN C C TERNINATE PR@G 1 RESTART WITH NEW TMAX C "PRINT 122 122 FORMAT ('iESTART WITH G'REATER TAX *) READ 114,0P I'F(.PNE1Y4ES>) G Ti 7 PRINT 125 1 7 PRINT 123 123 FORMAT(OENTER NEW TWAXIP9~IW4AXIQ0P*CFL,& DTDF*) READ I 3., TMAX l P I AX,IQ0PCFL~ TDF 11t F8RNAT ( Fl 5 e.6L5, F5, S,69I 5, 5F 1 5 *6 ) P1RINT 105 105 FlRMAT(''OCHANGE J(P ) SR J(CFL) VALUES?') READ 114 0P IF(P.NESYES) GO T1 $ PRINT 1.06 106 FVR'AT( "ENTER JPSI JP02,AN JCFL( 1-3) ) READ 117 JP I lJP 2, (JCFL(I),I =t,3) 9 IF( XMAXEQ.XIAX A D l.QP oPEQ IQSPX) RETURN I QPX =I QP CALL PR@PS CALL PRIPS (R T,2,&7) RETURN 6 PRINT 109,SAVREG(I ),SAVREG (2),XTI 109 F RMAT(~OP RG INTRPT(Ib) PSW. *,Z8lX,iZ'TIME:',F*94) CALL PGINTP(SAVREG,&6) 7 TIM(2000)=NREC ASSIGN 19 T$ WHERE CALL SETDSN(89,*T* "FDB8) CALL WRITE(TI~(l) LENl 90,LNM3,FDB8) WRITE(S) JJ9I RUN9 GANTT,NPAR,(GE0 (N),N=I.,NPAR),

238 & (AW(J),J=,JJI ), (AARC(J),J=I,JJ'),(AS(J),J:l1,JJI ) CALL SETBSN(8*,*MSINK* *,FDB8) PRINT I 2,vTIM,NRECTIM(NREC) 112 FSRMAT(*OFINAL DTIM =' -El I 4, 3X,'NREC,I4,3X, &'TIM(NREC) =*,9*5) 18 CALL SETD8(9,'*-RID',FDB9) WRITE(9) JJI,ISET,GAM,TT I,NPAR (GEBM(N),Ns-1,NPAR), & (AV(J),J,JJ I), (AARC(J),J=1,JJI ), (AS.(J),J. 1,JJI), & (R(J),J:,JJI),t(U(J)J.l,J JJI),(T(J),J=l,JJ ) CALL SETDSN(9,'*ISINK*',FDB9) G0 TO WHERE,(17,19) 19 PRINT 124 124 FORMAT( OLAST TW. REC0RDS & RESTART FILE(-R1D) WRITTEN'/ & IX,'RETURNING TO MTS; T0 START NEW RUN(1l) EMPTY -RID,'/ & IX,'(2) CHANGE -INPUT, (3)WTM, (4)DSN, (5)RES') CALL MTS RETURN I 204 PAUSE'FIRST SKIP(4)' 208 PAUSE'FIRST SKIP(8)' 212 PAUSE'FIRST SKIP(12)' 304 PAUSE'SEC0ND SKIP(4)' 508 PAUSE'SECOND SKIP(8)' 312 PAUSE'SECOND SKIP(12)' RETURN END SUBROUTINE PROPS REAL IMAXI 1,I2,LEN,L COMMlN /REG/ SAVREG( 8) COMMON /PROP/ P(5 ),Q(51),E(51),ACFL(3) C0MMON /If2P/ IMAX, IQ0P,GESM(20),PREF,TREF,CRR,IDFTIMDF C0MMON /I 02ALL/ JJ,JJ 1,DELX,AWV(51),AARC(51 ),L1 GAM,GMI, & DTG I,TT 1,XTIM,DTIM,CFL,DT (3),JCFL (3) CBMM0N /SCRACH/ RS,US,TS,T I,T2,T3,T4 DIMENSION CQ(51),R(51),T(51),C0N(51) EQUIVALENCE (LEN,GE0M(2)) DATA PI2,CQRCVR,FREQ/6.28518,3.22E7,3. 8E2,60./ CALL PGMTP(SAVREG,&6) C INITIALIZE VALUES eF ELECTR FIELD & Q C E(1)=0.0 IDF -0 IF(IQfP.NE.l) G0 TO 1 PRINT 102 102 F0RMAT("OENTER JHF') READ 103,JHF 103 FRMAT(I 5)

DO 10 J=2,JJ E(J ) =0.0 Q (J)=.00 IF(J.LToJHF) G0 TO 10 IF(INAXGTOo0) GO TO 11 Q (J).-I MAX Ge' TO I O 11 Q (J):(J -HF)*IMAX/(JJI-JHF) 10 CONTINUE RETURN C C CALC REF VALUES & C0EFFS FIXED THRU0UT INTEGRATI0N I RREF CRR*PREF/TREF VREFS SQRT (PRE /RRE ) QREF =CQR*PREF* VREF /LEN IF(IQ0P.EQ.3) G0 T9 2 I 1 MAX I2I1 1*1 I G9 TO 3 2 VREFs:CVR* VEF P12F=PI2*FREQ*LEN/VREF 11:0,0 I2:0,0 3 DO 20 J2,JJI 20 CQ (J ) = 1 /(QREF*A ARC (J )*A RC(J) 4 RETURN C C INTERM1ED ENTRIES: CALC P,Q,& E C ENTRY PR@PSI (RT,IS,*) IF(IS.EQ,2) GO T' 12 DI 30 I:1,3 30 ACFL(I )SQRT(GAM*T(JCFL I))) IF(XTIM.GE.TIMDF) IDOF= DI 40 J ~,JJI JX:J RS=R(J) TS T (J) 41 IF(RS.LT.O,,0RTS LTO.0) GO TI 5 40 P(J) RS*TS 12 IF(IQBPEQ.1) RETURN IF(IQP, EQ 2,2RISoEQ e) GO TO 14 I 1=I MAX*SIN (PI2F(g XTIM+DTIM) ) 12:= 1*I 1 14 CALL SIGI(TCON) DQ 50 J=2,JJI

240 Q(J) = CQ (J)I 2/CC(J) IF(IDF.EQ*0) G0 TO 50 E(J )=I I/(C N(J)*AARC(J)) 50 CONTINUE RETURN 5 PRINT 100,JXXTIM,IS,R(JX),T(JX) 100 FORMAT(ONEG DENSITY 0R TEMP -J =',t3,5X, TIME -',F9.5/ 3X,*STEP f#,I2,5X, *R & T:',2E15.5) RETURN I 6 PRINT 101,SVREG(I),SAVREG(2),XTIM 101 F0RMAT('OPR0G INTRPT PSW =:,iS IX,1XZ8,'Xt'TIME:',F9.4) CALL PGMTP(SAVREG,&6) RETURN I END SUBROUTINE SIGMA(TREF,JJI) REAL TT(12)/4E35,5,E3,6.E3,7E3, 8.-E,9.E3,I.E4,. 1E4,1,2E4, & I.3E4,.4E4,1.SE4/,S(12)/.75,35.56,236.5,650*.,1232., & 1920.,2656*,3397.,4 t O.,4780.,5440*,6 40./ DIMENSISN TND(I),C6N(I) C C INITIALIZE CURVE FITTING REUTINE C CALL ENTPTS( 12,TT,S) RETURN C C CALC ELECTRICAL CONDUCTIVITY C ENTRY SIGI(TSDCeN) D0 10 J=2,JJI TX TND (J) *TREF IF(TX.GT.5.E3) GO TO 2 I C~N(J)7. 112E-3*TX Ge TB 10 2 IF(TXLT,,l6E4) G9 TN 3 PRINT 100,J,TX 100 FiRMAT(OFR0M SUBPR0G SIGMAt T(',I2,') =',E12.5) 3 CALL GETPTS (l,TX,C0N(J)) 1 O CONTINUE RETURN END SUBROUTINE BDRYS REAL LI C0BMiN /REG/ SAVREG (18) C0MMSN /CVARS/ R(51),U(51),T(51) CbMMltN /PVARS/ RB(51),UB(51),TB(51) CtMMON /PRP/ P(51 ),Q( 51 ),E(51),ACFL (5)

241 C9INMN /SCRACH/ RSUSTS,T,T2 T3,T4 C0MM@N /I 2ALL/ JJ, JJI OELXAW (SI ),AARC(51 ),LI GAM,GM1, & DTGQ1 TTlXTIM1DTiMCFLDT(53) JCFL(5) DATA YES/'Y"/ CALL PGMTP(SAVREG,&6) C C INITIALIZE CONSTANTS C P(1)sR( 1 )* T(S) K2 JJ 1 KI K2-1 GM2 =.5*GMI /G AM Gi3 =GAf /G(t RETURN C C CALC PREDICTED INLET & EXIT VALUES F DEP VARIABLES C ENTRY BDBRY(*) IS=1 T2=U(2) oU( I ) T3=T(2)*T(1) UB( 1 ) U( I )-L l*(U( I )*T2+T3+T(l)*(R(2) -( ))/R( 1 ) ) T4-=U( l)* BL( l ) IF(UB(I),LT0,oOO.lR.RTTlIGTl, oO) Go To 12 I T4 GM2* T4 T8(I ).s1T4 T4|1,/TB(l) G60 T 14 12 TBi): T(I)-L * (U(I)*TS +Gf *T 1)*T2)'DTG1*T ( )*U(I)*AW( ) T4 S 1 s+GS2* T4/TB ( ) TS =Tb( l )*T4 IF((TS,3To1eO) GO TO 15 TI.l o 0 IF(UB(1)oG3EOO) G0 T1 14 T4 UB( I )*UB(1-) G9 T l I 15 TT.=TS 14 IF(T4oLT.OoO) GO TO 5 P(l ) T4**GM3 RB(S):P(I )/TB(l) RS (JJ) =RS(K )3o., (RBS (K2)-T8(JJ ) ) TB(JJI!)TB(KI)-3.*(TB(K2)-TB(JJ)) RETURN C CC CA CORRECTD IiNLET & EXIT VALUES.F DEP VARIABLES C ENTRY BORY2(*) IS 2 RS 53*(RB( )-R(C2))+Rt8(3)

242 US * -3 (UB( I )-UB(2) )+UB(3) TS 3.*(TB( )-T8() )+TB(3) T2 =UB( )-US T3WT3( i )-TS U( -) =,5*(U( )+UB( l )-LI*(UB( 1 )*T2T$+TB(1 )*(RB( I )-RS) & /RB(1))) T4 U(1)*U( 1) IF(U(. ).LT.O.0R o.T TI.NE..O) GB T 22 21 T4:GM2*T4 T(1)=:1.-T4 T4: 1./T (I) GO TO 24 22 T(l)=,S*(T(I +TB( )-Ll*<UB(i)*T3+GMl*TB( )*T2)-DTGI*T3(I) & * U(1 )*AW ( )') T4:s.+G t 2* T4/T ( ) TS =T( I )*T4 IF(TS.GT.BI.O) GO T0 25 TTIl 1,0. IF(U(I)-.GE.O.O) G0 T 24 T4=U( )* U( 1 ) GO T0 21 25 TT ITS 24 IF(T4.LT.OO) G0 TB 5 P ( I ) T4** GM3 R(l)=P(l)/T(l) R (JJdI) =R(K 1 )- -3*(R (K2) -R (JJ) ) U(JJ!) U(K )-3,*(U(K2)-U(JJ)) T(JJ ) T (K )-3.* (T(K2)T (JJ )) RETURN 5 PRINT l00,T4,XTIMIS,R(I),U>l),T(I),RB(I),UB(I),TB(l) 100 F0RMAT(COT4'E1253X, *AT T:'F9.5,3X,'& STEP',I2/ & (3E1 55)) RETURN I 6 PRINT 1 0lSAVREG(l ),SAVREG(2),XTIM OI FORMATCOPR 0 INTRPT PSW S: I,X,ZX,'83X"TIME:=,F9,4) CALL PGMTP(SAVREG,&6) RETURM I END

REFERENCES 1. Lee, T. H., "Plasma Physics and the Interruption of an Electric Arc," Proc. IEEE, Vol. 57, No. 3, p. 307, 1969. 2. Rieder, W., "Circuit Breakers-Physical and Engineering Problems, Pts. I, II, and III, " IEEE Spectrum, July, Aug., Sept., 1970. 3. Hertz, W., Motschmann, H., and Wittel, H., "Investigations of the Properties of SF6 as an Arc Quenching Medium," Proc. IEEE, Vol. 59, No. 4, p. 485, 1971. 4. Frost, L. S. and Liebermann, R.W., "Composition and Transport Properties of SF6 and Their Use in a Simplified Enthalpy Flow Model, " Proc. IEEE, Vol. 59, No. 4, p. 474, 1971. 5. Cassie, A. M., "Arc Rupture and Circuit Severity: A New Theory, " Conference Internationale des Grands Reseaux Electriques a Haute Tension, (Paris, France) Tech. Rept. 102, 1939. 6. Anderson, R.W., "Nonstationary Electric Arcs, " Thesis, University of Michigan, 1972. 7. Phillips, R. L., "Induced Radial Velocity in Nonstationary Electric Arcs," Proc. IEEE, Vol. 59, No. 4, p. 466, 1971. 8. Stine, H.A. and Watson, V.R., "The Theoretical Enthalpy Distribution of Air in Steady Flow Along the Axis of a DirectCurrent Electric Arc," NASA TN D-1331, 1962. 9. Watson, V.R. and Pegot, E.B., "Numerical Calculations for the Characteristics of a Gas Flowing Axially Through a Constricted Arc, " NASA TN D-4042, 1967. 10. Ragaller, K., Schneider, W.R., and Hermann, W., "A Special Transformation of the Differential Equations Describing Blown Arcs," ZAMP, Vol. 22, p. 920, 1971. 11. Kircos, L., "Piecewise Uniform One-Dimensional Supersonic Flow with Heat Addition, " Thesis, University of Michigan, 1970. 243

244 12. Swanson, B.W. and Roidt, P.M., "Some Numerical Solutions of the Boundary Layer Equations for an SF6 Arc," Proc. IEEE, Vol. 59, No. 4, p. 493, 1971. 13. Sichel, M., "Computer Experiments Related to Chemical Propulsion, " AIAA J., Vol. 5, p. 1537, 1967. 14. Emmons, H.W., "Critique of Numerical Modeling of FluidMechanics Phenomena, " Annual Rev. of Fluid Mechanics, Vol. 2, Annual Revs. Inc., Palo Alto, California, 1970. 15. Cheng, S. I., "Numerical Integration of Navier-Stokes Equations, " AIAA Paper No. 70-2, 1970. 16. Richardson, D.J., "The Solution of Two-Dimensional Hydrodynamic Equations by the Method of Characteristics, " Methods in Computational Physics, Vol. 3, p. 295, 1964. 17. Richtmyer, R.D. and Morton, K.W., Difference Methods for Initial-Value Problems, Interscience, New York, 1967. 18. Forsythe, G. E. and Wasow, W.R., Finite-Difference Methods for Partial Differential Equations, John Wiley and Sons, New York, 1960. 19. Moretti, G., "The Choice of a Time-Dependent Technique in Gas Dynamics, " PIBAL Rept. No. 69-26, 1969. 20. O'Brien, G.G., Hymon, M.A., and Kaplan, S., "A Study of the Numerical Solution of Partial Differential Equations," J. of Mathematics and Physics, p. 223, January 1951. 21. von Neumann, J. and Richtmyer, R.D., "A Method for the Numerical Calculation of Hydrodynamic Shocks," J. Appl. Physics, Vol. 21, p. 232, 1950. 22. Moretti, G., "A Critical Analysis of Numerical Techniques: The Piston-Driven Inviscid Flow," PIBAL Rept. No. 69-25, 1969. 23. Chorin, A.J., "On the Convergence of Discrete Approximations to the Navier-Stokes Equations, " AEC Res. Dev. Rept. No. NYO 1480-106, 1968.

245 24. Emery, A. F., "An Evaluation of Several Differencing Methods for Inviscid Fluid Flow Problems, " J. Comp. Phys., Vol. 2, p. 306, 1968. 25. Lax, P. and Wendroff, B., "Systems of Conservation Laws, Comm. Pure and Appl. Math., Vol. XIII, p. 217, 1960. 26. Burstein, S.Z., "Numerical Methods in Multidimensional Shocked Flows," AIAA J., Vol. 2, No. 12, p. 2111, 1964. 27. Serra, R.A., "The Determination of Internal Gas Flows By a Transient Numerical Technique, " Thesis, Rensselaer Polytechnic Institute, 1970. 28. Gourlay, A.R. and Morris, J., "Finite-Difference Methods for Nonlinear Hyperbolic Systems, " Math. Comp., Vol. 22, No. 102, p. 28, 1968. 29. Gourlay, A.R. and Morris, J., "Finite-Difference Methods for Nonlinear Hyperbolic Systems. II, " Math. Comp., Vol. 22, No. 103, p. 549, 1968. 30. Gourlay, A.R. and Morris, J., "A Multistep Formulation of the Optimized Lax-Wendroff Method for Nonlinear Hyperbolic Systems in Two Space Variables," Math. Comp., Vol. 22, No. 104, p. 715, 1968. 31. Gourlay, A.R. and Morris, J., "On the Comparison of Multistep Formulations of the Optimized Lax-Wendroff Method for Nonlinear Hyperbolic Systems in Two Space Variables," J. Comp. Phys., Vol. 5, p. 229, 1970. 32. Crocco, L., "A Suggestion for the Numerical Solution of the Steady Navier-Stokes Equations, " AIAA J., Vol. 3, p. 1824, 1965. 33. Callens, E. E., "A Time-Dependent Approach to Fluid Mechanical Phenomenology," AIAA Paper No. 70-46, 1970. 34. Crocco, L., "One-Dimensional Treatment of Steady Gas Dynamics," High Speed Aerodynamics and Jet Propulsion, Vol. 3, Princeton Univ. Press, Princeton, New Jersey, pp. 64-349., 1958.

246 35. Shapiro, A.H., The Dynamics and Thermodynamics of Compressible Fluid Flow, Vols. I and I, The Ronald Press Co., New York, 1958. 36. Sentman, L. H., Fundamentals of Quasi-One-Dimensional Flow, University of Illinois Dept. of Aeronautical and Astronautical Engrg., Urbana, Illinois, 1969. 37. Mac Cormack, R.W., "The Effect of Viscosity in Hypervelocity Impact Cratering," AIAA Paper No. 69-354, 1969. 38. Anderson, J.D., "Numerical Experiments Associated With Gas Dynamic Lasers," NOLTR 70-198, 1970. 39. Burstein, S.Z., "Finite-Difference Calculations for Hydrodynamic Flows Containing Discontinuities, " J. Comp. Physics, Vol. 1, p. 198, 1967. 40. Moretti, G., "Transient and Asymptotically Steady Flow of an Inviscid, Compressible Gas Past a Circular Cylinder, " PIBAL Rept. No. 70-20, 1970. 41. Moretti, G. and Pandolfi, M., "Entropy Layers, " PIBAL Rept. No. 71-33, 1971. 42. Moretti, G. and Pandolfi, M., "Analysis of the Inviscid Flow About a Yawed Cone. Preliminary Studies," PIBAL Rept. No. 72-18, 1972. 43. Moretti, G., Grossman, B., and Marconi, F., "A Complete Numerical Technique for the Calculation of Three-Dimensional Inviscid Supersonic Flows," AIAA Paper No. 72-192, 1972. 44. Kutler, P. and Lomax, H., "Shock-Capturing, Finite-Difference Approach to Supersonic Flows," J. Spacecraft, Vol. 8, No. 12, p. 1175, 1971. 45. Thomas, P.D., Vinokur, M., Bastionon, R.A., and Conti, R. J., "Numerical Solution for Three-Dimensional Inviscid Supersonic Flow," AIAA J., Vol. 10, No. 7, 1972. 46. MacCormack, R.W., "Numerical Solution of the Interaction of a Shock Wave with a Laminar Boundary Layer," Lecture Notes in Physics, Vol. 8, Springer-Verlag, p. 151, 1971.

247 47. MacCormack, R.W. and Paullay, A.J., "Computational Efficiency Achieved By Time Splitting of Finite-Difference Operators, "AIAA Paper No. 72-154, 1972. 48. Moretti, G., "Importance of Boundary Conditions in the Numerical Treatment of Hyperbolic Equations," Phys. Fluids Suppl., Vol. 12, No. 12, p. 13, 1969. 49. Kutler, P., Lomax, H., and Warming, R. F., "Computation of Space Shuttle Flow Fields Using Noncentered Finite-Difference Schemes," AIAA Paper No. 72-193, 1972. 50. Migdal, D., Klein, K., and Moretti, G., "Time-Dependent Calculations for Transonic Nozzle Flow, " AIAA J., Vol. 7, No. 2, p. 372, 1969. 51. Anderson, J.D., "A Time-Dependent Analysis for QuasiOne-Dimensional Nozzle Flows with Vibrational and Chemical Nonequilibrium," NOLTR 69-52, 1969. 52. Courant, R., Friedrichs, K., and Lewy, H., "On the Partial Difference Equations of Mathematical Physics, " IBM J., March 1967. 53. Rubin, E.L. and Burstein, S.Z., "Difference Methods for the Inviscid and Viscous Equations of a Compressible Gas," J. Comp. Phys., Vol. 2, p. 178, 1967. 54. Isaacson, E. and Keller, H. B., Analysis of Numeral Methods, John Wiley and Sons, New York, 1968. 55. Ames, W. F., Numerical Methods for Partial Differential Equations, Barnes and Noble, New York, 1969. 56. Liepmann, H.W. and Roshko, A., Elements of Gasdynamics, John Wiley and Sons, New York, 1957. 57. Akima, H., "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures," J. of Assoc. for Computing Machinery, Vol. 17, No. 4, p. 589, 1970. 58. Zueckler, K., "The Effect on the Gas Flow of an Electric Arc Burning in a Nozzle, " Royal Aircraft Estab. Translation No. 973, Ministry of Aviation (London), 1961.

248 59. Zueckler, K., "Uber wechselseitige Beeinflussung von GasstroSmungen und Bogen in Gasstronungs schaltern, Z. Angew. Phys., Vol. 22, p. 155, 1967. 60. Zueckler, K., "Beitrag zum Rfichstauproblem in Hochspannungs-Leistungsschaltein, " Elektrotech. Z. Ausg. A., Vol. 90, p. 711, 1969. 61. Frie, W., "Berechung der Gaszusammensetzung und der Materialfunktionen von SF6, " Z. fur Physik, Vol. 201, p. 269, 1967. 62. Phillips, R. L., "One-Dimensional Flow in Variable Area Nozzles With Heat Addition, " Dept. of Aerospace Engrg. Note, University of Michigan, 1971. 63. Penner, S.S. and Davidor, W., "Parametric Solution of the One-Dimensional Flow Equations With Heat Addition," Proc. of the 1967 Heat Transfer and Fluid Mechanics Institute, Stanford Univ. Press, Stanford, California, 1967. 64. Kogelschatz, U. and Schade, E., "Experimental Investigation of a High Pressure Arc in a Strong Axial Gas Flow, " Proc. 10th Ionization Conf., Oxford, England, 1971. 65. Urbanek, J., "The Time Constant of High Voltage Circuit Breaker Arcs Before Current Zero," Proc. IEEE, Vol. 59, No. 4, p. 502, 1971. 66. Siddons, D.J. and Heron, K., "High Speed Photographic Study of the Extinction of an SF6 Arc," presented at the Nottingham Power Engineering Conf., England, 1969. 67. Phillips, R. L., "Time Constants for Nonstationary Arcs," Z. fir Physik, Vol. 211, p. 113, 1968. 68. Kurzrock, J.W., "Exact Numerical Solutions of the TimeDependent Compressible Navier-Stokes Equations," CAL Rept. No. AG-2026-W-1, 1966. 69. Saunders, L. M., "Numerical Solution of the Flowfield in the Throat Region of a Nozzle, " NASA CR 82601, 1966.

249 70. Armitage, J. V., "The Lax-Wendroff Method Applied to Axial Symmetric Swirl Flow," AD 657 323, Blanch Anniversary Volume, p. 3, 1967. 71. Lapidus, A., "Computation of Radially Symmetric Shocked Flows," J. Comp. Physics, Vol. 8, p. 106, 1971. 72. Phillips, R. L., "The Behavior of Dynamic Electric Arcs, " Thesis, University of Michigan, 1964. 73. Cuffel, R. F.,Back, L.H., and Massier, P.F., "Transonic Flowfield in a Supersonic Nozzle with Small Throat Radius of Curvature, " AIAA J., Vol. 7, No. 7, p. 1364, 1969. 74. Back, L. H., Massier, P. F., and Cuffel, R. F., "Flow Phenomena and Convective Heat Transfer in a Conical Supersonic Nozzle," J. Spacecraft, Vol. 4, No. 8, p. 1040, 1967. 75. Back, L.H. and Cuffel, R. F., "Detection of Oblique Shocks in a Conical Nozzle with a Circular-Arc Throat," AIAA J., Vol. 4, No. 12, p. 2219, 1966. 76. Hopkins, D. F. and Hill, D. E., "Effect of Small Radius of Curvature on Transonic Flow in Axisymmetric Nozzles," AIAA J., Vol. 4, No. 8, p. 1337, 1966. 77. Kliegel, J. R. and Levine, J. N., "Transonic Flow in Small Throat Radius of Curvature Nozzles, " AIAA J., Vol. 7, No. 7, p. 1375, 1969. 78. Shelton, S.V., work done at Jet Propulsion Laboratory, Pasadena, California; see Reference 73 above. 79. Migdal., D and Kosson, R., "Shock Predictions in Conical Nozzles," AIAA J., Vol. 3, No. 8, p. 1554, 1965. 80. Hermann, W., as related to R.W. Anderson, private communication.

UNIVERSITY OF MICHIGAN 3 9015 03527 3286