Bagher Abareshi

Hassan Mohagheghi

Khosronia Kalalabadi

Investigate auto-ignition of stoichiometric methane-air mixture

Keywords:

Auto-ignition process of stoichiometric mixture of methane-air is investigated using detailed chemical kinetics in a single-zone combustion chamber. Effect of initial temperature on start of combustion (SOC). The Arrhenius expression for the specific reaction rate are calculated and auto-ignition was evaluated based on the species fractions and sensitivity analysis. Our results suggest that the SOC is directly related to initial temperature and the auto-ignition will not occur if the initial temperature low enough.

Introduction

Study of auto-ignition process in premixed air/fuel systems will help us understand the effect of different parameters on SOC. This will be very helpful for identifying and control of SOC in premixed fueled engines like Homogenous Charge Compression Ignition (HCCI) combustion engines [START_REF] Nobakht | A parametric study on natural gas fueled hcci combustion engine using a multi-zone combustion model[END_REF][START_REF] Andwari | Experimental investigation of the influence of internal and external egr on the combustion characteristics of a controlled auto-ignition two-stroke cycle engine[END_REF][START_REF] Tahsini | Numerical study of solid fuel evaporation and auto-ignition in a dump combustor[END_REF][START_REF] Maurya | Numerical investigation of ethanol fuelled hcci engine using stochastic reactor model. part 2: Parametric study of performance and emissions characteristics using new reduced ethanol oxidation mechanism[END_REF][START_REF] Nobakht | Optimization of gri-mech 3.0 mechanism using hcci combustion models and genetic algorithm[END_REF][START_REF] Mahabadipour | lost available imep: A second-law-based performance parameter for ic engines[END_REF]. Control of SOC is a very challenging task to commercialize the HCCI combustion engines. [START_REF] Razzaghi | Feasibility study of auto ignition in in-situ combustion process[END_REF][START_REF] Mahabadipour | Development and comparison of two expander cycles used in refrigeration system of olefin plant based on exergy analysis[END_REF].

In this study, we examined the effect of initial temperature on SOC and studied the mole fraction behavior of essential radicals for auto-ignition initialization using detailed chemical kinetics.

Methodology

To GRI-mech 3.0 chemical kinetic mechanism [START_REF] Gardiner | The gri-mechtm model for natural gas combustion and no formation and removal chemistry[END_REF] was used to model detailed chemical kinetics.

Thermodynamic properties are uniform and it is assumed that the gas is an ideal gas. We also assumed that the mass inside the cylinder is homogenous at the beginning of the simulation:

dm k dt = 0 (1) Nz k=1 m k = m tot (2)
Energy conservation equation is as follows:

m k ck v dT k dt = -m k Ns i=1 U i dY i,k dt -P dV k t + dQ k dt + ṁi-1→i h i-1 - ṁi→i+1 h i (3)
The net rate of species production can be calculated using:

dY i,k dt = ωi,k M W i ρ k (4)
Woschnis correlation [START_REF] Chang | New heat transfer correlation for an hcci engine derived from measurements of instantaneous surface heat flux[END_REF] is used to simulate convective thermal transport inside the cylinder to transfer heat to cylinder wall.

h c (t) = β × Height(t) -0.2 P (t) 0.8 T (t) -0.73 ν(t) 0.8 (5)
Thermal transport modeled by a mechanism which is like conduction. So, heat flux of thermal transport will be dependent on the temperature difference of the mixture and the wall which are neighbor and on the distance as follows:

q = -K tot dT dy (6)
In Eq. 6, the K tot is a conductivity of total flow which is found by Yang and Martin method [START_REF] Yang | Approximate solutionone-dimensional energy equation for transient, compressible, low mach number turbulent boundary layer flows[END_REF][START_REF] Nobakht | Cfd simulation of tesla-type microvalves[END_REF] and is the summation the turbulent and laminar component.

k tot = K l + K t (7)
The ratio of turbulent to laminar component is calculated by the following equation:

K t K l = P r l P r t µ t µ l (8)
Eq. 8 presupposes flows. The viscosity ratio can by calculated by:

