# The Numerical Simulation of two-dimension Flow with Strong Shocks.pdf

JOURNAL OF COMPUTATIONAL PHYSICS 54, 115-173 (1984) Review Article The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks PAUL WOODWARD Lawrence Livermore National Laboratory, University of Calgornia, Livermore, California 94550 AND PHILLIP COLELLA Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720 Received July 15, 1982; revised August 25, 1983 Results of an extensive comparison of numerical methods for simulating hydrodynamics are presented and discussed. This study focuses on the simulation of fluid flows with strong shocks in two dimensions. By “strong shocks,” we here refer to shocks in which there is substantial entropy production. For the case of shocks in air, we therefore refer to Mach numbers of three and greater. For flows containing such strong shocks we find that a careful treatment of flow discontinuities is of greatest importance in obtaining accurate numerical results. Three approaches to treating discontinuities in the flow are discussed-artificial viscosity, blending of low- and high-order-accurate fluxes, and the use of nonlinear solutions to Riemann’s problem. The advantages and disadvantages of each approach are discussed and illustrated by computed results for three test problems. In this comparison we have focused our attention entirely upon the performance of schemes for differencing the hydrodynamic equations. We have regarded the nature of the grid upon which such differencing schemes are applied as an independent issue outside the scope of this work. Therefore we have restricted our study to the case of uniform, square computational zones in Cartesian coordinates. For simplicity we have further restricted our attention to two-dimensional difference schemes which are built out of symmetrized products of one-dimensional difference operators. I. 1NTRoDucT10~ Over the last 35 years, a great number of numerical schemes have been devised for the simulation of compressible gas dynamics on digital computers. A major difficulty which has enlivened research in this area is the problem of representing the shock and contact discontinuities which arise in these simulations. As early as 1950 a solution to this problem was proposed by von Neumann and Richtmyer [ 11. With minor modifications this first solution to the problem is still in general use today. However, 115 0021-9991/84 $3.00 Copyright 0 1984 by Academic Press, Inc. saLp4p 9 All rights of reproduction in any form reserved. 116 WOODWAR‘DAND COLELLA certain disadvantages of von Neumann and Richtmyer’s approach have led to continued efforts at more convenient, more accurate, and/or more elegant solutions. In this article we will discuss three different approaches to the representation of discontinuities which have evolved over the years. We will limit our consideration to approaches to “shock capturing,” in which shocks are smeared out on the grid somewhat so that they may be treated in a relatively simple and convenient fashion. Our goal is to present a comprehensive comparison of these approaches, so that their relative merits can be accurately assessed. The motivation for this work is the great need for accurate simulations of flows with strong shocks which exists in many fields of physics. Much experience indicates that the overall accuracy of such simulations is very closely related to the accuracy with which flow discontinuities are represented. Several algorithms have been proposed which perform well when applied to one-dimensional flow problems but which encounter major difficulties in two dimensions. Because two-dimensional calculations dominate present applications of gas dynamic simulation, we will focus on these more difficult simulations here. Recently techniques for adapting a computational grid to aid in the resoluion of flow discontinuities have received a great deal of attention. Although we will make some remarks about this subject, our main interest in this work is to compare the difference schemes which may be used on a particular grid and not the methods for choosing the grid itself. We regard the issues of choosing a difference scheme and of choosing a grid as largely independent, so that it makes sense to study them separately. We will also leave aside all issues relating to the treatment of the geometric source terms which arise when curvilinear coordinates are employed. Thus we will work exclusively in Cartesian coordinates, and we will use only uniform grids. In the same spirit we will generate all our two-dimensional difference schemes by forming symmetrized products of one-dimensional difference operators. All flows computed here will also employ an especially simple equation of state, namely, that of a gamma-law gas with gamma equal to 1.4. Even with all these restrictions, many computational difficulties remain, and there is a wide spread in the performance of the various methods considered below. II. THREE MAJOR APPROACHES TO THE REPRESENTATION OF DISCONTINUITIES a. Artljkial Viscosity This earliest and most commonly used approach to representing discontinuities was originally suggested by von Neumann and Richtmyer [ 11. In the flow of real fluids we believe that there are no discontinuities. There are instead very thin regions of very steep gradients. If terms representing viscosity and heat conduction were included in the usual equations of hydrodynamics, the discontinuities in the solutions of these equations would no longer develop. However, we would have to resolve smooth structures on the very small distance scales characteristic of viscous NUMERICAL SIMULATION OF FLUID FLOW 117 momentum transport and molecular heat conduction. Without changing the flow very much we can increase the physical coefficients of viscosity and heat condution so that discontinuities are spread over distance scales which are negligible but still resolvable on a practical computational mesh. von Neumann and Richtmyer suggested that an artificial viscous pressure be used to smear shocks over a few computational cells, or zones. Later Lapidus [45] suggested that terms describing diffusion of mass, momentum, and energy be used for the same purpose. In both methods the terms added to the differential equations cause dissipation of kinetic energy into heat. b. Linear Hybridization In this approach the results of two difference schemes are blended together. A high- order difference scheme which is very accurate in smooth flow but badly behaved at discontinuities is blended with a low-order scheme. The low-order scheme should have sufficiently large dissipative truncation errors that it yields monotone profiles for flow discontinuities. In smooth flow, the high-order scheme is used, but near a discontinuity the low-order scheme is blended with it to an extent sufficient to guarantee monotone representations of the jumps at the discontinuity. We will refer to this criterion for determining the blending weight factor as a monotonicity constraint. In order to preserve the exact conservation of mass, momentum, and energy by the composite scheme when the blending factors vary over the grid, it is the fluxes of these conserved quantities at zone interfaces as computed by the two schemes which are blended. These combined fluxes are then differenced in order to update the mass, momentum, and energy of each zone. This approach of linear hybridization is similar to the artificial viscosity approach in one respect. To a high-order method which oscillates near discontinuities some terms are added which are negligible in smooth flow and which are strongly dissipative near discontinuities. However, these terms are not motivated by analogy with a more realistic physical model for the flow than is given by the inviscid flow equations. Instead, the added terms are designed specifically to give the sharpest possible discontinuity profiles which are also monotone. They are able to perform this task better than the artificial viscosity terms; but under certain circumstances the monotonicity constraint proves not to be the physically appropriate condition, and then difficulties can arise. c. Godunou’s Approach A third means of treating discontinuities, suggested by Godunov [2], is to introduce explicit nonlinearity into the difference method. In the two other approaches, difference schemes are derived from Taylor series expansions of the terms in the differential equations. This technique is fundamentally based upon the assumption that the solution is smooth. At a discontinuity this assumption is inap- propriate; hence the need to force a well-behaved solution by introducing an 118 WOODWARDAND COLELLA unphysically large viscosity or an unphysical constraint of monotonicity. In 1959 Godunov pioneered a new approach to this problem [2]. Instead of building up a full solution to the hydrodynamic equations by piecing together smooth, small-amplitude solutions, he built up a solution by piecing together discontinuous solutions. These discontinuous solutions closely approximate the smooth ones where those are appror- piate, but they have the great additional advantage of approximating the true solution reasonably well even when that solution is not smooth. Godunov made use of a nonlinear flow problem which is simple enough to permit exact solution-Riemann’s shock tube problem. This simple solution describes the nonlinear flow which develops from a discontinuous jump separating two constant states. In general, the jump develops into two nonlinear waves, either shocks or rarefactions, with a contact discontinuity in between (cf. Ref. [3]). Godunov’s approach was then to approximate a hydrodynamic flow by a large number of constant states, compute their interactions exactly, and average the results in a conservative fashion. This procedure leads to an accurate and very well-behaved treatment of shock discontinuities. Godunov’s approach has also been extended to higher-order methods [4-91. In all of these extensions, the Riemann solver is an essential element in allowing narrow discontinuities without unphysical oscillations. A different way of using Riemann problems was developed by Glimm, Chorin, and others [ 10-151. This method, the random choice method, represents the flow within each zone by the detailed solution to the Reimann problem sampled at a represen- tative point within the zone. Diffusive errors from spatial averaging are avoided, but errors in the time dimension take their place. Unfortunately the very desirable properties of this method in one dimension do not persist in its multidimensional formulations to date (cf. [ 13]), so that we will not discuss the method further here. III. ADVANTAGES AND DISADVANTAGES OF THE THREE APPROACHES a. Artificial Viscosity The major advantage of this approach is its simplicity. If a large enough artificial viscosity is added to the equations of hydrodynamics, their solutions will always be smooth and the design of effective difference methods can be done in a systematic, straightforward way. This procedure is especially successful when the artificial viscous effects are important on distance scales of a fixed length. When practical considerations force this length scale to vary with the size of the computational zones of a strongly nonuniform grid, large unphysical effects can result. This difficulty, which was first noticed in 1968 by Axelrod [ 161, has been discussed at length in Ref. [ 171 and will receive no further attention here. The straightforward use of an artificial viscosity method on practical problems involves a trade-off of computation time and computer memory resources for program simplicity. Because shocks are smeared out over three or more zones by NUMERICAL SIMULATION OF FLUID FLOW 119 these methods, fine grids are necessary if the grids are to be uniform to allow for the possibility of shock motion to any point in the problem domain. This necessity for fine grids is caused by the direct relationship between the accuracy of the computation and the thickness of the shocks in many problems of practical interest. The results we will show in a later section demonstrate that if a fully converged solution is desired, this straightforward approach is impractical. Even for the relatively simple 2-D test problems presented below full convergence requires 400,000 words of memory and rather large amounts of computer time. A more practical approach toward the use of these methods has been suggested by Oliger, and has been developed by Gropp [ 181, Bolstad [ 191, and Berger and Oliger [20]. In this approach the mesh is locally refined in both space and time in the neighborhood of a discontinuity. Of course, this idea of local grid refinement is not entirely new, but Oliger and his collaborators have carried it out in a very general and systematic fashion. Local mesh relinement saves the computer memory required for a global fine grid at the cost of introducing substantial program complexity. The computer time needed to achieve a given level of accuracy can also be reduced, but savings here are not as great as one might expect. Considerable time is required to decide where and by how much the grid should be refined. To prevent sudden jumps in the solution at the locations of discontinuous changes in the mesh width, the regions of mesh refinement must be wider than one would suspect at first. Finally, the comutation on the refined mesh regions must be specially organized for vector computers. This special organization involves both extra computer time and program complexity. A discussion of these issues can be found in [2 11. A detailed assessment of the costs and benefits of the local mesh refinement approach is beyond the scope of this article. An investigation of the usefulness of this approach when combined with some of the difference schemes discussed in this article is presently under way. An alternative approach is the use of a continuously adaptive mesh. In this approach the mesh is continuously deformed so that structures in the flow which have disparate characteristic length scales are all resolved most effectively. This sort of method is indispensible in problems which contain very thin features whose internal structures must be computed in detail. The method was pioneered in astrophysical and combustion physics simulations [22-291, in which very thin features whose internal structures are essential to the problem must be carefully tracked. In the work of Winkler [24-261 an additional equation relating the distribution of grid points to the flow solution is solved implicitly together with the flow equations. This is a very effective way to deal with the extreme grid distortions forced on the method in order to resolve very sharp features in the flow. Similar approaches have been developed by Dwyer et al. [27], and later by Miller and Miller [28], and Gelinas et al. [29]. In this last work the mesh equations are developed so that an estimate of the error