Theoretical solutions of near wake flow have been obtained by application of the particle-in-cell numerical technique to hypersonic flow over the base of a cylindrical afterbody.
Using the inviscid equations of motion and a uniform upstream flow, the flow angle after separation was found to be much smaller than experimental slender-body results suggested. However, the inclusion of a vertical layer representative of a boundary layer on the body resulted in a realistic wake flow pattern, although the lack of true viscosity was apparent in the profiles of flow properties.
It is concluded that satisfactory near wake solutions by this method might be obtained by the use of realistic viscosity terms throughout the flow field and of a much finer representation of the approaching boundary layer, but the required computing power would be presently impracticable.
The last line on some pages are partially obscured and have had to be guessed.
1-1 Structure of the hypersonic near wake
2-1 P-I-C cell indexing system
2-2 Area weighting system for effective particle velocities
4-1 Particle configuration - run 0. t=17.4
4-2 Particle configuration - run 0. t=43.4
4-3 Particle configuration - run 0. t=93.4
4-4 Particle configuration - run 0. t=132.6
4-5 Particle configuration - run 0. t=195.0
4-6 Effect of local particle number density on velocity profiles
4-7 Variation of total system mass and energy
4-8 Velocity change on passage of starting shock
4-9 Temperature change on passage of starting shock
4-10 Particle configuration - run 1. t=17.4
4-11 Particle configuration - run 1. t=60.2
4-12 Particle configuration - run 1. t=111.8
4-13 Particle configuration - run 1. t=172.4
4-14 Particle configuration - run 1. t=213.4
4-15 Effect on cycle time of number of particles in system
4-17 Temperature profiles - run 1
4-18 Pressure profiles - run 1
4-19 Mach number profiles - run 1
4-20 Axial variation of temperature & Mach number run 1
4-21 Scale effect on wake properties
4-22 Particle configuration - run 2. t = 30.8
4-23 Particle configuration - run 2. t = 102.2
4-24 Particle configuration - run 2. t = 176.0
4-25 Particle configuration - run 2. t = 232.4
4-26 Particle configuration - run 2. t = 301.4
4-27 Velocity profiles - run 2
4-29 Temperature profiles - run 2
4-30 Pressure profiles - run 2
4-31 Mach number profiles - run 2
4-32 Change in velocity profile with time - run 2
4-33 Particle configuration run 3. t = 28.0
4-34 Particle configuration run 3. t = 185.2
4-35 Particle configuration run 3. t = 300.6
4-36 Particle configuration run 3. t = 382.4
4-37 Particle configuration run 3. t = 477.2
4-38 Velocity profiles - run 3
4-40 Temperature profiles run 3
4-41 Pressure profiles - run 3
4-42 Mach number profiles run 3
4-43 Reversed velocity in the base region
4-44 Tracks of initial particles in the base region
4-45 Tracks of inner boundary layer particles
The near wake region behind a hypersonic vehicle forms a complex flow field. Local Mach numbers vary across the free shear layer from perhaps ten or twenty to subsonic values, and Reynolds numbers may change by two or three orders of magnitude, from near freestream level in the outer, inviscid part of the wake to low values in the viscous core, Fig. 1-1 shows the main features of the hypersonic laminar near wake, typical of that behind a slender cone at Mach numbers of around 10.
Determination of laminar wake flow is an interesting fluid mechanical problem in its own right, but one so complex that only the recent need to understand atmospheric re-entry phenomena has prompted the considerable effort, both theoretical and experimental, which has been made during the last five years or so. Calculation of far wake properties depends on accurate knowledge of flow properties at the wake neck, as stressed in review papers by Lees (1) and Lykoudis (2); this can only be obtained from a complete solution of the near wake.
Weiss (3) has recently reported a complete solution involving iterative coupling of separate solutions for three different sub-regions of the near wake, each resulting from some degree of mathematical approximation. However, at the time when the present study was begun, theoretical analyses, such as those described in ref. 4, related only to the free shear layer, and their coupling to crude models of the inviscid flow produced unsatisfactory results except for very limited ranges of flow conditions. A promising new computational technique for compressible flows was applied successfully to wakes at low supersonic Mach numbers by Amsden and Harlow (5), and the present report describes an attempt to use this method, the particle-in-cell (or P-I-C) method, for axisymmetrc hypersonic wakes.
The P-I-C technique, described in section 2, is a combination of traditional Eulerian and Lagrangian numerical methods which solves the complete partial differential eouations of compressible fluid f1ow, with no simplifying approximations. A fuller description than is possible in this report may be found in ref. 6 and ref 7.
True viscosity and heat conduction may be included, but, for reasons explained in section 2, these effects have been omitted in the present calculations; the artificial viscous and heat conduction terms implicit in the method remain. The method suffers from two disadvantages from the point of view of this study. First, no part of the flow must be subsonic, since interactions are propagated only from cell to cell in the Eulerian mesh, so that the recirculating flow behind the base of the body is unlikely to be correctly computed. Second, there should be no small region within the flow field for which detailed resolution is required. Now the flow very close to the separation point at the rear of the body, within the subsonic part of the boundary layer, is thought to have an appreciable effect on the entire base flow, as suggested by Weinbaum (8). Such an effect could not be calculated in a P-I-C mesh of manageable size.
The attractiveness of the method lies in its simplicity of application, and the ease with which cylindrical geometries may be treated. In addition, transient effects may be studied (e.g. the starting of a shock tube flow round the base of a model) since solutions are developed in time steps from a specified initial state.
The present computations were aimed not so much at the prediction of near wake flow properties for a particular case, but were intended as an evaluation of the capabilities of the P-I-C method when applied to this complex flow problem.
Four complete P-I-C steady-state solutions have been obtained, representing a Mach 8.6 airflow over the base of a cylindrical afterbody. Effects studied include a change in the ratio of mesh size to body size, and the inclusion of a mock boundary layer represented by an appropriate amount of vorticity. Results are in the form of profiles of flow properties at various locations in the wake, and computer-generated pictures of the flow configuration. The program, was written in FORTRAN for the I.C.T. Atlas ccmputer, and in anticipation of output such as Fig. 4-11, was named Jabberwock after Lewis Carroll's creature which burbled as it came.
Section 2 of this report describes the P-I-C method, and its application to the current problem is dealt with in section 3. The results of the computations are discussed in section 4.
The flow field under consideration is subdivided into a number of cells, which, for ease of computation, are normally rectangular and equal in size. This basic layout may represent either two-dimensional or axisymmetric geometry; in the latter case, the mesh becomes one half of a thin radial slice, with the axis of symmetry along one edge of the mesh.
Two or more different fluids may be represented in the same problem, but since this is unnecessary for wake calculations where no dissociation occurs, the method will be described for a single fluid only. The fluid is composed of a number of discrete mass points called particles; in cylindrical geometry, the particles are actually rings, and the cells are toroids of revolution. Different masses are assigned to the particles in a cylindrical problem, the values being proportional to the original radius of the particle from the axis, in the same way that the volume of a cell will be proportional to its distance from the axis. Thus in either geometry, the number of particles in each cell, the particle density, is proportional to the true density of the fluid, at the beginning of a problem. However, the initial mass of a particle remains fixed, regardless of its later positions in the mesh, so that at a given time, the density of the fluid in a particular cell is calculated as the sum of the masses of all particles currently within the cell, divided by the volume of the cell. Only in plane geometry will this be proportional to the particle density.
Properties of the fluid are represented by a single value for each cell at a given time, these being the two components of velocity, the internal energy, density and its mass.
From a chosen initial configuration, where all cell and particle quantities are specified, the calculation proceeds in cycles. On each cycle, new values are computed for each cell quantity and particle position, these representing a slightly later time than the previous set of values. This is not an iteration procedure; the solution after each cycle is a stage in the transient behaviour of the fluid. The object of running a P-I-C calculation might be the determination of the unsteady dynamics of the fluid, though in the present investigation, the steady state solution was of major interest. It is assumed that a steady state has been reached when further calculation cycles produce, on average, no further changes in any of the properties of the system.
The calculations in one cycle may be conveniently split into two phases, each of which will be described below, after a brief outline of the scheme.
Phase I begins with the calculation of new cell pressures from values of density and specific internal energy, via an appropriate equation of state. For this phase, the particles remain fixed in position, and tentative new values of each velocity component are calculated for each cell, from the gradients of these pressures. In other words, the equations of motion are solved, with a finite-difference.approximation, omitting all transport terms from the equations. A finite-difference form, omitting transport terms, is used to obtain tentative new values of specific internal energy.
Omission of the transport terms is rectified in phase II of the cycle, when particle movement takes place. First, the phase I values of tentative velocity components and specific internal energy are converted into values, for each cell, of total momentum in each direction, and total energy (internal + kinetic). An effective velocity is then calculated for every particle, based on the velocity of the cell containing it and its position within the cell; the particles are moved to new positions determined by the effective velocity components and the cycle time increment. Should a particle cross a cell boundary, it carries with it its share of the total momentum and energy of its original cell, as well as its mass, which are subtracted from the totals for the old cell and added to those for the new cell. Finally, the new cell totals of mass, momentum, and energy are converted back to values of velocity components and specific internal energy. The elapsed time is increased by the cycle increment, and the necessary quantities are now available to begin the next cycle.
This section describes the calculation procedure used in the Jabberwock program. A cylindrical coordinate system is used, with the x-axis as the axis of symmetry; the radial coordinate is denoted by y, and cells are identified by indices i and j, increasing in the x and y directions respectively. The x and y velocity components are denoted by u and v, as usual. Specific internal energy of the fluid has been taken as I = cvT, and the perfect gas equation of state has been used.
The inclusion of real viscosity in the P-I-C finite difference equations, although simple in principle, introduces complexities in the difference equations and boundary conditions which would have involved much extra programming effort; this was considered unwise during a first attempt at writing a P-I-C code. However, the chief reason for its omission is simply one of computing power. Accurate simulation of the boundary layer at the rear of a hypersonic vehicle would require a ratio of boundary layer thickness to base radius of the order of one tenth, and the number of cells required for a reasonable representation of the variation of properties through the layer would perhaps be about five. Thus a flow field covering the region 0 ≤ y/r ≤ 2, 0 ≤ x/d ≤ 3, would require of the order of 50000 cells, and a computer which could handle such a problem, producing an answer within three or four hours, would need to have 8 times the store and 40 times the speed of Atlas! An alternative would be just to let a boundary layer develop automatically along the wall, but this would be no more likely to represent the true features of the boundary layer separation than altogether omitting the layer.
There remains, however, the artificial viscosity implicit in the P-I-C method, which would first have to be subtracted if the inclusion of real viscosity were to be realistic. The terms of the finite-difference equations may be expanded about some central space and time in a Taylor series, and this results in the original differential equations of fluid motion plus some additional terms about a cell centre, and retaining the lowest order terms in the spatial finite-difference increments δx and δy, yields the following momentum and energy equations:
The effective artificial viscosity and heat conduction coefficients, ½ρ|u|δx and ½ρ|u|δy, are non-isotropic and are proportional to the fluid speed and cell size. In regions of stagnation, therefore, the viscosity becomes negligible, but P-I-C calculations may often be run successfully with no additional explicit artificial viscosity terms. This was the case with the wake flow calculation reported in ref. 5, and no additional terms were used in the present problem. The effective artificial viscosity is adequate as a dissipative mechanism for calculating shock waves.
In the free shear layer region, the important viscosity coefficient is expected to be ½ρ|v|δy, so that for a P-I-C base flow model using a base dimension equivalent to 10 cell lengths, a free shear layer Reynolds number may be constructed, equal to
For typical shear layer angles, this number will be around 50 to 100, whereas a very rough estimate for the same figure from a typical hypersonic shock tunnel test on a slender body gives a Reynolds number of order 104. However, the flow model in this problem has no viscosity upstream of the base of the body, and it becomes meaningless to attempt to predict how viscous a flow the P-I-C solution really represents.
The dependence of the artificial viscous terms on fluid speed is one reason why far subsonic flow cannot be calculated. As the speed decreases, the viscosity approaches zero, and the difference equations become unconditionally unstable; fluctuations of the variables will then grow until the velocities are large enough to produce sufficient dissipation.
Fig. 2-1 shows the cylindrical coordinate system used, with the method of cell indexing. Conversion to plane geometry requires no alteration of the equations, but is effected by putting the radial coordinate y equal to 1, and n equal to ½. The difference equations are basically those given by Amsden (7), re-written for the different cell indexing system used here. Steps in the computation procedure will now be described in more detail.
(a) Cell densities are calculated:
(b) Since the particle model automatically conserves mass, the continuity equation need not be considered. Tentative new cell velocities are now calculated from the compressible flow momentum equations:-
Tentative quantities in this phase of the calculation are denoted by an overhead tilde (~), and will be referred to as tilde quantities. The appropriate finite-difference approximations are:-
where yj is the radius to centre of cell (i,j).
Quantities with the halves appearing in the cell indices are cell boundary values, defined as the average of the values in the cells on either side, i.e.
If a neighbouring cell is empty, the boundary pressure is set to zero, in order that the only momentum changes in the system arise from forces on the system boundaries. Where the boundary pressure is equated to that of the cell.
Should a certain tilde velocity be too great, it becomes possible in Phase II of the cycle for particles from that cell to remove more energy than the cell possesses, resulting in negative internal energies. To minimise this possibility, each new tilde velocity is compared with a maximum allowable value, and if it is too large, the time increment δt is reduced, and the whole tilde velocity calculation is restarted. A stage may be reached where δt becomes so small that computer real-time is excessive; two features of the present program counteract this possibility. First, the time increment is automatically restored to its original value five cycles after each reduction, and second, a minimum allowable value of δt may be specified. If excessive tilde velocities still exist at this minimum ot value, then they are set equal to the maximum allowed velocity. This results in a reduction of kinetic energy, and increase in internal energy in the cells concerned, and thus acts as a kind of artificial viscosity. Momentum is not conserved when this happens, but in practice, occurrences of this velocity restriction were found to be sufficiently rare for the accuracy of the calculation not to be seriously affected.
(c) Having calculated u~ and v~ for all cells, the tentative new specific internal energies are computed from the energy equation:
where E' = I + ½ (u2 + v 2).
Since ∂/∂t(½(u2 + v2) ≡ u∂u/∂t + v∂v/∂t, the energy equation becomes, after omission of the transport terms,
For greater stability, the velocity components used in the finite-difference form of this equation are an average of the old and the tilde values. The difference equation is:
where
and
If a neighbouring cell is empty, e.g. cell (i + 1,j), the following form is used
and similarly for other bar quantities. At a cell boundary which forms part of a rigid wall or the axis of symmetry, the bar quantities are put equal to zero.
(a) Specific quantities are transformed into cell totals. From values of u~, v~ and I~, the total momentum in each of the two directions and the total energy for each cell are calculated. Being tentative values, the momenta X and Y, and energy E are also labelled with tildes, and are found from the equations:
and similarly for E~.
(b) Movement of the particles now takes place. For each particle, effective velocity components are determined by a weighting procedure, dependent on the particle's position within the cell. An imaginary cell is considered, centred on the particle, as shown in Fig. 2-2. A1, A2, A3 and A4 are the areas of overlap of this cell on that containing the particle and the three cells nearest to its position. If A, = δxδx, is the area of a cell, then the effective velocity components for the particle are calculated as:
and similarly for veff. u~1 is the tilde x-direction velocity component in the cell denoted as 1 in Fig. 2-2, etc.
The new particle coordinates after movement are calculated as:
Boundary conditions have been dealt with as follows. If the imaginary particle-centred cell overlaps an empty cell, the velocity of that cell is put equal to the velocity in the particle's own cell, for the purposes of the weighting formula. Should the imaginary cell be partly outside the computing mesh, the fictitious outside cell is given the velocity of the adjacent interior cell, and this may result in a particle finding itself outside the mesh after movement. If the boundary at this point is a rigid wall or the axis of symmetry, it is then necessary to reflect the particle back inside the mesh; imagining the wall to be a mirror, the final new position of the particle, inside the mesh, is the mirror image of its location immediately after movement with the effective velocity. Its velocity component normal to the wall is now reversed in sign, and the necessary change in the total momentum component of the cell concerned is made at this point in the calculation.
After movement, a particle may remain in the original cell, in which case no alteration to any cell quantities is required, but some of the particles will have crossed boundaries into new cells, and must transfer their mass, momentum, and energy as follows. If a particle of mass m moves from cell 1 to cell 2, then
Cell values of mass, momentum, and energy have their final value for the cycle once all particles in the system have been moved and the necessary adjustments made to all the cell-wise quantities in use at this stage.
(c) Phase II of the cycle is completed by converting all cell quantities back to their specific form. Final momentum components are divided by final cell masses to obtain final velocities, and final kinetic energies are subtracted from the final total energies to give the final specific internal energies.
All cell quantities required to begin the next cycle of the calculation are now available. Before beginning the next cycle, however, it is useful to perform some auxiliary calculations. For example, the total mass and energy of the system may be computed by summing the flux of these quantities across the system boundaries. During the development stage of the Jabberwock program, several programming errors were detected by noting anomalous values of various properties of the system as a whole.
The present problem required a continuous input of fluid particles above the base of the body, and a continuous output at the downstream end of the mesh.
A fictitious column of cells just upstream of the input bonndary, possessing the prescribed input flow conditions, was used in order to treat the first real column of cells in the mesh as true interior cells. A new column of input particles, each contributing the required input value of momentum and energy was generated whenever the fluid already in the system had moved a sufficient distance downstream. The test for inclusion of a new column was made after completion of the particle movement on each cycle, but before cell quantities were transformed to specific values.
The last column of cells before reaching the output boundary was also treated as a true interior column, but this time the cells in the fictitious column just outside the boundary were given properties identical at the particular time with those of the adjacent interior cells. Thus no gradients within the system were caused by this boundary. On crossing the output boundary, a particle's mass, momentum components and energy were subtracted from those of the cell it had just left, and arrangement made for the computer storage locations relevant to that particle to be made available for a new input particle. In this way, the number of particles in the system is constrained during the calculation regardless of the total number of particles that will have entered the system by completion of the problem.
The finite and relatively small number of particles in the average cell of the mesh results in rather severely quantized cell densities, and the pattern formed by the incoming particles must be arranged so that a large number of particles does not cross the same cell boundary in a given cycle; this would lead to violent spatial density fluctuations in a region where uniform density may be required. In the Jabberwock program, particles in a single input column were staggered in the x-direction (the direction of initial motion), but uniformly distributed in the y-direction. Alternate columns were inserted with the same pattern, but shifted in the y-direction, so that several distinct density steps were available in each direction.
In a problem requiring a continuous throughput of particles, it would seem that the contents of the cells at time zero has little or no effect on the steady state solution, since the majority of the initially stationary particles filling the mesh will have been swept away by the incoming flow. A staggered one-particle-per-cell arrangement was used in this problem to fill the mesh initially, except for the last run, where two particles occupied each cell.
Generally speaking, it is desirable to have the cycle time increment δt as small as possible, but a compromise has to be reached in order to avoid using excessive amounts of computer time on the calculation.
The equations used in Phase I of a P-I-C calculation are basically unstable, giving rise to a statistical fluctuation in the cell variables. The success of the method depends on the existence of a statistical averaging effect, which is enhanced by a small value of δt, so that the particle movement vectors experience many different orientations relative to the mesh in a given elapsed problem time. If umax is the greatest fluid speed in the system, and δx is the cell dimension, then in order that many steps are required to move a particle across a cell,
(umax δt) / δx ≪ 1
In all previous P-I-C calculations, catastrophe has resulted from even the slightest violation of the condition that pressure signals must not travel more than one cell width in one cycle, i.e.
a δt/δx < 1
where a is the local speed of sound.
A further restriction on δt may be made in order to minimise the possible occurrence of negative internal energies. Amsden (7) gives the following criterion for avoiding negative energies when the simple equation of state ρ = RT is used:
2(γ-1)umax δt(1/δx + 1/δy) < 1
This criterion was used to select a suitable δt value, but occasions can arise when isolated velocities occur which are greater than the estimated maximum velocity. Automatic reduction of the time increment may cure the trouble, as described earlier in this chapter.
There exists no analysis of the P-I-C method which shows in detail the nature of the errors to be expected, and it has been reported that a feeling for the accuracy and applicability of the methodhas been obtained only by a trial and error process, experimenting with different forms of difference equation, for example. It is fortunate that the most satisfactory conservative finite-difference forms had already been determined by the Los Alamos team, well before the present calculations were planned.
In a sense, a P-I-C calculation resembles a physical experiment in that the computer develops a solution at a sequence of later times after being given initial and boundary conditions together with a set of differential equations of motion. Such a computation is quite different from those often performed by computer in fluid flow problems, involving the precise solution of a set of ordinary differential equations which result from analytical approximation procedures. In contrast, the solutions of partial differential equations by the P-I-C technique are always approximate solutions.
One can say very little concerning the accuracy of the method other than to mention that it has been observed, though not proved, that in many circumstances of interest the approximate solutions are good, and can be improved by decreasing the size of the cells in the Eulerian mesh.
Accurate simulation of the flow over the rear of a typical slender body at hypersonic speeds is virtually impossible with the P-I-C method, unless the whole of the body is included in the computing mesh. This arises from the need to specify completely the boundary conditions at the upstream end of the mesh, so that the nose shock wave cannot be included if the nose of the body is not itself within the mesh.
The method can become very complex if the rigid walls representing the body do not lie along cell boundaries. It was thought advisable to simulate only parallel-sided two-dimensional bodies, or the equivalent in cylindrical geometry, cylinders aligned axially with the flow. All nose effects are absent, so the configuration is equivalent to a wind tunnel model which extends upstream into the throat of the nozzle.
Omission of true viscosity implies the absence of a boundary layer on the wall just upstream of the base, since the effective artificial viscosity term in the relevant direction is proportional to the cross-stream velocity component, v, and is therefore equal to zero. In relating such a flow model to any physical situation, one must consider that the body is sufficiently long for nose effects to be absent, yet short enough for any boundary layer to have a negligible thickness compared with the cell dimension. In ref. 7, Amsden gives details of a P-I-C wake calculation with real viscosity incorporated, showing a strong shock wave layer begins to form. However, the ideal requirement for a wake calculation is a fully developed boundary layer already in existence at the upstream boundary of the mesh, with no sudden interaction between this layer and the undisturbed flow above.
After two runs with no boundary layer representation, the Jabberwock program was developed to incorporate a mock boundary layer meeting the requirements mentioned above. This was achieved by introducing the vorticity representative of a boundary layer, while using the same inviscid computation scheme. Ideally, this layer should be several cells thick in order to simulate an initial velocity profile, but the limitation on computer store means that this thickness would have to be of the same order as the base radius. To obtain a realistic ratio of boundary layer thickness to base radius, the plan adopted was to restrict the vorticity to the row of cells adjacent to the wall, by specifying an input velocity for this row equal to half that of the free stream above. The flow in this mock boundary layer, after suitable adjustments to the fluid properties (described later in this section), is still supersonic.
Justification of this development may be found in the recent theoretical work on hypersonic boundary layer separation by Weiss and Weinbaum, described in reference 9. This work indicates that pressure gradients are far more important than viscous stresses in determining the fluid motion in the supersonic portion of the freshly separated flow. In the context of the present problem, it would seem that the inclusion of real viscosity would have little influence on the separation process, but it should be remembered that the development of the free shear layer does depend strongly on the viscosity, and it is in this region of the flow that the effective artificial viscosity of the P-I-C method can become large.
Only that portion of the flow field on one side of the plane or axis of symmetry need be calculated, and the axis will then form, in the case of this problem, the lower boundary of the mesh. The half height, or radius, of the base of the body should be not less than about five cell dimensions if a reasonable resolution of flow properties is to be obtained, and the upper boundary of the mesh must be sufficiently far from the axis for no significant reflection of the corner expansion wave to take place. Thus, the longer the mesh in the streamwise direction, the greater must be its cross-stream dimension. Having chosen a ratio of cell size to base dimension, the downstream distance (in dimensionless form) to which the wake ray be calculated is limited by the available computer memory.
For the first calculation (run 0), a mesh 50 cells long and 20 high was chosen, the cells being square for simplicity, with a base radius equal to 10 cell dimensions. A short length of the body, equal to 2 cells, was included in the mesh, so that the downstream boundary was located 2.4 diameters from the base; Fig. 3-1 (a) shows the configuration. On all subsequent runs, the base radius was reduced to the equivalent of 5 cells. Runs 1 and 2 used a mesh of 55 ×15 cells, placing the output boundary 5.3 diameters from the base, while on run 3 a 40×12 mesh was used, the maximum downstream distance becoming 3.8 base diameters. The mesh for runs 1 and 2 is shown by the solid lines in Fig. 3-1(b) and for run 3 by the broken lines. Fig. 3-1(c) shows the relation between the meshes and the overall geometry.
Properties of the uniform flow over the body are somewhat arbitrary, since the solution can be compared with experimental measurements only in a qualitative manner. The chosen flow actually represents a freestream flow theoretically attainable in the R.A.E. 7 in. hypersonic tunnel, and has the following properties.
Mach number M1 = 8.60 Static pressure p1 = 4.60 lbf/ft2 (220.2 N/m2) Static density ρ1 = 2.69×10-5 slug/ft3 (0.0139 kg/m3 Static temperature T1 = 55.3°K Speed u1 = 4210 ft/sec (1284 m/sec)
(This corresponds to stagnation conditions 500 p.s.i.a. and 600°C.)
These properties were scaled in such a way that the freestream density and velocity were both numerically equal to one. The base diameter was taken to represent 2 ins., and the length scale taken such that a cell was of unit length. By scaling length, velocity and density in this way, the scaling for all other properties was uniquely determined. For those runs using a base radius equivalent to only 5 cells, the flow properties were left unaltered in problem units, so a real base diameter of 1 in. in the same flow was then being simulated.
The experience of the Los Alamos team in running particle-in-cell problems suggests that an average over the whole mesh of about 4 particles per cell is desirable in the steady state configuration, though of course, the more particles per cell one provides, the smoother will be the property profiles of the flow. Again, computer storage space places a limit on the number of particles, and allowance must be made for the fact that for an impulsively started problem of this type, the maximum number in the system (as the starting shock wave is about to leave the mesh) can be almost half as many again as the number in the steady state. Also, the amount of computer time required for one cycle of the calculation is almost proportional to the number of particles in the system at the time, since the number of arithmetic operations required to update all cell quantities is much smaller than that required for operations on the particles.
An upstream particle density of 6 per cell was initially chosen for the current problem, in the hope that this would result in a reasonable resolution of flow properties in the wake, and at the same time produce a steady state configuration in a real time of two or three hours.
A total of 12 different cell quantities is required for a P-I-C code of the type used here, namely: M, ρ, p, I, I~, E~, u,u~, v, V~, X~, Y~. However, not all of these are required continuously throughout the cycle, and a sharing scheme may be devised to reduce the number of words of store to 7 per cell. In the Jabberwock program, 8 words per cell are actually used . Pressure and tilde internal energy could have shared store, but it was convenient to keep values of p for printing at the end of a cycle.
Each particle requires three words of store: two for its coordinates, and one for its mass (though not necessary for a calculation in plane geometry). Also required is a label so that particles may be indexed; this is normally provided simply as the subscript of the storage array element referring to that particle, but a scheme was envisaged at one time in which a scalar variable in a subroutine of the program would have the value of the particle's mass. The identity of that particle would then have been lost, so to avoid this, the quantity stored in the array was not simply the mass of the particle, but a number whose integer part consisted of the indexing label, and fractional part of one hundredth of the mass. (The quantity mass/100 was always less than one, yet large enough to enable all digits to be stored even for the lightest particles).
The major part of the store required for variables was therefore 8 times the number of cells, plus 3 times the maximum allowed number of particles. A large array was also kept to hold the labels of particles which had left the mesh, for use by fresh input particles. This totalled approximately 31000 words for run 0, 30000 for runs 1 and 2, and 36000 for run 3. (See the following subsection for the meaning of these run numbers). The program itself occupied about 9000 words of store, and a further 5000 words approximately were used for buffering input/output streams and the magnetic tapes.
Altogether, 44K words were used on run 0, 43K on runs 1 and 2, and 49K on run 3. (K = 1024 words).
This section describes the variations on the basic flow configuration which were made from run to run.
Though the program was written to allow either plane or cylindrical geometry, the latter was used in all cases; the properties of the uniform input flow were also identical for each run.
Run 0:-
This was the development run (hence the numbering), though the final solution is a valid one, having been reached from a fresh start while the only remaining errors concerned the system for retrieval of labels from used particles. A 50×20 cell mesh was used, covering the wake as far as 2.4 base diameters downstream.
Run 1:-
The results of run 0 showed that it was necessary to simulate the wake to a greater downstream distance, so a 55×15 mesh was introduced, and the base radius reduced from 10 to 5 cell dimensions, placing the output boundary 5.3 diameters from the base. Again, the input flow contained 6 particles per cell.
Run 2:-
The steady state flow pattern resulting from run 1 differed markedly from the expected form of wake flow for such a configuration, as discussed more fully in section 4. An attempt was therefore made to incorporate a representation of a boundary layer on the body, by introducing vorticity in the input flow. The row of cells adjacent to the wall was given a velocity component equal to half that of the free stream cells, and the density and temperature were specified to give the mock boundary layer the same static pressure and total temperature as the main portion of the input flow.
The mesh and input particle number density were as for run 1.
Run 3:-
This was a repeat of run 2, using a much larger number of particles, 16 per cell, in the input flow in an attempt to eliminate some of the voids appearing in low density regions of the flow. Superfluous portions of the mesh were removed, and the wake calculated as far as x/d=3.8. Certain other modifications were also made.
The Jabberwock program was written for the I.C.T. Atlas computer at the Science Research Council's Atlas Computer Laboratory, Chilton. This is a large high-speed machine with sufficient memory to enable the program to be run entirely with core and magnetic drum store. Transfers to and from the drum are completely automatic, and the program has been written as if for core store only.
Previous P-I-C codes reported in the literature have been written in machine codes which, while producing a more efficient program, are far more complex to use than a high-level language. It was decided to use FORTRAN for this program, for ease of writing and development, and to take advantage of the powerful and flexible operating system, HARTRAN, used at the Atlas laboratory.
Since each calculation was expected to take of the order of 2 or 3 hours computing time, and since frequent breaks were desirable, to assess the progress of the solution, provision had to be made to run the calculation as several separate jobs. The program incorporates a timing device which allows each subrun to be as long as desired in terms of computer time; having completed an integral number of P-I-C cycles, all the information necessary to restart the calculation is written to magnetic tape. This information consists of 4 arrays of cell variables (M, I, u, v) and 3 of particle quantities, plus all constants, etc. The program itself was kept in binary punched card form.
Where a subrun was required to take more than ten minutes computer time, the restart data was written to magnetic tape every ten minutes, so that in the event of a machine failure, no more than ten minutes calculation was lost. As a further safety precaution, two magnetic tapes were used for every subrun, being used alternately for the restart data. Then in the event of an accident to one of the tapes, again no more than ten minutes of caputing time has been lost.
Between subruns, it was usual to copy the restart data to another part of one of the tapes, thus forming a library of various stages in the solution. This was extremely useful in the development stage, when the effect of alterations to the program could be seen without having to make a completely fresh start to the calculation each time.
It should be mentioned here that although the maximum. total store required was in the region of 50000 words, ie, about 49K, including input/output buffers, etc., as much as 100K is normally available to users of Atlas. The Jabberwock program could therefore have been more ambitious in tenrms of mesh size and number of particles, and the decision to limit its size was based on purely practical reasons. Limiting the store request helped to some extent in speeding the throughput of work (Atlas has a multi-programming facility), but the main advantage was in the overall time required to obtain the final steady state solution, remembering that six or more long subruns may be required for each calculation made. Doubling the number of particles in a calculation would necessitate twice as many such subruns
The package of routines for use with the SC4020 microfilm recorder, required another 12000 words of store approximately. On early test subruns, these were included each time, and graphical output was generated as part of each subrun, However, this scheme required 56K store, and was soon replaced by running separate jobs to obtain this output, reading the necessary data from one of the tapes used to hold the restart information. These jobs needed about 34K store, but enabled the actual calculation programs to be run with the basic 44K store request, and with two magnetic tape decks instead of three. Turn-round times were considerably improved after implementation of the new scheme
The most useful form of output has been a plot showing the positions of all particles at a given stage in the problem. This was produced by writing all particle coordinates, and the data required for generating a grid, onto magnetio tape, using e set of FORTRAN routines supplied by the Atlas Laboratory for use with a Stromberg Carlson SC4020 microfilm recorder. The magnetic tape is used off-line to drive the recorder, which produces a permanent copy of the plot on paper, microfilm, or both. Some of the early plots were made on a Benson Lehner model J electroplotter (a conventional automatic graph plotter), and at other times, use was made of the Benson Lehner 120 microfilm recorder at the U.KA.E.A, Culham Laboratory.
Provision was made in the program for storing on magnetic tape the labels of those particles entering the mesh at certain selected points on the input boundary. With this information, it is possible to produce plots showing streaklines in the flow, though the facility was never used, for two reasons: observation of the particle position plots gave a good idea of the streakline pattern, due to the initial pattern of the input particles, and also the procedure would have increased the total computing time.
Further checks on the progress of the problem were made by printing values of significant quantities, such as total system energy, number of particles in the mesh, etc., every few cycles. Events such as the input of fresh particles, the occurrence of negative internal energy values, etc., were also recorded on the output stream, which thus formed a running commentary on the problem. A second output stream from the lineprinter was used to print profiles of cell properties at selected stations after the last cycle of each subrun. These were a useful supplement to particle position plots, especially as the particle density as seen on the plot is not proportional to the actual fluid density for a calculation in cylindrical coordinates, Because of the statistical fluctuations in cell properties, a profile of instantaneous values may not represent the short-time average state of the fluid at that stage of the solution, but general trends are well indicated. Final steady state profiles were always produced using a special program which kept running totals of all cell property values and printed an average value taken over 50 cycles; this was considered long enough to smooth out the large proportional changes from cell to cell caused by low particle number densities in some regions of the mesh.
In this section, the results of each of the four runs made using the Jabberwock program will be described, together with more general observations relating to all the calculations.
When interpreting the results, the following points should be noted. Since all runs were made in cylindrical geometry, the particle number density seen on the configuration plots is not proportional to the fluid density. An empty cell implies that the fluid in that region has a density less than that represented by one particle (of lowest mass) per cell, but the calculation will yield zero values for all properties in that cell. On graphs of property profiles, points referring to empty cells have been plotted as zero, except for Mach number, which is then of indeterminate form. The steady state profiles have been plotted as the average over 10 problem time units, resulting in smoother profiles in low density regions where a cell may contain particles for some of that time, but remain empty for the rest of that period.
Particle configuration plots showing the transient behaviour are given in Figs. 4-1 to 4-5. At t=17.4, the starting shock is about one third of the way down the mesh, and has already begun to sweep particles away from the base region. The pattern of curved lines at the front of the incoming flow is caused by bunching of the initially stationary fluid particles which occupied the upstream end of the mesh, and the division between this fluid and the input flow is marked by the line of highest particle number density, lying somewhat behind the shock front. By time t = 43.4, the shock has spread to the axis, and the base region is emptying further.
Between X = 10 and X = 20, the input flow is seen to have a more random pattern than that upstream or downstream of the region; this is due to an experiment where the time increment δt was allowed to be reduced indefinitely. For most of the time during this period, δt assumed values so low as to render the calculation quite impracticable, owing to a programming fault whereby tilde velocity components already set at the limiting value caused further reduction of δt on the following cycle.
The total number of particles in the system reaches a maximum as the starting shock is about to cross the downstream boundary, the number being about 4500 on this run. At this point the computer time required for one cycle also reaches a maximum of 7.5 seconds, then falls as the number of particles becomes less.
At time 93.4, the plot shows a uniform flow along the whole of the mesh for Y > 10, while those particles in the region Y < 10 are mostly those which formed part of the initial stationary fluid. By looking in the direction of the flow from the upstream end, at a shallow angle to the paper, streaklines may be seen, caused by preservation of the specified particle input pattern, which is not quite uniform in the Y direction. Bunching of the innermost lines at Y = 10, downstream of X = 35, may be interpreted as a shock wave, formed as the expanding flow is turned away from the axis on meeting the higher pressures near the axis at the downstream end of the flow field.
This shock has virtually disappeared at time 132.6, since many particles near the axis have now been swept out of the mesh, resulting in lower pressures in the wake. The low pressures cause the innermost portions of the input flow to expand through a larger angle, and the flow pattern at this time has almost reached its steady state configuration. (The symbols ∂ and Q which appear on the plot for t = 132.6 were caused by faulty operation of the Benson-Lehner 120 plotter).
By time 195.0, no significant change in flow pattern had occurred during the previous 60 time units; the total energy of the system had fallen by only 0.1% since t = 160, and the total mass in the system by the same amount since t = 180. It appeared that the steady state had almost been reached, and since the expansion angle of the flow was too small for the wake neck to be included in the calculated flow, the problem was stopped at this time.
The number of particles in the base region is far too small to allow a proper representation of the recirculating flow, so property profiles averaged over five time units provide no information on the flow close to the axis, y/r less than about 0.5. An example is shown in Fig. 4-6, which compares velocity profiles at an axial location x/d = 1.9 at two different stages in the solution. On a single cycle at t = 93, a smooth profile exists down to about y/r = 0.3, whereas the transition from free stream velocity, u1 to zero takes place over a much shorter distance at the end of the problem, and all cells below y/r =0 .7 have zero velocity, being empty.
It will be seen later that the scaled-down version of the same problem in run number 1 produced almost identical results at the same values of x/d, so no further profiles will be shown for run 0.
Fig. 4-7 shows the behaviour of the total mass ΣM and total energy ΣE of the system from t = 44 to t = 108. Until t = 58, no particles have left the flow field, and the rise of ΣM and ΣE is exactly linear with t; this as expected, since the properties of the input flow are constant, and the particle model conserves mass exactly while the difference equations are in a form which exactly conserves total energy. Since the specified free stream velocity is scaled as 1, fluid which is not decelerated should travel the length of the mesh (50 cells) in time 50; however, the first particles to enter the mesh at the beginning of the problem encounter a stationary fluid, so are retarded. The time taken for this first wave of particles to reach the output boundary is therefore rather greater than 50 units.
Shock waves calculated by the P-I-C method are smeared over several cell widths, as is the case with all finite difference methods possessing effective artificial viscosity. This is illustrated in Fig. 4-8, which shows how the velocity u in cell (2, 16) builds up as the starting shock wave passes through the cell early in the problem. About 8 time units elapse before the velocity reaches the free stream value of 1. The behaviour of the specific internal energy of that cell is similar, and is plotted (as temperature ratio) in Fig. 4-9. On this run, the specified temperature of the initially stationary fluid, occupying the mesh at t = 0, was equal to the free stream static temperature T1, but on the three subsequent runs a higher value was used (273°K), equal to 4.95T1. The peak temperature of a cell in the column i = 2 was then in the region of 7T, compared with about 6T, in Fig. 4-9, but the time taken for the internal energy to settle to free stream level was almost exactly the same as on run 0.
Although it was apparent from the results of run 0 that a simple P-I-C model, without boundary layer representation, would not produce a steady-state flow having a wake neck within one or two diameters of the base, it was decided to extend the flow field in order to calculate such a wake neck as might occur in the region x/d = 4 or 5.
Run 1 was virtually a half-scale repeat of run 0, as described in the previous chapter; this is equivalent to halving the size of a wind tunnel model in order to produce a wake twice the previous length in terms of base diameter. Particle configuration plots for this run are shown in Figs. 4-10 to 4-15, and it will be seen that some degree of randomness develops in the stationary fluid ahead of the starting shock. This is due to a particle placing pattern which did not repeat from one cell to the next; initial densities and pressures therefore differed slightly from cell to cell, causing a low-speed particle movement tending to equalise the pressures. The steady-state flow pattern should not be affected by this movement.
The plot at t = 60.2 shows that a wake shock is forming, with origin at about X = 20, Y = 5, before the starting shock has left the mesh. At t = 73.0, the number of particles in the system reaches a maximum value of 4590, and the computer time required for one cycle is now 7.7 seconds. The way in which the number of particles affects the rate of computing is clearly shown by 4-15 which plots the Atlas time, Tc, required for one P-I-C cycle, and the total number of particles in the mesh against the problem time t.
The next plot, at t = 111.8, shows the base region virtually empty, while the origin of the wake shock has moved further downstream. All the original particles, with one or two exceptions, have been swept out of the mesh by time 172.4, when the plot shows only the very innermost streamlines of the flow being deflected at the wake neck. At this time, the total energy and mass in the system are decreasing very slowly, and a steady state has almost been reached.
Computation was stopped at t = 213.4, with the flow pattern, as seen on the plot at this time, almost identical with that at t = 172.4. During the final 20 time units of the run, the total energy of the system varied randomly by ±0.047%, and the total mass by ±0.050%, further proof of the existence of steady flow. By the end of the run, a total of nearly 13000 particles had been input at X = 0, and the number in the system, in the steady state, was 3280, about 71% of the peak value; real time per cycle settled at 5.75 seconds. The total computer time used for the run was about 1¾ hours.
Reduction of the time increment δt occurred about ten times up to t = 150, two applications of the tilde velocity limit being made on six of these occasions. For the remainder of the run, no tilde velocity component ever reached the specified limiting value, allowing δt to maintain its maximum value. This situation was an improvement over than in run 0, when δt was still being reduced quite frequently at the end of the run; perhaps the reason lies in the fact that the empty base region on run 1 represents only about half the number of cells contained in the same region on run 0.
Property profiles in the steady flow were this time averaged over ten time units, and some of these are presented except where the cells near the axis at the particular x/d location are empty (u = v = 0), all cell values of u were very close to 1, the free stream value. Radial components, v, are all effectively zero outside the leading Mach wave of the corner expansion (actual values being of the order of 10-5 ) and reach maximum negative values of around 0.17 for x/d = 3, 0.09 for x/d ≈ 5. In the row j = 1 (closest to the axis), v is negative where the flow first meets the axis, becoming zero then positive as x/d increases; in cell (55, 1), next to the axis on the output boundary, v is about + 0.015.
A formal definition of the wake neck location may be provided by the axial position where v passes through the value zero in row j = 1. On run 1, this gives (x/d)neck = 4.2.
Density profiles are seen to be fairly smooth except at x/d = 3.3 (cell column i = 35); reference to the particle configuration plot for t = 213.4 (Fig. 4-14) shows that this is probably due to the way in which the streamlines of particles become separated, producing undesirable voids in the flow. The position of these voids in relation to the cell boundaries will clearly play a large part in determining the actual cell values of density. Close to the axis, high densities occur, ρ/ρ1, being about 3.4 in row j = 1 for x/d values greater than about 4; this rather unrealistic result is hardly surprising, since all the particles in this region began their existence at radii y/r > 1, and retain masses appropriate to those radii on moving into a region where cell volumes become much smaller. In order to keep a reasonable scale, these high densities have been omitted from Fig. 4-16.
Much less scatter is present in the temperature profiles, drawn from cell specific internal energy values. A peak value of T/T1, occurs, in row j = l, at x/d = 4.3, i.e. approximately at the neck of the wake. Fig. 4-20 shows the axial distribution of T/T1, on a single cycle at t= 213.4, together with the Mach number distribution at the same instant. Since the x-components of velocity on the axis are all very close to 1, and the y-components are small, Mach numbers in this region are determined almost entirely by the temperature. In common with other flow properties, Mach number shows very little variation from its free stream value of 8.6 in all non-empty regions of the mesh, as can be seen from Fig. 4-19.
The pressure profiles bear a very close resemblance to those of density, which is to be expected since variations of temperature with y/r are relatively small. Values of p/p1 of about 3.2 are reached in row j = 1 at the output end of the mesh; again, these high values have been omitted from the graphs. Strong normal pressure gradients exist in this region, acting to produce a flow away from the axis, and hence the increase in v in row j = 1 from the wake neck to the outflow boundary.
Two major points emerged from the results of run 1. First, it seems that realistic values of flow properties in the wake cannot be obtained when the input flow has uniform properties right down to the body wall, and second, an initial particle number density of 6 per cell is not sufficient to prevent the occurrence of undesirable voids in low density regions of the wake.
Comparison of the results of runs 0 and 1 will show the effect, if any, of changes in the mesh size relative to the base dimension. Strictly speaking, if runs 0 and 1 are to represent solutions of flows having the same free stream time t = t1, on run 0 should be compared with that at t = ½t1, on run 1, because the flow speed relative to the mesh was the same in both cases. However, once a steady state has been reached, this no longer matters. One might expect, from this argument, that the final flow configuration on run 0, at t = 195.0, should correspond with the flow at t = 97.5 on run 1, but comparison is not valid here because the transient solutions are strongly influenced by the density of the fluid filling the mesh at t = 0, and also by the particle number density used to represent this density.
Because the effective artificial viscosity of the P-I-C method is proportional to the cell size, the flows of run 0 and run 1 cannot be considered similar merely by changing the time scale, as discussed above. A Reynolds number (based on base diameter) may be constructed using the effective artificial viscosity, and will behave like d/δx; provided that the density and velocity in a particular region do not differ from run to run, this number will be approximately twice as large on run 0 as on run 1.
In comparing the steady flow of run 1 with the near-steady flow of run 0, the effect of the transient flow patterns will probably be insignificant, and any differences are likely to be due only to the change in the ratio d/δx. However, further consideration of such effects was found to be unnecessary, since averaged property profiles at the end of run 1 were found to be almost identical with those of run 0. Two of the profiles, static pressure and temperature, are shown in Fig. 4-21, for x/d = 1.9. The scatter of cell values is about the same on both runs, and noticeable differences occur only where small concentrations of the original have become empty, as on run 1, if calculation had been continued.
On the basis of this evidence, it was concluded that for this type of flow pattern, doubling the ratio of mesh size to reference length has little effect on the steady flow. The tendency for local instabilities to occur, in the form of excessive tilde velocities, does seem to depend on the mesh size, but this is thought to be largely a function of the actual size of empty regions of the mesh, in terms of the number of cells involved. Since better results in this respect were obtained with the coarser mesh, the effect is probably peculiar to this particular flow field.
This was the first attempt to incorporate a mock boundary layer, as described in section 3, the mesh being identical with that of run 1. The same number of particles per cell in the input flow was also used, because although a higher particle number density was known to be desirable, the effects of the mock boundary layer could have been obscured by such an extra modification.
Fig. 4-22 to Fig. 4-26 show particle configuration plots at various stages in the calculation. The first of these, at t = 30.8, appears to show spatial oscillations in the boundary layer flow after separation; this cannot represent the beginning of a vortex-street wake, since the flow is forced to have symmetry about the axis, but is likely to be caused by the rapid variations in cell pressures in the base region at this time. A cell may be quickly emptied of the initial fluid, and then reoccupied either by initial fluid from other cells, or by boundary layer fluid which turns much more sharply at the corner than the flow in previous runs. Although the boundary layer is supersonic, subsonic velocities are calculated from time to time in cells just behind the base, adjacent to the freshly separated layer. These may also be responsible for the fluctuations.
By time 102.2, the oscillations are less marked. Particles began to leave the mesh at t = 62, and the maximum number in the system, 4260, was reached at t = 73, as for run 1. A wake shock forms downstream of X = 40, but the behaviour of the separated boundary layer is not clear from the plot. The next plot, at t = 176.0, shows that the innermost particles of the input flow are now turning sharply at the corner, but are still meeting, near the axis, particles which formed part of the initial, stationary fluid behind the base. The scavenging effect of the input flow on the original fluid is far smaller than on run 1; by a similar time on that run, the higher velocities of the innermost streamlines of the input flow had resulted in an almost complete clearance from the mesh of the original particles.
A similar situation is seen at t = 232.4, but none of the original fluid remains upstream of X = 30. Expansion of the flow round the base corner results in a large divergence of streamlines, and the input flow contains too few particles per cell to cope with this situation; the void seen in the region X = 10, Y ≈ 4 is caused by the low particle number density, but could perhaps have been eliminated by a more random pattern in the input flow.
The wake shock has become weaker since time 102.2, and has also moved closer to the axis. Applications of the limit on tilde velocities, fairly frequent early in the calculation, now occur only rarely, and δt maintains its maximum value most of the time. A steady flow configuration was considered to have been reached by t = 301.4, with the number of particles fairly steady at 3350, and the total energy of the system varying randomly by ±0.03%. Average computer time per cycle at this stage was 5.82 seconds, and the total time expended on the run was 2 hours 40 minutes.
All original particles have disappeared, causing the shock near the output boundary also to vanish. A wake shock in the proper sense, originating near the neck of the wake, cannot be calculated with such a low particle number density as found in this region. The position of the neck, at x/d = 1, contrasts strongly with that of run 1.
The large increase in wake angle could not at this stage be attributed confidently to the vorticity in the mock boundary layer, because of the behaviour of the fluid in the layer before separation. In the steady state, the velocity in the layer was equal to 0.55u1, in cell (1, 6), and 0.59 u1, in cel1 (2, 6), compared with the specified value just outside the input boundary of 0.50u1. This is a result of the velocity weighting procedure, whereby particles near the bottom of a boundary layer cell will be given effective velocities of about 0.5u1, but those near the top (outer edge of the layer) will acquire values near 0.75u1. More momentum will thus be transported to the next cell downstream than if all particles moved at a speed 0.5u1 so the velocity value of the next cell will be calculated as rather more than that upstream. This process is perpetuated, which accounts for the increasing velocity along the wall. The acceleration decreases the density, and hence the pressure, with the result that the pressure in cell (2, 6) was only 0.73p1.
The large wake angle could have been caused by this normal pressure gradient over the wall, and a further computation was made to determine the importance of this effect. When assigning masses to the input particles, the required fluid density will only be obtained if the particle number density is estimated correctly; since the frequency of input is fixed, this figure depends only on the particle velocities. Starting at time 280, a short run of 25 time units was made, using 0.5u1, as the specified boundary layer cell velocity, as before, but basing the masses assigned to fresh boundary layer particles on a nominal particle velocity of 0.6u1. By the end of this run, the velocity in cell (2, 6) had fallen from 0.59u1, to 0.53u1, and the pressure had risen from 0.73p1, to 0.96p1. The artificial normal pressure gradient over the wall was thus almost eliminated, and yet the normal velocity components in cells just downstream of the corner were almost exactly the same as before (slightly larger, in a few cases).
These results suggested that the differences in flow pattern between runs 1 and 2 are due merely to the layer of vorticity in the input flow.
Figs. 4-27 to 4-31 show property profiles in the steady state, averaged over ten time units. Axial values of u are seen to increase with downstream distance, reaching almost 0.9u1, at x/d = 5, whereas on run 1, no variation from u1 existed for non-empty cells downstream of the neck. In the transient stage of the solution, it was again noted that the presence of large quantities of original particles near the axis resulted in smoother profiles than those in the steady state, as can be seen in Fig. 4-32, which compares average velocity profiles at times t = 232 - 242 and t = 302 - 312. The relatively low cell velocities in row j = 4 (y/r = 0.7) are a consequence of the separation between streams of particles, resulting in low particle number densities.
Profiles of density, ρ/ρ1, differ from those of run 1 mainly in that the large peaks close to the axis have been eliminated as a result of the much lower density of the input fluid in the mock boundary layer. Coalescence of particle streams causes the peak at y/r = 0.7 on the x/d = 5 profile.
Temperature profiles show peaks close to, or on, the axis except at x/d = 1.9, and the magnitude of the peak temperature falls steadily with increasing x/d. The vastly increased values of T/T1, on this run are again due to the properties of the mock boundary layer, which has the same total temperature as the flow outside.
Near the output boundary, axial static pressure is of the same order as p1 compared with about 3p1, on run 1, but little significance can be attached to the rather random axial distribution, again due to the scarcity of particles.
Mach numbers in the flow now have more realistic values, axial values increasing with x/d. In the region of the neck, M is equal to about 1.5 to 2.0.
From the results of run 2, it was apparent that a layer of vorticity, one tenth of the base diameter in thickness, has a large effect on the calculated flow, producing more realistic wake geometry and fluid properties. As the next stage of this investigation, it was decided to repeat the calculation with a much larger number of particles per cell in the input flow, in order to avoid the gaps caused by separating streams of particles.
For this run, a particle number density of 16 per cell was chosen for the input flow, in place of 6 on all previous runs. In addition, the prescribed pattern for these particles was adjusted in order to eliminate as far as possible the periodicity in the y direction; this reduces further the tendency for voids to appear between streams of particles, though it makes streaklines less easy to detect by looking along the particle configuration plots at a small angle to the paper. The assignment of masses to fresh boundary layer particles was based on an average effective velocity for particles equal to the mean of the prescribed boundary layer cell velocity and the interpolated cell velocity at the outer edge of the layer. This alteration should, in theory, prevent the increase in velocity and reduction of pressure in the cells of the layer which occurred on run 2.
Arrays for particle quantities were increased to 9000 elements, but the overall size of the flow field was reduced to 40 cells by 12, eliminating a strip of redundant uniform flow along the top of the mesh, and all of the region x/d > 3.8. The reason for the long mesh introduced on run 1 no longer exists.
Particle configuration plots at various stages of the run are shown in Figs. 4-33 to 4-37. On that at time t = 28.0, the pattern between the shock front and the line dividing original fluid from the input flow is seen to differ from that formed on previous runs. This is due to a new arrangement of the original particles, which also eliminates the cell to cell pressure differences in the stationary fluid. The density of the initial fluid on this run was ρ1/8, compared with ρ1/6 on previous calculations; since the temperatures are the same as before, the pressure difference between the initial and the input fluid was more favourable.
At t = 30.8 on run two shows the effect of altering the pressure difference. The average speed of the shock front, up to the respective time was 0.91u1 on run 2, and 0.95u1 on run 3, while the average speed of the first input column was 0.75u1 for run 2, 0.80 for run 3.
At t = 52, the number of particles reached a maximum of 6036, and the computer time for one cycle rose to 9.5 seconds. It is interesting to note that, since time is scaled such that one problem time unit corresponds to 1.98 microseconds in the actual flow being simulated, the computer has slowed down the rate of flow development by a factor of 24 million (for δt = 0.2).
Applications of the tilde velocity limit occurred far less frequently than on run 2, and after t = 160 no further limiting was necessary.
Comparison of the particle configuration plot at t = 185.2 with that at t = 176.0 on run 2 shows that particles of the initial fluid near the base are causing the innermost streams of boundary layer particles to be turned parallel to the axis at a radius y/r ∼ 0.6, whereas on run 2, the empty base region allowed these streams to continue as far as the axis. The wake angle seems to be the same as on the previous run, but the slower development of a steady state gave some indication at this stage that it might take considerably more than about 300 time units to complete the calculation.
Acceleration of fluid in the boundary layer was considerably reduced by the new scheme for assignment of particle masses. Cell velocities averaged 0.507u1, and 0.530u1, in cells (1, 6) and (2, 6) respectively, while the pressure in cell (2, 6) averaged 99.5% of p1; the resulting normal velocity component was only 0.004u1. Thus although the specified boundary layer velocity has not quite been achieved, the effects of the spurious normal pressure gradient are negligible.
After 300 time units had elapsed, equal to the total time used for run 2, a set of averaged profiles exhibited axial velocities much lower than the steady values of the previous run. The corresponding static temperatures were much greater, suggesting a state rather like that at about t = 200 on run 2; the continuing decrease in the number of particles in the system also suggested that much more computation would be required before a steady flow pattern emerged.
During the next 130 time units, a further reduction of 100 in the total number of particles took place, while velocities and temperatures near the axis altered significantly, moving closer to the steady state values of run 2. However, the final 30 minutes of computing time, which took the solution to t= 477, made virtually no difference to the flow properties, and the number of particles in the mesh remained at 4650. Over the last 20 time units, the total energy of the system varied randomly by less than ± 0.009%, about one quarter of the size of the fluctuations at the end of run 2. Total computer time for the run was 5 hours.
From the particle pattern at t = 300.6, it can be seen that the flow close to the base resembles that on run 2 at the same time, but the large voids in the free shear layer have been eliminated. Downstream of the neck, a line of high particle concentration indicates a wake shock, but as the region near the axis empties, this line breaks up, as the plot for time 382.4 shows. From this time onwards, little change occurs in the flow geometry, and the final plot at t = 477.2 corresponds closely with the steady flow of run 2; the concentrated particle lines match similar but less clearly defined lines on that run.
The wake neck is situated at x/d ≈ 0.9, and there appears to be a wake shock originating near this point. However, in view of its subsequent behaviour, it is doubtful whether this line may be interpreted as a wake shock; the particle number density is still too low, and the mesh too coarse to allow a shock of this strength to manifest itself in the property profiles.
Profiles of flow properties are shown in Figs. 4-38 to 4-42, averaged over the period t = 477 to 487; direct comparison with run 2 is possible at x/d = 1.9 and 3.3. The velocity profiles are much smoother than before, and values of axial velocity on a single cycle show also a slight reduction from the equivalent figures for run 2, e.g. on the axis in the region of x/d = 3.5, u/u1 was about 0.75 on run 2, but only 0.70 on run 3.
Profiles of density are very similar to those of the previous run, but the smoothing resulting from the increased number of particles is not as marked as one might have expected. The worst scatter occurred at x/d = 5 on the previous run, where unfortunately no comparison can be made. A small peak is seen on the axis at x/d = 1.0, where the neck is situated.
Static temperature profiles show a peak on the axis at x/d = 1.0, but further downstream the peak moves to y/r = 0.5, as before. The magnitude of the peak temperature at the neck is greater than on run 2, and now values of T/T1 ≈ 10 persist as far as x/d = 3.3; on run 2, for the axial value of T/T1, energy appears to be directly related to the difference in axial velocities between the two runs; in cell (35, 1) (next to the axis at x/d = 3.3), the total energy, effectively equal to I + ½u2, is virtually the same for both runs. This change in the distribution of total energy may be the result of shock heating at the lines of high particle concentration downstream of the neck.
Values of static pressure show as much scatter as on the previous run, but are similar in magnitude. Close to the axis, of course, the higher temperatures cause the pressures also to be greater than before. At x/d = 1.0, p/p1, reaches approximately 2 on the axis, as was the case on run 2, but beyond the neck, axial values less than 1 occur until x/d = 3.
Mach number profiles at x/d = 1.9 and 3.3 are smoother than the corresponding ones on run 2, but otherwise are very similar. On the axis, values are slightly lower than before owing to the changes in temperature and velocity.
The region of low speed recirculating flow found in all hypersonic wakes, bounded by the axis of symmetry, the base wall, and the dividing streamline, has not been included in the steady flow solutions. It seems that the base pressure for this particular flow is too low for representation even by one of the lightest input particles per cell, and the mixing of the input flow with the low density initial fluid empties the base region before the wake neck has properly formed. In the Mach 1.89 wake flow calculation made by Amsden and Harlow (5), the neck formed at a much earlier stage in the transient solution, resulting in a balance of fluid in particles and reverse flow from the rear stagnation point. The particle number density was sufficient to represent the low pressure recirculating flow, assisted by the fact that all particles had the same mass, since that calculation was in plane geometry.
However, evidence that a recirculating flow would exist in the present solution, if the conditions mentioned above could be satisfied, is provided by velocity profiles and particle tracks at early stages in the solution. An extensive region of reversed flow exists at x/d = 0.25 at time 69 on run 3, as shown by the velocity profile in Fig. 4-43, whereas by time 290 the three innermost cells at this location are empty on most cycles. Fig. 4-44 shows the paths followed by selected particles originating near the axis at x/d 0.3, on runs 1 and 3. On run 3, the chosen particle took much longer to traverse a similar though less tortuous path to that on run 1, an indication that the flow on run 3 might take considerably longer, in problem time units as well as computer time, to attain a steady state.
Exact comparison of the results obtained here with other theoretical and experimental results is not possible, since experiments have generally been performed with cones or blunt bodies, and the published applications of theory refer to wakes behind such bodies. However, some general comments may be made on both the wake geometry and the variation of flow properties.
The distinct difference in flow pattern resulting from the inclusion of a layer of vorticity has already been mentioned. A wake closure angle of 7° was obtained on run 1, and 29° on runs 2 and 3, defining the angle by the free stream direction and the line joining the base corner to the point where the innermost stream of particles meets the axis of symmetry. An idea of the streakline pattern close to the base may be obtained from Fig. 4-45, which shows the paths followed by three particles which were tracked towards the end of run 3.
The extensions to the Chapman theory by Denison and Baum, and by Baum (ref. 10 and ref. 11) predict very small wake closure angles. For a 5° semi-angle cone at a Mach number of about 10, these solutions give a wake angle of 2° if a Blasius initial profile is assumed, or about 4° allowing for distortion of the profile by the corner expansion; reducing the cone apex angle also reduces the wake angle. The inviscid boundary layer expansion solution of Weiss and Weinbaum (9) has at present been used only for planar flows. However, for the particular case of a 10° semi-angle wedge under typical shock tunnel conditions, their solution places the rear stagnation point about one diameter from the base, where the models of ref. 10 and ref. 11 would predict the distance to be more than ten times this value.
In ref. 12, Waldbusser presents an analysis of ballistic range shadowgraph pictures of pointed and blunt slender cones, over a wide range of freestream Mach and Reynolds numbers. This data indicates that the distance of the wake neck from the base of the cone is a roughly linear function of Mach number, with no large effects due to cone angle, nose bluntness, or Reynolds number (over a range 8 < M∞ < 16, 6.5×104<Re∞ = 8.6, the neck position is indicated as x/d=1.25 or 1.50 for a value of M∞ which would give an inviscid surface Mach number of 8.6 on a 5° semi-angle cone. The neck, as defined for Waldbusser's study, is the point where the wake shock intersects the edge of the viscous core, but the shock extends further upstream than this, and its origin indicates roughly the position of the rear stagnation point. The shock origin position from this correlation agrees well with the estimate of the stagnation point behind a 10° semi-angle cone at M∞ = 12.8 made by Muntz and Softley (13).
Corresponding to neck locations x/d = 1.25 and 1.50, Waldbusser's correlation gives the shock origin positions as x/d = 0.65 and 0.75 respectively. The definition of wake neck used in describing the steady-state P-I-C configurations probably corresponds more closely with the rear stagnation point than a true wake neck, since it refers to the intersection of the inner edge of the shear layer with the axis, and not the outer edge with the recompression shock. On this basis, then, the P-I-C model indicates a stagnation point location, x/d = 0.87, in fairly close agreement with experimental results for slender cones at equivalent Mach numbers (referring to Jabberwock run 3). On a single cycle at the end of run 3, a local static temperature peak on the axis occurred in cell (11, 1), i.e. x/d = 0.85, further evidence that the neck of the particle plot represents roughly the rear stagnation point.
On the schlieren wake photographs taken as part of the experimental investigation of ref. 14, the wake shock behind a 10° semi-angle cone can be traced upstream to about x/d = 1.0. Because the appearance of the photographs depends as much on the optical system as on the strength of the compression waves, one can say that the shock origin, and also the stagnation point, lie within one diameter of the base, again agreeing approximately with the Jabberwock result.
A further analysis of ballistic range results for laminar cone wakes is reported by Levensteins and Krumins (15). Data relating to a 6.3° semi-angle cone at M∞ = 9 shows a slight Reynolds number effect, the distances of both shock origin and wake neck from the base increasing as Re∞,d = 105. A variation of this kind has little relevance to the P-I-C calculations, because of the impossibility of assigning a Reynolds number, especially when a mock boundary layer is incorporated.
It seems, therefore, that the particle-in-cell model, with a crude inviscid boundary layer representation, will predict the shape of the base flow region with reasonable accuracy, if it can be assumed that experimental results from slender cones would not alter significantly if the cone semi-angle were reduced to zero, and all nose effects were absent. The model compares, in this respect, much more closely with an inviscid expansion model than with the earlier models which were essentially extensions of Chapman's original attempt at a free shear layer solution (ref. 16). This is perhaps to be expected, since the particle-in-cell solution contains none of the simplifying assumptions of the extended Chapman solutions concerning the shape and angle of the dividing streamline, and does essentially perform an inviscid expansion of the input flow.
Particularly noticeable on the final plot for run 3 is the cavity next to the axis, downstream of the neck. This is probably a consequence of the boundary condition applied at the axis during the particle movement stage of each cycle of the calculation, whereby a fictitious cell just beyond the axis is given the same velocity as the adjacent interior cell; particles are reflected back into the mesh, with a momentum adjustment as though they had entered from a cell with reflected velocity. Amsden, in ref. 7, mentions that for cases where the fluid motion is primarily along the boundary, the midline is likely to open up and produce a cavity when this type of boundary condition is used, but that the alternative procedure, of giving the fictitious outside cells reflected velocities during the weighting procedure, is undesirable when fluid is moving rapidly towards the boundary. Now in the free shear layer, a rapid movement of fluid towards the lower boundary of the mesh does occur, so the former method of treating this boundary was chosen, in order to avoid a so-called boundary catastrophe in the region of the wake neck. Perhaps the latter method would have resulted in a more realistic flow close to the axis, but the calculation could easily have become unstable while the wake neck was forming.
It is perhaps worth mentioning an apparent anomaly in the relative sizes of the wake closure angles on runs 1 and 2. Comparing the free shear layer solution of Denison and Baum with that of Chapman, the inclusion of a boundary layer of finite thickness is seen to decrease the wake angle from that of the zero thickness case, the opposite effect to that calculated here. The smaller theoretical wake angle is, however, a result of an increased base pressure, whereas in the P-I-C calculations the base region has zero pressure on both runs 1 and 2. Just after the corner, the normal pressure gradient will be similar for both runs, so the flow which has the smaller x-component of velocity close to the wall will attain a larger wake angle.
The artificially induced normal velocity components near the wall which occurred on run 2, and were eliminated, virtually, on run 3, raise the question of the possible effect on the wake angle of artificial viscosity. Close to the wall on both these runs, ∂u/∂y is large, and the relevant artificial viscosity coefficient is ½ρ|v|δy using the average value of v in cell (2, 6), the coefficient has a value equivalent to about 11.0 × 107 slug/ft.sec. for run 2, and 1.6 x 107 slug/ft.sec. for run 3. (Sutherland's formula results in a value for μ of 6.7 × 107 slug/ft.sec., using the temperature in the boundary layer). Now the wake angles on the two runs are almost identical, so it may be concluded that the artificial viscosity is not, in this problem, of much importance in fixing the wake angle. In any case, since the mock boundary layer occupies only one row of cells, the artificial viscosity cannot be acting in the same way as true viscosity in a real boundary layer.
Very little information has been published on the distributions of fluid properties which result from various theoretical treatments of the near wake. This state of affairs is to be expected, since no satisfactory overall solution of slender body near wakes was available before that of Weiss (3) . Without performing a solution of this type for the Jabberwock flow configuration, the calculated properties in the steady state must be judged by a comparison with experimental results for slender cone wakes. Because of the large changes which result from altering the freestream, properties (non-dimensionalised by the freestream values) will not be meaningful. For example, increasing Re∞,d from 105 to 106 typically halves wake centreline values of ρ/ρ∞ . Also, it is doubtful whether the uniform input flow of the P-I-C model would correspond more with the free stream than with the flow behind the bow shock in an experiment on a cone. However, general trends in the Jabberwock results may be compared with experiments on laminar wakes.
Zakkay and Cresci (17), experimenting with a sharp 5° semi-angle cone at M∞ = 8.0 and 11.8, found the centreline value of density to increase monotonically with x/d. The results of Muntz and Softley (13), for a slightly blunted 10° semi-angle cone at M∞ = 12.8, show a peak near the rear stagnation point (x/d ≈= 0.8), after which the density levels out at x/d ≈ 2.5. These measurements were made using an electron beam technique in a shock tunnel having a conical nozzle, whereas Zakkay and Cresci used conventional methods in a blowdown tunnel. A peak in ρ/ρ1, is present in the P-I-C results of run 3, as seen on the profile at x/d = 1.0, but further downstream, a gradual rise takes place in the centre-line value, which may be seen more easily by imagining smooth curves drawn on the profiles for run 2. Axial values of ρ/ρ1 are in fact of the same order as experimental values of ρ/ρ∞ typically 0.2 - 0.3 at x/d ≈ 4.
Static temperatures in the near wake are sensitive functions of the body wall temperature, particularly the peak value which occurs at the rear stagnation point. Experimental data indicates, however, that for downstream distances greater than about x/d = 2, the effect of wall temperature on the local stagnation temperature is unimportant. Zakkay and Cresci (17) present data on centreline variation of static temperature from several different investigations where wall temperature ratios were similar; downstream of the stagnation point, little variation in T/T∞ occurs until x/d = 3 or 4, when the values begin to fall off from about 0.3. The results of Schmidt (18), for a sharp 10° semi-angle cone at M=∞ = 8.0, confirm this behaviour. On run 2, the P-I-C results for x/d > 3 show a similar pattern, though run 3 is somewhat inconsistent in the calculated value of T/T1 on the axis at x/d = 3.3, as mentioned previously.
A more interesting feature of the results of both runs 2 and 3 is the occurrence of off-axis maxima at x/d = 1.9. An examination of available experimental data for sharp cones shows that the maximum in the radial direction occurs on the axis, for all x/d, and that T/T0∞ typically falls to about half its axial value at y/r = 0.5. It is difficult to detect the wake shock on experimental temperature profiles; however, the peak in the P-I-C profiles of T/T1, at x/d = 1.9, coincides in radial position with the particle concentration that represents the wake shock.
An interesting comparison is possible here with the theoretical solution of Weiss (3). The temperature distributions in the outer inviscid flow of that solution exhibited off-axis maxima in the expansion fan, but the magnitudes of these peak temperatures were smaller than the on-axis maxima produced by viscous effects in the shear layer and recirculation region. Weiss did not carry his outer flow solution beyond the position of the rear stagnation point, but it is reasonable to suppose that the off-axis maximum temperatures would have persisted downstream.
The peak in the P-I-C profiles is probably due both to the effect just mentioned and to the axial cavity, where lack of particles could result in a spurious on-axis minimum temperature. At x/d = 3.3, the cavity is narrower than at x/d = 1.9, and the cell next to the axis contains much more fluid; here the axial minimum temperature is not apparent. On run 2, the cavity had closed by x/d = 3.3, and at this position the maximum temperature does occur on the axis, in contrast with the profile at x/d = 4.0, where, once more, the cell next to the axis is almost empty. An off-axis maximum also occurs at x/d = 1.0 on run 3, though this peak is closer to the axis than those at x/d = 1.9 and 2.5. Here, the effect cannot be explained by the axial cavity, which is therefore likely to be responsible only in part for the similar effect further downstream.
Anomalies in the P-I-C density and temperature profiles cause the static pressure profiles to have a rather random appearance for y/r < 1, although in the expansion fan outside this region, the pressure variation on run 3 is smooth and realistic. The axial variation deduced in ref. 13, from electron beam measurements, has a peak near the stagnation point, but again this is not in agreement with the data of Zakkay and Cresci, whose values of p/p ∞ rise typically from about 0.5 at x/d = 1 to 2.5 at x/d = 4. Schmidt's measurements indicate that p/p∞ does not begin to fall again until x/d reaches 5 or 6, while Todisco and Pallone (19) come to similar conclusions for the wake of a 5° semi-angle cone at M∞ = 12.5. Interpretation of Muntz and Softley's data is again complicated by the falling freestream pressure produced by their conical nozzle.
The P-I-C pressure profiles of run 2 show, on average, no rise in axial value from x/d = 3.3 onwards, though from x/d = 1.9 up to this point, a considerable rise occurs on the run 3 profiles. However, the cavity probably has a large effect on these centreline values. Experimental radial profiles show a fairly constant pressure inside the wake shock, with a sharp drop on moving outwards through the shock, even at its origin near the stagnation point. The shock-like particle concentration of run 3 has its origin at x/d ≈ 1.0, y/r ≈ 0.4, and could be responsible for the sharp fall in p/p1, with y/r which occurs on the x/d = 1.0 profile, though the centreline pressure seems excessively high, in common with the density at this point. The type of axial variation found by Muntz and Softley, however, does include a peak pressure at the stagnation point (x/d = 0.8) equal to twice that at x/d = 2.0. If the axial boundary treatment during particle movement had been carried out by the alternative method, the pressure and density in this region might have been lower, since the flow in the inner edge of the shear layer would probably have been turned back parallel to the axis rather less suddenly.
The velocity profiles of run 3 show a rather small increase in the downstream direction, which is in general agreement with values calculated from Zakkay and Cresci's Mach number and static temperature data, although the actual magnitudes are different, e.g. at x/d = 3.3, u/u∞= 0.46 compared with the P-I-C value u/u1, = 0.65. Todisco and Pallone present a single velocity profile for x/d = 2, showing a roughly constant u/u∞ (equal to 0.3) for y/r < 0.5, whereas the P-I-C velocities tend to increase more rapidly with y/r. A comparison is possible here with the overall near wake solution by Weiss (3), who calculates the velocity profile at the rear stagnation point in the laminar wake of a 10° semi-angle wedge at M∞ = 16. The velocity rises almost linearly with distance from the axis until y = 0.6 times the base half-height, and is thereafter constant (approximately), equal to the value at the edge of the shear layer. At x/d = 1.0, the profile of run 3 is of very similar appearance, though it is somewhat downstream of the stagnation point .
Axial Mach numbers deduced from pitot and cone static probe measurements for a 10° semi-angle cone at M∞ = 11.8 (ref.17) are very similar to those behind a 5° semi-angle cone at M∞ = 12.0 (ref. 19). Values rise from about 1.3 at x/d = 1.9 to around 1.8 at x/d = 3.3, reaching 2.3 at about 5 diameters downstream. AtM∞ = 8.0, the results of ref. 4 show a rather faster rise, from 1.8 at x/d = 1.9 to 2.8 at x/d = 5.0. From the P-I-C profiles of run 3, Mach numbers near the centreline are all approximately equal to 2; the small variation found experimentally cannot be expected to appear in the P-I-C results when the nature of the flow near the axis is taken into account. The expected trend to higher values of M with increasing x/d does show up when the wake is calculated as far as x/d = 5, on run 2, though the value of about 4 reached here at the downstream end of the flow field is slightly higher than experiments indicate for cone wakes. Except at x/d = 3.3, the profiles of run 3 show a region of approximately constant Mach number from the axis to y/r = 0.5, exactly as on the radial profile at x/d = 2.0 presented by Todisco and Pallone (19). This is probably not significant, since the constant experimental value arises from constant values of both u and T, whereas both these properties are increasing with y/r in the P-I-C solution.
The P-I-C method has been used to compute the near wake now behind a cylindricaJ. afterbody in a Mach 8.6 air stream. Where the input flow over the base of the body was completely uniform, the steady flow pattern did not correspond to that expected for such geometry and flow conditions. The neck of the wake was located at around four times the distance fran the body which one would expect by comparison with experiments on slender cones, coarse representation of a boundary layer was then included by introducing an appropriate amount of vorticity next to the all upstream of the corner, the result of this modification was to move the neck to a realistic position, supporting the recent theories that vorticity in a hypersonic boundary layer plays a more important part than viscous effects in determining the rear stagnation point location. Too much should not be made of this conclusion, however, since the computing mesh was too coarse for a proper simulation of reflected expansion wave phenomena.
Flow in the recirculation region was absent in the steady solutions because the density in that region was too low to be represented by the particle model of the fluid, even then the particle number density was increased on the fourth calculation. Properties of the flow downstream of the neck agreed reasonably well with experimental values for slender body wakes in sane cases, but the lack of real viscosity was evident, especially in the temperature profiles. The recompression shock wave was too weak for adequate representation by the particle number density in use.
Alteration by a factor of two of the ratio of mesh size to body diameter produced no noticeable effect on the solutions when the mock boundary layer was absent.
Stable near wake solutions by the P-I-C method have been shown to be possible, but in the form used here, the steady state profiles would not be of any use as initial conditions for a far wake calculation. The overall geometry of the wake was simulated with reasonable accuracy, but to obtain sufficiently accurate values of flow properties it would be necessary to include real viscous terms in the equations, and probably subtraction of the artificial viscous terms in some regions of the mesh. A correct description of the free shear layer and the recirculating flow, and their effect on the flow downstream, would require a much finer representation of the boundary layer and also a much greater density of particles in the input flow. The consequent increase in the required computing power would be enormous. For the flow over a cone rather than an axially-aligned cylinder, triangular cells would be necessary upstream of the base, introducing further complication, unless the input boundary were at the base position, in the latter case, conical flow could be simulated by giving the input flow a radial velocity component.
The author wishes to thank Dr. D. A. Spence for supervising this work, and also the staff of the Atlas Computer Laboratory for their assistance.
cv specific heat at constant volume
d body base diameter
E a total energy of cell [= M(I+ ½(U2+v2))]
I specific internal energy (=cvT)
i cell index in x-direction
j cell index in y-direction
M cell mass; Mach number
m particle mass
p static pressure
R gas constant
Re Reynolds number
r body base radius
T static temperature
t scaled time
u velocity component in x-direction
v velocity component in y-direction
X cell momentum in x-direction (section 2 only);
mesh coordinate in x-direction measured from input boundary
x coordinate in freestream direction measure from base of body
Y cell momentum in y-direction (section 2 only);
mesh coordinate in y-direction measured from lower boundary
y coordinate in radial direction measure from axis of symmetry
γ ratio of specific heats
δ prefix for P-I-C finite difference increment
μ dynamic viscosity coefficient
ρ density
τ cell volume
Subscripts
d based on body base diameter
i,j value in cell (i,j)
1 value in uniform input flow
∞ value in free stream
Superscript
~ tentative value calculated in phase I of P-I-C cycle
1. Lees, L, Hypersonic wakes and trails, AIAA J. 2, 417 (Mar. 1964).
2. Lykoudis, P, A review of hypersonic wake studies, Rand Corpn Memo. RM-4493-ARPA (May 1965); also AIAA J. 4, 557, (Apr 1966).
3. Weiss, R F, A new theoretical solution of the laminar hypersonic near wake, AIAA J. 5, 2142 (Dec 1967).
4. Baum, E, King, H H, Denison, M R , Recent studies of the laminar base-flow region, AIAA 2, 1527 (Sep 1964).
5. Amsden, A A, Harlow, F H, Numerical calculation of supersonic wake flow, AIAA J. 3, 2081 (Nov 1965).
6. Harlow, F H, The particle-in-cell computing method for fluid dynamics, in Methods in Computational Physics, Ed. Alder, B. et al. (Academic Press, New York, 1964), Vol. 3, p.319.
7. Amsden, A A, The particle-in-cell method for the calculation of the dynamics of compressible fluids, Los Alamos Sci. Lab. Rep. LA-3466 (1966).
8. Weinbaum, S, On the singular points in the laminar two-dimensional near wake flow field, J. Fluid Mech. 33(1), 39 (Jul 1968).
9. Weiss, R F, Weinbaum, S, Hypersonic boundary layer separation and the base flow problem, AIAA J. 1, 1321 (Aug 1961).
10. Denison, M R, Baum, E, Compressible free shear layer with finite initial thickness, AIAA J. l, 342 (Feb. 1963).
11. Baum, E, Effect of boundary layer distortion at separation on the laminar base flow , ElectroOptical Systems Inc. RN-l6 (Oct 1963).
12. Waldbusser, E, Geometry of the near wake of pointed and blunt hypersonic cones, AIAA J. 4, 1874 (Oct 1966).
13. Muntz, E P, Softley, E J, A study of laminar near wakes, AIAA J. A, 961 (Jun 1966).
14. Crane, R I, Studies of near wakes in hypersonic flow, D. Phil. Thesis, Oxford University (Aug 1968),
15. Levensteins, Z J, Krumins, M V, Aerodynamic characteristics of hypersonic wakes, AIAA J. 5, 1596 (Sep 1967).
16. Chapman, D R, Laminar mixing of a compressible fluid", NACA Rep. 958 (1950).
17. Zakkay, V, Cresci, R J, An experimental investigation of the near wake of a slender cone at M = 8 and 12, AIAA J. 4, 41 (Jan 1966).
18. Schmidt, E M, An investigation of hypersonic flow around a slender cone, AIAA Student J 3, 129 (Dec 1965).
19. Todisco, A, Pallone, A J, Measurements in laminar and turbulent near wakes, AIAA Paper 67-30 (Jan 1967).