µ t µ l = ky + n [1 -exp(-2αky + n ] (9)
and

y + n = u * µ w u 0 ρdy n (10)
here k = 0.44 which is the Von Karman Constant, α = 0.08 and y + n is the distance to the wall. The velocity is proportional to engine rpm.

Results and Discussion

The effect of initial temperature on auto-ignition time is investigated in Fig 1 Increasing the initial temperature will increase the rate of essential radials production in the mixture and will help the SOC to happen faster. We can see as the initial temperature increases the SOC happens sooner and as the initial temperature drops below 1200K the auto ignition happens very late. Fig. 2 shows the zoomed graph for the earlier times to see the effect of temperature for higher temperatures in more details. We have calculated the auto-ignition delay time using this simulations. Considering 300 K delta T for auto-ignition time, the auto-ignition delay times are:

For sensitivity analysis when the mixture temperature has increased by 100 K, 200 K, and 300K, first corresponding time for 1600K, 1700K and 1800K are obtained:

Figure 1 :

 1 Figure 1: Effect of initial temperature on auto-ignition time.

Figure 2 :

 2 Figure 2: Effect of initial temperature on auto-ignition time (zoomed).

  that with increasing the initial temperature, 1/T term reduces and auto-ignition term decreased as well. It is due to exponent term of Arrhenius equation which has negative sign which results to increase of rate of reaction. Therefore, reactions happen faster which results in faster production of radicals and decreasing of auto ignition delay. Fig.3demonstrates this auto-ignition delay time.

Figure 3 :

 3 Figure 3: Auto-ignition delay time.

Figure 4 :

 4 Figure 4: Mole fractions of OH, H, O, and CH radicals.

Figure 5 :

 5 Figure 5: Sensitivity analysis of O.

Figure 6 :

 6 Figure 6: Sensitivity analysis of H.

Figure 7 :

 7 Figure 7: Sensitivity analysis of CH.

Figure 8 :

 8 Figure 8: Sensitivity analysis of OH.

Figure 9 :

 9 Figure 9: Sensitivity analysis of O.

Figure 10 :

 10 Figure 10: Sensitivity analysis of H.

Figure 11 :

 11 Figure 11: Sensitivity analysis of CH.

Figure 12 :

 12 Figure 12: Sensitivity analysis of OH.

Figure 13 :

 13 Figure 13: Sensitivity analysis of O.

Figure 14 :

 14 Figure 14: Sensitivity analysis of H.

Table 1 :

 1 Auto-ignition delay times

Table 2 :

 2 Sensitivity analysis of the mixture

	Time(s) T (K)
	1.10E-03 1.61E+03
	1.10E-03 1.61E+03
	1.11E-03 1.62E+03
	1.12E-03 1.63E+03
	1.12E-03 1.65E+03
	1.13E-03 1.67E+03
	1.14E-03 1.70E+03
	1.15E-03 1.73E+03
	1.15E-03 1.75E+03
	1.15E-03 1.77E+03
	1.16E-03 1.80E+03
	At time 0.001099 s corresponding to 1600 K the results of sensitivity analysis

for O, H, CH and OH are as follow:

Conclusions

The effects of initial temperature on auto-ignition behavior of premixed methane-air fuel is investigated using detailed chemical kinetics method. Increasing temperature from 1000 K to 1500 K reduced the Auto-ignition delay significantly. The highest reduction rate is between 1000K and 1100K and between 1100K and 1200K. The Arrhenius analysis shows that the AD is increases


Makoto Mita

Manabu Ataka

Hiroyuki Fujita

Hiroshi Toshiyoshi

An Inertia Driven Micro-Actuator for Space Applications

Keywords:

Introduction

Primarily motors and piezo elements have been used conventionally as drive devices for space probe mechanisms. However, planets located close to the Sun, such as Mercury and Venus, have high temperature environments of several hundred degrees Celsius, while on the other hand, cryogenic environments of tens to hundreds of degrees below 0 • C may exist in deep space or on the dark side of astronomical bodies located far from the sun. In such high temperature or cryogenic environments, motors and piezo elements show drastic degradations in performance, or even fail in the worst case. Electrostatically driven actuators are suitable for such extreme environments. Such actuators use air-gaps as dielectrics, and their operating characteristics theoretically do not depend on temperature. This means that drive systems with superior environmental tolerance can be implemented. However, the output of electrostatically driven actuators is low compared to that of magneticallydriven or piezoelectrically driven units. In order to re-solve this problem, we developed an inertia-driven actuator featuring a high output [START_REF] Mita | A micromachined impact microactuator driven by electrostatic force[END_REF][START_REF] Mita | Highly mobile 2D micro impact actuator for space applications[END_REF]. However, in the previous research the output was not sufficiently high to implement the actuator on a printed board rather than on a single chip. In this study we report improvements that result in a greater driving force.

Driving Principle

The working principle of the proposed device is illustrated in Fig. 1. Here the driving voltage is applied to the drive electrode and the stopper is set to the same potential as the mass (grounded) so as to avoid short circuiting.

(a) First, the driving voltage is applied between the movable mass and the stopper (actually, the drive electrode). The mass is accelerated toward the stopper by an electrostatic attractive force. (b) When the mass collides with the stopper, the whole chip moves due to the impact force. (c) When the voltage is removed, the mass returns to its original position due to the restoring force of the spring (suspension). The chip moves forward as operating sequence (a) to (c) is repeated. The travel distance due to one collision is small, but as the collisions are repeated, the chip can move over a long distance similarly to an inchworm.

Operation Analysis and Simulation

Aiming at operation analysis for the purpose of design, we determine the collision velocity of the movable mass, which has a strong effect on the driving force generated by impact. A simplified model of the device used in analysis is shown in Fig. 2. In this model, a collision phenomenon is represented by a stiff virtual spring applying pressure to a wall [START_REF] Faik | Modeling of impact dynamics: A literature survey[END_REF]. Here k is the spring constant of the suspension supporting the mass, k w is the spring constant of the virtual spring, m is the mass, g is the initial gap, and V d is the driving voltage. In the initial state (time t = 0, displacement x = 0), the system energy W 0 includes only electrostatic energy, that is,

W 0 = 1 2 C 0 V 2 d . ( 1
)
Here C 0 is the electrostatic capacitance of the parallel plate actuator in the initial state. Immediately before collision, the system energy W f can be expressed as follows:

W f = 1 2 C f V 2 d + 1 2 mv 2 f + 1 2 kx 2 f , (2)
Fig. 2. Model of actuator.

where C f is the electrostatic capacitance of the parallel plate actuator immediately before collision, x f is the displacement of the mass at collision, and m is the mass. In addition, according to the law of conservation of energy, the work W of the power source equals the difference between W f and W 0 , and therefore,

W =ΔQV d =ΔCV 2 d =(C f -C 0 )V 2 d = W f -W 0 , ( 3
) where ΔC is the difference in electrostatic capacitance between the initial state and that immediately before collision. Now the following can be obtained by substituting Eqs. ( 1) and (2) into Eq. ( 3):

1 2 mv 2 f = 1 2 (C f -C 0 )V 2 d - 1 2 kx 2 f . ( 4
)
Thus the velocity v f immediately before collision is derived as follows:

v f = √ 1 m (C f -C 0 )V 2 d -kx 2 f . ( 5
)
In addition the electrostatic capacitance C 0 in the initial state and C f immediately before collision can be expressed as follows:

C 0 = ε 0 S g C f = ε 0 S g -x f . ( 6
)
Here ε 0 is the permittivity of a vacuum and S is the area of the parallel plates. Since C f represents the capacitance immediately before collision, deformation of the spring k w (Fig. 2) after collision is not considered. The following is obtained by substituting the above equations into Eq. ( 5):

v f = √ 1 m ε 0 S ( 1 g -x f - 1 g ) V 2 d -kx 2 f . ( 7
)
In the real device, the driving voltage is V d = 150 V, the mass is m = 20 mg, the initial gap is g = 20 μm, the area of the parallel plate actuator is S = 2.4 × 10 -7 m 2 . and the gap immediately before collision is (gx f ) = 2 μm. Thus the displacement x f immediately before collision is 18 μm. Therefore, the energy kx 2 f accumulated in the spring is smaller than the first term by about three orders of magnitude and thus can be ignored. In addition x f is usually made as close to the gap g as possible, so that 1/g is negligibly small compared to 1∕(gx f ). Thus Eq. ( 7) can be simplified as follows:

v f = √ √ √ √ 1 m ε 0 SV 2 d g -x f . ( 8
)
Substituting the real device parameters into the above equation, we obtain v f = 0.035 m/s. In order to test the operation of the device, we carried out simulations using an analysis of electrostatic actuators based on a circuit simulator that was proposed earlier by our group [START_REF] Mita | An equivalent-circuit model for MEMS electrostatic actuator using open-source software Qucs[END_REF][START_REF] Mita | Multi-Physics analysis for micro electromechanical systems based on electrical circuit simulator[END_REF]. The device model in Fig. 2 was modified for the simulations as shown in Fig. 3. In the diagram, K1 is suspension module supporting the movable mass, MIII1 is the equation-of-motion solver, and PPpi1 is the parallel plate actuator module, including the virtual spring that simulates collision. As explained above, collision is modeled by application of pressure to a wall by a spring with an extremely large spring constant (10 6 N/m in the simulations). Using this model, we confirmed the conservation of momentum prior to the simulations. However, these simulations did not involve forward movement due to the impact force.

The simulation results are presented in Fig. diagrams, the collision phenomenon is realistically analyzed. In addition velocity at collision was 0.02 m/s, which is close to the value calculated by Eq. ( 8). The instantaneous impact force produced by collision was 15 mN to 35 mN, which is about 1000 times the force of several μN or several tens of micro Newtons generated by conventional electrostatic microactuators. In these simulations, we investigated how the impact force and velocity were affected by various parameters, including the mass m, spring constant k, and initial gap g. We found that a greater mass m contributes strongly to the impact force. As follows from Eq. ( 8), v f decreases with increasing mass; however, this decrease can be compensated by a greater driving force. That is, the driving voltage and the parallel plate area should be increased. These results agree with the above analysis of the velocity before collision.

Figure 5 shows the simulation results for operation at the resonant frequency obtained with the same model. The diagrams represent (a) the mass displacement, (b) the velocity, (c) the impact force, and (d) the displacement and velocity on an expanded time axis. In diagram (d), the dashed line represents the displacement, shown on the left axis, and the solid line represents the velocity, shown on the right axis. When the device is operated at the resonant frequency, the amplitude gradually increases and eventually the mass collides with the stopper. At collision, the waveforms are disturbed but essentially remain sinusoidal. The real device parameters were used in these simulations, and the obtained collision velocity (0.03 m/s) was close to that in Fig. 4 (0.02 m/s). When the spring constant is set higher, the vibration velocity increases and collision velocity likewise grows. Therefore, in the case of resonant driving, the vibration velocity can be greater than that obtained by Eq. (8) if the spring constant k is sufficiently large, that is, if the resonant frequency is sufficiently high. In addition the collision rate per unit time grows, resulting in a higher travel speed. Thus, a greater driving force and a higher travel speed than in single collisions can be obtained by driving at the resonant frequency. The simulation results showed that in case of resonant operation, a sinusoidal approximation is acceptable, and therefore we describe motion of the mass as follows:

x(t) = A sin(ωt). (9)
Here x(t) is the mass displacement, A is the amplitude, and ω is the angular frequency. The mass velocity v(t) is obtained by time differentiation of x(t):

v(t) = Aω cos(ωt). (10)
Because operation at the resonant frequency is assumed, the angular velocity ω is can be found from the spring constant k and mass m as follows:

ω = √ k m . ( 11
)
The vibration is limited by the stopper, and hence A #∼ g. Substituting this relation into Eq. ( 10), the maximum mass velocity v max at resonance is

v max = g √ k m . ( 12
)
When the driving voltage is sufficiently high, the collision velocity approaches v max ; thus this value can be used in design.

From the above results it follows that the collision velocity can be increased by setting the initial gap, spring constant, and driving speed as large as possible. The mass velocity is greatest at a displacement x = 0, and therefore the maximum collision velocity is not equal to v max calculated above. However, an increase in v max results in a higher collision velocity, and thus we set v max high in design. Now the impact force F c produced at collision can be expressed in terms of the collision time Δt and the velocity change ΔV as follows:

F c = m ΔV Δt . ( 13
)
That is, the impact force varies directly with the mass m.

The driving voltage is limited by the dielectric strength of the device, which depends on the material and the structure.

The initial gap and the drive electrode area are limited by the device dimensions. With these limitations in mind, we designed the device as explained below.

Design and Structure

The structure of the actuator is shown in Fig. 6. The driven mass is supported by an outer frame consisting of Fig. 6. Structure of device.

four springs. This device is driven by a parallel plate electrostatic actuator. The bottom layer of the SOI substrate provides a mechanical structure including the suspension and mass, and the top layer acts as a parallel plate actuator. The electrodes placed on the mass are grounded via cantilevers to the bottom level [START_REF] Takahashi | Switched-layer design of SOI bulk micromachined XYZ stage using stiction bar for interlayer electrical connection[END_REF].

As explained above, a large mass and a high collision velocity, that is, high v max , are required to obtain a strong impact force. For this purpose, the spring constant and gap should be set large, as suggested by Eq. ( 12). In addition, not only the impact force at collision but also the collision rate per unit time should be increased to obtain a higher travel speed of the whole device. In this study, a large portion of the chip area is allocated to the mass in order to implement a large mass. In addition there is a limit to the dielectric strength of the device, and therefore, after determination of the driving voltage, the drive electrode area is set as large as possible so as to obtain a stronger electrostatic attractive force at the given voltage. The spring constant too should be set as high as possible in order to increase the collision velocity and the collision rate per unit time; on the other hand, operation must be restricted to the abovementioned range of the driving force. In the case of driving at the resonant frequency, the amplitude grows by a factor equal to the Q-factor compared to DC driving, and Fig. 7. Fabrication process. therefore the gap g is set within Q times the displacement determined by the spring constant and the driving force. From the above considerations, we determined the mass, driving voltage, and spring constant and fabricated a prototype device.

Fabrication Process

The device was fabricated on an Silicon On Insulator (SOI) substrate. The SOI substrate had an active layer 80 μm thick, a silicon oxide film layer (BOX layer) 2 μm thick, and a handle layer 300 μm thick. Stable results were assured by using a simple fabrication process, as shown in Fig. 7. First, (a) aluminum was vacuum-deposited on both sides of the substrate to form the etching mask. Then, (b) the pattern was transferred to the resist using the photomask, the deposited aluminum was etched, and a mask pattern for silicon etching was created. After that, (c) silicon was deepetched from the backside using aluminum as the etching mask, thus forming the outer frame, mass, and suspension on the bottom layer (handle layer). Then, (d) deep etching with the aluminum mask was performed in the same way from the front side to form the electrostatic actuator on the top layer (active layer). Finally, (e) the silicon oxide film was removed from between the silicon layers.

Fabrication Results

The fabricated device is shown in Fig. 8. The device size was 5 mm × 5 mm and the thickness was about 0.4 mm. The suspension spring was 15 μm wide, 1.5 mm long, and 300 μm high. The size of the mass was 3 mm × 3 mm and the width of the outer frame was 0.5 mm. The mass and the outer frame had the same thickness, and therefore the mass ratio was equal to the area ratio. The mass constituted 50% of the entire device. The size of the PC board was 8 mm × 8 mm; copper wires 0.05 mm in diameter were soldered to the board in order to supply the driving voltage. The layout of the printed board is shown in Fig. 9. In order to connect the device to the power supply, the fabricated chip was mounted on a special printed circuit board, and power supply from the board was implemented by wire bonding. (In Fig. 9, the mass module does not contact the printed circuit board.) Fig. 10. Operation photographs.

Driving Results

In operation tests with a sinusoidal driving voltage, the resonant frequency of the device was 311 Hz. The driving results of the whole chip are illustrated in Fig. 10. The top, middle, and bottom photographs were taken, respectively, before actuation began, 10 s after actuation began, and 20 s after actuation began. Using thin copper wires, a sinusoidal voltage with an amplitude of 150 V, a bias of 75 V, and a frequency of 311 Hz was applied to the device mounted on the printed board. We found that the whole printed board including the chip traveled a distance of about 1 mm in 10 s while resisting frictional and restoring forces. Thus, the device traveled about 100 μm per second; considering the driving frequency of 311 Hz, the average travel distance per collision can be estimated as about 320 nm.

The total weight including the printed board was about 0.2 g, and the weight of the chip was about 20 mg. Thus the device is capable of driving about 10 times its own weight. Inclining the stand on which the device was placed, we found that the device started off at an angle of about 22 • . Thus, the coefficient of static friction between the device and the stand can be estimated as about 0.37. Considering the frictional force calculated with this coefficient as well as the restoring force of the copper wire, the driving force was of mN order. This is an order of magnitude smaller than the simulation value; however, we can expect to achieve more adequate simulations by making more accurate measurement of the collision time Δt and the actual collision velocity.

Conclusions

We considered an inertia-driven micro-actuator for space applications. In particular, we succeeded in impact driving of an entire chip including the printed circuit board, which was difficult with the previous MEMS actuator that produced driving force of about 10 μN. In the future, we plan to continue development aiming at a further increase in the driving force. In addition, we are planning to carry out operation tests in high-temperature and cryogenic environments in order to improve the for space applications.

Fig. 1 .

 1 Fig. 1. Principle of actuator.

Fig. 3 .

 3 Fig. 3. Simulation model of actuator.

  [START_REF] Mita | An equivalent-circuit model for MEMS electrostatic actuator using open-source software Qucs[END_REF]

  Fig. 4. Simulation results.

Fig. 5 .

 5 Fig. 5. Simulation results.

Fig. 8 .

 8 Fig. 8. Photograph of device.

Fig. 9 .

 9 Fig. 9. Cross-section view of device.


Phillip Colella

GLIMM'S METHOD FOR GAS DYNAMICS

Keywords: random choice method, gas dynamics, Glimm 's method

We investigate Glimm's method, a method for constructing approximate solutions to systems of hyperbolic conservation laws in one space variable by sampling explicit wave solutions. It is extended to several space variables by operator splitting. We consider two problems. 1) We propose a highly accurate form of the sampling procedure, in one space variable, based on the van der Corput sampling sequence. We test the improved sampling procedure numerically in the case of inviscid compressible flow in one space dimension and find that it gives high resolution results both in the smooth parts of the solution, as well as at discontinuities. 2) We investigate the operator splitting procedure by means of which the multidimensional method is constructed. An 0(1) error stemming from the use of this procedure near shocks oblique to the spatial grid is analyzed numerically in the case of the equations for in viscid compressible flow in two space dimensions. We present a hybrid method which eliminates this error, consisting of Glimm's method, used in continuous parts of the flow, and the nonlinear Godunov method, used in regions where large pressure jumps are generated. The resulting method is seen to be a substantial improvement over either of the component methods for multidimensional calculations.

1. Introduction. The problem which motivates this study is the numerical calculation of time-dependent, discontinuous solutions to compressible fluid flow problems in one or more space variables. There are three criteria which such approximate solutions must simultaneously satisfy.

1) The approximate solution must be reasonably accurate in regions where the flow is smooth. Continuous waves should move at the correct speed, have the correct shape, steepen or spread at the correct rate.

