2427
Numerical Simulation of Complicated Flow Fields
t
IO 0 ,
I
I
I
I
I
I
TEST PROBLEM ~
Concerning accuracy, we feel that 4 less than lo4, in the current test problem, gives no useful additional information and could be a waste of valuable computer time in the context of a coupled chemical kinetic-hydrodynamic reactive flow problem. The hydrodynamics are not sufficiently accurate to deserve more than a 1% calculation of the reaction rates. This is reinforced, in any case, by rate coefficient uncertainties. References and Notes
aFigure 1. A graph that gives a comparison of efficiency for various numerical integration schemes used to solve the test problem. The graph gives normalized run time vs. (the error criterion).
CHEMEQ over
classical methods a couple of points from a second-order predictor-corrector scheme without asymptotics has been plotted. CHEMEQ is approximately 25 times faster. There is clearly room for a substantial further improvement in CHEMEQ since the calculations were performed by our library integrator and no hand “tuning” has taken place. For example, improvement could be gained by performing the asymptotics on the equation(s) known to be stiff without testing.
P. R. Byron, Simulation, 11, 219 (1968). T. J. Keneshea, Research Report AF CRL-0221 (April, 1967). C. E. Treanor, Math. Comput., 20, 39 (1966). M. E. Fowler and R. N. Warten, I5M J . , 11, 537 (1967). C. W. Gear, Commun. ACM, 14, 176 (1971); DIFSUB for Solution of Ordinary Differential Equations (D2), bid., 14, 185 (1971). M. D. Kregel, J . Atmos. Terr. Phys., 34, 1601 (1972). J. P. Boris, 6. E. McDonald, T. P. Coffey, and T. R. Young, Proceedings of the DNA High Akitude Nuclear Effects Symposium, Vol. 2, DASIAC SR-130 (Dec, 1971). J. P. Borls, Proceedings of the 2nd European Conference on Computational Physics (also NRL Memorandum Report 3327 and references therein). A. W. Ali, T. Coffey, P. Kepple, S.Piacsek, A. Warn-Varnas, and T. R. Young, Trans., Am. Geophys. Union, 53, 458 (1972). R. Bulirsch and J. Stoer, Numer. Math., 8, 1 (1966). J. P. Boris and N. K. Winsor, MATT-652, Princeton Plasma Physics Laboratory Report (Nov, 1970). M. Kregel, private communication. D. Edelson, J. Chem. Ed., 52, 642 (1975). Proceedings of the DNA 1973 Atmospheric Effects Symposium (U), Vol. 4, DNA 3131P-9 (VI 75-442), pp 571-591.
Parallel Computation of Unsteady, Three-Dimensional, Chemically Reacting, Nonequilibrium Flow Using a Time-Split Finite-Volume Method on the Illiac I V Walter A. Reinhardt Ames Research Center, National Aeronautics and Space Administration, Moffett Field, California 94035 (Received May 12, 1977) Publication costs asslsted by the Ames Research Center, National Aeronautics and Space Administration
The system of unsteady, three-dimensional, partial differential equations used to simulate the inviscid flow of air in chemical nonequilibrium is approximated by a set of factored, finite-volume difference operators where the effect of chemical production is also contained in the factorization. The method is similar to that of Rizzi and Bailey, except for the emphasis on vector-matrix re-formulation designed to be suitable for the special architecture of modern advanced computers (e.g., the Illiac IV, CDC 7600, or Star). The systematic application of the operators yields a second-order accurate numerical algorithm. The method is programmed in the vector Fortran-like language called CFD developed in the CFD Branch at the NASA Ames Research Center. Translators are available for this language which allows debugging on serial computers, but all results for the examples given were obtained from the Illiac. The problem described is a numerical simulation of the flow in the high temperature stagnation region of a reentering space shuttle orbiter flying at large angles of attack (-40O). Capability for treating arbitrary geometry in a flow containing subsonic, sonic, and supersonic regions is demonstrated by this method. The air chemistry is described by a five-reaction model which includes the three dissociation reactions for Nz,02, and NO and the two rearrangement reactions involving NO. The vector-matrix formulation and the unique disk-memory mapping results in extremely efficient data management for the architecture of the Illiac and makes maximum use of the Illiac’s “data crunching” capability. Comparative running times are given for the Illiac IV and the CDC 7600.
Introduction The new generation of very fast, special purpose vector computers (e.g., Illiac Iv, CDC Star, CDC 7600, TI ASC, Gray 1, and IBM 370/195) has made possible the numerical simulation of complicated flow fields, including chemical reactions, about geometrically complex b0dies.l The need for these solutions results partly from the
continuing interest and usefulness of more sophisticated atmospheric entry vehicles such as the space shuttle. TO obtain such results, the split finite-volume method discussed in this paper is a viable numerical method. The equations that are approximated using this scheme are quite general and, with the exception of the easily modifiable chemical reaction model, are applicable for studies The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
Waiter A. Relnhardi
Figure 1.
Shock structure about
the
shuttle orbiter
of combustion, pollution, and other chemically reacting flow phenomena, where convective transport effects dominate the influence of radiative, viscous, and other transport mechanisms. The resulting numerical simulations are particularly valuable to the vehicle d e s i g n e F as a source of information for estimating heat-transfer rates, boundary-layer effects6 (e.g., the influence of flow separation and entropy layer “swallowing”), surface-material corrosion, as well as the aerodynamic loads acting on the spacecraft during atmospheric entry. Wind-tunnel tests alone cannot provide such information.%‘ The effect of chemical reactions greatly complicates the scaling of such data to what happens in full scale actual flight. The shock perturbed flow about a shuttle orbiter flying at a large angle-of-attack during atmospheric entry is illustrated in Figure l. The three-dimensional flow is interlaced with embedded discontinuities that inclose nonreacting or reacting regions, depending on the altitude and velocity along the flight trajectory. The flow field itself contains a large variety of possible flow phenomena. To numerically simulate these flows requires several varied methods. Within the nose region (circled in Figure 1)there exist subsonic flow (in the stagnation region), transonic flow, and supersonic flow. Here the numerical simulation^'-^' generally involve marching the unsteady fluid flow equations in time, starting with an initially specified estimated flow field. The marching continues until unsteady effects are no longer observed. The flow field on the exit boundary of this solution (i.e., on the “data surface”) serves as the initial condition for numerical methods that approximates the steady representations of the flow equations. As long as the exit boundary lies in supersonic flow, the problem is hyperbolic in the direction of the flow. The coordinate in this direction is timelike; the data surface can then be marched stepwise down the body either as a generalized co~rdinate’~-~‘ surface, as a lane,'^.'^ or by method-ofcharacteristics.” Still other methods solve this flow using the unsteady flow equations similarly as in the nose reg i ~ n . ’ ~ ,Canopy ’~ shocks, induced by body curvature, and cross-flow shocks that result from strong cross flow at large angles of attack are also found in the numerical simulation~.’~,~,*’ The intersection of the bow and wing shocks, besides introducing subsonic flow pockets at the wing leading edge, yields a variety of complicated flow discontinuity effects such as multiple shocks and slip surfaces. These have also been investigated in numerical simulation~.’~.~~ The subject of this paper is the finite-volume method first proposed by MacCormack, Warming, and Pa~llay2”~ and generalized by Rizzi, Inouye, and Bailey>* Schiff?‘ The Jovmsl of F h y s b l Ctmnfstry. Val. 81, No. 25, 7977
D e i ~ e r t , 2and ~ Hung and MacCormack.26 The method is quite flexible and has been employed by Rizzi, Klavins, MacCormack, and Bailey’zJ3 in the calculation of the supersonic flow as well as the subsonic flow regions about the space shuttle. The method has also been used in inviscid studies of jet counterflowz‘ as well as the twodimensional viscous studies of separated transonic flow over an airfoilz5and of separated supersonic flow over a compression The method is based on the integral conservation-law representations of the fluid flow equations. “Finitevolume” denotes the partitioning of the entire flow region into arrays of topological hexahedra that are the computational elements. The calculation procedure of Rizzi and Bailef involves calculating the fluxes through the hexahedra faces and the chemical production within these elements. In applying the method to special purpose advanced computers, a valuable adjunct is “time-splitting”, that is, factoring the three-dimensional spatial differencing operator into three one-dimensional operators. This method remains second-order accurate, but improves operational efficiency on conventional computers,7a8and has especially profound effects on Illiac IV efficiencies. This occurs because either data arrays or their transposes are equally accessible within the IUiac’s main memory disk storage and thus data transfer is equally optimal, regardless of data order requirements of the operator being executed. Another “splitting” discussed here is that of separating the species convection from chemical production. In this case the production terms are contained within a separate operator that is also one-dimensional. Several advantages occur. The chemical effects can be “advanced with a smaller time step than that used for convection (several applications of the chemical production operator are still required, however, so that the aggregate step is that of the convection); and, depending on whether the chemistry is ’ or explicit numerical algorithms can “&iff‘ ,27.28 implicit be used without penalty to the accuracy of the overall method. The computation of chemically-reacting, three-dimensional flows, even with the simplest chemical models, seriously strains the capability of other than the new generation of vector computers. These computers achieve their rapidity principally through special hardware features (overlap, parallelism, or pipeline), but to take greatest advantage of their computational capability requires careful use of vector-matrix formalism and programming in a vector language. An easily learned language for programmers who have used Fortran is CFD.29 This language, developed within the Computational Fluid Dynamics Branch at the NASA Ames Research Center,30 leads through the use of translators and compilers to very efficient code for both the Illiac IV parallel computer and the CDC 7600 pipeline computer. Mathematical Formulation In this paper the basic equations will be introduced first. Then the procedure for approximating these equations using the finite-volume method will be described. The generality of the equations and discussion will be relaxed when the coordinate system, which has proved valuable for solving the flow in the nose region of the shuttle at large angles of attack, is introduced. The discussion becomes more specific when the Illiac IV architecture is described and we point out the procedure for solving the flow on the IUiac. Finally, several results me presented to demonstrate the viability of the method as well as of parallel processing. Conservation Equations. The unsteady equations of fluid dynamics, which govern the flow of a multicomponent
Numerical Simulation of Complicated Flow Fields
2429
TABLE I: Chemical Reaction Model
AS) + I
0,+M==20+M N, + M * 2N + M NO + M = N-i- 0 -i- M NO+O=N+O, NO + N = 0 + N,
(1) (2)
(3) (4 (5)
1,
I + 1, k I
reacting mixture of gases, are described in vector integral conservation-law form by the representation I,
where the column vectors U and Q and the second-order tensor H , whose elements are flux yectorg, are defined in eq A for flow velocity tj = ui, + ui, + wi,, total specific
1
y
I, k +
+ I, 1, k + I
I
v
Figure 2. A typical computational cell i , j , k, with volume
-I U
and the atomic species (oxygen (0)and nitrogen (N)). The production terms w, are based on a relatively simple chemical reaction model given in Table I. It includes those reactions that most significantly affect the enthalpy. The reaction rates used in this study, which appear explicitly in the chemical production terms w1, eq 1,are the same as used by Davy and Reinhardt (ref 16). Using the reaction model given in Table I, the specific heat ratio may now be written explicitly:
u=p
1
7 = [4(cO,
[3(co2 1 = 1, 2,
* e ' )
s
(A)
internal energy eT = e + q 2 / 2 ,and total enthalpy h T = eT p / p , pressure p, density p, concentration (mass fraction) c,, and chemical production w1. The explicit formulation for the w,will not be given here. The relation used for this study, which depends on a chemical reaction model to be defined, is a conventional expression which can be found in a number of references (e.g., see ref 31 or 32). The equations given above refer, respectively, to conservation of mass, of the three components of momentum, of energy, and of species within an unsteady volume-region enclosed by surfaces which move with a velocity A. The above equations are made complete with the addition of a state relation for pressure p ( e , p , c l ) . The Lighthill model is introduced where the translational and molecular rotational modes are assumed to be fully equilibrated while the molecular vibrational mode is assumed half-excited.16 The following equation for pressure results
+
+ C N , + CNO) + 5 / 2 ( c 0 + c N ) l / + CN, + CNO)+ 3 / 2 ( C ~+ C N ) ~
(2~)
Split Finite-Difference Operators. Finite-difference approximations to the gas dynamic conservation-law equations described in the previous section are used to advance the flow in time from specified initial data. The finite-difference operators to be defined here approximate eq 1 for the labeled computational cell illustrated in Figure 2. If the solution is known at time t (= Cl,lnAt,) inside the topological hexahedron i, j , k with volume ArLJ,kand a s k , and bounded by the six sides As,, As,,,, As,, Ask+l,then it can be determined at time t + At from the time-split sequence denoted by unt1/3 =L,U" (3a) I
un+2/3
un+l
=
= L,
un+1/3
(3b)
~ ~ U n + Z f 3
(3c) The symbol U is used to denote that we are temporarily assuming for this illustration that wl = 0 (i.e,, chemical effects are frozen). The fractional powers imply that three fractional steps are required to advance one time step. The operator representation Un+lj3= L,U" denotes At Ujn+l/3 = U." - -( Hjn*Asj+1 H ~ - , ~ . A < ) (4a) Arj Y
4
Ujn+lI3= 0.5[UJn
+
+ 61n+1i3 -
S
P = (7 - l ) p ( e - 1=1c h O )
(2a)
where e = eT - (u2+ u2 + w 2 ) / 2is the internal energy per unit mass and hp is the heat-of-formation corresponding to the species c,. The ratio for specific heats is y = c p / c u , and the sound velocity, needed in subsequent expressions, is given by a = (YP/P)1'2
(2b)
The model air mixture is assumed to contain the molecular species (oxygen (02), nitrogen (NJ, and nitric oxide (NO))
0
Mean values of the flow variables in the cell are used in the above representation as defined by (5) where Ar,j,k is the small but finite vol_ume of the cell _at that time step. Also, the flux vector Hjn+lI3denotes H( Ujn+'i3),In eq 4 the subscripts i and k,which do not vary, are implied but are not written, to simplify notation. The The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
Walter A. Reinhardt
2430
operational relation for Lk appears identical with that for Lj except for the replacement of k with j and appropriate modification of the fractional power denoting suhstep, and similarly for Lj. This notation exemplifies the one-dimensional character of the operators, and it is this essence that characterizes “splitting”. The condition on At necessary for the stahility of the above method is that the numerical domain of dependence must include the physical 0 n e . 8 ~ Stability ~ ~ ~ conditions ~ ~ ~ ~ can be determined analytically for each operator. For Lj we have
Similar relations are used to obtain Atk and Ati. The operator sequence denoted by eq 3 is then stable if
Flgure 3. Mesh
geometry determined by Series of nested wnes.
At < min ( A t , , Ati, A t k )
(7) This discussion of time splitting has been brief; additional detail can be obtained by referring to the original sources (see ref 8, 22, or 23). Chemical production effects were not considered in the above development as a result the symbol 0 does not appear in these relations. The species concentrations may, however, still he contained in the vector U (see eq 1). Therefore, given an initial nonuniform distribution of the species cl, the above relations allow an accurate simulation of the convection of cI through the flow field. To account for chemical production another operator, denoted Lp! is introduced. Consistent with the notation above we define Un+’ = L ,y+‘ to denote the set of operations Figure 4.
where, similarily as before 0 ss Q(@+’). Accuracy requirements may necessitate that the time step At, for chemistry be different from the At found by inequality 7. For the simulation of shuttle flows Atr = 0.5At gave satisfactory results, but smaller values can he used. The Lp is successively applied (Le., Lehem= Lp,;N = At/Atd until the aggregate chemistry step matches that used for the convection, that is, At = xn=lNAtf. One advantage of splitting the chemistry is that implicit operators may he exchanged with the explicit method given above. This may he necessary if the system of equations becomes “stiff’ (ref 2’7,281. Then methods similar to those described by Lomax and Baileyz7 may be applied. The sequence of operations U”+‘= Lc,&jLkLjun represents a complete time step that properly accounts for chemical production as well as convection. This sequence of operations, however, is only accurate to first order in At. Second-order accuracy is achieved by reversing the operator order during the next time step. The proper second-order sequence, therefore, is given by the two time-step sequence u n + l = L.L J
L.
nL
)
, L ~ L ~ L ~ U ~
(9)
Computational Cell Network. To apply the finitedifference operators requires that the entire flow region be divided into a network of small topological hexahedra. For the nose region flow field of the shuttle orbiter discussed here, the coordinate surfaces are cones, shells, and planes. The cones are arbitrarily positioned and translated The Journal of phvsical Chemlshy. Vol. 81, No. 25. 1977
Partitioning of the shock layer into finite volums.
in the manner displayed in Figure 3. Translation of the cone is accounted for by the coordinate x of its apex measured along the body axis, rotation by the angle $ between its axis and the free stream and, lastly, dilation by its vertex angle w. Each of these conical surfaces is then divided by rays from the apex into equally spaced angular increments. The planes formed by the ray of one cone and the corresponding ray of the next (see Figure 4) delineate a system of contiguous pyramidal columns. All that is needed to specify the ray i, k are its two angles 0 , and &k, made with the z and x axes, and its intersection xi with the body axis. The columns are partitioned into small hexahedra by a sequence of shells that coincide with the body and shock and divide the distance E along each ray into J equal segments. The cells compose a nonorthogonal mesh network floating in time that fills the time-dependent shock layer. The other boundaries are the pitch plane of symmetry and a downstream boundary immersed in supersonic flow. The mesh network is quite general and allows a wide range of flow regions and computational spaces to be studied. Initial and Boundary Conditions. Because the governing equations are hyperbolic and the subsonic region is bounded hy supersonic flow, the time-dependent method is well posed as an initial-boundary-value problem. To commence the calculation, an initial approximation must be specified for the complete field. Our initial flow field is built up as follows.‘,8 A shock surface, axisymmetric about the wind direction and positioned a t an estimated standoff distance, is generated by a quadratic function of the latitudinal angle 0. The slope of this surface can thus he determined at any given point, and the flow properties there are then calculated from the free-stream conditions by use of the Rankine-
Numerical Simulation
of Complicated Flow
2431
Fields
Hugoniot shock relations. On the body, pressure is derived from a Newtonian formula, and the entropy there is set to the same value as that of the streamline which has passed through a normal shock (the stagnation streamline). With these two properties, the density and velocity components can be found by using the equation-of-state and integrated steady-energy relation (the species are assumed constant at their free-stream values throughout the field). Finally, the flow properties within the shock layer are specified by linearly interpolating between the shock and body values along each pyramidal column of cells. This procedure yields a satisfactory set of initial conditions for perfect gas flow over a broad range of Mach numbers and for angles of attach approaching 45”. For nonequilibrium flow, however, the transients generated from the impulse start of the estimated flow field with frozen species composition can cause difficulties. A more satisfactory initial condition in this case is a perfect gas solution after most unsteady effects have decayed. For the inviscid calculations presented here, three distinct types of boundaries are encountered at the edges of the overall mesh: entrance, exit, and streamline boundaries. Along the entrance boundary the dependent variables U (in eq l),are held constant at their supersonic free-stream values, while at the exit they are calculated using one-sided differences. Across cell faces coincident with a streamline boundary, such as an impervious body, no transport is allowed. The only variable actually needed at such a cell face is the pressure, which can be expressed in terms of the interior mesh values of pressure and the derivative of pressure normal to the face. This derivative, ap/a&,,+,, is obtained from the momentum equation normal to the streamline ( u ~ F +~ u~ ~ F , ,+ w ~ F , ,+ 2uuFx, + 2uwFx, + 2uwFy,)/(Fx2+ Fyz + l$z)1’2
aplaQibody = P
(B)
where the body is the surface F ( x , y, z ) = 0 and the subscripts indicate partial differentiation with respect to that variable. The bow shock-wave itself is treated as an interior feature of the flow fields and is not assigned any special attention within the difference operators Li, Lj, and Lk. After every iteration the mesh is readjusted to maintain alignment with the shock. The conservation form of the difference operators will then implicitly satisfy the Rankine-Hugoniot shock-wave “jump” relations and, in addition, accurately determine the solution in the vicinity of the shock. To maintain alignment, the mesh surface coincident with the shock must move with the unsteady shock itself. This is accomplished within an operator called LBSHK.The velocity of each cell segment of this mesh surface is obtained from the simultaneous solution of the shock jump relations for a moving discontinuity and a local characteristic relation, which is valid in the plane defined by the free-stream velocity and the shock normal direction (see ref 8). An iteration procedure yields the shock velocity X for each ray i , k shown in Figure 4. The shock-mesh surface is then moved by the increment AAt computed for each ray. Coordinates are assigned new values to maintain equal spacing of the shells between the shock and body surface, and values of the variables, U, are then found by interpolation.
Illiac IV The current generation of advanced computers, which includes the CDC 7600, CDC Star, TI ASC, Cray 1, and
CONTROL UNIT
PROCESS ELEMEN
BAND
1
2
...
51
52
Figure 5. Illiac I V data routing network.
IBM 370/195 in addition to the Illiac IV, are unconventional digital computers. These computers rely on some form of replication rather than just on circuit speed to obtain their rapidity. Their designs evolved from the speed-of-light limitation of electronic circuits; light travels s). Present day computers about 1ft in a nanosecond operate close to this range. Conventional computers contain the functional components: input/output, control unit, arithmetic and logical unit, and memory. The control unit processes instructions fetched from memory and executes these instructions. The instructions may require memory data transfer to or from peripheral devices (input/output) or arithmetic or logical operations on data within memory (arithmetic and logical unit). The unconventional advanced computers replicate these functional components in varied ways depending on the computer. In the case of the Illiac the principal replication is the 64-fold increase in the number of arithmetic and logic units (called processing elements); this is generally referred to as “parallelism”. There is also an additional replication on the Illiac that is the simultaneous processing of several commands within the control unit (overlap). In the sections that follow some detail on the Illiac operation is described. Architecture of the Illiac IV. A brief understanding of the Illiac’s architecture is essential to ensure proper problem formulation for efficient utilization of the machine. A user’s point of view is that the Illiac IV computer system (see Figure 5) consists of a contrQl unit (CU), 64 processing elements (PE’s), and a main memory (MM) rotating disk system. The CU fetches instructions stored in the PE memories and decodes and broadcasts the decoded instructions to the PE’s for execution. The CU is itself a small unsophisticated “scalar” computer with a 64 64-bit word memory which, besides decoding instructions and serving as a storage device for scalar constants, can perform simple integer arithmetic (not multiplication or division) for such uses as program loop control. The processing elements are where the “data crunching’’ of the computer system takes place. Each PE simultaneously receives and executes, in a manner depending on “mode” settings, the same arithmetic or logical instruction The Journal of Physical Chemistry, Vo/. 81, No. 25. 1977
2432
broadcast from the CU. Data processed by a PE is either fetched from its own processing element memory (each PE contains a random access memory, PEM, with 2048 64-bit words) or “routed” to it from neighboring PE’s. P E operation is bimodal; it depends on mode settings which are either enabled or disabled. If enabled, a PE can store or modify data contained within its own memory; if disabled, its memory is protected and cannot be changed. Execution speed is roughly equivalent to a CDC 6600 computer, but rapid data processing capability occurs because of the 64-fold P E redundancy. The P E random access memory (PEM) is very limited for most problems applicable for the Illiac, thus careful data management to and from the MM is essential to obtain efficiency. Logically, the main memory may be thought of as a drum although, in fact, it is a set of 13 synchronized rotating disks (with angular velocity of 1500 rpm). Each disk contains four bands and each band contains 300 pages of 16 data rows per page. A row contains one word for each PE. The page is the smallest unit of data transferred to and from the disks, and data transfer occurs asynchronously, that is, computation continues during disk access. The data transfer rate is 0.5 X lo9 bits/s. The MM has a total storage capacity of 16 X lo6 64-bit words; its most significant feature is that it can be “formatted”. Formatting is a user controlled mechanism that allows skewed disk storage so that stored matrices or their transposes are addressable with similar efficiencies. This feature profoundly impacts efficient coding of multidimensional problems on the Illiac IV. Although bulk data transfer from PEM to disk is extremely rapid and can occur asynchronously, the mean disk access time is relatively slow compared to the PE compute speed. About 30000 arithmetic row operations (or about 2 X lo6 actual operations) can occur during the time required to position arbitrarily stored disk data under a reading head. Therefore, efficient Illiac programs maximize the number of arithmetic operations between minimal data transfers to or from the MM disk storage. Operations on the llliac IV. Discussed below are considerations unique to vector machines for selecting and programming a method as well as accessing the machine and running the code. The arithmetic units and replication features of advanced computers (Illiac, CDC 7600, CDC Star, Cray 1,TI ASC, and 370/195) have such high cycle frequencies that most often machine speed is controlled not by cycle time for an operation, but by time for data transfer to and from a massive bulk storage area. Thus a numerical method designed for minimizing arithmetic operations without considering data transfer may yield very inefficient vector computer programs. As described in the previous section, the MM and PEM have features that lead to very efficient array operations on the Illiac. The entire data base required for a problem is stored in the MM; selected portions of these data are then transferred to or from PEM, which can be considered as the “working” storage area where data are actually modified. (Problems not requiring the large data base of the blunt-body program discussed here may be designed to operate within PEM and hence use the MM only for data output, e.g., see ref 16.) We denote the massive data stored in MM by the symbol M(1, J). This symbol refers to the two-page blocks of data stored in the MM that are conceptually labeled I and J. As pointed out earlier, a page is the smallest unit of data transferred and assigns 16 words to each PE. Data assignment within these MM blocks is such that for each mesh point, the conservative dependent variables U (eq 1)are sequentially stored first The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
,
Walter A. Reinhardt
SHOCK RING
a J=
J,
J,-1
/
BODY RING
i
Figure 6. (a) The cone coordinate surface B(*, L, J). (b) The shell coordinate surface B(”, L, I).
(10 variables). In addition, also stored are the coordinates (3), surface areas (9), volume (l),and the shock velocity along a ray (1)for a total of 24 variables (eight additional locations are reserved for species variables in studies involving more complicated chemistry models), The mesh points in the meridional direction are stored along the PE’s, and the subscripts I and J of the array M(1, J) in MM refer, respectively, to cones and to shells. We denote by B the actual array of data transferred as it appears stored in the PEM. Ideally, B should be large because of the rapidity of data transfer on the Illiac (0.5 X 10’ bits/s). For example, we may have B(*, L, I) = M(1, 3) or just as easily, we might get B(*, L, J) = Mt(J, 4); that is (see Figure 6), B contains all of the variables for the computations involving the third shell or the fourth cone (Mt designates the matrix transpose of M(1, J)). The asterisk denotes the vector row alignment along the PE’s; L denotes variable type (e.g., p, u,u, w, etc.). The array M(1, J) need not be square. The advantage of the data transfer flexibility on the Illiac should now be apparent; PEM storage is really a “buffer” area and data stored depend only on operator requirements. The operator L,, eq 3a, involves data only on a cone; Lk, eq 3b, involves data on a ring which may be on a cone or shell (see Figure 6); and L, considers data on a shell. The chemistry operator Lchem, eq 8, contains no implicit geometrical preference. The operator Lremesh, designating the reinterpolation of data after the advancement of the unsteady shock-wave surface, has a cone preference. Finding the shock-wave velocity (LBSHK)involves only the single shell that is the shock surface. The actual sequence of operations implemented during each step is illustrated on the flow-chart in Figure 7. The looping about the operators shown in the block diagram denotes that the entire sequence of cones or shells is processed by the loop. Each two-step sequence requires three complete passes through the entire data base stored
2433
Numerical Simulation of Complicated Flow Fields un+2 =
1
---it LREMESH LjLk
L(LchernLchemLi L k L i "n
L
I
I SET INITIAL CONDITIONS
a
ENTHALPY ERROR
u
I-----1
s
I
> 10%
3
z
0
400
200
579
800
STEP NUMBER MAIN OPERATION LOOP
LBSHK
Figure 8. Enthalpy error measure. %.
6 o o L ~ 581.4
DELTA
\
$
Nok-----l
I
LEESIDE ON OUTER CONE
400-
-I
$E
1. WIND-SIDE ON OUTER CONE
Flgure 7. Program flow chart.
in the MM. The operator Lptappearing in the last cone processing loop in Figure 7 denotes the sequence of operations (inequalities 6 and 7) required to find the time increment At which is needed in the difference equations (eq 4 and 8) and in LBsHK(see discussion on Initial and Boundary Conditions).
Results Results from two entirely different calculations are discussed here to demonstrate the viability of the method as well as of parallel processing. The first is a perfect gas calculation (i.e., with frozen chemistry and specific heat ratio y = 1.4) for a flight Mach number of 22.0 and angle-of-attack of 40.2'. The second is a chemical nonequilibrium calculation for a trajectory point corresponding to a Mach number of 21.7, an altitude of about 65.1 km, and an angle of attack of 30' (free-stream conditions are p,(pressure) = 106.2 dyn/cm2,p,(density) = 1.55 X g/cm3, V,(velocity) = 6.544 km/s). Several parameters from the perfect gas calculation will be used to illustrate the convergence to steady state. Several contours of selected variables will be presented from the nonequilibrium solution to illustrate features about the flow. As discussed in the previous section, the calculation starts with a specified set of initial conditions. The entire flow field is then marched in time until the solution becomes steady. Steady state is determined by monitoring the fractional total enthalpy given by cHT = IH, - ( h + q2/2)1/Hmwhere H , = y m p m / ( y-ml)p,. This difference is a measure of deviation from steady flow and is computed for each point in the entire flow field. Displayed in Figure 8 are two curves; the upper curve denotes the total number of mesh points with €HTgreater than 1%, while the lower curve similarly denotes the number of points greater than 10%. The perfect gas calculation used a coarse grid network of 2295 points: 9 shells, 15 cones, and 17 meridional planes (Le., 17 enabled PE's). These results are preliminary and were for comparison with results obtained on the CDC 7600. The grid network, however, is easily refined. We see in Figure 8 that the number of points that
c
y Y
8
/-
L L. .--
97.0
7 64.05
1
,
v)
0
'-
/ 200
-
THE"STAND-OFF" ON INNERMOST CONE
~-
400
- L - 7
600
8.72 800
STEP NUMBER
Figure 9. Shock-wave location.
satisfy the 1 and 10% error criteria increases to a maximum and then decreases. After 579 steps all points have an error less than 10%; after 800 steps, 289 points still have an enthalpy error between 1 and 10%. The shock distance measured from the body on the lee-side ray and on the wind-side ray is shown in Figure 9. Also shown is the standoff distance on the innermost cone whose axis points in the wind direction (see Figure 3). On this cone there is negligible variation of the shock distance around the cone. This distance, therefore, approximates what is normally referred to as the stagnation stream line "shock stand-off distance". We observe that the shock-wave locations in the windward region, where pressures are highest, decay most rapidly to a constant value (i.e., within about 100 steps). In contrast the lee-side shock-wave location on the outermost cone, where the flow has expanded most greatly with considerably lower pressures, shows the slowest convergence rate. Here the shock-wave position oscillates with decaying amplitude to the final constant value. Even though the shockwave position is constant everywhere after about 600 steps, the flow within the shock may still have pockets with errors of between 1 and 10% (Figure 8). Displayed in Figure 10 are the shock-wave position relative to the body surface and contours of molecular oxygen, nitrogen, and nitric oxide and of temperature. Within each frame the inner closed curve represents the body and the outer curves are the shock. The left-half set of curves show the body and shock positions properly scaled relative to each other. In the right set of curves the shock perturbed region is expanded five times so that contour features can be seen. The plots are on the last coordinate surface after the cones are entirely opened to become an axis normal plane located 1.6 m from the shuttle orbiter nose. The grid network is also more refined, The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
Walter A. Reinhardi
2434
translator was intended to simulate CFD nliac code so that program debugging could take advantage of the Fortran aids that have been developed over the years. However, the Fortran generated by CFD translators was found to be very efficient for CDC 7600 operation, because of increased "instack usage of the registers. For example, a given nonequilibrium test generally required about half as much computer time as the "hand-coded Fortran (this optimization, however, can also be achieved with careful programming). The most interesting comparisons, however, are for equivalent runs of the CFD program on the 7600 and on the Mac. Results are listed as follows: MOLECULAR OXYGEN
Perfect gas Nonequilibriurn
N I T R I C OXIDE
MOLECULAR NITROGEN
Flgure IO. Contours on axis normal plane 1.6 m from nose = 65.1 km, Mach number = 21.7, a = 30").
(altitude
which explains the smooth contours. The molecular oxygen contours are easiest to explain; their positions are principally caused by the dissociation reaction, eq 1in Table I. Approximately 100% diasociation of O2occurs and on the windward symmetry plane this dissociation occurs near the shock wave. In the lee side, however, where the flow is cooler and the shock-wave strength is considerably less, the dissociation is not as abrupt and occurs somewhat more uniformly throughout the field. This explanation is simplistic because convection effects also play an important role in the leeside oxygen concentration field. The temperature relaxation observed on the windward symmetry plane is due principally to oxygen dissociation. (The contour closest to the windward shock has a temperature of about 8ooO K each adjacent contour represents an 8% change of this temperature.) The field away from this region is complex, depending on convection as well as on the chemical effects. The nitric oxide contours show a complex interplay of production, destruction, and convection of this species as evidenced by the closed contour lines. The maximum concentration of NO within the field, however, does not exceed around 5% of the mixture. The molecular nitrogen contours (each level corresponds to about a 3% change of the free-stream concentration) show very little dissociation but do show a coupling effect with the nitric oxide as expected (see chemical reaction model in Table I). The results described above can be obtained from either a Fortran language or from a CFD language program. The Fortran version is that used by Rizzi and Bailey (ref 8), and can be run on any computer system that supports Fortran. The CFD code is the principal language for running on the Illiac (less efficient programming languages are available, see ref 29). Also, a translator that generates Fortran from the CFD is available. Originally, the The Jwm-1 of myskal chemishy. vd. 81, No. 25, 1977
7600
Illiac
4.5 slstep
4.0 slstep 5.6 slstep
9.5 slstep
These results are for the very coarse mesh (9 shells, 15 cones, and 17 meridional planes) used to compare 7600 and Illiac Iv results. The Illiac was largely idle (Le., only 17/64, or a little more than one-quarter, of the available PESwere used). Therefore, if comparisons were made on a problem involving optimal mesh network for both machines (i.e., with at least 64 meridional planes) the Illiac would be about four times faster. It should be apparent on comparing the Illiac perfect gas and nonequilibrium time that, probably, some input/output overhead is embedded in the perfect gas time (the ratio is different between both machines). A factor of about 7 in speed-up is obtained for the nonequilibrium calculations. Conclusions The present numerical method permits an efficient and accurate calculation of three-dimensional reacting flow. The method itself was developed originally' using Fortran language on a serial computer (Le., IBM 360/67) and allowed relatively efficient studies of perfect gas blunt-body flows. The current design, using CFD language,29yields the most efficient computer codes for either the CDC 7600 or the Illiac IV. The Illiac results require the least computational time (i.e., about 115 the CDC 7600 time). The method uses a time splitting of the convection differencing operator to achieve efficient data management between random access and disk access storage on the Illiac. The efficient calculation of the effects of the chemical reactions is achieved by an additional splitting of chemical production from convection. The demonstration reported in this paper continues the successful series of applications (ref 8, 13,22,24, 25, and 26) of the finite-volume method for solving complicated multidimensional fluid flow problems. References a n d Notes (1) F. R. W a y . "Cunputaaaal Aerojyna-lLLIAC Nand Beyond". Presented at the Meeting of the IEEE Computer Society. Compwn
Spring. Feb 28-Mar 3. 1977. San Francisco, Calif. A. Ladi, R. J. Vdal. and C. B. Jahnson, NASA Tech. Note, D7189
(2) J.
(Ma,1973). (3)W.D.Goodlch,C.P.U.C.K.Houston.R.M.Meyerr,andL.Omedo.
NASA Spec. Publ.. SP-347 (Mar 1975). (4) W. D. Gocdrkh. C. P. Li, C. K. Houston. P. Chiu, and L. Olmedo. "Numerical cWnputat!mm of orbiier Flow Fields and Heating Rates". AIAA Paper No. 76359, July 1976. (5) J. V. Rakich and M. J. Lanfranco. '"NumericalComputation of Space Shuttle Heating and Surface Streamlines". AIAA Paper No. 76464. July 1976. (6) J. C. Adams. Jr.. W. R. Martindale. A. W. Mayne, Jr., and E. 0. Marchand, '"RealGas Eflects on Hype~nicLaminar Boundarylayer Parameters Including Effects of Entropy-Layer Swallowing". AIAA Paper No. 76358. July 1976. (7) A. W . Rizzl and M. Inouye. AIAA J., 11. 1478-1485 (1973). (8) A. W. Rizzi and H. E. Bailey. NASA Spec. M I . , SP-347 1327-1349 (1975). (9) C. P. Li. J. Spacecr. Rockets. 9, 571-572 (1972).
Checking of Chemical Equations for Mass Balance G. Morettl and G. Bleich, AIAA J . , 5 , 1557-1562 (1967). R. W. Barnwell, NASA Tech. Note, TN D-6283 (1971). A. W. Rizzi, A. Klavins, and R. W. MacCormack, Lecf. Notes Phys., 35, 341-346 (1975). A. W. Rizzi and H. E. Bailey, Proceedings of the AIAA 2nd Computational Fluid Dynamics Conference, June 1975, pp 38-46. A. Rizzi and H. Bailey, “Finite-Volume Solution of the Euler Equations for Steady Three-Dlmenslonal Transonic Flow”, presented at 5th Conference on Numerical Methods in Fluid Dynamics, Enschede, Holland, June 1976. P. Kutler, W. A. Reinhardt, and R. F. Warming, AIAA J., 11,657-664 (1973). W. C. Davy and W. A. Reinhardt, NASA Spec. Pub/., SP-347, 1351-1376 (1975). J. V. Rakich, AIAA J., 5 , 1906-1908 (1967). A. W. Rizzi, “Proceedings of the Symposium on Transsonicum”, Vol. 11, K. Oswatisch and D. Rues, Ed., Springer-Verlag, New York, N.Y., 1976, pp 567-574. J. E. Daywit! and D. A. Anderson, “Analysis of a Time-Dependent Finite-Difference Technique for Shock Interaction and Blunt-Body Flows”, Engineering Research Institute, Iowa State University, ERI Project 101, May 1974. P. Kutler, Lect. Notes Phys., 41 287-374 (1975). J. Daywitt. D. Anderson, and P. Kutler, “Supersonic Flow About Circular Cones at Large Angles of Attack: A Floating Discontinuity Approach”, AIAA Paper 77-86, Jan 1977.
2435 (22) R. W. MacCormack and A. J. Paullay, “Computational Efficiency Achieved by Time Splitting of Finite Difference Operators”, AIAA Paper 72-154, 1972. (23) R. W. MacCormack and R. F. Warming, “Survey of Computational Methods for Three-Dimensional Supersonic Inviscid Flows with Shocks“ in “Advances In Numerical Fluid Dynamics”, AGARD Lecture Series No. 84, Brussels, Belgium, Feb 1973. (24) L. B. Schiff, “The Axisymmetric Jet Counterflow Problem”, AIAA Paper No. 76-325, July 1976. AIAA 9th Fluid and Plasma Dynamics Conference, San Diego, Calif., July 14-16, 1976. (25) G. S. Deiwert, “Computation of Separated Transonic Turbulent Flows”, AIAA Paper No. 75-829, June 1975. (26) C. M. Hung and R. W. MacCormack, “Numerical Solutions of Supersonic and Hypersonic Laminar Flows over a Two-Dimensional Compression Corner”, AIAA Paper No. 75-2, Jan 1975. (27) H. Lomax and H. E. Bailey, NASA Tech. Note, TN D-4109 (1967). (28) R. J. Gelinas, J . Comp. Phys., 9, No. 2, 222-236 (1972). (29) K. G. Stevens, Jr., “CFD-A Fortran-like Language for the Illiac IV”, ACM SIGPLAN Notices, 10, No. 3, March 1975, pp 72-76. (30) Computational Fluid Dynamics Branch, “CFD A Fortran-Based Language for Mac IV”, CFB Branch, 202-1, NASA-Ames Research Center, Moffett Field, Calif. 94035. (31) W. G. Vincenti and C. H. Kruger, Jr., “Introduction to Physical Gas Dynamics”, Wiley, New York, N.Y., 1965. (32) J. F. Clarke and M. McChesney, “The Dynamics of Real Gases”, Butterworths, London, 1964.
A Systematic Method of Checking of Systems of Chemical Equations for Mass Balance Gweneth M. Ridler,” Philip F. Ridler, and John G. Sheppard University of Rhodesia, P.O. Box MP 167, Saiisbury, Rhodesia (Received May 12, 1977) Publication costs assisted by the University of Rhodesia
A method is presented for the systematic checking of systems of chemical equations for mass balance. Recent developments in numerical techniques have made possible the rapid and efficient solution of large systems of differential rate equations. Mass balance considerations need not necessarily be applied to chemical equations that summarize kinetically the effect of a number of chemical reactions. However, where chemically detailed mechanisms are involved, the accurate solution of the differential rate equations requires that mass balance be preserved for the complete system. The method described involves the writing of the system of equations in matrix form as the product of a matrix containing the stoichiometric coefficients of the reactions and a column vector of the chemical species involved. A Gauss elimination procedure immediately determines any mass balance inconsistency in the system of equations.
Introduction Detailed kinetic predictions based on the numerical integration of differential rate equations are playing a role of increasing importance in many areas of kinetics, particularly in aeronomy,l in combustion processes,2 and in the investigation of fundamental reaction mechani~rns.~ The computer simulation of chemical rate processes is most useful in those areas where detailed information about the constituent chemical reactions and their associated rate constants is available. In addition, there now exist well-established methods for the estimation of rate constants for postulated reactions, the detailed rates of which have either not been measured or are experimentally inacce~sible.~ The increasing use of computer simulation is also closely aligned to the development of numerical techniques for the rapid and accurate integration of “stiff” differential equation^.^^^ With the development of these techniques, the need to resort to oversimplified chemical mechanisms and mathematical approximations has largely disappeared.
Computer simulation of reaction mechanisms ranges from the use of models, of greater or lesser complexity, to the simulation of large systems of reactions in full detail. A good example of a model is the Oregonator7 which is a model of the Belousov oscillating reaction.8 The Oregonator and the Reversible Oregonatorg consist of shortened or incomplete chemical reactions and “rate constants” which are kinetically and mathematically necessary to describe the behavior of the intermediates X, Y, and Z, corresponding to HBr02, Br-, and Ce4+ in the detailed reaction scheme. Mass Balance It may be convenient at times to ignore mass balance considerations in limited model “chemical equations” that summarize kinetically the behavior of a small number of species in the overall system of reactions. However, the accurate solution of the differential rate equations for fully detailed chemical mechanisms requires that mass balance be preserved for the overall system. The set of equations The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
