Diffusion regulation for Euler solvers

Journal of Computational Physics 221 (2007) 577–599 www.elsevier.com/locate/jcp Diffusion regulation for Euler solvers S. Jaisankar, S.V. Raghurama Rao * CFD Centre, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, Karnataka 560012, India Received 10 March 2006; received in revised form 17 June 2006; accepted 20 June 2006 Available online 14 August 2006 Abstract A diffusion regulation parameter, which operates based on the jump in the Mach number, is presented for implementation in Euler solvers. This diffusion regulation parameter adjusts itself automatically in different regimes of the flow and leads to the exact capturing of steady contact discontinuities which are aligned with the grid-lines. This diffusion regulator parameter reduces numerical dissipation, is very simple and can be easily incorporated in any Euler solver. By coupling such a parameter with a simple numerical method like the Local Lax-Friedrichs (Rusanov) method, an accurate and yet simple numerical method is developed for the numerical simulation of inviscid compressible fluid flows. To demonstrate the applicability of this approach to any Euler solver, the diffusion regulation parameter is also applied in the framework of a Kinetic Scheme which is very diffusive and the improvements in the accuracy for both the methods are demonstrated through several bench-mark test problems.  2006 Elsevier Inc. All rights reserved. Keywords: Diffusion regulation parameter; Euler solvers; Exact capturing of contact discontinuities 1. Introduction The motivation for this work is to develop a numerical method for solving Euler equations, capturing discontinuities accurately without oscillations and with algorithmic simplicity. Of the upwind Euler solvers developed in the last few decades, the approximate Riemann solvers of Roe [1] and Osher and Solomon [2], the advective upstream splitting method (AUSM) and its several variants of Liou [3–5] and the convective upwind split pressure (CUSP) scheme of Jameson [6,7] are the most popular methods which can capture discontinuities with very good accuracy. Apart from these popular methods, the flux vector splitting methods of Steger and Warming [8] and van Leer [9], the Kinetic Schemes [10] and Relaxation Schemes [11] also provide interesting alternatives for solving compressible fluid flow equations. However, these later numerical methods are known to be overly diffusive in capturing discontinuities. Recently, there have been attempts to reduce the amount of numerical diffusion in Relaxation Schemes [12] and Kinetic Schemes [13–15]. Here, we present a new and simple approach to capture discontinuities accurately and demonstrate it with a simple Euler solver, the Local Lax-Friedrichs (Rusanov) method [16,17] and also with a Kinetic Scheme, the peculiar velocity * Corresponding author. Tel.: +91 80 2293 3031; fax: +91 80 2360 0134. E-mail address: raghu@aero.iisc.ernet.in (S.V. Raghurama Rao). 0021-9991/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2006.06.030 578 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 based upwind (PVU) method [18–21]. The coupling of the diffusion regulation parameter makes these otherwise diffusive schemes very accurate and leads to exact capturing of steady contact discontinuities when the discontinuities are aligned with the grid lines. Our approach of diffusion regulation is applicable to any Euler solver. We demonstrate its usefulness by applying it also to Roe’s method and Osher’s scheme and even these less diffusive schemes give improved results when coupled with the diffusion regulation parameter for multidimensional flows. Several bench-mark problems are solved with this approach and the results demonstrate the efficiency of this approach. 2. Diffusion regulator (DR) model Consider one-dimensional system of conservation laws for inviscid fluid flow oU oG þ ¼0 ot ox ð1Þ where U = (q, qu, qE)T is the vector of conserved variables and G = (qu, qu2 + p, quE + pu)T is the inviscid flux vector. An explicit scheme in conservation form on a three point stencil and with piecewise constant approximation for conserved variable is given by U jnþ1 ¼ U nj  Dt ðG 1  Gj12 Þ Dx jþ2 ð2Þ where Dt and Dx are respectively time and space steps, j represents the cell centroid and j  12 represent the cell interfaces. In the finite-volume formulation (2), the differences among all the numerical schemes lie essentially in the definition of numerical flux GI ¼ Gjþ1 evaluated at the cell interface. The upwind flux splitting schemes 2 define the interface flux in terms of split fluxes of left and right states as  GI ¼ Gþ L þ GR ð3Þ 1 GI ¼ ðGL þ GR Þ  DI 2 ð4Þ The interface flux in any scheme can be written in the form where L and R represent the left and right states, the first term on the right hand side (average flux) represents the central discretization of the flux terms and DI represents the flux corresponding to numerical diffusion. The basic idea of our diffusion regulator (DR) model [22,23] is to regulate the dissipative flux DI using a parameter U (DR parameter) so that optimal amount of dissipation will be added in the different flow regimes, for capturing discontinuities accurately without oscillations and to avoid unphysical expansion shocks. The parameter U is introduced to multiply the diffusive flux in the definition of the interface flux, as given below. 1 GI ¼ ðGL þ GR Þ  UDI ð5Þ 2 Note that the above format of the interface flux being a sum of an average flux (which leads to central discretization) and a dissipative flux is similar to the format used in the JST scheme [24], as a central discretization based scheme stabilized with an artificial dissipation term which is adapted to the local properties of flow. JST scheme, known for its efficiency, uses higher order Runge–Kutta time stepping and has a dissipation term which is a blend of second and fourth differences along with pressure gradient based shock sensor. Here, we regulate the dissipation term in a different manner, as described in the following section. 2.1. Diffusion regulator (DR) parameter The solution of hyperbolic conservation laws in a stable manner requires sufficient artificial viscosity, either as a direct addition to the averaged flux (which leads to central discretization) or as the implicitly present diffusive flux in an upwind based formulation. There are at least three requirements which the numerical viscosity should satisfy, for good shock capturing. Firstly, in most of the accurate schemes, including the methods of Roe or Osher, the dissipation flux DI is of the form aDU and assumes small values away from shock S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 579 discontinuities as DU is negligible. Therefore, a minimal dissipation for a stable solution is sufficient in such flow regimes. Secondly, near to a shock but not at the shock discontinuity itself, the dissipation must be sufficient enough and do the job of diminishing the influence of left state on right state and vice versa, for a greater accuracy. But, most schemes for vector conservation laws do not completely eliminate the above influence and hence shock is always smeared due to the contribution from both sides of the shock. Thirdly, the numerical dissipation at the discontinuities should be zero to have an exact capture. The diffusion regulator parameter described here is modeled to adopt the first two roles to a greater extent, thereby leading to a low dissipative and yet robust scheme. In addition, the numerical dissipation is also driven to zero at a steady contact discontinuity, with use of an exponential function of an average Mach number at the cell interfaces, leading to exact capturing of the steady contact discontinuities. We use a simple expression for the diffusion regulator parameter ‘U’ as a function of jump in the magnitude of Mach number across the interface. The maximum value of ‘U’ is restricted to unity so that maximum value of the numerical dissipation corresponds to the numerical dissipation present in the original scheme chosen for diffusion regulation. In an upwind scheme, U = 1 corresponds to the full upwind flux. The discrete approximation to spatial derivative takes values ranging from a central difference approximation (discounting the residual dissipation) to a full dissipative approximation allowed by the original scheme chosen, for example, an upwind scheme. The diffusion regulator parameter is expressed mathematically as (Ma is the arithmetic mean of left and right Mach numbers) ðDM 2 þ d2 Þ ð1  ejM a Þ when jDMj 6 d 2d U ¼ jDMj when d < jDMj 6 1:0 U ¼ 1:0 when jDMj > 1:0 U¼ ð6Þ ð7Þ ð8Þ DM ¼ M L  M R ð9Þ The diffusion regulation parameter without the exponential term is sketched in Fig. 1. The role of the exponential term, which is a function of the average Mach number, is to drive the numerical diffusion to zero when Mach numbers become small so that steady contact discontinuities are captured accurately. The parameter j ensures a steep and still smooth variation of U in the very low Mach number region. To have the steep exponential variation, a high value of j is needed as the Mach number is reduced. Empirically, for all the supersonic test cases, j = 10 is found to work very well and for subsonic cases, higher value of j is needed as Mach number keeps decreasing. Numerical dissipation in the vicinity of a shock wave is sufficiently high and takes its maximum value (U = 1) for capturing strong shocks. The dissipative flux varies smoothly when moving from smooth solution to discontinuities. A minimum residual dissipation is retained, except when the average Mach number is close to zero, to avoid expansion shocks and to prevent the violation of the entropy condition. The DR parameter is modeled as a function of |DM| and takes values from d/2 to 1.0 as follows. Φ 1.0 δ –1.0 0 δ 1.0 ΔΜ Fig. 1. Diffusion regulator parameter without the exponential term. 580 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599  Case (1) – At steady contact surfaces. Magnitude of both ML and MR will be close to zero. Ma, which is a arithmetic mean of left and right Mach numbers, is zero and hence U = 0.0.  Case (2) – At shocks. ML and MR will be of the same sign. |DM| will vary based on whether shock is weak or strong. For relatively weak shocks when |DM| 6 d, U varies between d/2 and d. For stronger shocks, as |DM| P d, U varies from d to 1.0.  Case (3) – In smooth regions. Generally, |DM| 6 d, 0.0 6 U 6 d. The value of d is taken as 0.5 and is found to work in all the test cases presented in this study, for both internal and external flows. The diffusion regulator parameter of the model takes a low value far away from a shock discontinuity but picks up gradually in the shock proximity. This helps in decreasing the influence of the right state (of the shock) on the left state and vice versa, thus leading to a crisper shock. Since the dissipation is non-zero near the shock, the segregation of the information from left and right sides of a shock wave is not completely achieved, but gets optimized and the points which retain influence from both the sides become few in number. The final effect is that a better solution emerges than the parent scheme. In the smooth regions, including the expansion zones, sufficient artificial viscosity as a function of DM(6d) is maintained. At stagnation sonic points U becomes zero and hence grid aligned contact discontinuities are captured exactly. This algorithm is implemented on one of the simplest algorithms, the Local Lax-Fredrichs (Rusanov) method, to obtain a low dissipative scheme and is presented in the following section. 3. Local Lax-Friedrichs (Rusanov) method with diffusion regulation A low dissipative and yet simple numerical method, the Diffusion Regulated Local Lax-Friedrichs (DRLLF) scheme, is presented here. The dissipation term in Rusanov’s scheme, also known as Local Lax-Fredrichs method [16,17], is written as a function of local maxima of the eigenvalues of the flux Jacobian matrix and the state variable difference as 1 DI ¼ jkjmax ðU R  U L Þ 2 Thus, the interface flux in LLF method is expressed as 1 GI ¼ ðGL þ GR  jkjmax ðU R  U L ÞÞ 2 ð10Þ ð11Þ We modify the diffusive flux by multiplying it with the diffusion regulation parameter and obtain the interface flux as 1 GI ¼ ðGL þ GR  Ujkjmax ðU R  U L ÞÞ ð12Þ 2 The results of the standard 1D shock tube test case [25] (for the steady contact discontinuity and Sod test cases) are presented in Fig. 2 and for the quasi-1D nozzle flow test case [26] in Fig. 3. It can be observed that there is a substantial improvement in the resolution of shocks in addition to exact resolution of the steady contact discontinuity. The shocks with DR-LLF scheme are captured with fewer number of points than with the LLF scheme. 3.1. Higher order extension The concept of uniformly accurate, arbitrarily high order, essentially non-oscillatory (UNO) type upwind interpolation of Chakravarthy et al. [27] is adopted with the LLF method to obtain a higher order accurate scheme. A form of linear reconstruction of split fluxes is used to obtain the interface flux and hence will be influenced by four points (two each on left and right sides of the interface). Thus, the flux at the right interface on an equally spaced five point stencil is written as Gjþ12 ¼ Gþ j þ Dx Dx  ðGþ  Gþ ðG  G j1 Þ þ Gjþ1 þ jþ2 Þ 2Dx j 2Dx jþ1 ð13Þ 581 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 1 EXACT LLF DR–LLF 1 LLF DR–LLF 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 (a) 0– –10 –8 –6 –4 –2 0 2 4 6 8 10 (b) –10 0 10 Fig. 2. Comparison of DR-LLF and LLF methods: (a) 1D steady contact discontinuity in a shock tube; and (b) 1D Sod’s shock tube results on 100 grid points. 1.4 1.4 (a) (b) EXACT LLF DR–LLF 1.2 1.2 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0 0.5 1 1.5 2 2.5 3 0.2 EXACT LLF DR–LLF 0 0.5 1 1.5 2 2.5 3 Fig. 3. Comparison of DR-LLF and LLF methods for quasi-1D Flow through a converging–diverging nozzle with 100 points: (a) first order and (b) second order solutions. The split fluxes are calculated by rewriting first order scheme (Eq. (12)) as 1  Gjþ12 ¼ ðGj þ Gjþ1  Ujkjmax ðU jþ1  U j ÞÞ ¼ Gþ j þ Gjþ1 2 ð14Þ and hence 1 1 Gþ and G ð15Þ j ¼ ðGj þ Ujkjmax U j Þ jþ1 ¼ ðGjþ1  Ujkjmax U jþ1 Þ 2 2 DR parameter ‘U’ is computed using the stencil maxima of the jump in Mach number (which is the difference of the minimum and maximum Mach numbers in the stencil). In order to suppress the wiggles which will appear in higher order schemes, the split fluxes in Eq. (13) are modified by a minmod limiter ‘g’ [23] as Gjþ12 ¼ Gþ j þ gða; bÞ ¼ ( Dx Dx þ  þ  þ   g½ðGþ g½ðG j  Gj1 Þ; aðGjþ1  Gj Þ þ Gjþ1 þ jþ1  Gjþ2 Þ; aðGj  Gjþ1 Þ 2Dx 2Dx a if jaj < jbj b if jbj < jaj ð16Þ ð17Þ 582 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 3.2. Results and discussions The proposed low dissipative scheme is tested on the following standard 2D bench-mark problems: (i) gridaligned shear flow [28]; (ii) shock reflection problem [29]; (iii) supersonic flow over a compression ramp of 15 in a channel [30]; (iv) supersonic flow over a step in wind tunnel [31]; (v) supersonic [32], transonic [32] and subsonic [33] flows over NACA0012 airfoil; and (vi) hypersonic flow over a semi-cylinder [34,35]. The diffusion regulator parameter for multi-dimensional problems is computed based on the normal Mach numbers across the finite volume interface. The results obtained with first order accurate and second order accurate DR-LLF method are compared with those obtained with first order accurate and second order accurate LLF Scheme in Figs. 4–9. While the first order accurate LLF method is highly diffusive, the first order accurate DR-LLF method captures the shocks with much less numerical dissipation and is uniformly more accurate. In comparison, the second order accurate LLF method captures the shock solutions more accurately than the first order accurate LLF method, but in some cases the solutions are still found to be more dissipative than the first order DR-LLF scheme. Whenever a contact discontinuity (or shear layer) is aligned with the grid, DR-LLF scheme captures it exactly, without any numerical dissipation (Fig. 4). The accurate capturing of slip surfaces, bow shock, Mach stem and expansion corner for the step-in-a-wind-tunnel problem (Fig. 9) demonstrates the higher accuracy of the DR-LLF method, compared to the LLF method. With the original LLF method, the slip surface at the downstream is not at all seen and the position of the beginning of the Lambda-shock is not right above the beginning of the step. These drawbacks are corrected by the diffusion regulation parameter in the DR-LLF scheme. The zero normal velocity across a steady contact discontinuity results in a zero normal mean Mach number and hence the numerical dissipation is completely turned off. Thus, the grid-aligned contact discontinuities are captured exactly. However, if the contact discontinuity is not aligned to the finite volume interface, then the non-zero normal mean Mach number can result in dissipative solutions of contact discontinuities, depending on their deviation from grid lines. This problem is basically due to the limitation of the finite volume method as only the fluxes normal to the cell-interfaces are considered. This problem can be cured only by genuinely multi-dimensional models, which is beyond the scope of the present work. However, the numerical dissipation in capturing discontinuities in multi-dimensional cases is still expected to be low with a diffusion regulated scheme compared to a scheme without such diffusion regulation. Solutions obtained for external supersonic, transonic and subsonic flows over NACA0012 airfoil presented in Figs. 10 and 11 corroborate the efficacy of the new diffusion regulation model in capturing discontinuities with high accuracy. The modified scheme brings the bow shock and fish tail shocks in supersonic flow into prominence which are hardly seen in the LLF solution. Pressure contours for transonic flow case show the neatly resolved body shocks which are unseen in the LLF solution. The solution for external hypersonic flow (Fig. 12) around a semi-cylinder presented for Mach number 6.0 [3] shows a dissipation reduction with DRLLF method. The result obtained for a tangentially stretched grid for a high free stream Mach number of 20.0 (typical test for carbuncle phenomenon [34,35]) is carbuncle free with the low dissipative scheme. Fig. 13 shows the higher values assumed by diffusion regulator parameter closer to the shock and a minimum value LLF DR–LLF(Β) DR–LLF(A) 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0.5 1 0 0 0.5 1 0 0 0.5 1 Fig. 4. Comparison of first order DR-LLF and LLF schemes: Mach number contours (2.0:0.1:3.0) for grid aligned shear flow with 40 · 40 grid; regulator parameter (A) without exponential term (B) with exponential term. 583 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 LLF(0.9:0.05:2.9) 1 DR–LLF(0.9:0.05:2.9) 1 (a) 0.5 0 –1 1 0.5 –0.5 0 0.5 1 1.5 2 (b) 0.5 0 –1 1 0 –1 1 –0.5 0 0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 (b) 0.5 –0.5 0 0.5 1 1.5 2 (c) 0.5 0 –1 (a) 0 –1 1 –0.5 0 (c) 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 0 Fig. 5. Comparison of first order DR-LLF and LLF schemes: pressure contours (0.7:0.1:2.9) for oblique shock reflection with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. LLF(0.9:0.05:2.9) 1 DR–LLF(0.9:0.05:2.9) 1 (a) 0.5 0 –1 1 0.5 –0.5 0 0.5 1 1.5 2 (b) 0.5 0 –1 1 0 –1 1 –0.5 0 0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 (b) 0.5 –0.5 0 0.5 1 1.5 2 (c) 0.5 0 –1 (a) 0 –1 1 –0.5 0 (c) 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 0 Fig. 6. Comparison of second order DR-LLF and LLF schemes: pressure contours (0.7:0.1:2.9) for oblique shock reflection with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. ‘d’ near to the solid body. The local dissipation correction (see section on diffusion regulator parameter) in the small zone upstream to the shock, which is explained as a cure for the carbuncle phenomenon [34,36], may be naturally available with this regulator although this understanding needs a detailed study. The convergence rate for the new scheme is found comparable to that of the LLF scheme (Figs. 14 and 15) up to medium sized meshes but slows down for very fine mesh. Thus, the DR model allows numerical viscosity to be computed according to flow conditions and its minimal value of as low as a quarter of the upwind dissipation is sufficient to maintain a workable scheme. 584 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 LLF(1.0:0.05:3.2) DR–LLF(1.0:0.05:3.2) 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 (b) 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 1 (c) (c) 0.5 0 –1 0.5 (b) 0.5 0 –1 0 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 Fig. 7. Comparison of first order DR-LLF and LLF schemes: pressure contours (1.1:0.05:3.8) for ramp in channel flow with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. LLF(1.0:0.05:3.2) DR–LLF(1.0:0.05:3.2) 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 (b) 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 1 (c) (c) 0.5 0 –1 0 (b) 0.5 0 –1 –0.5 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 Fig. 8. Comparison of second order DR-LLF and LLF schemes: pressure contours (1.1:0.05:3.8) for ramp in channel flow with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. 4. Diffusion regulation in a Kinetic Scheme A low diffusion Kinetic Scheme, the diffusion regulated peculiar velocity based upwind (DR-PVU) method, is presented in this section, using the concept of diffusion regulation mentioned in the previous sections. The basic Kinetic Scheme, peculiar velocity based upwind (PVU) method, is based on the strategy of introducing an upwind scheme for the Boltzmann equation of the kinetic theory of gases and then taking moments to obtain an upwind scheme for the Euler equations. The relative velocity of the molecules with respect to the fluid velocity, called as the peculiar velocity, is used for introducing upwinding. The details of 585 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 LLF(1.0:0.15:6.5) DR–LLF(1.0:0.15:6.5) 1 1 (a) (a) 0.5 0.5 0 0 1 1 (b) (b) 0.5 0.5 0 0 1 1 (c) (c) 0.5 0.5 0 0 1 1 (d) (d) 0.5 0 0.5 0 0.5 1 1.5 2 2.5 3 0 0 0.5 1 1.5 2 2.5 3 Fig. 9. Comparison of first order DR-LLF and LLF schemes: density contours (1.0:0.15:6.5) for step in wind tunnel at (a) t = 2, (b) t = 3, (c) t = 4 with grid Dx = Dy = 1/40 and (d) t = 4 with grid Dx = Dy = 1/80. LLF(0.4:0.05:1.8) DR–LLF(0.4:0.05:1.8) 1.5 1.5 1 1 0.5 0.5 0 0 –0.5 –0.5 –1 –1 –1.5 –1 0 1 2 –1.5 –1 0 1 2 Fig. 10. Comparison of first order DR-LLF and LLF schemes for transonic flow across NACA0012 airfoil at Mach number 0.85 and angle of attack 1: pressure contours (0.4:0.05:1.8). the parent scheme are available in [18–21]. The Kinetic Schemes, like the flux vector splitting methods, are known to be robust but overly diffusive, compared to the Riemann solver of Roe. Our motivation here is to introduce low dissipation in the above mentioned Kinetic Scheme, using the diffusion regulation parameter. It is convenient to use the modified courant splitting of Dominic et al. [15] for upwinding of the peculiar velocity c = v  u, where v and u are the molecular and macroscopic fluid velocities, respectively. The new method can also be formulated in a different way directly by modifying the diffusion added to central discretization part of the interface fluxes. Since the introduction of peculiar velocity leads to the separation of the convective (transport) and pressure (acoustic) fluxes, each of these can be considered for upwinding separately, as shown below. 586 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 LLF(0.4:0.05:1.8) DR–LLF(0.4:0.05:1.8) 1.5 1.5 1 1 0.5 0.5 0 0 –0.5 –0.5 –1 –1 (a) –1.5 –1 0 1 (b) –1.5 –1 2 0 1 1 0.5 0.5 0 0 –0.5 –0.5 –1 –1 –1.5 –1.5 (c) –2 0 1 2 DR–LLF LLF 0.5 (d) –2 0 1 0.5 1 Fig. 11. Comparison of first order DR-LLF and LLF schemes for flow across NACA0012 airfoil: (a, b) pressure contours (0.4:0.05:2.8) at Mach number 1.2 and angle of attack 0; (c, d) Cp distribution at Mach number 0.63 and angle of attack 2. LLF DR–LLF 1.5 LLF 1.5 (a) DR–LLF 1.5 (a) 1.5 (b) (b) 1 1 1 1 0.5 0.5 0.5 0.5 0 0 0 0 –0.5 –0.5 –0.5 –0.5 –1 –1 –1 –1 –1.5 –1 –0.5 0 –1.5 –1 –0.5 0 –1.5 –1 –0.5 0 –1.5 –1 –0.5 0 Fig. 12. Comparison of first order DR-LLF and LLF schemes: density contours (2.0:0.2:5.0) for hypersonic flow across a semi-cylinder: (a) for Mach 6.0 with grid 45 · 45; and (b) for Mach 20.0 with tangentially stretched grid 320 · 20. 587 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 Φ(0.0:0.1:1.0) 1.5 0.6 1 1.0 0.8 0.5 0.4 0 0.6 0.2 1.0 0.4 0.2 –0.5 0.8 0.6 1.0 –1 0.8 –1.5 –1.5 –1 –0.5 0 0.5 Fig. 13. Contours of diffusion regulator parameter (0:0.2:1.0) for hypersonic flow across a semi-cylinder for Mach 20.0 with tangentially stretched grid 320 · 20. 0 0 10 10 –2 (c) 10 –5 10 (b) (a) –4 10 (c) –10 10 –6 (a) 10 (c) (b) (a) –8 10 (b) –15 0 0.5 1 1.5 2 2.5 4 x 10 10 0 0.5 1 1.5 2 2.5 3 4 x 10 Fig. 14. Convergence history of (left): first (- Æ - Æ -) and second (––) order accurate DR-LLF scheme for ramp in channel flow and (right): comparison of LLF (- Æ - Æ -) and DR-LLF (––) methods, with grid size: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. 10 10 10 LLF 0 10 –1 10 –2 10 (c) DR–LLF 0 –1 –2 (c) (b) 10 –3 10 (a) –3 (b) 10 –4 10 0 1 2 3 4 x 10 4 (a) –4 0 1 2 3 4 x 10 4 Fig. 15. Comparison of convergence history of second order solutions for ramp in channel flow (left): LLF scheme and (right): DR-LLF scheme with grid size: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. 588 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 4.1. Upwinding of the u-part The splitting for upwinding the transport flux oGt o o ¼ hWuþ F i þ hWu F i ox ox ot ð18Þ is reformulated as uþ ¼ u þ /t juj 2 u ¼ and u  /t juj 2 ð19Þ instead of uþ ¼ u þ juj 2 and u ¼ u  juj 2 ð20Þ Thus oGt oGtþ oGt ¼ þ ox ox ox ð21Þ where Gt ¼ u U ð22Þ 4.2. Upwinding of the c-part The splitting of the peculiar velocity (c) in the acoustic flux term oGa o o ¼ hWcþ F i þ hWc F i ox ox ot ð23Þ is reformulated as where cþ ¼ c þ /a jcj 2 and c ¼ c  /a jcj 2 ð24Þ Thus oGa oGaþ oGa ¼ þ ox ox ox where Ga ¼ hWc F i c jcj hWcþ F i ¼ hW F i þ /a hW F i 2 Z 12 Z 1 c cF hW F i ¼ dI dcW 2 2 0 1 Z 1 Z 1 jcj cF dI dcW hW F i ¼ 2 2 0 0 Therefore, the expression for split transport and acoustic fluxes becomes 3 2 /q  paffiffiffiffiffi 7 6 2 pb 7 6 7 6 / qu p t a  a 7 6 p ffi ffi ffi ffi ffi  G ¼ u U; G ¼ 6 7 2 2 pb 6 n o7 5 4 pu /q p  paffiffiffiffiffi þE 2 2 pb 2q ð25Þ ð26Þ ð27Þ ð28Þ ð29Þ ð30Þ 589 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 EXACT PVU 1 PVU DR–PVU 1 DR–PVU 0.8 0.6 0.5 0.4 0.2 (a) (b) 0 –10 0 10 –10 0 10 Fig. 16. Comparison of DR-PVU and PVU methods: (a) 1D steady contact discontinuity in a shock tube; and (b) 1D Sod’s shock tube results on 100 grid points. 1.5 1.5 EXACT EXACT PVU PVU DR–PVU DR–PVU 1 1 (b) (a) 0.5 0.5 0 0 0 1 2 3 0 1 2 3 Fig. 17. Comparison of DR-PVU and PVU methods for quasi-1D Flow through a converging–diverging nozzle with 100 points: (a) first order and (b) second order solutions. PVU DR–PVU 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.5 1 0 0 0.5 1 Fig. 18. Comparison of first order DR-PVU and PVU schemes: Mach number contours (2.0:0.1:3.0) for grid aligned shear flow with 40 · 40 grid. 590 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 PVU DR–PVU 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 0 (b) 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 0 1 (c) (c) 0.5 0 –1 1 (b) 0.5 0 –1 0.5 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 0 Fig. 19. Comparison of first Order DR-PVU and PVU schemes: pressure contours (0.7:0.1:2.9) for oblique shock reflection with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. PVU DR–PVU 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 0 (b) 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 0 1 (c) (c) 0.5 0 –1 1 (b) 0.5 0 –1 0.5 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 0 Fig. 20. Comparison of second order DR-PVU and PVU schemes: pressure contours (0.7:0.1:2.9) for oblique shock reflection with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. where u± are as defined in (19). The diffusion regulator parameters /a and /t take the same value of U as per Eqs. (6)–(9) with the minimum finite dissipation required for avoiding expansion shocks being used only in /t as /t ¼ ðDM 2 þ d2 Þ 2d if jDMj < d ð31Þ The following results presented for PVU and DR-PVU schemes for all the test cases (Figs. 16–26) demonstrate that the diffusion regulation parameter works very effectively with this Kinetic Scheme, reducing numerical 591 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 PVU DR–PVU 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 (b) 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 1 (c) (c) 0.5 0 –1 0.5 (b) 0.5 0 –1 0 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 Fig. 21. Comparison of first order DR-PVU and PVU schemes: pressure contours (1.1:0.05:3.8) for ramp in channel flow with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. PVU DR–PVU 1 1 (a) (a) 0.5 0 –1 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 (b) 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 0.5 –0.5 0 0.5 1 1.5 2 1 0 –1 –0.5 1 (c) (c) 0.5 0 –1 0 (b) 0.5 0 –1 –0.5 1 0.5 –0.5 0 0.5 1 1.5 2 0 –1 –0.5 Fig. 22. Comparison of second order DR-PVU and PVU schemes: pressure contours (1.1:0.05:3.8) for ramp in channel flow with grids: (a) 60 · 20; (b) 120 · 40; and (c) 240 · 80. dissipation in capturing the shocks and capturing steady contact discontinuities aligned with the grid lines exactly, while still maintaining the features of a robust scheme. 5. Diffusion regulation with other schemes and on unstructured meshes; demonstration with Roe and Osher schemes With the reduction of numerical dissipation in Local Lax-Friedrichs scheme and a Kinetic Scheme, which are known to be less accurate than the approximate Riemann solvers of Roe and Osher, a natural question 592 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 PVU DR–PVU 1 1 (a) (a) 0.5 0 0.5 0 0.5 1 1.5 2 2.5 3 1 0 0 (b) 0 0.5 1 1.5 2 2.5 3 0 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 1 (c) (c) 0.5 0.5 0 0.5 1 1.5 2 2.5 3 1 0 0 1 (d) (d) 0.5 0 1.5 0.5 1 0 1 (b) 0.5 0 0.5 1 0.5 0 0.5 1 1.5 2 2.5 3 0 0 Fig. 23. Comparison of first order DR-PVU and PVU schemes: density contours for step in wind tunnel at (a) t = 2, (b) t = 3, (c) t = 4 with grid Dx = Dy = 1/40 and (d) t = 4 with grid Dx = Dy = 1/80. DR–PVU PVU 1.5 1 1.5 (a) 1 0.5 0.5 0 0 –0.5 –0.5 –1 –1 –1.5 –1 0 1 2 –1.5 –1 PVU 1 0 –1 –1 0.5 1 1 (c) 0 0 2 DR–PVU 0 –2 (b) 1 (d) –2 0 0.5 1 Fig. 24. Comparison of first order PVU and DR-PVU schemes for flow across NACA0012 airfoil: (a, b) pressure contours (0.4:0.05:2.8) at Mach number 1.2 and angle of attack 0; (c, d) Cp distribution at Mach number 0.63 and angle of attack 2. S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 PVU 593 DR–PVU 1.5 1.5 (a) 1 (a) 1 0.5 0.5 0 0 –0.5 –0.5 –1 –1 –1.5 –1.5 1 0 1 2 1 1 0 1 2 1 (b) 0.5 (b) 0.5 0 0 –0.5 –0.5 –1 –1 –1.5 –1.5 –2 –0.5 0 –2 1 –0.5 0 1 Fig. 25. Comparison of first order DR-PVU and PVU schemes for transonic flow across NACA0012 airfoil at Mach number 0.85 and angle of attack 1: (a) pressure contours (0.4:0.05:1.8); and (b) Cp distribution. PVU DR–PVU 1.5 PVU 1.5 (a) DR–PVU 1.5 (a) 1.5 (b) (b) 1 1 1 1 0.5 0.5 0.5 0.5 0 0 0 0 –0.5 –0.5 –0.5 –0.5 –1 –1 –1 –1 –1.5 –0.5 0 –1.5 –0.5 0 –1.5 –0.5 0 –1.5 –0.5 0 Fig. 26. Comparison of first order DR-PVU and PVU schemes: density contours (2.0:0.2:5.0) for hypersonic flow across a semi-cylinder: (a) for Mach 6.0 with grid 45 · 45; and (b) for Mach 20.0 with tangentially stretched grid 320 · 20. would be how this new diffusion regulation method works with these later and more accurate methods. In one dimension, Roe’s approximate Riemann solver is designed to capture the steady discontinuities exactly and hence no further reduction in numerical dissipation is possible. However, in multi-dimensional cases, even 594 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 ROE(1.0:0.1:3.8) GRID=60x20 0.8 0.6 0.4 0.2 −0.5 0 0.5 1 1.5 DR–ROE(1.0:0.1:3.8) GRID=60x20 0.8 0.6 0.4 0.2 −0.5 0 0.5 1 ROE(2.0:0.2:5.0) 1.5 1 1 0.5 0.5 0 0 −0.5 −0.5 −1 −1 −1.5 −1.5 −1 −0.5 1.5 DR–ROE(2.0:0.2:5.0) 1.5 0 −1.5 −1.5 0.5 −1 −0.5 0 0.5 Fig. 27. Results with first order DR-Roe and Roe’s methods: pressure contours for supersonic flow over ramp on 60 · 20 grid (top) and hypersonic flow at Mach number 6.0 across a semi-cylinder on 45 · 43 grid (bottom). 1 0.8 0.6 0.4 0.2 0 -1 -0.5 0 0.5 1 1 1.5 2 1 PVU 0.8 OSHER 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -1 -0.5 0 0.5 1 1.5 2 1 -1 -0.5 0 0.5 1 1.5 2 0 0.5 1 1.5 2 1 DR-PVU 0.8 0.8 DR-OSHER 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -1 -0.5 0 0.5 1 1.5 2 -1 -0.5 Fig. 28. Supersonic flow with Mach number 2 across ramp in channel, pressure contours with DR-PVU, DR-OSHER (first order) compared to respective upwind schemes on unstructured mesh (triangular cells: 1286 and edges: 2413). S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 595 Roe scheme cannot capture the discontinuities oblique to the grid lines exactly using a finite volume method. This is due to the inherent dimensional splitting present in the finite volume method. When the diffusion regulation parameter is introduced in Roe scheme for 2D cases, a further reduction in numerical dissipation is seen (Fig. 27). Thus, the diffusion regulation parameter improves the results of even Roe’s approximate Riemann solver in multi-dimensions. Lastly, the diffusion regulation approach is also tried on unstructured meshes, with Osher’s approximate Riemann solver [2] and the Kinetic Scheme, peculiar velocity based upwind (PVU) method. The reduction in numerical dissipation is demonstrated also on unstructured meshes, with the above methods (Figs. 28 and 29). Mesh adaptation also works well with diffusion regulated methods, as can be seen in Figs. 30 and 31. 5.1. Working of diffusion regulator (DR) parameter The diffusion regulator parameter is formulated to take peak values in the shock zone while taking a definite positive value at all other flow regimes except at contact discontinuities. An insight into how the diffusion regulator works during time marching computations is elaborated here (see Fig. 32a). Consider the finite volume method of solution to a Riemann problem and focus on a shock positioned between points J and J + 1 at the end of Nth time level. Assume the shock to move by one grid spacing in the (N + 1)th time iteration and reposition itself at cell interface between (J + 1, J + 2). Since U is designed to take the value of local maximum at the shock, numerical diffusion will be a local maximum at the cell interface between points (J, J + 1) in the (N + 1)th iteration. Thus, the diffusion regulator enforces higher numerical viscosity at the nearest interface to the shock which is crucial to the stabilization of the scheme. 1.5 1.5 1.5 PVU DR-PVU 1 1 1 0.5 0.5 0.5 0 0 0 -0.5 -0.5 -0.5 -1 -1 -1 -1.5 -1.5 -1 -0.5 -1.5 -1.5 0 -1 1.5 -0.5 -1 -0.5 0 1.5 OSHER DR-OSHER 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -1 -1.5 -1.5 -1.5 -1.5 0 -1 -0.5 0 -1.5 -1.5 -1 -0.5 0 Fig. 29. Hypersonic flow with Mach number 6 across semi-cylinder, pressure contours of diffusion regulated upwind schemes (first order) compared to respective upwind schemes on unstructured mesh (triangular cells: 3987 and edges: 7481). 596 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -0.5 0 0.5 1 -1 -0.5 1.5 0 0.5 1 1.5 Fig. 30. Flow over NACA0012 airfoil at Mach number 0.8 and zero angle of attack, pressure contours with diffusion regulated upwind scheme (DR-PVU) with the starting mesh. 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -0.5 0 0.5 1 -1 -0.5 1.5 0 0.5 1 1.5 Fig. 31. Flow over NACA0012 airfoil at Mach 0.8 and zero angle of attack, pressure contours with diffusion regulated upwind scheme (DR-PVU) with adapted mesh. R R J J+1 J+2 J+3 ’LEVEL=N+1’ φ J+1 J+2 J+ J+3 J+4 J+4 J+5 J+5 φ J–2 J–1 J J+1 J+2 R J+3 J J+1 J+2 J+3 J–2 J–1 J J+1 J+3 J J+1 J+2 J+3 J+4 J+5 J J+1 J+2 J+3 J+4 J+5 J+4 J+5 R J+2 φ ’LEVEL=N+2’ J φ J–2 J–2 J–1 J–1 J J J+1 R J+2 J+1 J+2 J+3 R J+3 J J+1 J+2 J+3 ’LEVEL=S’ J–1 STEADY ’LEVEL=S+1,2’ J–2 (b) UNSTEADY ’LEVEL=S+3,4’ ’LEVEL=N’ (a) Fig. 32. Evolution of DR parameter (U) with time level for (a) unsteady and (b) steady cases (R-Riemann problem solution). S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 1.5 1.5 1.5 N=250 N=500 N=750 1 1 1 0.5 0.5 0.5 0 0 1 2 3 1.5 0 0 1 2 3 1.5 0 1 0.5 0.5 0.5 2 3 1.5 0 0 1 2 3 1.5 0 1 0.5 0.5 0.5 2 1 3 0 2 3 N=4000 1 1 0 N=2000 1 0 3 1.5 N=1750 0 2 N=1500 1 1 1 N=1250 1 0 0 1.5 N=1000 0 597 0 1 2 3 0 0 1 2 3 Fig. 33. Evolution of density (  ) and DR parameter (—–) with number of iterations (N) for quasi-1D converging–diverging nozzle flow. Once the residual becomes considerably low, the shock position becomes steady. With further time marching, the residual falls further and the shock becomes crisper (see Fig. 32b). Assume that the shock position becomes steady at time level S at the cell interface between points (J + 2, J + 3). Then, further marching to levels S + 1 and S + 2 results in further fall of residue and self-regulation of numerical diffusion around the shock zone takes place as described below. DM at the entrance to the shock zone becomes lower compared to previous time levels while the DM at the shock interface itself becomes higher. Thus, a two way process of reduction in numerical diffusion at the entrance to shock zone along with an increase in dissipation at the shock interface occurs. Since the maximum and minimum value of DR parameter is restricted to unity and d respectively, this two way process stabilizes with shock being captured with fewer points. The evolution of the DR parameter with number of iterations, while moving from unsteady to steady solution for quasi-1D flow through a nozzle, is given in Fig. 33. A constant value of parameter due to the residual dissipative flux is maintained till the shock discontinuity appears. When a shock discontinuity evolves, the parameter increases in the adjacent zones with its peak value at the shock. Once the position of shock starts becoming steady, the neighborhood points to the shock are pulled to the smooth solution regime in the subsequent iterations (shock becomes sharper), due to the dissipation regulation mechanism. This can be seen in the figure for iterations between 1750 and 2000. Ultimately, a steady solution sets in when the peak of DR parameter becomes steady (seen for iterations between 2000 and 4000) with fewer points across the shock than the parent scheme. It can be seen from the numerical experiments that the gradual variation of the diffusion regulation parameter based on the jump in Mach number strongly aids in self adjustment of numerical dissipation around the shock, ultimately resulting in a sharper shock solution. The model expressed as DM is found to be highly effective in a gradual increase of dissipative flux in shock vicinity while a minimum non-zero dissipation farther from shocks significantly contributes to stability of the scheme. Added to the mechanism described above, the exponential term in the DR model helps in capturing the steady contact discontinuities exactly. 6. Conclusions A diffusion regulation parameter is presented for Euler solvers. This parameter is based on jump in the Mach number and adjusts itself automatically in different flow regimes. With this diffusion regulation 598 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 parameter, any Euler solver can capture steady contact discontinuities aligned with the grid lines exactly and this approach leads to low diffusion schemes. By coupling this parameter with one of the simplest of the numerical methods, the Local Lax-Friedrichs (Rusanov) method, a very efficient Euler solver is developed. The efficacy of this approach is also demonstrated with a diffusive Kinetic Scheme and this approach improves multi-dimensional solutions of even the low diffusive methods like Roe and Osher schemes. Several benchmark problems are presented, demonstrating the exact capturing of steady contact discontinuities and reduction in numerical dissipation. This approach is currently being applied to viscous flows. Acknowledgment We thank the reviewer for his valuable suggestions which improved the quality of this manuscript. References [1] P.L. Roe, Approximate Riemann solvers, parameter vectors and difference schemes, Journal of Computational Physics 43 (1981) 357– 372. [2] S. Osher, F. Solomon, Upwind difference schemes for hyperbolic systems of equations, Mathematics of Computation 38 (1982) 339– 374. [3] Meng-Sing Liou, C.J. Steffen Jr., A new flux splitting scheme, Journal of Computational Physics 107 (1993) 23–39. [4] Meng-Sing Liou, Progress towards an improved CFD method: AUSM+, AIAA Paper, AIAA-1995-1701-CP, 1995. [5] Yasuhiro Wada, Meng-Sing Liou, An accurate and Robust flux splitting scheme for shock and contact discontinuities, SIAM Journal of Scientific Computing 18 (3) (1997) 633–657. [6] A. Jameson, Analysis and design of numerical schemes for gas dynamics 2: Artificial diffusion and discrete shock structure, International Journal of Computational Fluid Dynamics 5 (1995) 1–38. [7] A. Jameson, Analysis and design of numerical schemes for gas dynamics 1: Artificial diffusion, upwind biasing, limiters and their effect on accuracy and multigrid convergence, International Journal of Computational Fluid Dynamics 4 (1995) 171–218. [8] J.L. Steger, R.F. Warming, Flux vector splitting of the inviscid gas dynamic equations with applications to finite difference methods, Journal of Computational Physics 40 (1981) 263–293. [9] B. van Leer, Flux vector splitting for Euler equations, Lecture Notes in Physics 170 (1982) 507. [10] S.M. Deshpande, Kinetic flux splitting schemes, in: M.M. Hafez, K. Oshima (Eds.), Computational Fluid Dynamics Review 1995: A State-of-the-Art Reference to the Latest Developments in CFD, Wiley, Chichester, 1995. [11] S. Jin, Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Communications in Pure and Applied Mathematics 48 (1995) 235–276. [12] S.V. Raghurama Rao, K. Balakrishna, An Accurate Shock Capturing Algorithm with a Relaxation System (ASCARS) for Hyperbolic Conservation Laws, AIAA Paper, AIAA-2003-4115, 2003. [13] K. Prendergast, K. Xu, Numerical hydrodynamics from gas-kinetic theory, Journal of Computational Physics 109 (1993) 53–66. [14] K. Xu, K. Prendergast, Numerical Navier–Stokes solutions from gas-kinetic theory, Journal of Computational Physics 114 (1994) 9– 17. [15] Dominic D.J. Chandar, S.V. Raghurama Rao, S.M. Deshpande, A one point shock capturing kinetic scheme for hyperbolic conservation laws, in: Proceedings of the 3rd International Conference on CFD, 12–16 July, 2004, to be published by Springer. [16] V.V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, Journal of Computational Mathematical Physics, USSR 1 (1961) 267–279. [17] Randall J. LeVeque, Finite Volume Methods for Hyperbolic Problems (Cambridge Texts in Applied Mathematics), Cambridge University Press, Cambridge, 2002, p. 232. [18] S.V. Raghurama Rao, S.M. Deshpande, A Class of Efficient Kinetic Upwind Methods for Compressible Flows, 91 FM 11, Fluid Dynamics Reports, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India, 1991. [19] S.V. Raghurama Rao, New upwind methods based on kinetic theory for inviscid compressible flows, PhD Thesis, September 1994, Department of Mechanical Engineering, Indian Institute of Science, Bangalore, India. [20] S.V. Raghurama Rao, S.M. Deshpande, Peculiar velocity based upwind method for inviscid compressible flows, Computational Fluid Dynamics Journal 3 (1995) 415–432. [21] S.V. Raghurama Rao, Peculiar Velocity based Upwind Method for Inviscid Compressible Flows, Lecture Notes in Physics, vol. 453, Springer, Berlin, 1995, pp. 112–116. [22] S. Jaisankar, S.V. Raghurama Rao, A Low Dissipative Kinetic Scheme for Compressible Flows, 2004 FM 18, Fluid Mechanics Reports, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India, 2004. [23] S. Jaisankar, S.V. Raghurama Rao, A Diffusion Regulator (DR) Model for Low Dissipation in Upwind/Kinetic Schemes, AIAA Paper No. AIAA-2005-4700. [24] A. Jameson, W. Schmidt, E. Turkel, Numerical Solutions of the Euler Equations by Finite Volume Methods using Runge–Kutta Time-Stepping Schemes, AIAA Paper No. AIAA-1981-1259. [25] Culbert B. Laney, Computational Gasdynamics, Cambridge University Press, Cambridge, 1998 (Chapter 17–18). S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 221 (2007) 577–599 599 [26] John D. Anderson Jr., Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill Inc., 1995, p. 356. [27] S.R. Chakravarthy, A. Harten, S. Osher, Essentially Non-Oscillatory Shock Capturing Schemes of Arbitrarily-High Accuracy, AIAA Paper No. AIAA-86-0339, 1986. [28] M. Manna, A Three Dimensional High Resolution Upwind Finite Volume Euler Solver, von Karman Institute for Fluid Dynamics, Technical Note 180, April 1992. [29] H.C. Yee, R.F. Warming, A. Harten, in: Proceedings in Eighth International Conference on Numerical Methods in Fluid DynamicsLecture Notes in Physics, vol. 141, Springer, New York/Berlin, 1982, p. 547. [30] D. Levy, K.G. Powell, B. van Leer, Use of a rotated Riemann solver for two dimensional Euler equations, Journal of Computational Physics 106 (1993) 201–214. [31] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, Journal of Computational Physics 54 (1984) 115–173. [32] H. Viviand, Numerical Solutions of two-dimensional reference test cases; in test cases for inviscid flow filed methods, AGARD 211, 1985. [33] GAMM workshop on Numerical Solutions of Compressible Euler Flows, June 1986. [34] Meng-Sing Liou, Mass flux schemes and connection to shock instability, Journal of Computational Physics 160 (2000) 623–648. [35] M. Pandolfi, D. D’Ambrosio, Numerical instabilities in upwind methods: analysis and cures for the Carbuncle phenomenon, Journal of Computational Physics 166 (2001) 271–301. [36] M. Dumbser, J.M. Moschetta, J. Gressier, A matrix stability analysis of the carbuncle phenomenon, Journal of Computational Physics 197 (2004) 647–670.