2) Discontinuities which are transported along characteristics should be modeled in the approximate solution by sharp jumps which are transported at the correct speed. Examples of such discontinuities are: contact discontinuities (across which the density and temperature have jump discontinuities while the pressure and velocity remain continuous); the interface between two different materials, or between two different thermodynamic phases of the same material; lines or surfaces across which the solution is continuous, but some derivative of the solution is ,not.

3) Nonlinear discontinuities should be computed stably and accurately. Such discontinuities occur, for example, when there is mass transported across the discontinuity, as in the case of shock fronts in an ideal gas.

The main method used for computing such solutions has been to solve a set of finite difference equations which approximate the differential equations of motion. However, it is difficult to construct difference methods which satisfy all three of the above criteria simultaneously. For example, it is well known that a high order difference method may generate oscillations behind a shock. A first-order method will generally treat the same shock correctly, but numerical diffusion will cause it to give lowresolution results in continuous parts of the flow.

We will be examining an alternative approach to computing discontinuous flows, known variously as Glimm's method, the random choice method, or the piecewise sampling method. This method was first used by Glimm [START_REF] Glimm | A numerical method for two phase flow with an unstable interface[END_REF] as part of a constructive existence proof of existence of solutions to systems of nonlinear hyperbolic conservation laws. It was developed by Chorin [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], [START_REF]Random choice methods with application to reacting gas flow[END_REF] into an effective numerical method in the case of gas dynamics. In the first reference Chorin also introduced a multidimensional version of the scheme; in the second, he applied the method of reacting gas flow in one space variable. Since that time, the method had been used to compute compressible flow in cylindrical or spherical geometry (Sod [START_REF] Son | A numerical study of a converging cylindrical shock[END_REF], [START_REF]A numerical method for a slightly viscous axisymmetric flow with application to internal combustion engines[END_REF]), and in applications to some problems in petroleum engineering (Concus and Proskurowski [START_REF] Concus | Numerical solution of a nonlinear hyperbolic equation by a random choice method[END_REF], Albright, Concus and Proskurowski [START_REF] Albright | Numerical solution of the multidimensional Buckley-Leverett equation by a sampling method[END_REF], Glimm, Marchesin, and McBryan [START_REF] Glimm | A numerical method for two phase flow with an unstable interface[END_REF], Glimm, Marchesin, Isaacson, and McBryan [START_REF] Glimm | A shock tracking method for hyperbolic systems[END_REF]).

Although one computes solutions on a grid with Glimm's method, it is not a difference method. Rather than computing a weighted sum to arrive at the value of the solution at a grid point, one samples values from an explicit wave solution. Thus, the method has built into it an approximate form of wave transport and interaction, without the smoothing of such information inherent in averaging.

The introduction of such a sampling technique as a numerical method is quite recent, compared to the length of time difference methods have been in use, and has not been subject to the extensive scrutiny and application from which the latter has benefited. One of the purposes of this study is to indicate some of the features of Glimm's method which might make developing it worth the effort, as well as a few of the directions in which the development might go.

We consider in this study two fundamental problems. 1) We introduce a more accurate form of the sampling procedure for the onedimensional method than that used in [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], based on the van der Corput sampling sequence. We compare the performance of the van der Corput sampling and the previously used random sampling schemes.

2) We investigate the operator splitting procedure by which Chorin constructs a multidimensional scheme from the one-dimensional method. A source of error stemming from this procedure, not noticed in [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], is analyzed here and a method for eliminating it is proposed and tested.

In one dimension Glimm's method, with the appropriate sampling, is seen to be superior to any difference method in meeting the three criteria given above. The final method obtained for multidimensional calculations, although it does not share the special properties of the one-dimensional method, has a number of interesting features, and is worthy of further investigation for its own sake.

This paper is divided into three sections. In § 2 we discuss Glimm's method as applied to gas dynamics in one space variable. We define the van der Corput sampling sequences, and compare the van der Corput and random sampling strategies. In § 3 we describe the operator splitting technique. In § 4 we compare Glimm's method to some difference methods and give some conclusions, and suggestions for future work.

Gas dynamics in one space variable.

We want to construct approximate solutions to the initial value problem for Euler's equations for the motion of a one-dimensional, compressible, inviscid gas with a polytropic equation of state:

au + aF(U) o, at ax (2.1) U(x,t)=U:Rx[O,T]~R 3 , U(x,O)=u(x) gtven, U= (~), F(U)= ( m 2 ~+p ), E (E+p)m/p
where p is the density, m is the momentum, and E is the total energy per unit volume of the gas. We can express in terms of these variables the more familiar quantities v the velocity, and £, the internal energy of the gas:

m v=-, p
The pressure p which appears in the equation is a function of p, e: p = ( y -l)pe = Ap '\

