Computer experiments may be defined as theoretical analyses in which the computer plays the central role in evaluating the consequences of the mathematical model of the process or phenomenon under investigation. The use of such computer experiments to study problems related to chemical propulsion forms the subject of the present survey. Several different types of experiment are considered. In certain problems the physics is well known but the mathematical equations are intractable. Then high speed computers are used to evaluate the consequences of the equations. In a second class of problems the detailed physics is unknown and computer experiments are used to identify the most important physical factors. A third type of experiment consists of an extended parametric study of numerical solutions leading, finally, to a correlation formula for the numerical results. Computer experiments can also be classified according to computational technique, or according to subject area. With these classifications in mind the present survey begins with a brief discussion of computational methods. Then, computer experiments carried out in various subject areas related to chemical propulsion such as chemical kinetics, fluid dynamics, transport processes, shock and detonation phenomena, combustion, etc., are surveyed.
The use of high speed computers has played an increasingly important role in the analysis of problems relating to chemical propulsion. In many cases the computer is merely used for the straightforward evaluation of the analytical expressions resulting from theoretical analysis; however, in an increasing number of problems the computer has played a central and indispensable role being used as an apparatus to test theoretical models of the phenomenon under study. The use of such computer experiments" to study problems related to chemical propulsion forms the subject of the present survey.
It is proposed that a computer experiment be defined as a theoretical analysis in which the use of the computer plays the central role in evaluating the consequences of the mathematical model of the physical phenomenon or process under investigation. All such experiments involve first a formulation of the problem; second the development of a numerical scheme for solving the equations describing the process in question; then, the actual computations, which would generally be impossible without modern high speed computers; and finally an evaluation and interpretation of the results. The boundary between what can be classed as a computer experiment and a straightforward calculation will be somewhat arbitrary.
Computer experiments take several different forms. In certain cases the physics of the problem is completely under stood; however, the mathematical equations describing the phenomenon in question are intractable to analytical solution without the use of severe assumptions. Then high speed computers are used to numerically evaluate the consequences of the equations. The calculations of viscous flows over obstacles made by Fromm(1) and others, and resulting in movies of the computed flow solutions provide an outstanding example of this type of experiment.
In other problems the detailed physics underlying a given phenomenon are unknown, or the value of key physical constants may be uncertain. Thus, for example, reaction rate constant values needed to compute nonequilibrium gas flows may be known only to an accuracy of several orders of magnitude. Such problems are sometimes approached by developing a mathematical model including all possible factors of importance. Computer experiments or extensive calculations are then made to establish which factors are important and which can be ignored, Duff's (2) calculations of detonation structure which were carried out for several possible values of the key reaction rate constant provides an early example of this type of experiment.
A third type of experiment can also be distinquished, and consists of a parametric study of numerical solutions that finally leads to an empirical correlation formula describing the numerical results. The stagnation point heat transfer calculations of Fay and Ridde11 (3) may be considered an example of this type of experiment. This classification of computer experiments into three categories, while convenient for discussion, is highly arbitrary. It is no doubt possible to cite examples of computer experiments which fit into none of the above categories.
Computer experiments cannot replace either laboratory experiments or analytical solutions, rather these three approaches to the study of physical phenomena supplement each other. An analytical solution, where available, is usually far more valuable than computer solutions for analytical solutions contain extensive information about any given phenomenon in a very condensed form and show exactly how various factors affect the phenomenon under study. Analytical solutions can be used to calibrate a particular numerical scheme which then can be applied to more complicated problems for which analytical solutions cannot be found. Often the results of a computer experiment will provide a guide to the development of an analytical theory, for example, in indicating what simplifying assumptions can reasonably be made in a given situation. Both analytical solutions and computer experiments may suggest meaningful laboratory experiments. Since a computer experiment is only as good as the physical model on which it is based it must be emphasized that the study of new physical phenomena will always require experimentation in the laboratory, though computer results can provide a powerful guide to the form such experimentation should take and can result in a drastic decrease in the number of experiments required.
Although the present full blown use of computer experiments has had to await the development of the modern high speed computer, the use of computer or computational experiments is a technique which has been in use for some time. Thus for example the result that the recovery factor for a laminar boundary layer on flat plate is given by Pr1/2 represents the correlation of numerical computations of recovery factors over a range of Prandtl numbers(4), (5). This work dates back to Polhausen's 1921 paper (6) on boundary layer heat transfer. Crocco's (7) laborious computations of flat plate compressible boundary layer skin friction coefficient and the variation of enthalpy with velocity, for different relations between viscosity and temperature, although the computations were made by hand, is certainly a relatively early example of a computational experiment. While the computational experiment is thus not a new technique a completely new element has been introduced by the evolution of the modern computer with its high speed, large storage capacity and versatility. Thus the computational experiment has changed from a laborious brute force procedure to a powerful tool for the theoretical analysis of complex problems and phenomena.
As mentioned above, a key step in all computer experiments is the development of a suitable numerical scheme for solving the equations describing the phenomenon of interest. Particularly in the case of partial differential equations numerical analysis is still essentially an art rather than a science. The stability and convergence of finite difference schemes can only be rigorously established in a few special cases. The choice of the best finite difference scheme for a particular problem is very much a matter of intuition and trial and error. Thus, for instance, changing from forward to central finite differences or vise versa may change a completely unstable finite difference scheme to a stable one.
Extensive and sophisticated mathematical development will be required to place numerical analysis on more rational basis.
Computer experiments can also be classified according to the computational technique used or according to subject area. With these classifications in mind the present survey begins with a brief description of computational methods followed by a discussion of specific computer experiments carried out in various subject areas related to chemical propulsion.
A complete treatise on numerical methods is clearly beyond the scope of the present survey and an extensive discussion of this subject may be found in books such as those of Fox (8), Richtmeyor(9). Ralston and Wilf (10)), and Todd (11) to name only a few. However, some of the main methods which have been used in computer experiments will be outlined below.
Many problems require the numerical integration of one or a set of ordinary differential equations with initial conditions or two point boundary conditions. The Runge-Kutta procedure with the modification by Gill(10) or predictor corrector methods such as those described by Milne(12) and Fox(8), for example, are frequently used for this purpose.
In the case of partial differential equations the numerical procedure employed depends on whether the equations are elliptic, parabolic, or hyperbolic. The boundary conditions for elliptic equations must be specified about a closed contour surrounding the domain of interest. Consequently, when elliptic partial differential equations are replaced by finite difference equations it is in effect necessary to solve a set of simultaneous algebraic equations, equal in number to the number of mesh points in the domain of interest. Direct solution of this set of equations is practical only if the number of mesh points is relatively small, say less than about 100. Usually iterative methods are used to solve the finite difference equations representing an elliptic system. Perhaps the simplest iteration scheme is a variation of Liebmann's method, as used for example by Fromm (13). This method is readily illustrated in the case of the LaPlace equation, ∇2ψ=0 for which a commonly used finite difference scheme is
If superscripts are used to indicate the order of iteration then Liebmann's method consists in using the kth iterates to compute the values of the (k + 1) th iterate at each point so that
while the modification used by Fromm(13) and others is to use the most recent iterates at neighboring points so that, for example, if the net of points is scanned from top to bottom and left to right
Important considerations in the choice of an iteration method are the rate of convergence and the computational complexity of each step. Generally the simple iteration method described above converges toward the actual solution very slowly; however, the computations are very simple. A number of much more rapidly converging iteration schemes such as, for example, the under and over relaxation methods and the Paceman-Rachford alternating direction procedure (11) are also available, but of course each iteration step is more complicated. As indicated by Fromm(13) in certain cases a slower but simpler iteration scheme may be preferable. In applying any iteration scheme to elliptic finite difference equations or to, what is in effect the same thing, the inversion of a large matrix, there is no guarantee that the iteration will ever converge to the solution of the finite difference equation; however, in certain cases convergence can be proved, An extensive treatment of the convergence of iteration methods of solving finite difference equations is, for example, given by Varga(14).
In the case of parabolic equations the situation is somewhat simpler. Only initial and boundary conditions are required and therefore it is possible to use forward marching techniques to solve the parabolic finite difference equations. Explicit finite difference methods are perhaps the simplest to use. Thus, for example, the heat equation ψxx−ψt=0 can be written in the finite difference form
Then knowing the solution on line j the values on line j+1 can be found from the explicit formula
This explicit finite difference scheme will only be stable if k/h2 < 1/2, a condition which often results in exceedingly small timewise step sizes. Implicit finite difference schemes, though more complicated, avoid this difficulty. Thus, for example, the heat conduction equation can be replaced by the finite difference equation
and in this case values of ψ on the new line j+1 occur implicitly in a set of simultaneous algebraic equations for all the unknown values of ψ on the j+1 line. With a= 1/2 equation (6) reduces to the well known Crank-Nicolson formula (8) and the forward marching procedure will be unconditionally stable for all values of k/h2.
Hyperbolic equations have long been solved numerically using the method of characteristics. When carried out by hand, characteristic calculations tend to be long and laborious; however, recently it has become possible to program these calculations for digital computers. Hyperbolic equations can also be solved employing both explicit and implicit finite difference methods (8), (9).
While the above are perhaps the classical approaches to the numerical integration of partial differential equations many other schemes, some of which are described below, have been used.
A technique often used in two dimensional probems is to use finite differences in only one direction or in effect to break the field of interest into a finite number of strips. This procedure reduces the partial differential equation to a set of ordinary differential equations for each of the strips. A simple example of this method is provided by the heat conduction equation, ψxx−ψt=0. Taking finite differences only with respect to x this equation becomes
or considering ψ as a function t for each strip corresponding to a particular xi.
so that the heat equation reduces to a set of simultaneous ordinary differential equations for the ψi with initial conditions applied at t=0. Analogue computers can be used to integrate this system of equations if extreme accuracy is not required. The method of integral relations (15), which has been used extensively in the Soviet Union and is coming into increasing use in the United States, is another example of this technique. In obtaining numerical solutions of the boundary layer equations Smith and Clutter (16) used finite differences along the boundary and ordinary differential equations across the boundary layers. Two point boundary conditions made it necessary to solve the resultant set of ordinary differential equations iteratively.
Many problems of steady fluid mechanics require the numerical solution of elliptic partial differential equations with al! the difficulties that this entails. However, as pointed out by Crocco (17) steady flows may be considered the asymptotic limit of some time unsteady flow, and taking advantage of this fact steady problems can be approached by studying the long time behavior of a suitable unsteady flow. Since the equation for unsteady flow are usually of parabolic or hyperbolic character the much simpler finite difference marching procedures, described above, can then be applied to the investigation of steady problems. Crocco (17) used this technique in the study of shock wave development in converging diverging nozzles while Bohachevsky and Rubin (18) investigated blunt body flow using this time unsteady approach.
In certain problems, particularly when large distortions of a fluid or solid are involved, approaches which combine the advantages of the Eulerian and Lagrangian points of view are very effective. The Particle in Cell or PIC method developed by Evans and Harlow (19), (20) is one example of such a method. The fluid field is divided into a number of Eulerian cells with mean properties ascribed to each. In addition a Lagrangian mesh of particles is established. The computation proceeds in two steps. First the particles are presumed fixed and the conservation equations are used to compute new cell wise quantities, a procedure which is equivalent to neglecting the convective terms. Then, using suitably interpolated values of the velocity the particle motion is determined after which new average cell properties are computed. The output consists of the motion of the particles as a function of time. Evans and Harlow(21) applied this method to a study of the interaction between a shock wave and a bubble.
The difficulty of treating unsteady flows with shock waves has led to methods in which artificial viscosity terms are introduced into the equations of motion (8) either explicitly (22) or by the choice of a suitable finite difference method(23). Then shock waves, or rather shock layers across which the Rankine-Hugoniot conditions hold occur naturally in the solution. The artificial viscosity terms are proportional to the gradients of the various flow parameters and so have a negligible effect on the regions outside the shock layers; however, these methods do not provide an accurate description of the actual shock structure.
The computation of certain complex processes, as for example the determination of molecular dissociation rates from the detailed motion of the atoms constituting the molecules, can be reduced to tractable proportions by using a statistical method depending upon the random sampling of the results of a particular mathematical model. These statistical methods, which vary widely with the particular application are collectively known as the Monte Carlo method. This method is discussed by Fox (8), and extensive references are given by Todd (11). The use of the Monte Carlo method to evaluate a multiple integral provides a simple illustration of how this method works (24). Consider the evaluation of the integral
where x is the vector (ξ1,ξ2... ξn), and dv is a volume element in ξi space. It is supposed that the integral (9) is normalized such that 0 ≤ φ(x) ≤ 1, and 0 ≤ ξi ≤ 1. The integral (9) represents a hypervolume in an (n+1) dimensional space defined by the coordinates ξ1, ... ξn,η. If points are drawn at random from this (n+1) dimensional space, then Φ is the probability that any one of the random points has η satisfying η ≤ φ(x) where x=(ξ1..ξn) are the other n coordinates of this point. If N drawings are made and N' points satisfy the above inequality then N'/N provides an estimate of Φ.
Finite difference equations for initial value problems are sometimes unstable, that is, given certain initial conditions, values begin oscillating wildly at some stage of the calculation and bare no relation to the solution of the partial differential equation of interest. In all finite difference calculations the use of finite intervals also leads to a certain truncation error as compared to the solution of the continuous equation. Questions of stability and truncation error are discussed in some detail by Fox (8) and Richtmeyer (9), A simple method of developing stability criteria is to study the propagation of a single row of errors of the form eiβx. Then the error at any point of the mesh will propagate as eαjkeβih where h and k are the x and y grid spacings respectively. The requirement for stability, obtained by substituting this result in the finite difference equation, is then that eαk< 1. This scheme will only work in the case of linear finite difference equations. In dealing; with non linear equations a common procedure (13) is to study a linearized version of the equations in order to obtain at least a guide as to the requirements for stability.
The difficulties inherent in some differential equations can sometimes be removed by recasting the problem in question into an integral equation, and the numerical solution of integral equations is discussed by Fox (8).
It must also be mentioned that various output techniques which directly convert digital computer results into graphs or maps have played a key role in increasing the effectiveness of modern computers.
Computer experiments carried out in various subject areas related to chemical propulsion are discussed below. Chemical propulsion is such a broad field that it has not been possible to include all related computer experiments. However, it is not the objective of the present survey to provide a complete catalog of computer studies, but rather to show how various computational methods have been applied in a number of areas related to chemical propulsion and, perhaps, to provide a guide to the future use of such computer experiments.
Converging-diverging nozzles are used to convert thermal to directed kinetic energy in rockets and air breathing propulsion devices as well as in supersonic and hypersonic wind tunnels. With the high energy flows encountered in rocket and hypersonic wind tunnel nozzles the effect of chemical reactions occurring within the nozzle has become increasingly important. While results for the extreme cases of frozen and equilibrium chemistry have been available for some time these are not adequate in many cases. In the case of rocket and SCRAM (Supersonic Combustion Ram Jet) nozzles it is important to know the extent of recombination of dissociated species as this can have an important effect on the specific impulse attained, In hypersonic nozzles the free stream velocity and composition must be known for the proper interpretation of test data. Solution of these problems requires an analysis of nonequilibrium nozzle flow.
For nozzle flows involving only one reaction approximate nonequilibrium flow soluUons based on the sudden freezing approximation developed by Bray (25) and Penner(26) and discussed extensively by Cheng (27) can be used.
For more complex reaction systems numerical nonequilibrium computations are required, and due to their complexity the high speed computer plays a central role in these computations. It is therefore felt that such computations can be considered to be computer experiments in the sense of the definition given above. Since even these computations are expensive and time consuming a few carefully selected equilibrium computations or computer experiments have often been used to develop empirical equations for use in an approximate analysis of nonequilibrium flow.
In these computations the flow is generally assumed to be one dimensional and inviscid. The resultant system of ordinary differential equations is usually solved by marching forward from some point upstream or downstream of the nozzle throat using the Runge-Katta procedure for example. The rate equations become singular at the combustion chamber or nozzle reservoir where the system is in chemical equilibrium; therefore, special perturbation expansions as developed by Hall and Russ (28) and Vincenti (29) and Emanuel (30) must be used in this singular region. Another problem is encountered in continuing the solution through the saddle point singularity at the throat. This problem has been discussed by Emanuel (31). An indirect approach to this difficulty has been used by Eschenroeder et al (32) who specified nozzle mass flow and then computed the nozzle area variation A(x) to some point beyond the throat while Sarli et al (33) used extrapolation to carry their computations through the throat region. In nonequilibrium flow calculations, particularly near equilibrium, very small step sizes must be used in the numerical integration in order to avoid instability. A numerical technique based upon local linearization of the rate equations which overcomes this difficulty has been proposed by Moretti (34) and resulted in appreciable reduction in computation time for the test cases considered.
Extensive computations of the nonequilibrium nozzle flow of air have been carried out by Eschenroeder et a1 (32), Vincenti (35), Lordi and Mates (36), (37), and Emanuel and Vincenti (74). Eschenroeder et al. considered the species O2, N2, Ar, O, N, and NO governed by fourteen reactions. Nonequilibrium expansions were compared to infinite reaction rate expansions in order to assess the influence of nonequilibrium on various flow parameters and species concentrations. An interesting result of these computations is that N atoms do not freeze out as rapidly as O atoms because of the shuffle reaction N + O2= NO+O, O+2=NO + N. Lordi and Mates (36), (37) extended these nozzle computations to reservoir temperatures up to 15,000K° and reservoir densities up to 100 times atmospheric. The effect of reservoir conditions upon nonequilibrium effects was clearly shown by these computations. The results of the calculations were used to determine empirical relations for nonequilibrium effects.
By studying the results of exact calculations for the nonequilibrium expansion of O2 Hall and Russo(38) developed a freezing length criterion based on the local value of the chemical relaxation length. This freezing criterion was then applied to the approximate computation of nonequilibrium H2 expansions. In particular the effect of the recombination rate constant upon the vacuum specific impulse was explored. These calculations were extended by Hall, Eschenroeder, and Klein (39), and provide an excellent illustration of how computer experiments can be used to develop rapid approximate methods to deal with complex problems. The possibility of using carbon to catalyze atomic H recombination in nozzles was later explored by Eschenroeder and Lordi (40) who made an extensive numerical study of the effects of rate constants, carbon hydrogen ratios and reservoir conditions upon recombination.
The expansion of hydrogen air combustion products was investigated by Sarli et al (38) since the degree of recombination at the exit of SCRAM jet nozzle will have an important influence upon the thrust. Equilibrium, frozen and nonequilibrium solutions were compared, typical results being shown in Fig. 1. Westenborg and Farvin (41) carried out extensive investigations of nonequilibrium expansions involving the species combinations H2, H22O, CO, CO2, H, OH and H2, H2O, O2, H, OH, O. On the basis of the exact computations various approximation schemes sometimes used in nozzle calculations were assessed, and the effects of changes in some of the critical rate constants were examined. This typical computer experiment also made it possible to assess the importance of the various reactions occurring in the nozzle. A similar comparison between approximate nonequilibrium nozzle flow theories and exact computations was carried out by Wilde (42) for the chemical systems UDM H/N2H4-N2O4 and CH2-O2. With increasing length nonequilibrium losses in a rocket or ram jet nozzle should decrease; however, frictional losses will increase. On the basis of numerical calculations (or a computer experiment) Sarli et a1 (43) established that the increase in friction will generally outweigh any effects of increased recombination.
The combustion products of many solid propellants contain solid particles; hence, it is extremely important to determine the effect of the particles upon nozzle specific impulse. As indicated in the review by Hoglund analytical results for gas-particle flow through nozzles have been obtained only for very special cases as for example for zero particle lag, or for a fixed ratio of particle to fluid velocity. Inclusion of realistic particle drag and heat transfer laws as well as real air properties has required computer solutions so that computer experiments have played a key role in the study of gas-particle flows. One dimensional flows with particles of uniform size have, for example, been investigated by Kuegel and Nickerson (45) and Bailey et a1.(46) who have computed particle velocity and temperature lags as well as nozzle specific impulse for different particle sizes and different nozzle contours, As in the nonequilibrium computations, considerable finesse is required to continue the calculations through the nozzle throat (44). Kliegel (47) has also computed nozzle flows with nonuniform particles having a logarithmic size distribution and with drag coefficients and Nusselt numbers including rarefaction corrections. The method of characteristics was used to compute axially symmetric gas-particle nozzle flows taking the inability of the particles to follow the streamlines into account. Results were compared to the analytical constant fractional velocity and temperature lag predictions. Extensive parametric computer studies of gas-particle nozzle flow to study the effects of nozzle cone angle, expansion ratio, throat radius, particle size and particle flow rate ratio upon specific impulse were made by Hoffman and Lorenc (48) for aluminum oxide particles and a combustion chamber temperature ranging from 5500°R to 6500°R.
Condensation phenomena are of considerable importance in hypersonic wind tunnels and in expansions of metallic vapors involved in various nuclear space propulsion devices. Griffin and Sherman(49) used a computer experiment to study the expansions of pure vapors in two dimensional nozzles. The equations of spontaneous nucleation theory were combined with the equations of steady one dimensional adiabatic flow in order to study the onset and rate of condensation. These computations made it possible to follow in detail nucleation rates and particle sizes in the condensation region. Although there is considerable uncertainty about the nucleation theory, the excellent agreement between the computed results and measurements for pure nitrogen were quite encouraging. By varying the accomodation coefficient, heat of vaporization and surface tension it was possible to assess the importance of these parameters, at least with respect to the nucleation model employed. It was found that small variations in surface tension could have a large effect upon the onset of condensation. An important conclusion was that the lack of accuracy in the value of surface tension leads to the largest source of error in the calculations. It was also possible to study the effect of the nozzle expansion rate upon the condensation process.
A key source of uncertainty in nonequilibrium flow computations such as those described above are the values of the reaction rate constants. Often these constants are known only to an order of magnitude. An interesting recent development in the study of reaction rates has been the use of computer experiments to compute rate constants and to evaluate various reaction rate theories (49).
Bunker's study of the dissociation of triatomic molecules (50) is a typical example of such an experiment. A triatomic molecule with sufficient energy will dissociate when one of the bonds gains sufficient energy to break. When and whether or not this happens will depend on the initial configuration of the molecule, the total energy and the nature of the interatomic forces. By studying the detailed atomic motions for many initial configurations Bunker (50) was able to statistically compute a normalized distribution of molecular lifetimes for molecules in an equilibrium distribution as a function of energy. The first step in this analysis is to derive the equations of motion of the system consisting of the three atoms connected by the two bonds r1 and r2 at angle α to each other as sketched below:
It is of course possible to use various interatomic potential functions in specifying the nature of these bonds. Bunker for example investigated harmonic and Morse type potential functions. For given initial condition and a total energy H these equations are integrated using a Runge-Kutta procedure and the atomic motions are followed to see whether either of the two bonds is stretched sufficiently to consider the molecule dissociated. A key part of the analysis is the choice of the initial conditions used in these calculations, and a Monte Carlo procedure is employed in that initial distributions are chosen at random from the hyper surface H = constant in the phase space. A Monte Carlo procedure is also used to determine S(H), the relative volume of the phase space between H and H + dH. The rate constants obtained by determining what proportion of the initial conditions used result in dissociation were averaged over the phase space by using the energy weighting factor S(H) exp (-H/kT), Bunker (50) integrated 90,026 initial configurations of which 10,559 resulted in dissociation. On the basis of this experiment Bunker was able to discuss his results in the light of various reaction rate theories and also was able to determine how various assumptions regarding the molecular bonds affected the results. The crux of Bunker's experiment and others mentioned below is to obtain statistical averages from the study a large number of randomly chosen patterns of atomic motion. Bunker's calculations were carried out on the Los Alamos STRETCH computer which is extremely fast, but even then required lengthy computation times.
Perhaps the earliest computer experiment of this type was the study of dissociation of H2 due to collisions with H atoms by Wall et a1 (51). The collinear H + H2 system was investigated under a great variety of conditions both with and without collisions. Thiele and Wilson (52) studied the dissociation of triatomic linear molecules considering the effects various models of the bond-stretching potential. Benson et al (53) investigated colinear collisions of diatomic molecules with atoms over a broad range of dynamical parameters as a model of the collisional excitation of molecular vibrations. Two dimensional kinematic computations have been carried out for exothermic reactions by Kuntz et a1 (54).
Computer experiments have also been applied to the difficult problem of particle-surface interaction. Rogers (55) studied a simple model of the interaction between a gas particle and a surface atom. The coupled surface atom-gas particle equations of motion were integrated on an analogue computer for many different values of the incoming particle velocity. By determining the energy change in each case and using an average weighted with the Boltzmann distribution it was possible to compute accomodation coefficients for this simple interact ion model. Particle-surface interaction experiments are also discussed by Gentry et al(56).
The study of various reaction kinetic schemes is another area where computer experiments have been extremely useful. Some of this work is discussed in the section on Combustion below.
The H2-O2 reaction is of interest both in rocket propulsion and recently with respect to supersonic combustion in SCRAM jets. The ignition delay which occurs when mixtures of hydrogen and oxygen are suddenly heated beyond their ignition temperature is of particular importance in SCRAM jets where excessive delays may make supersonic combustion impossible. As in the case of nonequilibrium nozzle flow the complexity of the reaction equations has necessitated considerable computer experimentation. Thus Belles and Lauver (57) for example, integrated the H2-O2 rate equations numerically to determine the ignition delay period. The results of these computations (or computer experiment) were used by Brokaw (58) to develop approximate analytical expressions for the ignition delay time, thus providing an example of how computer experiments can provide support for theoretical analyses.
Momtchiloffe et al (59) carried out a detailed computer study of ignition delays in hydrogen air flow systems in one dimensional supersonic flow. The reaction rate equations coupled with the one dimensional gas dynamic equations were integrated using the Runge-Kutta procedure with the modifications by Gill (10) starting from an arbitrarily specified initial temperature, equivalence ratio and flow velocity. The species O,H,O2, H2, OH, H2O, N2, N, and O governed by ten reaction rate equations were considered. Including the gas dynamic equations, fifteen simultaneous ordinary differential equations had to be integrated. Double precision computations were used in order to avoid instability near chemical equilibrium. Detailed species concentration profiles were obtained. Computed ignition delays were compared to measured values and the effects of initial temperature, pressure, equivalence ratio, and velocity upon ignition delay time and ignition delay length were determined. It is doubtful whether Momtchiloff et al (59) could have obtained the above information without extensive computer experimentation, and their results have also provided a guide to the development of approximate expressions for ignition delay.
Slutsky et al (60) carried out an even more ambitious supersonic combustion computer experiment for they studied the supersonic mixing and combustion of an axisymmetric premixed hydrogen-air mixture ignited by a pilot flame. The rate equation had to be solved simultaneously with boundary layer type equations for the species, momentum, and energy conservation with in the jet. Empirical values were used for the turbulent viscosity. Detailed temperature, concentration, and velocity profiles were obtained for the axisymmetric combustion region. Experimental and computed total temperature profiles were found to be in good agreement, a result which provides ground for some confidence in the computed results. This complicated computer experiment was made possible by using the local linearization technique of Moretti (34) to integrate the rate equations. Use of the more usual Runge- Kutta technique would undoubtedly have resulted in excessive computation time.
Computer experiments also have made advances possible in the understanding of solid propellant combustion. Hermance et a1 (61) studied the gas phase ignition of a condensed fuel by hot oxidizing gases. A simplified mathematical model using a single overall reaction rate and neglecting non linear convective terms in the conservation equation was used. The three simultaneous heat conduction type equations were integrated using a forward difference method similar to those described in the section on numerical methods above. The effects of fuel heat release, fuel oxidizer ratio and initial gas temperature on ignition time delay were examined and simplified solutions of the governing equations were compared to the exact numerical results. In another paper Hermance (62) carried out a computer study of the combustion of composite solid propellants. The mathematical model employed considered fuel pyrolysis, oxidizer decomposition, heterogeneous reaction between the fuel and the oxidizer in small fissures surrounding the oxidizer particles and gas phase combustion of the final decomposition products. The resultant set of implicit algebraic expressions were solved numerically for propellant burning rate for a wide range of parameters.
The overall performance of a solid propellant rocket during the starting process was investigated by DeSoto and Friedman(73). A finite difference calculation using the Peaceman-Rochoral (11) alternating direction method was used to compute the two dimensional heat conduction within the propellant grain, and by introducing an ignition temperature it was possible to determine the rate of increase of the burning surface. These propellant grain calculations were coupled to a quasi analytical solution of the gas flow within the propellant grain and in the combustion chamber in order to determine the chamber pressure variation. Effects of varying heat transfer rates across the boundary layer at the grain surface, and of various values of the radial and axial thermal conductivity of the nonisotropic propellant grain upon the rate of pressure increase were investigated. The effect of a burst diaphragm in the nozzle throat also was determined. In effect the idealized solid propellant rocket can be operated on the computer subject to varying material properties and external conditions.
Gilbert and Scala (63) have made an extended numerical study of the combustion and sublimation of graphite cones, spheres, and wedges. The boundary layer equations including the effects of chemical reactions and gas-solid interactions were integrated over a range of conditions subject to the assumption of local similarity and local thermochemical equilibrium. Empirical correlation. formulas were then derived to describe the numerical results. Gilbert and Scala's study is a typical and advanced example of a parametric computer experiment yielding empirical correlations of exact numerical results.
Computer experiments have yielded some remarkable results in the study of detonations in gases and solids; therefore, even though detonation theory may not have direct application to chemical propulsion it is felt that the interesting nature of these experiments warrants their discussion below. The experiments of White (64) and others showing that self supported detonations tend to be turbulent has led to studies of the stability of detonations such as that of Erpenbeck (65). In order to test this theory Fickett and Wood (68) used a computer experiment to study the initiation of a detonation due to the motion of a piston through a combustible gas. The numerical solution of the time dependent one dimensional reactive flow equations was obtained using the method of characteristics and treating shock waves as Rankine-Hugoniot discontinuities. Only a single step reaction of the form A-B, with activation energy E was considered. The long time behavior of overdriven detonations was studied for different values of E corresponding to unstable, neutrally stable, and stable regimes in the Erpenbeck theory (65). Typical shock pressure variations corresponding to a stable and an unstable detonation are shown in Fig. 2. The results of this experiment were in good agreement with Erpenbeck's theory, though of course the applicability of one dimensional calculations such as these to real detonations, where transverse waves may be quite important, is subject to question. Numerical computation of two dimensional instabilities is surely a subject for future research. Gilbert and Strehlow (67) studied the initiation of detonations behind reflected shock waves numerically using several models based on different laws of induction time and heat addition. This computer experiment was used to evaluate the ability of various models of the initiation process to represent experimental results.
Initiation of detonations in solid explosives due to an incident shock wave was studied numerically by Warner (68). The one dimensional unsteady flow equations were integrated using Lax's finite difference scheme (23). A modified Abel type equation of state and a pressure dependent burning law were used. Incident shock waves of various strengths were considered and were found either to decelerate or to accelerate to form a Chapman-Jouguet detonation. Typical pressure profile development as a function of time is shown in Fig. 3. Warner concluded from his experiment that use of a pressure dependent burning law was important in obtaining an adequate representation of explosive behavior. Shock waves interacting with bubbles or inhomogeneities in explosives give rise to hot spots which, if of sufficient size, can also initiate detonation. Mader (69) used a computer experiment to study both shock and hot spot initiation of detonation in solid explosives. The computations were performed by writing the equations of motion in a Lagrangian finite difference form. The hot spot model which is in part based on the computer experiment of Evans et al (70) on shock-bubble interactions was idealized as a spherically disturbed region with uniform temperature and density. Mader's computations made it possible to follow the complex processes which occur during hot spot initiation in detail. The critical hot spot radii determined by Mader(69) were in remarkable agreement with experimentally determined values. In a more ambitious study Mader (71) used the particle in cell method to simultaneously determine the interaction of a shock with a bubble and an aluminum cylinder in nitromethane and the accompanying initiation of detonation in the explosive. Hot spot initiation was also studied numerically by Friedman (72) in order to determine explosive impact sensitivity.
The study of steady state detonation structure has also been a subject of extensive theoretical and experimental investigations since the coupling of theory with experiment can yield extensive chemical kinetic information. Duff's numerical calculation of detonation structure (2) is an early and classic example of a computer experiment. The Von Neumann Zeldovich model of a shock followed by an inviscid reaction zone was chosen and the reaction rate equations were integrated using a Runge-Kutta procedure. The effect of various possible values of the rate constants for the key reactions H2+M =H+OH+M and O2+H=OH-O upon the density and species profiles through the detonation was investigated. Thus, for example, Duff found the ignition delay to be extremely sensitive to the chain branching reaction O2+H=OH+O. Typical density profiles found by Duff are shown in Fig. 4. Duff's numerical results have been an important touchstone for numerous theoretical and experimental studies of H2-O2 detonation structure. Numerical experimentation has also been used to investigate the coupling between the shock structure and the reaction zone by using the full Navier-Stokes equations to compute the detonation structure. This problem was first studied by Hirschfelder and Curtiss (75), (76) for irreversible unimolecular reactions, and their results are extended to the thermal decomposition of Ozone by Oppenheim and Rosciszewski (77). Complex singularities make integration of the ordinary differential equation describing coupled detonation structure quite difficult.
Campbell et al (78) used numerical experiments to study the effects of various reaction schemes, rate constants, and diffusion coefficients on the mass flow eigenvalue of a laminar flame and compared their results to various laminar flame theories. A computer experiment was also used by Kuehl (79) in analysis of the burning of aluminum and beryllium droplets. The model used assumed vapor phase burning at a diffusion flame front some distance from the particle surface. One of the results of the computations was to define the domain of oxidizer concentration and particle diameter in which the vapor phase burning model was valid. Burning rates over a wide range of ambient temperature oxidizer mole fractions, and particle diameter were obtained. Kuehl's results, which could not have been obtained without high speed computation, should provide a valuable guide for experiments designed to validate the vapor phase combustion model.
A close interaction between chemical and fluid dynamic processes is an essential element of many problems related to chemical propulsion, and it is in the field of fluid mechanics that some of the most exciting computer experiments have been developed, The study of the vortex development behind a blunt obstacle by Fromm and Harlow (80), (81) is perhaps one of the most remarkable of these computer experiments. A finite difference solution of the incompressible Navier-Stokes equations in the form
was obtained for the two dimensional viscous flow over a rectangular obstacle. A forward marching procedure was used to solve the vorticity equation for each succeeding time interval, but at each time interval the variation of the Liebmann method described in Section II above was used to iteratively solve the Poisson equation for the stream function. The fact that the computer output was converted directly into graphical form is an extremely important feature of this experiment and others like it, for detailed interpretation of the numerical data produced by such a calculation would be a hopeless task. In Fig. 5, taken from Ref. 80, the experimentally observed dye trails behind a circular cylinder are compared to the computed streaklines behind a flat plate. This amazing picture might perhaps be considered a monument to the computer experiment: The results of the computations were also presented in the form of a movie film showing the detailed development of the vortex street behind a flat plate. By including the heat conduction equation Harlow and Fromm (82) are also able to compute isotherms and the heat transfer from a heated cylinder to the surrounding fluid. Correlation curves of Nusselt versus Reynolds number were obtained from the calculations. A finite difference solution for the motion of a circular disc through an incompressible viscous fluid has been obtained by Michael (83) for low Reynolds numbers. Drag coefficients determined in this computer experiment were in excellent agreement with experimental values. Amsden and Harlow (84) used the Particle in Cell method to study the development of the wake behind the flat base of a projectile moving supersonically through air. Wake experiments, such as those described above coupled with exothermic chemical reactions might shed some interesting light on the behavior of flame holders and might be an interesting subject for future study.
Another significant computer experiment is the calculation of the interaction of a shock or rarefaction with a bubble by Evans et al (70). The Particle in Cell method was used to compute bubble shapes and temperature distributions due to the interactions between shocks and rarefactions and helium bubbles in neon and neon bubbles in helium. The impetus for this study was in part provided by a desire to study the role of inhomogeneities in hot spot formation in explosives. Typical particle configurations obtained when a shock wave in neon hits a helium bubble are shown in Fig. 6, and a corresponding temperature profile is shown in Fig. 7. As mentioned above, the results of this experiment were used by Mader (69) in his investigation of hot spot initiation of explosives. Computer experiments similar to this one might also shed some light on the mechanism of droplet breakup, a problem of central importance in the formation and combustion of liquid sprays as well as in the propagation of detonations through fields of liquid fuel and oxidizer droplets. Amsden and Harlow(85) calculations of the instability and growth of disturbances at the slip interface between two viscous fluids moving with different velocities also bears some relation to the flame holder problem.
External flows with chemical nonequilibrium such as might be encountered in hypersonic inlets or external flow type rocket nozzles have also been the subject of extensive numerical experimentation. Thus for example Spurk et al (86) used the method of characteristics as described by Sedney and Gerber (87) to obtain exact solutions for high temperature air flowing over pointed bodies. A realistic air model was used including the species N2, O2, NO, N, and O; however, vibrational, rotational, and translational equilibrium was assumed. Calculations were made for various free stream densities and it was found that the results could be correlated with ρ∞L, the product of free stream density and characteristic length, as predicted by Gibson's nonequilibrium scaling law (88). The results were also compared to other calculations using idealized air models. South (89) applied the method of integral solutions to the computation of nonequilibrium flow over cones and wedges. The method of characteristics was used by Lee (90) to study the jet structure at the exit of a plug nozzle. Because the calculations were programmed for a computer it was possible to explore the structure of the jet for various secondary flow pressures at the base of the nozzle, and for various ambient pressures. In effect the plume shapes and shock locations could be determined experimentally on the computer. There have been, of course, many other numerical computations of external flows, particularly in relation to the re-entry problem, but these are considered beyond the scope of the present review.
This discussion of fluid dynamics experiments is by no means complete, though the general character of these experiments has been illustrated. It has also been possible to compute the free surface behavior of fluids in tanks subject to acceleration (91). Similarly the computer experiments of Brode and others have made possible the effective design of high velocity launchers for ballistic ranges (92).
The study of boundary layers invariably requires the solution of non linear ordinary or partial differential equations so that almost all boundary layer analyes starting with the pioneering work of Blasius (93) have required numerical or computer experimentation. Generally numerical procedures must be used to solve the non linear ordinary differential equations which result from a similarity or momentum integral theory. In chemical propulsion boundary layers play an important role in determining inlet and nozzle flows and nozzle and combustion chamber heat transfer. The subject of boundary layer theory has been discussed extensively in a number of reviews and books (94), (95) and so will not be discussed in detail here.
It does appear appropriate to mention the work of Smith et al (96) who made a direct numerical attack upon the partial differential boundary layer equations instead of using the more usual and somewhat restrictive, similarity or momentum integral approach. In this numerical approach arbitrary initial boundary layer velocity profiles as well as free stream conditions could be specified. In the compressible case a Howarth-Dorodnitsyn type transformation was applied to the coordinate normal to the wall. Finite difference equations were taken with respect to the x direction thus reducing the partial differential boundary layer equations to a system of ordinary differential equations. Since these equations had two point boundary conditions they had to be solved iteratively. The parabolic nature of the boundary layer equations, however, made it possible to solve these equations one by one rather than simultaneously. The method was tested against several known exact solutions. In the case of heat transfer to a flat nosed cylinder excellent agreement with measurements and with calculations based on the local similarity assumption was found. The exact method is able to treat laminar compressible boundary layers with nonequilibrium chemical reactions and provides a computer experiment which can be used to test various approximate boundary layer theories. The earlier numerical boundary layer computations of Flugge-Lotz et al (97), (98) should also be mentioned, in regard to numerical integration of the full boundary layer equations.
Use of extensive numerical computation or computer experiments has led to considerable progress in the understanding of shock wave phenomena. Most of these computations either deal with unsteady problems such as shock tube flow or use Crocco's (17) suggestion of changing steady flow problems to a study of the asymptotic behavior of an unsteady flow. This procedure was used by Crocco(17) to study the one dimensional formation of a standing shock wave in a diverging duct with initially supersonic flow when the back pressure is suddenty increased. An explicit finite difference scheme was used to solve the Navier-Stokes equations, and the von Neumann procedure was used on a linearized-form of the difference equation in order to establish approximate stability criteria, It was found possible to compute the final steady shock configuration by this method. An attractive feature of Crocco's scheme is that the details of the unsteady behavior appear to be immaterial to the final steady flow. Bohachevsky and Rubin (18) also used the time unsteady technique to compute nonequilibrium flows about blunt bodies with detached shocks. Lax's finite difference scheme, which, as discussed in Section II, implicitly contains an artificial viscosity term, was used, A stability analysis of the linearized finite difference equations was supplemented by stability and convergence tests on the computer. Bohachevsky and Rubin (18) studied shock reflection at the closed end of a shock tube, plane flow over a rectangular step, axisymmetric flow over a sphere and supersonic flow in the afterbody region. The gas chemistry was idealized by assuming a Lighthill dissociating gas. Burstein(100) studied the regular and Mach reflections of oblique shock waves using artificial viscosity and studying asymptotic t-∞ behavior of an initial unsteady flow. The good agreement between the analytical and numerical results for regular reflection provides an excellent check of Burstein's (100) results.
Both real and artificial viscosity have been usld in the study of a number unsteady shock wave problems. Butler (99) used this technique to study the transient loading of blunt objects due to the passage of a shock wave. Of course a finite difference technique in which artificial viscosity is used cannot yield adequate information about shock wave structure. Kurzrock and Mates (101) used real viscosity in their extensive calculation of shock tube flow development. These calculations yielded the development of both the shock structure and expansion wave. The detailed behavior of the shock structure during reflection from the end of a tube has been studied numerically using real viscosity by Scala and Gordon (102) and by Petty (103). An examination of temperature and velocity gradients within the shock wave made it possible to assess the validity of this reflection in detail. These calculations are of course also of considerable value in interpreting various shock tube end wall measurements.
Since about 1955-58 the high speed computer has in many ways revolutionized the study of the complex physical problems related to chemical propulsion as should be evident from the survey above. Essentially typical examples of computer experiments have been presented though it has been impossible to take all recent computer experiments into account.
Although the development of computer experiments during the last decade has been remarkable many important and difficult problems remain to be solved, Thus, for example, a detailed study of the important problem of droplet breakup would appear worthwhile. Computer experiments no doubt will make contributions to the study of turbulence and combustion instability, and important advances in the computer determination of reaction rate constants are to be expected. Many of the present computer experiments are limited in scope by computer speed and storage capacity, so that the numerical analysis of some of these problems must await the further development of digital and hybrid digital-analog computers, It is further to be hoped that further improvements in the numerical analysis of partial differential equations will put the use of finite difference methods on a less empirical level.
1. Fromm, J E, Harlow, F H, Numerical Solution of the Problem of Vortex Street Development, Phys Fl, 6, 975, (1963).
2. Duff, R E, Calculation of Reaction Profiles Behind Steady State Shock Waves. I. Application to Detonation Waves, J. Chem Phys, 28, 1193-1197, (1958).
3. Fay, J A, Riddell, F R, Theory of Stagnation Point Heat Transfer in Dissociated Air, J. Aero Sci, 28, 73, (1958).
4. Eckert E, Drewitz O, Der Warmeubergang an eine mit groser Geschwindigkeit langsangestromte Platte, NACA Tech. Memo, No 1045, (1943).
5. Schlichting, H, Boundary Layer Theory, McGraw-Hill, New York, (1955), p 266.
5. Polhausen, E, Der Warmeaustunsch zwischen festen Korpern und Flussigkeiten mit kleiner Reibung und kleiner Warmeleitung, ZAMM, 1, 115, (1921).
7. Howarth, L, Modern Developments in Fluid Dynamics, High Speed Flow, Vol. I, Oxford at the Clarendon Press (1953) (see chapter on boundary layer theory).
8. Fox, L, Numerical Solution of Ordinary and Partial Differential Equations, Pergamon Press, (1962).
9. Richtmeyer, R D, Difference Methods for Initial Value Problems, Interscience, New York, (1957).
10. Ralston, A, Wilf, H S, Mathematical Methods for Digital Computers, Wiley, New York, (1960).
11. Todd, J, ed, Survey of Numerical Analysis, McGraw-Hill, New York, (1962).
12. Milne, W E, Numerical Solution of Differential Equations, Wiley, New York, (1953).
13. Fromm, J E, A Method for Computing Nonsteady, Incompressible Viscous Fluid Flows, LASL Rept-2910. (Sept. 1963).
14. Varga, R S, Matrix Iterative Analysis, Prentice Hall Inc, Englewood Cliffs, N.J, (1962).
15. Belotserkovskii, O M, Chushking, P L, The Numerical Solution of Problems in Gas Dynamics, in Basic Developments in Fluid Dynamics, Vol 1, M. Holt, editor, Academic Press, New York, (1965).
16. Smith, A M O, Clutter, D W., Machine Calculation of Compressible Laminar Boundary Layers, AIAA J, 3, 639, (1965).
17. Crocco, L, A Suggestion for the Numerical Solution of the Steady Navier-Stokes Equations, AIAA J, 3, 1824, (1965).
18. Bohachevsky, I O, Rubin, E L, NonEquilibrium Flows with Detached Shock Waves,AIAA J, 4, 600, (1966).
19. Evans, M W, Harlow, F H, The Particle-in-Cell Method for Hydrodynamic Calculation, LASL Rpt No. LA 2139, (1957).
20. Harlow, F H, The Particle in Cell Computing Method for Fluid Dynamics, in Methods in Computational Physics, Vol. 3, Fundamental Methods in Hydrodynamics, B. Alder, S. Fernbach, and M. Rotenberg, editors, Academic Press, (1964).
21. Evans, M W, Harlow, F H, Interaction of a Shock or Rarefaction with a Bubble, Phys Fl, 5, 651, (1962).
22. von Neumann, J, Richtmeyer, R D, A Method for Numerical Calculations of Hydrodynamic Shocks, J Appl Phys, 21, 232-237, 353-359, (1950).
23. Lax, P D, Weak Solutions of Non-Linear Hyperbolic Equations and their Numerical Computations, Commun Pure Appl Math, 7, 159-193, (1954).
24. Householder, A L, Principles of Numerical Analysis, McGraw Hill (1953) Chap. 8.
25. Bray, K N C, Atomic Recombination in a Hypersonic Wind Tunnel Nozzle, J Fluid Mech, pt 1, 1.
26, Penner, S J, Chemical Reactions in Flow Systems, Butterworth's, 1955.
27. Cheng, H K, Lee, R S, The Freezing of Dissociation in Supersonic Nozzle Flows, AIAA Preprint No. 66-1 (1966).
28. Hall, J C, Russo, A L, Cornell Aero. Lab Report AD-1118-A-6, (1959).
29. Vincenti, W C, Arnold Engineering and Development Center Report AEDC-TN-61-65, May 1961.
30. Emanuel, G, Problems Underlying the Numerical Integration of the Chemical and Vibrational Rate Equations in a Near Equilibrium Flow, Arnold Engineering and Development Center Report AEDC-TDR- 63-82.
31. Emanuel, G, A General Method for Numerical Integration through a Saddle Point Singularity with Application to One Dimensional Nonequilibrium Nozzle Flow, Arnold Engineering and Development Center Report AEDC-TDR-64-29, 1964.
32. Eschenroeder, A Q, Boyer, D W, Hall, J G, Non Equilibrium Expansions of Air with Coupled Chemical Reactions, Phys. Fl. 5, 615-624, (1962).
33. Sarli, V J, Blackman, A W, Buswell, R J, Kinetics of Hydrogen-Air Flow Systems II. Calculations of Nozzle Flows for Ramjets, Ninth Symposium (International) on Combustion, Academic Press, New York, 231-240, (1963).
34. Moretti, G, A New Technique for the Numerical Analysis of Non Equilibrium Flows, ALAA J, 3, 223 (1965).
35. Vincenti, W G, Calculation of the One Dimensional Nonequilibrium Flow of Air through a Hypersonic Nozzle-Interim Report, Stanford Univ., Dept. of Aero Engineering, Report No. 101, January 1961.
36. Lordi, J A, Mates, R E, Nonequilibrium Effects on High Enthalpy Expansions of Air, AIAA J, 3, 1972, (1965).
37. Mates, R E , Lordi, J A, Technique for Solving Non-Equilibrium Expanding Flow Problems, Aerospace Res. Labs Report 65-2, (1965).
38. Hall, J G, Russo, A L, Studies of Chemical Nonequilibrium in Hypersonic Nozzle Flows, Kinetics Equilibria and Performance of High Temperature Systems, G, Bahn and E, Zukoski, editors, Butterworths, Washington, DC, p. 219, (1960).
39. Hall, J G, Eschenroeder, A Q, Klein, J T, Chemical Nonequilibrium Effects on Hydrogen Rocket Impulse at Low Pressures, ARS J , 30, 188, (1960).
40. Eschenroeder, A Q, Lordi, J A,"Catalysis of Recombination in Nonequilibrium Nozzle Flows, Ninth Symposium (International) on Combustion, Academic Press, New York, 241, (1963).
41. Westenberg, A A, Favin, S, Complex Chemical Kinetics in Supersonic Nozzle Flow, Ninth Symposium (International) on Combustion, Academic Press, New York, 785, (1963)
42. Wilde, K A, Complex Kinetics in Adiabatic Flow: C-H-O(N) Systems, AIAA J, 3, 1847 (1965),
43. Sarli, V J, Blackman, A W, Migdal, D, Chemical and Nozzle Flow Losses in Hypersonic Ramjets, AIAA J, 3, 139, (1965).
44. Hoglund, R F, Recent Advances in GasParticle Nozzle Flows, ARS J, 32, 662, (1962).
45. Kliegel, J R, Nickerson, G R, Flow of Gas-Particle Mixtures in Axially Symmetric Nozzles, ARS Preprint 1713-61, (April 1961).
46. Bailey, W S, Nilson, E N, Serra, R A, Zupnik, T F, Gas Particle Flow in an Axisymmetric Nozzle, ARS J, 31, 793, (1961).
47. Kliegel, J R, Gas-Particle Nozzle Flows, Ninth Symposium (International) on Combustion, Academic Press, New York, 811, (1963).
48. Hoffman, J D. Lorenc, S A, A Parametric Study of Gas-Particle Flows in Conical Nozzles, AIAA J, 3, 103, (1965).
49. Bunker, D L, Computer Experiments in Chemistry, Scientific American, July 1964.
50. Bunker, D L, Monte Carlo Calculation of Triatomic Dissociation Rates. I. N2O and O3, J Chem Phys, 37, 393, (1962).
51. Wall, F J, Hiller, L A, Mazur, J, Statistical Computation of Reaction Probabilities, J, Chem. Phys, 29, 259, (1958).
52. Thiele, E, Wilson, D J, Anharmonicity in Unimolecular Reactions, J. Chem. Phys, 35, 1256, (1961).
53. Benson, S W, Berend, G S, Wu, J C, Classical Model for Vibrational and Rotational Excitation of Diatomic Molecules by Collision. I. Hard Sphere Collision, J Chem Phys, 38, 25, (1963).
54. Kuntz, P J, Energy Distribution Among Products of Exothermic Reactions. II. Repulsive Mixed and Attractive Energy Release, J Chem Phys, 44, 1168, (1966).
55. Rogers, M, Analog Computer Studies of Particle Surface Interactions, Rarefied Gas Dynamics Symposium, Toronto, Vol. 2, (Adv in Appl. Mech. Supplement), Academic Press, New York, 429, (1966).
56. Gentry, R A, Harlow, F H, Martin, R E, Computer Experiments for Molecular Dynamics Problems, Methods in Computational Physics, B. Alder, S. Fernbach, and M. Rotenberg, editors, Academic Press, (1965), pp 211-245.
57. Belles, F E, Lauver, M R, Origin of OH Chemiluminescence During the Induction Period of the Hg-O9 Reaction Behind Shock Waves, J. Chem Phys, 40, 415, (1964).
58. Brokaw, R L, Analytical Solutions to the Ignition Kinetics of the Hydrogen-Oxygen Reaction, NASA TN D-2542, (1964).
59. Momtchiloff, I N, Taback, E D, Buswell, R J, Kinetics in Hydrogen-Air Flow Systems. I. Calculation of Ignition Delays for Hypersonic Ramjets, Ninth Symposium (International) on Combustion, Academic Press, New York, (1963), pp. 220-230.
60. Slutsky, S, Tamagno, J, Trentacoste, N, Supersonic Combustion in Premixed Hydrogen-Air Flows, AIAA J, 3, 1599, (1965). -
61. Hermance, C E, Shinnar, R, Summerfield, M, Ignition of an Evaporating Fuel in a Hot, Stagnant Gas Containing an Oxidizer, AIAA J, 3, 1584, (1965).
62. Hermance, C E, A Model of Composite Propellant Combustion Including Surface Heterogeneity and Heat Generation, AIAA J, 4, 1629, (1966).
63. Gilbert, L M, Scala, S M, Combustion and Sublimation of Cones, Spheres, and Wedges at Hypersonic Speeds, AIAA J, 3, 2124, (1965). -
64. White, D R, Turbulent Structure of Gaseous Detonation, Phys Fl, 4, 465, (1961).
65. Erpenbeck, J J, Stability of Idealized One Reaction Detonations: Zero Activation Energy, Phys. Fl, 8, 1192, (1965).
66. Fickett, W, Wood, W W, Flow Calculations for Pulsating One-Dimensional Detonations, Phys. Fl, 9, 903, (1966).
67. Gilbert, R B , Strehlow, R A, Theory of Detonation Initiation Behind Reflected Shock Waves, AAA J, 4, 1777, (1966).
68. Warner, F J, The Initiation of Detonation in Solid Explosives, Ninth Symposium (International) on Combustion, Academic Press, New York, (1963), pp. 536-543.
69. Mader, C L, Shock and Hot Spot Initiation of Homogeneous Explosives, Phys Fl, 6, 375, (1963).
70. Evans, M W, Harlow, F H, Meixner, B D, Interaction of Shock or Rarefaction with a Bubble, Phys. Fl, 5, 651, (1962).
71. Mader, C L, Initiation of Detonation by the Interaction of Shocks with Density Discontinuities, Phys. Fl, 8, 1811, (1965).
72. Friedman, M H, A Correlation of Impact Sensitivities by Means of the Hot Spot Method, Ninth Symposium (International) on Combustion, Academic Press, New York, (1963), pp 294-300.
73. DeSoto, S, Friedmann, A A, Flame Spreading and Ignition Transients in Solid Grain Propellants, AIAA J, 3, 405, (1965).
74. Emanuel, G, Vincenti, W G, Method for the Calculation of the One-Dimensional Nonequilibrium Flow of a General Gas Mixture through a Hypersonic Nozzle, Arnold Engineering and Development Center, AEDC-TDR-62-131, (1962).
75. Hirschfelder, J O, Curtiss, C F, Theory of Detonations. I. Irreversible Unimolecular Reaction, J. Chem Phys, 28, 1130, (1958).
76. Hirschfelder, J O, Curtiss, C F, Propagation of Flames and Detonations, Advances in Chemical Physics, Vol. III, pp. 59-129, Interscience Publishers, (1961).
77. Oppenheim, A K, Rosciszewski J, Determination of the Detonation Wave Structure, Ninth Symposium (International) on Combustion, Academic Press, New York, 1963, pp 424-434.
78. Campbell, E S, Heinen, F J, Schalit, L M, A Theoretical Study of Some Properties of Laminar Steady State Flames as a Function of Properties of their Chemical Components, Ninth Symposium (International) on Combustion, Academic Press, New York, (1963), pp 72-80.
79. Kuehl, D K, Ignition and Combustion of Aluminum and Beryllium, AIAA J, 3, 2239, (1965).
80. Fromm, J E, Harlow, F H, Numerical Solution of the Problem of Vortex Street Development, Phys Fl, 6, 975, (1963).
81. Harlow, F H, Fromm, J E, Computer Experiments in Fluid Dynamics, Scientific American, Mar 1965,
82. Harlow, F H, Fromm, J E, Dynamics and Heat Transfer in the von Karman Wake of a Rectangular Cylinder, Phys Fl, 7, 1147, (1964).
83. Michael, P, Steady Motion of a Disk in a Viscous Fluid, Preprint of paper submitted to Physics of Fluids, October 1965, Brookhaven National Laboratory, Upton, New York, October 1965.
84. Amsden, A A, Harlow, F H, Numerical Calculation of Supersonic Wake Flow, AIAA J, 3, 2081, (1965).
85. Amsden, A A, Harlow, F H, Slip Instability, Phys Fl, 1, 327, (1964).
86. Spurk, J H, Gerber, N, and Sedney, R. , Characteristic Calculations of Flow Fields with Chemical Reactions, AIAA J, 4, 30, (1966).
87. Sidney, R, Gerber, N, Nonequilibrium Flow over a Cone, AIAA J, 1, 2482, (1963).
88. Gibson, W E, Dissociation Scaling for Nonequilibrium Blunt Nose Flows, ARS J, 32, 285, (1962).
89. South, J C, Application of the Method of Integral Relations to Supersonic Nonequilibrium Flow Past Wedges and Cones, NASA TR R-205, (1964).
90. Lee, C C, Gasdynamic Structure of Jets from Plug Nozzles, AIAA J, 4, 1114, (1966).
91. Griffith, W C , invited talk on Problems in Manned Atmospheric Flight, Amer Phys Soc, Fluid Dynamics Div meeting, Palo Alto,, November 21-23, 1966.
92. Glass, I I , Univ. of Toronto, Institute of Aerophysics, private communcation, November, 1966.
93. Blasius, H, Grenzschichten in Flassigkeiten mit kleiner Reibung, z, Math U, Phys, 56, 1, (1908), also NACA TM 1256.
94. Moore, F K, Laminar Flow Theory, Vol. IV, High Speed Aerodynamics and Jet Propulsion, Princeton University Press, (1964).
95. Bartz, D R, Turbulent Boundary-Layer Heat Transfer from Rapidly Accelerating Flow of Rocket Combustion Gases and of Heated air,, and Chung, P. M., Chemically Reacting Nonequilibrium Boundary Layers, both review articles in Advances in Heat Transfer, Hartnett, J.P. and Irvine, T. F., Jr., editors, Academic Press, New York, (1965).
96. Smith, A M O, Jaffe, N A, General Method for Solving the Laminar Nonequilibrium Boundary-Layer Equations of a Dissociating Gas, AIAA J, 4, 611, (1966).
Smith, A M O, Clutter, D W, Machine Calculation of Compressible Laminar Boundary Layers, AIAA J, 3, 639, (1965). Smith, A M O, Clutter, D W, Solution of the Incompressible Laminar Boundary Layer Equations, AIAA J, 1, 2062, (1963). -97. Flugge-Lotz, I, Baxter, D C , The Solution of Compressible Laminar Boundary Layer Problems by a Finite Difference Method, Part I and II, Tech Rept No 103 and 110, Div of Eng Mech, Stanford Univ, September 1956, and September 1957.
98. Flugge-Lotz, I, Blottner, F,G, Computation of Compressible Laminar Boundary Layer Flow Including Displacement Thickness Interaction using Finite Difference Methods, Tech Rept No 131, Div of Eng Mech, Stanford Univ, January 1962.
99. Butler, T D, Numerical Calculations of the Transient Loading of Blunt Obstacles by Shocks in Air, AIAA J, 4, 460, (1966).
100. Burstein, S Z, Numerical Methods in Multidimensional Shocked Flows, AIAA J, 2, 2111, (1964).
101. Kurzrock, J W. Mates, R E, Exact Numerical Solutions of the Time Dependent Compressible Navier-Stokes Equations, AIAA Preprint 66-30, Third Aerospace Sciences Meeting, New York, January 1966.
102. Scala, S M, Gordon, P, Reflection of a Shock Wave at a Surface, Phys. Fl, 9 1158, (1966).
103. Petty, J, paper presented at APS, Div. of Fluid Dynamics meeting, Cleveland, Ohio, November 1965,