where the constant y > 1 is the ratio of specific heats for the gas. Another quantity of interest is the entropy, S, defined (up to an additive and a multiplicative constant) S =log (pp -y) =log A.

The system (2.1) is a first-order, hyperbolic system of conservation laws; i.e., the 3 x 3 matrix A ( U), the J acobian ofF, has three real eigenvalues,

A1(U)=u-c, Az(U)=u, A 3 (U)=u+c,
where c = J yp/ p is the adiabatic sound speed. Ai, i = 1, 2, 3 are the characteristic velocities associated with the three modes of wave propagation for (2.1).

Since we are dealing with piecewise smooth solutions, we interpret (2.1) in the sense of distributions. That is, if U(x, t) is discontinuous along a piecewise smooth curve (/(!), t), then dl/ dt = s(t) must satisfy

(2.2) s(t)(UL(/(t), t)-UR(l(t), t)) = F(UL(l(t), t))-F(UR(l(t), t)),
where

UL,R(l(t), t) = lim U(/(t) + e, t). € tO,€ LO
The discontinuity must also satisfy the entropy conditions (Courant and Friedrichs [START_REF] Friedrichs | Supersonic Flow and Shock Waves[END_REF]). Discontinuities calculated using any of the numerical methods discussed in this paper satisfy the entropy conditions if they satisfy (2.2).

The simplest initial value problem for which discontinuities appear is one for which the initial data is constant on either side of the origin, where it has a jump discon tin ui ty:

where we denote X <0, X >0, This problem is known as the Riemann problem; its solution is a fundamental component of Glimm's method. The special case of the Riemann problem in which UL = UR = 0 is often referred to as the shock tube problem. The solution of the Riemann problem is discussed extensively in Chorin [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], Courant and Friedrichs [START_REF] Friedrichs | Supersonic Flow and Shock Waves[END_REF], Godunov [START_REF] Godunov | Difference methods for the numerical calculations of discontinuous solutions of the equations of fluid dynamics[END_REF] and Sod [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] and van Leer. [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF] In [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF] and [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] detailed instructions for constructing solutions numerically are given; thus we will describe only qualitatively the structure of the solution.

Two general properties of the solution to a Riemann problem are that it is self-similar, i.e., U(x, t) = h(x/t) for some piecewise continuous h: R ~ R 3 , and that it has the following additivity property: if U(~, t~) = UM E R 3 for some R, then the function

{ U(x, t), Ut(X, t) = UM, X -<I: t ~' X ->~ t
is the solution to the Riemann problem with left and right states UL, UM. Similarly, the function X -<I:

t ~' X ->~ t
is the solution to the Riemann problem with left and right states UM, UR. Geometrically, this says that the solutions Ut, U2 fit together to form U. One can divide the (x, t)-plane into four regions, I, II, Ill, IV, where U(x, t) is constant (Fig. 1). These four regions are connected by three waves, each associated with one of the characteristic speeds. These are: a backward facing sonic wave (associated with u-c = At(U), between l 1 ,b and l2,b; a contact discontinuity (associated with u = ,\2( U), occurring across the line Is; and a forward facing sonic wave (associated with u + c = A 3 ( U) between /1.1 and /2.1• The pressure and velocity are continuous across the line Is so they are equal to some fixed values p*, u* in II and Ill. Only the density changes across As is discussed in [START_REF] Friedrichs | Supersonic Flow and Shock Waves[END_REF], the hydrodynamic waves are uniquely determined by knowing the state of the gas on one side of the wave, and only the pressure on the other. For the backward facing wave, for example, there are two possibilities. If p* > PL then u * < uL, pI > PL, lt,b = l2,b and the wave is a shock associated with the characteristic velocity u -c. If p * < PL, then we have a backward facing centered rarefaction wave:

Is (t) = u * t, from p f to p ~ . U(x,t)=UL u(x,t)=u* p(x,t)=p* p(x,
l 1 ,b < l 2 ,
b, p(x, t) and u(x, t) are continuous strictly monotone decreasing functions of x/ t, and u (x, t) a continuous strictly monotone increasing function of x/ t, for (x, t) between lt,b and l2,b• The description of the forward facing wave is the same, replacing UL by UR, u by -u, and u+c by u-c.

In Fig. 2 we show the solution at a fixed time to the shock tube problem

(2.3) PL = 1.0, PL = 1.0, PR= 0.1, PR= 0.125, 1' = 1.4.
The waves which occur are a backward facing rarefaction wave (A), a forward facing shock (B), and a contact discontinuity (C).

+-(B)

We can now describe Glimm's method for solving approximately the initial value problem (1.1). Let ax be a spatial increment, llt a time increment. We assume that, at time n at, the approximate solution is constant on intervals of length ax:

(2.4) u<~x) (x, nat) = Uj E R 3 (j-!)llx <X< (j + !)llx, j = 0, ± 1, ±2, • • • .
We wish to compute an approximate solution which at time nllt has the same property:

U(~x) (x, (n + 1)/lt) = U'j+ 1 •
The procedure is given as follows:

(1) Define U~ (x, t)nat < t < (n + l)at to be the exact solution to the initial value problem for ( 1.1) with initial data given by (2 .4). The initial data consist of intervals where the solution is constant, separated by jump discontinuities; i.e., we have a succession of Riemann problems. If at is sufficiently small, then by finite propagation speed the waves from adjacent discontinuities do not intersect each other and the solutions to the adjacent Riemann problems fit together to given U~ (Fig. 3). A condition on llt which guarantees that the waves do not intersect is When doing calculations, one usually uses the more easily verified (2.5)

( j-3/2)L1x n u. 1 J- llt 1 -=A<- Ilx 2 (j-l/2)
~t =A< u sup lu(~x\x, nAt)l + c(~x) (x, ndt) = u sup (lull+ cj), aX xeR j
where u is a constant, 0 < u < !.

(2) Choose a n+l E [0, 1) and take u;+l = u~ ((j-!+a n+l) ax, (n + 1) at).

See Fig. 4.

Thus we obtain a solution at time (n + 1) which depends on a sequence a= a 1 , a 2 , • • • ; much of the remainder of this section will be devoted to determining the best choice for the sequence a.

At first glance, this method might look complicated, but in fact it requires the evaluation of the solution to a Riemann problem once per zone per timestep. Let ------------------------------~L------

h. 1 (n-(j -i)ax) J-2 ,n t-nat

Then

The procedure given here is slightly different from that used previously, in that the mesh is fixed, rather than shifting by ax/2 every timestep. Sampling back to a fixed grid shows that the relation to Godunov's method is immediate: in Godunov's method, u;+l is taken to be

1 J (j+~)~x -;- U~(x, at) dx; ax u-t)~x
in Glimm's method, one chooses a representative point value of the local exact solution.

The mechanism by which Glimm's method models wave propagation in a gas is most easily demonstrated by the following example. Let UL, UR be the left and right states of a Riemann problem whose solution consists of a single discontinuity propagating at speed s > 0. The exact solution for this problem is x<sT, x>sT.

We will solve this initial value problem using Glimm's method. First, it is obvious that, for any time step n there is an 1 (n) =jo-t, jo an integer, such that j<l(n), j>l(n). I (n) is the location of the shock in the approximate solution, and satisfies l(n+1)={11((n) 1) ~ffa:::>As, 1(0)=!.

n + 1 a <As, so that

I ( n) = I ( 0) + N {j = 1, • • • , n ; a i E [ 0, As)},
where N{j = n1 + 1, • • •, n2, ai E I} denotes the number of j, n1 <j < n2 such that ai is contained in I. We say that a sequence ii is equidistributed if the proportion of times that ai is contained in I is asympotically equal to I I I, the length of I; i.e., if we define then ii is equidistributed if limn 2 ~oo, 5(ii, nt, n 2 , I)= 0 for each fixed nt, I. Given this notation we write

l(n)~x = I(O)~x + ~xAsn + B(ii, 0, n, 0, As))n ~x = l(O)ax +nat( 1 + ~ 8(ii, o, n[(O, As))). If ii is equidistributed, then I (n) ~ sT in the limit n ~ oo, In ~t-Tl < ~t, ~t/ ~x =A > 0.
Thus the shock in the approximate solution at each time step either moves by ~x or does not move at all. Over many time steps, the cumulative displacement is close to that of the exact solution, the leading term in the .error being proportional to 5(ii, 0, n, [0, As)). In general, a piecewise continuous flow will be represented by 0(1/ ~x) waves of strength O(~x ), all having differing speeds, as well as an 0(1) number of discontinuities of strength 0(1). Furthermore, the speeds and strengths of the waves will be changing in a piecewise continuous fashion as a function of time. In order to model such a flow correctly by the above mechanism, one needs to choose ii such that 5 is as small as possible, uniformly in I, n 1 for n2-n 1 large relative to 1, but (n 2 -n 1 )~t small relative to the characteristic times in which the wave speeds change. The sampling procedure given below seems to be optimal from the point of view of these requirements.

The simplest form of this sampling sequence is due to van der Corput (see [START_REF] Hammersley | Monte Carlo Methods[END_REF]).

Let be the binary expansion of

n = 1, 2, • • • . Then m n _ ~ • 2 -(k+1) a -i.J lk • k=O
The easiest way to see how the sequence is constructed is to write down the first few elements in it:

k i k+1 . 4< a < - 4 -tf i = j(k) mod 4, k = 0, 1, 2, 3
, where j(O) = 0, j(1) = 2, j(2) = 1, j(3) = 3. In general, if one divides the unit interval into the subintervals (r2-s, (r + 1)2-s)r = 0, • • • , 2s -1, then for each r there is exactly one q for which q 0 < q < q 0 + 2s such that aq E [r2-s, (r + 1)2-s).

We will have need of a variant of this procedure for use in multidimensional problems. Let kt, k2 > 0 be integers, kt > k 2 relatively prime. The (kt, k 2 ) van der Corput sampling sequence ii is given by where

m a n_ ~ k-{l+l) -i..J ql . 1 ' l=O m qt = k2i1 mod kt and L itki = n l=O
is the base kt expansion of n. Thus the binary van der Corput sampling sequence given above is the special case kt = 2, k2 = 1.

All the van der Corput sampling sequences are equidistributed. In fact 8(ii,nt,n2,I)<(Ctlogkt(n2-nt)+C2)/(n2-nt), where Ct, C2 are constants depending on kt, k2 but not on nt, n 2 or I. For the binary van der Corput sequence Ct = 3, C2 = 1. In the example given above, this gives an error bound of O(axjlog axj).

In previous computational work for gas dynamics using Glimm's method, random sampling was used; i.e., the values were drawn from a random number generator implemented on the computer, usually with some variance ~duction technique, such as stratification for random sampling, for which 8 = 0(1/~n), giving an error bound in our simple shock example of 0(~.

Lax [START_REF] Lax | Hyperbolic systems of equations and computing[END_REF] proposed the use of a nonrandom equidistributed sequence due to Richtmyer and Ostrowski, defined by an =~~mod 1, where r is an integer which is not the square of another integer.

We shall not discuss the Richtmyer-Ostrowski sampling sequence in detail here, save to note that in numerical experiments, and in simple analytical examples, one obtains results using the Richtmyer-Ostrowski sequence similar to those obtained using van der Corput sampling. However, the bound on 8 is stronger for the van der Corput sequence than that obtained for the Richtmyer-Ostrowski sequence, as well as being explicitly uniform in I and n1; uniformity in I and n1, does not hold explicitly for the Richtmyer-Ostrowski sequence. Also, van der Corput sampling has some special properties which guarantee that certain qualitative features of the continuous part of the solution preserved in the approximate solution, at least for simple waves (see [START_REF] Colella | An analysis of the effect of operator splitting and of the sampling procedure on the accuracy of Glimm's method[END_REF], [START_REF]Error analysis of Glimm's method[END_REF]). Finally, van der Corput sampling has several straightforward extensions to two or more dimensions which guarantee good distribution properties in the square, even for finite sample sizes. In contrast, it has been pointed out by Maltz and Hitzl [START_REF] Maltz | Variance reduction in Monte Carlo computations using multidimensional Hermite polynomials[END_REF] that such an extension of the Richtmyer-Ostrowski sequence can give rise to poor distribution in the square for finite sample sizes.

In an effort to understand the errors introduced by the interaction of the sampling and variations in time in the wave speeds, we consider the following class of test problems. The initial data consist of two discontinuities located at x and x,, separated by constant states:

x <xt, x 1 <x <x, x >x, ( Pl,m,r ) u = p u l,m,r l,m,r l,m,r • 1 +1 u2 Pt,m,r/ ( 'Y-) 2Pl,m,r l,rn,r
We choose U~, Um, U, such that U 1 and Urn can be connected by a forward facing shock and that Urn and U, can be connected by a forward facing centered rarefaction wave (Fig. 5).

The shock overtakes the rarefaction, the cancellation between them weakening both (Fig. 6,(A)). The nonlinear coupling between the tnodes produces waves of the other two families in back of the shock and moving to the left, away from the shock. These are, a backward facing compression wave (Fig. 6, (B)), and a strong entropy I density wave (Fig. 6,(C)) advected passively by the velocity field u (x, t).

In Figs. 789we show the calculation of such a problem using Glimm's method with, respectively, a random sampling sequence, a stratified random sampling, and the binary van der Corput sampling sequence. The initial data are: Also plotted on each of the graphs, with a solid line, is a solution obtained using Glimm's method, with van der Corput and Llx = .0025. Having compared the latter solution with a similar one done for Llx = .005 we found that the two results differed by less than .5°/o, so that the method has converged for Llx = .0025. For the purposes of comparing the various Llx = .01 solutions, we treat the Llx = .0025 solutions as exact, against which the Llx = .01 solutions can be compared. The sampling governs the rate at which the shock and rarefaction interact. If s~+ 1 ; 2 is the speed of the shock, located between zones q and q + 1 at time step n, and A~= u~+l +c~+b then the shock will cancel with a piece of the rarefaction wave, and produce more wave of the other two families, at time

Pt = 28.68,
step n + 1, if and only if n+l [ Llt ( n ) Llt n ) a E Llx max A+,o , Llx Sq+! •
Thus the loss of gradient information observed in the randomly sampled solution (Fig. 7) is a result of random fluctuations in the rate of interaction between the shock and rarefaction which is producing the wave. The use of stratified random sampling (Fig. 8) produces smoother profiles than those obtained with the unmodified random sequence, but the shape of the entropy wave is incorrect; in particular, there is a sizable deviation in the density profile, a failure to conserve mass. The profile obtained using van der Corput sampling (Fig. 9) is in much closer agreement with the ax = .0025 result, the rate of wave production being modeled much better than in the other two cases. In fact, if one uses van der Corput sampling, one can use a much coarser mesh and still get good results for this problem. In Fig. 10 we present the results obtained on this problem with binary van der Corput sampling, and ax = lo. The absolute locations of the waves, and their locations relative to each other, are correct to within ax; more important, the size and shape of the waves, which are more sensitive to the cumulative error introduced by the sampling, are in very close agreement with the ax = .0025 result. In all the calculations, the shock discontinuity is sharp, as guaranteed by Glimm's method.

3. Operator splitting. In [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], Chorin proposed a method for computing multidimensional unsteady compressible flow using Glimm's method by means of operator splitting. We can write the equations of motion for an ideal gas in two space

dimensions as

U= p m n E au a a -+-(F(U))+-(G(U)) = 0,
at ax ay

U(x, y, t)= U: R 2 x[O, T]~R 4 ,
U(x, y, 0)=</J(x, y)<fJ: Here pis the density, m is the x-component of momentum, n is they-component of momentum, and E is the total energy. We can express the velocity i5 and the internal energy e in terms of the above variables: Vx =m/ p is the x-component of the velocity, Vy = n/ p is the y-component of the velocity, and e = E/ p -!(v~ + v;).

The pressure p is a function of p and e: p = ( y -l)pe, where y, the ratio of specific .

o.o-•• CO lS) lS)
heats, is a constant assumed to be greater than 1. Thus, as was the case for one space variable, the value U at a given point is uniquely determined by the values of p, p and v at that point.

We wish to construct approximate solutions

U Ax,Ay ( A ) un R4 x, y, n at = i,i E , (i -!) dx <x < (i +!) dx, (j -!) dy < y < (j +!) dy,
where x, y are spatial increments, t is a time increment, and i, j, n are integers, n > 0.

Assume we know UZi and want to find UZt 1

; the procedure is as follows. 1) For each j perform one time step of Glimm's method for the equation av a

at+ ax (F(V)) = o. k . . .. 1 d V 0 un s h 1 V 1 un+l ( d h' d
ta tng as tnttta ata i = i,i• et t e resu ts i = i,i 2 we enote t ts proce ure by (L~tUn)i,i = uzt!).

2) For each fixed i perform one time step of Glimm's method for the equations av a at+ ay (G(V))=O, taking as initial data vJ = uzt!, time step t. Set the result V] = UZt 1 (we denote this procedure by (LXtun+~)i,j = uzt 1

).

The solution thus derived at time (n + 1) ~t is interpreted as being the piecewise constant function

UAx,Ay (X, y, (n + 1) ~f)= uztl, (i -!)~X <x < (i +!)~X, (j -!) ~y <y <(j +!) ~y.

A necessary condition on the time step t is that it must satisfy (1.4) for each of the one-dimensional calculations ~t < u Il'l:ax (\v~.i.i\ + cZi),

'-lX z,J (3.1)

O<CT<!.

The above procedure is formally the same as is done to construct multidimensional difference methods from one-dimensional ones. However, the mechanism by which Glimm's• method propagates the solution to the equation in one dimension is rather different from that of different methods, as it requires many time steps for the cumulative effect of the sampling to give the correct wave speeds; therefore the actual justification of the splitting procedure, currently unknown, is likely to be quite different than the usual truncation error analysis for difference methods.

The Riemann problems in question are easily solved, given the solution for one-dimensional gas dynamics. For example, to solve Riemann's problem for

av a -+-(F(V)) = 0, at ax V ( ----- p +(-2+-2)~ = p, PVx, pvy, ' }' _ 1 V x V y 2J, take the solution p (x, t), p(x, t), u (x, t) in § 2 with p(x, t) = p(x, t), p(x, t) = p(x, t), Vx(x, t) = u(x, t),
and Vy (x, t) = Vy,L if (x, t) is to the left of the contact discontinuity Is, Vy (x, t) = Vy,R if (x, t) is to the right of the contact discontinuity Is. Thus in the x-sweep, we have ordinary one dimensional gas dynamics, with the discontinuity in Vy passively advected.

To solve the Riemann problem for avjat+ajay (G(V))=O, interchange the roles of Vx and Vy.

To test the validity of this procedure, we looked at the simplest two-dimensional test problem possible. We took our computational domain to be the unit square with the computational mesh aligned with the x-and y-axes, and took the initial conditions to be

X <y, X> y, PL PLVx,L PLVy,L PLf(y-1)+(pL/2)(v;,L +v~.L) PR PRVx,R PRVy,R PR!( 'Y -1) + (pR/2)(v;,R + v~.R)
This is the Riemann problem, for which we have an analytic solution. Computationally, it is a two-dimensional problem, since the initial discontinuity is at a 45° angle to the mesh directions.

We denote by v" the component of the velocity normal to x = y, Vt the component parallel to x = y:

Vx-Vy Vn = J 2 , Vx,R-Vy,R Vn,R = Vx,L-Vy,L Vn L = 1 - , ' -v2 Vx,R + Vy,R Vt,R = J 2 , Vx,L + Vy,L VtL = 1 - • ' -v2
Throughout these test calculations we will set Vr,L = Vt,R = 0; i.e., we will be looking at problems for which there is no slip line in the exact solution. Unless otherwise indicated, the calculations shown were done on a 50 x 50 grid: ~x = ~y = .02. The results of the calculations are displayed by plotting the profiles of various quantities along the line y = 1x, and comparing them with the exact solution. In these plots, the computed values at the mesh points are graphed as circles, interpolated by a dotted line: the exact solution is plotted as a solid line. When boundary conditions are required, we assume the solution is constant on lines parallel to the initial jump. This was quite effective in preserving the symmetry of the solution, and enabled us to run for long times without noise from the boundary affecting the results. The one-dimensional calculations using Glimm's method in the x and y directions require sampling sequences ax, ay which we took to be two independent van der Corput sampling sequences: dx was the (3, 2) van der Corput sequence, and dy was the (5, 3) van der Corput sequence. This insured optimal distribution in the square [0, 1) X [0, 1).

(3.2)

In Fig. 11, we show the results for the following problem: The exact solution is a strong, right facing shock. It is almost stagnant (after 175 time steps, the exact shock point has moved only two zones). By this time, the oscillations (80°/o of the exact post-shock value in the pressure) have begun to make themselves known by a three-zone error in the shock location, the shock moving a distance more than two times greater than it should have. We see substantial values (60°/o of lvn,L-Vn,R j) for vr (x, y, t), the tangential component of the velocity appearing.

Finally, the density profile shows a substantial deviation from conservation of mass.

The fundamental reason why large errors occur in this problem is that, although each half-step L~r, L~r models the resulting one-dimensional gas dynamics well, the problem it is modeling is 0(1) incorrect from the point of view of the two-dimensional flow. For example, consider the problem one solves (one for each value of j) in the first x-pass in the test problem (3.1). They are each the same Riemann problem for a one-dimensional gas flow, with the jump taking place along the diagonal. The left and right states

PL,R PL,RUL,R PL,R(Vy)L,R PL,R!(, -1) + cai,R + cv;)L,R)(PL,R/2)
for the one-dimensional problem are

Vn,L Vy,L = ~2'

The jump in the velocity, uL-uR, is less than Vn,L-Vn,R so a weaker forward facing shock than that of the original two-dimensional problem is produced, as well as a backward facing rarefaction wave. If we sample anywhere in the fan other than the left or right states, we get ( Vx) zt! > UL, UR. The new values

( n+! n+! n+! V X ) i,j ' p i,j ' p i,j
depend only on the sampling value a! and the ratio At/ Ax but not on il.t and il.x separately. So the difference between these and the exact answer is an 0(1) quantity relative to the mesh spacing. In particular, there is an 0(1) contribution to the tangential component of the velocity. Since there has been an 0(1) change in the thermodynamic variables p and p, there is no reason for the y-pass to produce a tangential velocity to cancel the one produced by the x-pass, and in fact it does not. Similar phenomena occur for a shock tube, (Fig. 12) or a Riemann problem whose solution consists of two centered rarefaction waves.

The above failures in the splitting procedure in situations when there are discontinuities in p, i5 can be viewed as a consequence of an invalid interchange of limiting procedures. Analytically, shock solutions are obtained as limits of viscous solutions as some set of diffusion coefficients go to zero. One might try to obtain the shocked solutions by using an operator splitting method to solve the viscous equations; the splitting procedure is then known to converge as at~ 0. Then, in the in viscid limit, the viscous solutions converge to the physically correct shocked solutions. In a difference method, the two limiting procedures take place simultaneously, with the coefficients multiplying the numerical diffusion approaching zero with at. The use of operator splitting with Glimm's method corresponds to letting the diffusion coefficients vanish for nonzero at. This interchange of limits is valid for continuous solutions, or near contact discontinuities, but near discontinuities in p or v the two limiting procedures are singular with respect to each other, and cannot be interchanged freely.

In order to correct this problem, we replace Glimm's method at discontinuities in p and v with a conservative finite difference method; the method we use in the nonlinear Godunov method (Godunov [START_REF] Boris | Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works[END_REF], Richtmyer and Morton [START_REF] Morton | Difference Methods for Initial-Value Problems[END_REF]) adapted for Eulerian coordinates (Godunov et al. [START_REF] Goounov | A comutational scheme for twodimensional nonstationary problems of gas dynamics and calculations of the flow from a shock wave approaching a stationary state[END_REF]). We describe the procedure for advancing this hybrid Glimm-Godunov method by one timestep, in one space dimension; the extension to two space dimensions, is achieved by an operator splitting procedure like the one described above. First, one calculates the exact solution to the initial value problem to (1.1), as before. At those mesh points where one uses Glimm's method, one samples as before. At those mesh points (jax, (n + 1) at) where one wishes to use Godunov's method, one sets 1

1 U+! )ax ( Uj+l )Godunov =A ue (x, (n + 1) at) dx. L.lX U-! )ax

By integrating the conservation law over the rectangle [{j-!) ax, (j+!) ax]x[n at,

(n + 1) at] we obtain, using the notation of § 2,

( Uj+l )Godunov = Uj + (F(h i-i,n (0))-F(h j+i,n (0))) :;,
where F is the vector of fluxes for the conservation law being integrated. Finally, we need a prescription for deciding whether to use Glimm or Godunov. Let -ko<k-j<ko+1, -ko<k-j<ko+1, where p/~! is the pressure in the region separating the two sonic waves which come from the Riemann problem centered at ((j -!) ax, n at) (see [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF], and the appendix to this paper). Then the prescription for choosing Glimm and Godunov is max min

{(u n+l) "fPi -pi C n+l _ j Godunov 1 J?in > o, ( u j ) hybrid - p 1 ( u;+l )alimm otherwise.
Here eo, ko are constants to be set at the beginning of the calculation. Roughly, eo is a measure of the strength of the weakest sonic wave in the problem that must be treated as a discontinuity, and k 0 + 1 is the effective width of a discontinuity. For weak problems (excess pressure ratios< 5), it suffices to set k 0 = 1. For stronger shocks, it appears to be necessary to set k 0 = 2. In all the calculations presented here, .05 < C 0 < .2. The consequence of the nonoptimal choice of parameters is a loss of accuracy, not of stability: failing to detect a pressure jump results in noise; using Godunov's method unnecessarily results in the smoothing of relevant wave structures. One can make several minor modifications of the method described above. The fact that we are using Godunov's method near strong pressure jumps makes it possible to introduce some approximations into the solution of the Riemann problem. Another consequence of sampling only if the sonic waves are weak is that, by accounting approximately for the interaction of those waves, one can use a larger time step, allowing u < 1 in (3.1), which is the time step restruction for Godunov's method. We have implemented both of these changes in the method for the examples computed here; in an appendix to this paper, we describe the details of the algorithms used. Finally, we noticed that the first-order splitting algorithm described above can lead to errors near very strong shocks (excess pressure ratios greater than lOO) computed using Godunov's method. In particular, large tangential components of velocity are generated behind the shock. We found that the use of the Strang splitting algorithm U n+2 LX LYLY LX un = ~t ~t ~t ~t reduced this error to the level found in regions where the flow is continuous.

In Fig. 13, we show the results for the problem (3.2) using the hybrid method. Since the solution is a shock discontinuity separating two constant states, this calculation is mostly a test of how well the nonlinear Godunov method computes a strong 15 ---~----~----~--~~--- shock. The dip in the density behind the shock is a starting error, common to most conservative difference methods; it comes from starting a strong shock as a jump discontinuity. Since there is no numerical viscosity away from the shock the oscillation is not damped, but flows downstream unchanged.

Figure 14 shows the results obtained for the shock tube problem (1.3) calculated as a diagonal Riemann problem; Fig. 15 shows the result for the same problem using Godunov's method alone. The hybrid method treats the shock correctly, as opposed to the Glimm's method alone (Fig. 12). The hybrid method is also an improvement over Godunov's method alone: the three waves are clearly resolved; in particular, the contact discontinuity is spread over only three zones. We have found that, in general, the hybrid method spreads any discontinuity over a small (1-4) number of zones, independent of the zone size, regardless of whether the discontinuity is a shock, contact discontinuity, or slip surface.

In order to test this method on a more complex problem, we computed a two-dimensional Cartesian shock reflection problem used by van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF] as a test problem; see also Woodward and Colella [START_REF] Woodward | High resolution difference schemes for compressible gas dynamics[END_REF]. The computational domain is a channel of unit length, open at both ends. For x < .2, the channel has width ~; at x = .2, the lower side of the channel is constricted, so that the width of the channel is 1~ for x > .2. Reflecting boundary conditions are imposed on the upper and lower sides of the channel, and on the segment x = .2, 0 < y < 1 1 s. The solution is assumed to be continuous at both ends. The initial conditions for this problem are those of uniform flow throughout the tube: p(X, y, 0) = 1, p(X, y, 0) = 1.4, Vx(X, y, 0) = 3, Vy(X, y, 0) = 0, y = 1.4; with these initial conditions, a detached shock reflects off the constriction, and reflects off the upper side of the channel, having formed by time t =! a three-shock Mach reflection configuration. According to Woodward [START_REF] Woodward | High resolution difference schemes for compressible gas dynamics[END_REF], the correct location of the Mach stem along the side of the channel is right above the constriction at x = .2; the Mach stem should extend about one-fourth the distance across the left end of the channel.

' Q •n ••o '0 0.4 0 0.4 QQ 0.2 Q 0.2 b 0 IS) N :j"i (.() CO lS) IS) N :j"i (.() CO lS) CS) CS) IS) IS) ..-t IS) IS) IS) IS) ..-t PRESSURE DENSITY 1. 0 3.0 0.8 0.6 2.5 •oo 0 0.4 •o •o 2.0 •o 0 '0 0.2 •o 0 0 0 0 oooon 0 0 ao 1.
We show results obtained using Godunov's method .. by itself (Fig. 16), and the hybrid method, with two different zone sites (Fig. 17 stem about half the correct length. The slip line emerging to the right of the triple point is spread over four or more zones, except right at the shock. The solutions obtained with the hybrid method both have the shock in the correct position to within one zone length, and the length of the stem differs from the cofrect length by two zone lengths. Both of the hybrid calculations have the slip line spread across two zones for its entire length.

4. Discussion and conclusions. In one space variable, Glimm's method has built into it an approximate form of linear and nonlinear wave propagation along characteristics, without the smoothing of such information, as occurs in difference methods, and without any complicated bookkeeping; the sampling procedure determining the weakest wave or wave interaction to be resolved. The motivation for using van der Corput sampling is that one obtains the best possible representation of the wave propagation in Glimm's method, independent of the speed of the waves. This is essential for the correct representation of continuous waves, particularly those produced by nonlinear wave interactions.

We would like to compare the performance of Glimm's method to that of difference methods. Sod [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] performed such a comparison, using a one-dimensional shock tube problem. The results obtained there were not the best possible, due to the use of stratified random sampling. On the other hand, comparing difference methods to Glimm's method on this problem is not entirely fair, either. As is pointed out in [START_REF]Random choice methods with application to reacting gas flow[END_REF], it follows from the additivity property for solutions to the Riemann problem that the only values taken on by the computed solution are ones taken on by the exact solution, as well. In any case, we present in Fig. 19 the calculation done with Glimm's method, but using van der Corput sampling. The result obtained here is clearly superior to any of those in [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF].

We compared the performance of Glimm's method to that of two difference methods on a shock and rarefaction interaction problem (Fig. 20) like the one described in § 3, but with the waves an order of magnitude stronger: Pt = 473.9, Pm= 1.077, Pr = 100, PI= 23.27, Pm= 3.930, Pr = 100, Ut= 6.0, Urn= -4.0, u, = -1.181, X1 = .3, X,= .9, y=1.4.

The solution has the same qualitative features as those of the weaker problem, except that the backward facing compression wave produced by the shock-rarefaction interaction has itself steepened into a shock at the time the solutions are compared. Otherwise, the waves are all much stronger; in particular, the passively advected density wave is a spike, two zones in width for the ax = .01 cases. The two difference methods compared are the version of Godunov's method (Fig. 21) discussed in the previous section, and the MUSCL code written by Paul Woodward of Lawrence Livermore National Laboratory, based on the scheme of van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF] (Fig. 22). These two methods represent, respectively, one of the most accurate of the first-order methods, and a state-of-the-art representative of the adaptive or hybrid difference methods.

(For other examples, see Boris and Book [START_REF] Boris | Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works[END_REF], Harten [START_REF] Harten | The method of artificial compression: I. Shocks and contact discontinuities[END_REF], Harten and Zwas [START_REF] Harten | Self-adjusting hybrid schemes for shock computations[END_REF], Zalesak [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF].) As before, we compare all three results with the answer obtained using Glimm's method, van der Corput sampling with ax = .0025. All three methods obtain reasonably good answers for the pressure and velocity profiles, modulo the varying widths for the shock transition region. However, neither of the difference methods are able to get the correct peak value of the density, nor the correct width for the density spike; Glimm's method, by virtue of its direct simulation of the wave interaction process without averaging, gets the correct answer. The original proposal in [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF] for using Glimm's method with operator splitting was seen to give incorrect results for flows in which there occur large jumps in the pressure and velocity along surfaces oblique to the mesh directions. By coupling Glimm's method with Godunov's method, we lose many of the special properties of the Glimm's method with respect to its treatment of shock interactions. However, the resulting method has a number of attractive properties. Of all the first-order difference methods, Godunov's method produces the narrowest shocks (2-3 zones wide). Both Glimm's method and Godunov's method are extremely stable, even in the strongly nonlinear region (the problem with Glimm's method at shocks is a loss of accuracy, not of stability). Finally, Glimm's method has no numerical diffusion, so that the hybrid method has no numerical diffusion away from regions where large pressure gradients are generated. which have narrower shocks, in computing problems such as the Mach reflection problem discussed in the previous section. The number of timesteps required after the time of reflection for the Mach stem to form increases as does the width of the shock; consequently, the MUSCL code having shocks which are 1-2 zones wide obtains, on coarser meshes, results comparable to those obtained here. There are several directions in which further work is indicated. For onedimensional flows, Glimm's method with van der Corput sampling is quite effective in modelling the interaction of discontinuities with smooth parts of the flow, without introducing unacceptable errors in the latter. The fact that the solutions to the Riemann problem we use in the numerical scheme satisfy exactly the conservation laws is probably not essential to the accuracy of the method, since much of that information is lost in the sampling. What is essential is that the solution which is sampled has built into it the physically correct waves and wave speeds to some reasonable order of accuracy. Thus it is feasible to try to model with Glimm's method the dynamics of other media than an ideal gas in Cartesian coordinates: for example, gas dynamics with source terms or unusual equations of state, or elastic-plastic flow. The central advantage of the hybrid Glimm-Godunov method is that it has the simplicity and stability of a first-order method, with substantially less numerical diffusion than is usually seen in first-order methods. As the method is currently engineered, it seems to be more accurate than.the nonadaptive first-and second-order methods, but not as accurate as some of the adaptive methods. The main question to be answered is the determination of a set of optimum engineering decisions. One problem is that we have seen that the criterion for whether to use Glimm's method or Godunov's method at a point is different depending on the strength of the waves; the distinction between strong and weak waves should be made locally, by the algorithm. More generally, although the general principle for switching between the two methods is clear, the actual details of the procedure are still determined in a fairly ad hoc, problem-dependent fashion, and a more systematic algorithm is needed.

Appendix. Calculation of approximate solutions to the Riemann problem. In this appendix, we present a detailed description of the procedure used to calculate approximate values to the solution to the Riemann Problem in the two-dimensional hybrid Glimm-Godunov calculations described in § 2. The approximations introduced here Godunov method, C 0 = .05, k 0 = 1, u = .9.

are designed to give sufficiently accurate answers for the minimum computational effort in the two situations which arise in those calculations: 1) the sonic waves are weak, or 2) the sonic waves are strong, but the values calculated are used only for computing ftuxes in Godunov's method. The algorithm given here is also better suited for efficient implementation on a vector processor, such as the Cray-1, than those given previously. The first step is an iteration to calculate p*, the pressure of the gas between the two sonic waves. We use the Newton's method algorithm given in van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF], with one important modification: we assume that the formulae for WL,R, the mean Lagrangian wave speeds, are the same for both shocks and rarefactions: y-1co+c which does not require evaluation of a rational power.

W ( *) = C /1 + y+ 1 (p* -pL,R) L,R p L,R y
If the iteration is carried out to convergence, all the approximations introduced here are correct to third order in the pressure jump across the strongest rarefaction wave present. If one of the waves is a strong rarefaction, we will be using the results to calculate fluxes for Godunov's method, and the error committed by using the rarefaction shock formula in the iteration is lost in the averaging. However, it is essential to evaluate the solution inside the rarefaction fans, rather than treating the waves as jump discontinuities, as is done in [START_REF] Goounov | A comutational scheme for twodimensional nonstationary problems of gas dynamics and calculations of the flow from a shock wave approaching a stationary state[END_REF]. Otherwise, sampling would not spread weak rarefactions; nor would the averaging in Godunov's method spread strong rarefaction shocks into rarefaction waves, if the speed of the rarefaction shock is close to zero.

Using approximations similar to those used above in computing the solution to the Riemann problem, we extend the sampling procedure to the case where the time step satisfies (3 .1) for ! < u < 1, assuming that the sonic waves in the solution being sampled are weak. The procedure described below accounts correctly for the possible interpenetration of waves from successive Riemann problems to first order in the strengths of the sonic waves, and reduces to the previous sampling procedure when the waves do not intersect.

To update the solution at zone j, we first evaluate at ((j -!+an+ 1 ) ax, (n + 1) at) the solution to the Riemann problems on either side of the zone: If U~~1 ¥-Uj, Ui-1 and U':_~1 ¥-Uj, Uj+b then we assume that the waves contained in the jump ( U~~1 , Uj) have reached and interacted with the waves contained in the jump (Uj, U~~1 ). Using again the assumption that the sonic waves are weak, we see that the waves which intersect are: a part of a backward facing sonic wave, R,i-h u i+ 2 are t e mean agrangtan wave spee s an t e velocity of the gas between the two sonic waves for Riemann problems centered at ((j + !) ax, n dt).

The above scheme can be implemented in such a way that almost all the calculations are vectorizable. For the test problems discussed in § 4, a program run on the Cray-I at the LLNLCC, compiled using the CFT compiler, took about 14 JJ..S/zone/time step/space dimension, or about 35,000 zones/second for a two-dimensional problem. The iteration scheme for computing p*, u* is done once per zone, independent of whether one is sampling or averaging in that zone, and takes ( 1. 7 + 1.15 x la) f.J.-S = 4 JJ.-S for la= 2, or less than one third of the time per timestep. A good deal of redundant work is performed because of the limited number of vectorized logical operations available at the Fortran level using the CFT compiler. As more of the Cray-l's capabilities become accessible, such as vectorized gather/scatter operations and bitvector logic, the timing should improve substantially.

FIG. 1 .

 1 FIG.1. The Riemann problem for gas dynamics.

FIG. 2 .

 2 FIG.2. Solution at fixed time to shock tube problem (2.3).

uj+ 1 =

 1 u~ ( ( j-1 I 2+a n+ 1 ) 6x, ( n+ 1 ) 6 t)((j-1/2+an+ 1 )6x,(n+1)6t) ~ ---+------...----..--___.~-T---------......_--) ( X , t) = U~ ( X , t)

FIG. 5 .

 5 FIG. 5. Waves initially generated in lD test problem.

FIG. 6 .

 6 FIG. 6. Computed solution to lD test problem, van der Corput sampling, ~x = .0025.

  Pm= 1.39, Pr = 10.0, Pt = .6878, Pm= .146, Pr = .6, Ut= .0181, Um=-11.9, Ur = -5.98, Xt = .4, Xr = .9, y= 1.4. All calculations were done on the spatial interval [0, 1], with boundary conditions at 0 and 1 obtained by assuming the solutions satisfy au;axlx=o,l = 0. The various solutions being compared were computed with ax = .01, and are represented graphi- cally by circles for the computed values at mesh points, interpolated by a dotted line.

FIG. 7 .

 7 FIG.[START_REF] Concus | Numerical solution of a nonlinear hyperbolic equation by a random choice method[END_REF]. Computed solution to 1D test problem, random sampling, ~x = .01.

FIG. 8 .

 8 FIG.[START_REF] Friedrichs | Supersonic Flow and Shock Waves[END_REF]. Computed solution to lD test problem, (7, 3) stratified random sampling, ~x = .01.

FIG. 9 .

 9 FIG.[START_REF] Glimm | Solution in the large for nonlinear hyperlinear systems of equations[END_REF]. Computed solution to test problem, van der Corput sampling, ~x = .01.

FIG. 10 .

 10 FIG.[START_REF] Glimm | A numerical method for two phase flow with an unstable interface[END_REF]. Computed solution to test problem, van der Corput sampling, ax = fo.

PR= . 5 ,

 5 VN,L = -1.78, VN,R = 11.6, 'Y = 1.667.

FIG. 11 .

 11 FIG.[START_REF] Glimm | A shock tracking method for hyperbolic systems[END_REF]. Diagonal Reimann problem (3.2) computed using Glimm's method, u = .5.

FIG. 12 .

 12 FIG. 12. Shock tube problem (2.3) computed as diagonal Riemann problem using Glimm's method, u = .5.

FIG. 13 .

 13 FIG.[START_REF] Goounov | A comutational scheme for twodimensional nonstationary problems of gas dynamics and calculations of the flow from a shock wave approaching a stationary state[END_REF]. Diagonal Riemann problem (3.2) computed using hybrid Glimm-Godunov method, u=.S, C 0 =.1, k 0 =2.

FIG. 14 .

 14 FIG.[START_REF] Hammersley | Monte Carlo Methods[END_REF]. Shock tube problem computed as a diagonal Riemann problem using hybrid Glimm-Godunov, method u = .5, C 0 = .1, ko = 1.

Q

FIG. 16 .FIG. 17 .FIG. 18 .

 161718 FIG. 16.Wind tunnel problem at timet=~ computed using Godunov method, ~x = t~o, u = .9.

FIG. 19 .

 19 FIG.[START_REF] Morton | Difference Methods for Initial-Value Problems[END_REF]. Solution to shock tube problem (2.3) in one dimension, computed using Glimm's method with van der Corput sampling, ~x = .01.

FIG. 21 .

 21 FIG.[START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]. Solution to strong lD interaction problem, computed using Godunov's method.

FIG. 22 .

 22 FIG.[START_REF]A numerical method for a slightly viscous axisymmetric flow with application to internal combustion engines[END_REF]. Solution to strong lD interaction problem, computed using MUSCL.

  We evaluate p, p, u as follows:{ p*,p*,u* p, p, u = Po, Po, Uo if (f <A*, if (f > Ao.If Ao > (f >A*, then the solution is being evaluated inside a centered rarefaction wave, 1-u*+----, u=s(l/1-c), p=-po, p=-.One can replace the expression for p inside a rarefaction fan by the formula

un+. 1 i

 1 =h . .l (a n+1 ax) or i+1• Then we can assume, since the sonic waves are weak, that the waves from the Riemann problem at ((j-!) ax, n at) and ((j + !) ax, n at) do not intersect, and setU n+1 "f un+1 -un ( n + 1 { -,j l + ,j -j ,Ui )Giimm = U~~l Otherwise.

FIG. A 1 .

 1 FIG. A 1. Interaction of waves from adjacent Riemann problems.

Acknowledgments. Part of this work is contained in a Ph.D. dissertation written at the University of California, Berkeley. The author wishes to thank his thesis advisor, Alexandre Chorin, for his encouragement and support. The author also thanks G. Sod and P. Woodward for many helpful discussions, and the Mathematics Group at Lawrence Berkeley Laboratory and the Theoretical Physics Division at Lawrence Livermore National Laboratory for their generous financial support.

This work was performed under the auspices of the U .S. Department of Energy at the Lawrence Livermore National Laboratory and the Lawrence Berkeley Laboratory, and supported in part by the Director, Office of Basic Energy Sciences, Engineering, Mathematical, and Geosciences Division, under contract W-7405-ENG-48.

 -----. -------, -----, ------For the purpose of comparison with the results of Sod [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF], obtained by the various difference methods, we computed the shock tube problem (2.3) as a diagonal Riemann problem (Fig. 23), but on a lOO x 100 grid (ax =ay= .01); results for this problem compted using the MUSCL code are also given in [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF]. In principle, the problem solved here is more difficult than the one solved in [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF], since the latter is solved as a one-dimensional problem. But the answer is the same for both, and the results are worth comparing.

The calculation of the rarefaction, and the width of the shock transition in the results obtained with the hybrid Glimm-Godunov method, compare favorably with those obtained by any of the difference methods. The treatment of the contact discontinuity is clearly superior to that given by any of the difference methods in [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF], and comparable to that obtained in [START_REF] Van Leer | Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov's methods[END_REF]. The difference methods in [START_REF]A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] either spread the contact discontinuity over 6-10 zones, with the number of zones increasing as a function of time, or introduce substantial oscillations near the contact discontinuity.

The major weakness of the hybrid method is that it computes shocks which are 2-3 zones wide. This puts the method at a disadvantage compared to the methods Then the iteration proceeds as follows:

One criterion for terminating the iteration is to terminate if

where e > 0 is some predetermined tolerance, and set p* = p*' 1 • In programming this procedure for the Cray-1, we iterated a fixed number of times lo, i.e., set p* = p*' 10 , independent of the left and right states. We obtained more than adequate accuracy using / 0 = 4, for even the strongest problems; it appears to be sufficient, for a wide class of problems, to set lo to be 1 or 2.

Having obtained p*, we calculate the other quantities we will need:

u* =PL -pR +uLWL +uRWR WL+WR .

If, at any point in the iteration, p*• 1 < 0, we reset p* to be equal to some floor value 1 >> Pmin > 0. If p*' 1 < 0 for two iterations in a row, we terminate the iteration (or ignore the results), setting p* = Pmin•

The second part of the procedure is to calculate the value of the solution at some given point (x, t); we denote the values of the pressure, density, velocity, and passive component of the velocity at that point by p, p, u, v. We follow the procedure given in [START_REF] Chorin | Random choice solution of hyperbolic systems[END_REF], but use explicitly the reflection symmetry of the equations to consolidate some of the formulae.

Let t/1 = x / t and s = sgn ( t/1 -u *). Then we define ( ,,,

We then compute * _ ( -1 ( P * -Po)) - wave from the Riemann problem at ((j-~) ax, n at) and at most one contact discon- tinuity, which might come from either one of the Riemann problems (Fig. A1). In h 1 1 ( n + 1 n + 1 n + 1 n + 1) f II t atcase,wecacuate Pi ,pi ,ui ,vi

Glimmas o ows: ~((j-l/2+an+l)~x.(n+l)At)