A Central Rankine-Hugoniot Solver for Hyperbolic Conservation Laws

Journal of Computational Physics 228 (2009) 770–798 Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp A central Rankine–Hugoniot solver for hyperbolic conservation laws S. Jaisankar, S.V. Raghurama Rao * CFD Centre, Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560 012, India a r t i c l e i n f o Article history: Received 1 December 2007 Received in revised form 26 September 2008 Accepted 1 October 2008 Available online 11 October 2008 Keywords: Rankine–Hugoniot jump condition Central scheme Low numerical diffusion Exact capturing of shocks and contact discontinuities a b s t r a c t A numerical method in which the Rankine–Hugoniot condition is enforced at the discrete level is developed. The simple format of central discretization in a finite volume method is used together with the jump condition to develop a simple and yet accurate numerical method free of Riemann solvers and complicated flux splittings. The steady discontinuities are captured accurately by this numerical method. The basic idea is to fix the coefficient of numerical dissipation based on the Rankine–Hugoniot (jump) condition. Several numerical examples for scalar and vector hyperbolic conservation laws representing the inviscid Burgers equation, the Euler equations of gas dynamics, shallow water equations and ideal MHD equations in one and two dimensions are presented which demonstrate the efficiency and accuracy of this numerical method in capturing the flow features. Ó 2008 Elsevier Inc. All rights reserved. 1. Introduction Numerical simulation of hyperbolic conservation equations, which model the physics of gas dynamics, magnetohydrodynamics (MHD) or shallow water flows, is non-trivial and is still a topic of intense research, despite being attempted for the last five decades, especially by the researchers in Computational Fluid Dynamics (CFD). The challenge is to device a numerical method that is stable, conservative, respects the hyperbolicity, satisfies the entropy condition, is robust and yet accurate. Since the non-linear hyperbolic partial differential equations give rise to discontinuities like shock waves, contact discontinuities, slip surfaces and also develop expansion waves with possible sonic points, any numerical method developed for such systems must be able to capture these features with sufficient accuracy without losing robustness. The numerical methods developed for solving hyperbolic conservation laws can be broadly classified into central discretization methods and upwind methods. While the central discretization methods are simpler than the upwind methods, they often contain tuning parameters which are problem dependent. Because of this drawback, upwind methods became more popular from 1970s. Reviews of these methods are available in several books, for example, by Hirsch [1,2], Toro [3], Laney [4] and Leveque [5,6]. Of all the upwind methods, the approximate Riemann solver of Roe [7] is the most popular, as it can capture steady discontinuities exactly, without any numerical dissipation. However, this low dissipation comes at the cost of the loss of robustness, leading to a list of problems like carbuncle shocks, kinked Mach stems and odd–even decoupling [8], apart from the violation of the entropy condition. Many of these failures seem to be common to most of the Riemann solvers and, therefore, some of the researchers started looking for numerical methods that are free of Riemann solvers. The reader is referred to the following literature for some of the numerical methods developed as interesting alternatives: Nessyahu and Tadmor [9] and Swanson and Turkel [10] (see also [11–13]). * Corresponding author. Tel.: +91 80 2293 3031; fax: +91 80 2360 0134. E-mail addresses: jshankar.s@gmail.com (S. Jaisankar), raghu@aero.iisc.ernet.in, svraghuramarao@yahoo.com (S.V. Raghurama Rao). 0021-9991/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2008.10.002 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 771 Of all the central discretization methods, Lax–Friedrichs method is the simplest (it is also one of the first successful numerical methods for hyperbolic conservation laws) [14,15]. This simple central discretization method has been the building block for several improved algorithms developed later. In a finite volume method, with a staggered grid interpretation of the Lax–Friedrichs method, the cell interface flux consists of an average of the fluxes on the left and right states, augmented by a diffusive flux which represents the numerical dissipation. The coefficient of numerical diffusion in the Lax–Friedrichs method has a fixed magnitude. In Rusanov’s method [16,6], which is also called as Local Lax–Friedrichs (LLF) method, the dissipation coefficient in the Lax–Friedrichs method is modified according to the local information (maximum of the left and right states of a cell interface) and consists of the spectral radius of the flux Jacobian matrix. Since the numerical dissipation in either the basic Lax–Friedrichs method or the LLF method is large and unacceptable for practical problems, higher order accurate formulations in this format were introduced by Nessyahu and Tadmor [9] which are free of Riemann solvers and complicated flux splittings and yet simple and accurate. For further developments of this method, the reader is referred to the work of Kurganov and Tadmor [13] and to the references therein. Based on a central discretization strategy, a diffusion regulation mechanism, which makes any numerical method capture steady contact discontinuities exactly, was introduced by Jaisankar and Raghurama Rao [17]. In this paper, a numerical method, which retains the simplicity of the central discretization methods and is free of Riemann solvers or complicated flux splittings but yet is very accurate and captures the steady discontinuities exactly, is presented (see also [18]). The basic idea of the present work is to enforce the Rankine–Hugoniot condition at the discrete level, based on a simple central type discretization (modifying the LLF method). For the history of the other classical alternative, the Lax–Wendroff method [19], and various related second order central schemes, the reader is refer to the book of Hirsch [2]. One popular central scheme of interest is due to Jameson, Schmidt and Turkel (JST) [20] in which spatial discretization is a blend of second and fourth central differences. A shock sensor is used to switch over from one type of difference to another. In smooth flows, a linear fourth difference will operate, damping the high frequencies, while near the large gradients, the second difference is switched on. Another popular scheme is of Swanson and Turkel [10] (see also [11]), in which the scalar coefficient of numerical dissipation is replaced by a matrix, mimicking the upwind schemes, still retaining the simpler format of central discretization. The method presented in this work is closer in spirit to the first group of central schemes mentioned above, as the basic building block is the Local Lax–Friedrichs method and the higher order accuracy is obtained by reconstruction at the cell interfaces, as in the Kurganov and Tadmor [13] scheme. The motive for the present work is to develop a simple central discretization method, building on the Local Lax– Friedrichs (LLF) method, capturing steady discontinuities aligned with the grid lines exactly and enhancing the resolution of the discontinuities otherwise. It is also worth mentioning here that modifying the scalar dissipation is also found in the popular AUSM family of schemes of Meng-Sing Liou (see [21,22] and the references therein), and the spirit of those works is closer to that of upwind schemes. The exact shock capturing feature of Roe scheme [7] (in steady state) is due to the enforcement of the conservation condition, which is related to the Rankine–Hugoniot (R–H) condition. Another scheme which can capture the steady shock waves with a single interior point is due to Jameson [24] (a sequel to [23]). In this method too, the R–H condition is enforced in the discretization process, due to the special discrete shock structure being forced in such a way that entire shock is captured in one cell. The low dissipative nature of these two methods is owing to the jump condition being used in the discretization process, though not enforced directly. Motivated by the idea of using the R–H directly in the discretization process, Raghurama Rao and Balakrishna [25] developed an accurate Relaxation Scheme which can capture steady discontinuities exactly. In this work, the jump condition is enforced directly in the discretization of the hyperbolic conservation laws. In any conservative numerical method (usually a finite volume method), the interface flux can be written as a sum of an average of the fluxes on the left and right sides of the interface and a diffusive flux. The average flux, if present alone, leads to a central discretization. The diffusive flux represents the numerical diffusion or the artificial viscosity. In the central discretization methods developed earlier, the diffusive fluxes were added explicitly and since the amount of artificial viscosity to be added was not known a priori, tuning parameters became necessary. It was in this background that upwind methods became popular, as they had inherent numerical dissipation and thus did not require tuning parameters. However, the right amount of artificial viscosity, required to capture discontinuities accurately and yet make the numerical method satisfy the entropy condition for avoiding unphysical expansion shocks, is not easy to obtain. These two features actually demand opposite requirements on the numerical dissipation: the accurate capturing of discontinuities requires as less numerical dissipation as possible, whereas capturing expansion waves with sonic points without violating the entropy condition requires a finite amount of it. Many of the upwind methods developed for vector conservation laws have either too much of artificial viscosity (which makes them robust but inaccurate) or too little (which makes them accurate but less robust). Thus, these methods are either robust or accurate and it seems difficult to obtain both these features simultaneously. Another feature of most of the upwind methods is rather sophisticated procedures of separating the information coming from different sides, especially for vector conservation laws, leading to quite complicated flux splittings. In the method proposed here, the simplicity of the central discretization is retained and yet the resolution of the shocks and contact discontinuities is enhanced, by using the Rankine–Hugoniot (jump) condition as a basis for deriving the diffusive flux. Thus, the numerical method presented in this work can be considered as a natural extension of both the central and upwind discretization methods, retaining the advantages of both. Since the Rankine–Hugoniot condition is used in discretization, and thus the optimum value of artificial viscosity is enforced, we term this method as the method of optimal viscosity for enhanced resolution of shocks (MOVERS). Note that though the Rankine–Hugoniot conditions were derived first in gas dynamics for shock waves, these jump conditions hold across both shock waves and contact discontinuities [26]. In this pa- 772 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 per, we develop an accurate shock capturing method but the advantages of the method carry over to the capturing of contact discontinuities or slip surfaces also. The basic idea of the new method for the scalar Burgers equation is demonstrated in the next section. In the following sections, this procedure is extended to the vector conservation laws (Euler equations of gas dynamics, shallow water equations and ideal MHD equations). The idea of enforcing the jump condition at the discrete level to the vector case is presented with three different variations. Two of these variants of MOVERS require an entropy fix while the third variant does not need it as it is based on a novel use of a limiter in switching over from MOVERS to LLF method. The results for several bench-mark test problems for scalar and vector conservation laws are presented in the later sections, followed by some remarks about the MOVERS and its relation to upwinding strategy, with the conclusions presented in the last section. 2. Method of optimal viscosity for enhanced resolution of shocks (MOVERS) for a scalar hyperbolic conservation law Consider the one-dimensional scalar hyperbolic conservation law ou of ðuÞ þ ¼0 ot ox ð1Þ An explicit scheme in conservation form on a three point stencil and with a piecewise polynomial approximation for the conserved variable is given by unþ1 ¼ unj  j  Dt  n n fjþ1  fj 1 2 2 Dx ð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 fI ¼ fjþ1 evaluated at the cell interface (represented by the subscript I). The upwind flux splitting schemes define the 2 interface flux in terms of split fluxes of left and right states as fI ¼ fLþ þ fR ð3Þ The interface flux in a central scheme is written in the form fI ¼ 1 1 ðfL þ fR Þ  aDu 2 2 ð4Þ where L and R represent the left and right states of the interface, the first term on the right hand side (average flux) represents the central discretization of the flux terms and the second term is the diffusive flux, with a being the coefficient of numerical diffusion. The interface flux for the upwind method (3) can also be recast in the form (4) by using the definitions of the split fluxes f  . Then the coefficient of numerical dissipation a will be fixed as per the chosen scheme. We can see that the central discretization method is simpler than the upwind method, but the amount of a required to produce a stable, robust and accurate scheme is unknown a priori. The basic idea in this work is to fix a based on a relevant criterion to enhance the resolution of shocks (and also contact discontinuities). Note that a must have the dimension of the wave speed, aðuÞ ¼ ofouðuÞ (e.g. for the case of Burgers equation, 2 f ðuÞ ¼ u2 and aðuÞ ¼ u) in (4) as aDu must represent a diffusive flux. In Local Lax–Friedrichs method [16], the value of a is taken as the maximum value of the wave speed of the left and right states of a cell interface. This makes the LLF method very simple but not so accurate. A better choice will be to use some physically relevant speed for choosing a. Since the motivation in this work is to enhance the resolution of the shocks (and other discontinuities), the most relevant physical speed is the shock speed (the speed with which a discontinuity travels) which is given by the Rankine–Hugoniot (jump) condition. Here, the Rankine–Hugoniot condition is enforced in the discretization process, by choosing a relevant expression for a. 2.1. Rankine–Hugoniot condition and cell interface flux For a shock wave with a left state and right state given by uL and uR , the Rankine–Hugoniot (jump) condition is given by Df ¼ sDu where Df ¼ fR  fL and Du ¼ uR  uL ð5Þ To enforce this jump condition in the numerical method, the best way is to obtain the interface flux from the same condition. The jump condition does not directly give the interface flux, as it describes only the relation between left and right states of the shock and its speed. In the finite volume method, because of the piece-wise polynomial approximation, a discontinuity in the conserved variable is present at the cell interface. Let us apply the above R–H (jump) condition at the cell interface, j þ 12 fjþ1  fj ¼ sjþ1 ðujþ1  uj Þ 2 ð6Þ To obtain the interface flux, we split the above R–H condition into two parts (see [27]). Suppose there is a shock wave between j and j þ 1 moving with a speed sjþ1 as shown in Fig. 1. 2 This shock wave can move to the left or right. The jump condition does not indicate the direction of the shock movement. To distinguish the positive and negative speeds, let us split the shock speed into two parts: a positive one and a negative one S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 s– j+1/2 j–1 s+ j+1/2 j j–1/2 773 j+1 j+1/2 Fig. 1. Shock moving to left or right. sjþ1 ¼ sþjþ1 þ sjþ1 2 ð7Þ 2 2 where sþjþ1 ¼ 2 We have     sjþ1 þ sjþ1  2 2 2 sþjþ1  sjþ1 ¼ 2 2 and sjþ1 ¼ 2     sjþ1  sjþ1  2 2 2 1 1   1 1     sjþ1 þ sjþ1   sjþ1 þ sjþ1  ¼ sjþ1  2 2 2 2 2 2 2 2 2 ð8Þ ð9Þ Since the left and right moving shocks are separated, we can now split the Rankine–Hugoniot condition into two parts as fjþ1  fjþ1 ¼ sþjþ1 ðujþ1  uj Þ 2 2 ð10Þ fjþ1  fj ¼ sjþ1 ðujþ1  uj Þ 2 2 Subtracting the second from the first in the above Eq. (10), we obtain fjþ1 ¼ 2 1 1   ðfj þ fjþ1 Þ  sjþ1 ðujþ1  uj Þ 2 2 2 ð11Þ This is the interface flux obtained from the Rankine–Hugoniot condition. Let us now consider the cell interface flux (4) obtained from the finite volume method, given by fjþ1 ¼ 2 1 1 ðfj þ fjþ1 Þ  ajþ1 ðujþ1  uj Þ 2 2 2 ð12Þ Comparing this flux (12) with the interface flux obtained from the Rankine–Hugoniot condition (11), we obtain   ajþ12 ¼ sjþ12  ð13Þ 2.2. Wave speed evaluation and update formula We can evaluate the value of sjþ1 numerically as sjþ1 2 8 f fj > < ujþ1 jþ1 uj  ¼ of  > : ou j 2 if uj – ujþ1 if uj ¼ ujþ1 ð14Þ Similarly, the flux at the cell interface j  12 can be computed. In the numerical evaluation of the shock speed as shown above, the denominator involves the differences of the conserved variables across a cell interface. When the left and right states of the cell interface are very close, this denominator can be very small and the shock speed can become unphysically large. To avoid this problem, we restrict its numerical value within the minimum and maximum of the wave speed in the computational domain if sjþ1 > amax then sjþ1 ¼ amax 2 2 if sjþ1 < amin then sjþ1 ¼ amin 2 2 of ðuÞ is the wave speed where a ¼ ou ð15Þ 774 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 The cell interface fluxes derived above can now be used in the finite volume update formula unþ1 ¼ unj  j i Dt h n n fjþ1  fj 1 2 2 Dx ð16Þ Since the coefficient of numerical dissipation is now fixed in such a way that the Rankine–Hugoniot condition is satisfied in the expressions for the cell interface fluxes, steady shocks are captured exactly by this numerical method. The final update formula for the above central scheme, which we name as the method of optimal viscosity for enhanced resolution of shocks (MOVERS), for a scalar conservation law is almost the same as that of the approximate Riemann solver of Roe [7], despite the derivation and implementation being different and the shock speed limiting within the maximum and minimum of the wave speeds being added. Though the approximate Riemann solver of Roe was derived first for vector conservation laws, it can be simplified for scalar conservation laws (see [2,3]) and the above comparison illuminates how the present method is both similar to and yet different from the popular upwind method of Roe. Note that this method is not based on Roe linearization and will be advantageous for use in the case of vector conservation laws, as will be seen in the next section, where we obtain a novel numerical method which is simple and yet accurate and neither the characteristic decomposition used in the approximate Riemann solvers nor the Roe averages will be needed. The coefficient of numerical dissipation, aj1 is linearly proportional to the numerical wave speed and thus the dissipa2 tion goes to zero when the wave speed goes to zero. This happens near the sonic points where the wave speed changes sign. Zero dissipation results in pure central differencing at these points, resulting in unphysical solutions such as expansion shocks which violate the entropy condition. To overcome this problem, Harten’s entropy fix [28], like in the case of Roe’s scheme, was used in [25]. Alternatively, for a scalar equation we can sense the diverging of waves at the expansion and locally avoid the R–H condition, taking for a the numerical value of the local maximum of the wave speeds. This strategy ensures avoiding of unphysical solutions, while capturing the steady shocks exactly. In the later sections, we also present a variation of this method in which the entropy fix is avoided by a novel use of a limiter, which helps in switching over from MOVERS to LLF method. In the next section, we extend the MOVERS to the hyperbolic vector conservation laws. 3. Method of optimal viscosity for enhanced resolution of shocks (MOVERS) for a hyperbolic system of equations Consider the hyperbolic system of conservation laws in one dimension, as oU oF þ ¼0 ot ox ð17Þ where U is vector of conserved variables, U i ; i ¼ 1; . . . ; n and F is the flux vector which may have a non-linear functional relationship with U. The hyperbolicity can be seen in the quasi-linear form of the above equations, given by oU oU þA ¼0 ot ox ð18Þ oF where A is the flux Jacobian matrix (A ¼ oU : a n  n matrix). The system is hyperbolic if the eigenvalues of the flux Jacobian matrix A, denoted by ki ; i ¼ 1; . . . ; n are real and the corresponding eigenvectors are linearly independent. If the eigenvalues are also distinct, the system is strictly hyperbolic. For a hyperbolic system, we can diagonalize the matrix as A ¼ RDR1 , where the diagonal elements of D are the characteristic values or eigenvalues of A. However, the numerical method presented in this work is not dependent on the eigenvectors and neither the evaluation of the eigenvectors nor the projection of the solution on to the space of eigenvectors is required for this scheme, unlike in the case of Riemann solvers. An explicit scheme in conservation form on a three point stencil and with piecewise polynomial approximation for conserved variable is given by U jnþ1 ¼ U nj   Dt  n F 1  F nj1 2 Dx jþ2 ð19Þ where Dt and Dx are respectively time step and the grid-spacing, j represents the cell centroid and j  12 represent the cell interfaces. The interface flux is written as in a central scheme in the form FI ¼ 1 1 ðF L þ F R Þ  aðU R  U L Þ 2 2 ð20Þ where L and R represent the left and right states of the cell interface. Note that a is the coefficient of numerical dissipation which will be chosen later based on the Rankine–Hugoniot condition. The Rankine–Hugoniot condition is DF ¼ sDU ð21Þ Note that while conserved variable vector U and the flux vector F are vectors, the shock speed s is a scalar. Proceeding in the same way as was done for the scalar case, we can obtain the coefficient of numerical dissipation, based on enforcing the Rankine–Hugoniot jump condition in the discretization, as S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798   ajþ12 ¼ sjþ12  775 ð22Þ So far, the situation is similar to that of the scalar case. Now, evaluating s from the Rankine–Hugoniot condition, DF ¼ sDU, is not straightforward, as U and F are vectors but s is a scalar. We present three different possibilities of evaluating the shock speed, leading to three different variations of the MOVERS, which we name as (1) MOVERS-n – with a nwave approximation for the wave speed, (2) MOVERS-1 – with a one wave approximation for the wave speed and (3) MOVERS-L – a limiter based central solver which switches over from R–H condition based dissipation of n (or 1)-wave approximation method to that of LLF method. The first two variants will require an entropy fix, while the last one does not need it. All the three variants will capture steady shocks (and also contact discontinuities or slip surfaces) with enhanced resolution, resolving them exactly wherever appropriate (in 1D and when they are aligned with grid-lines in multi-dimensions). 3.1. MOVERS-n One possibility is to choose s and, therefore, a as diagonal matrices. Thus, in the relation DF ¼ sDU, since U and F are n  1 vectors, we can choose s as a n  n matrix. The simplest of such n  n matrices is a diagonal matrix with n diagonal elements 3 2 s1;jþ1 DF 1 2 6 DF 7 6 0 6 27 6 6 . 7¼6 . 6 . 7 6 4 . 5 6 4 .. DF n 0 2 0 32 3 DU 1 7  0 76 DU 2 7 7 76 6 . 7 7 .. 76 . 7  . 54 . 5 DU n    sn;jþ1 0  s2;jþ1 2 .. . 0 ð23Þ 2 We can now evaluate the shock speeds as sm;jþ1 ¼ 2 DF m ; DU m m ¼ 1; 2; . . . ; n ð24Þ As an example, let us consider the case n ¼ 3. Then, we obtain 3 2s 1 1;jþ2 DF 1 6 7 6 0 4 DF 2 5 ¼ 6 4 DF 3 0 2 0 32 3 DU 1 76 7 0 7 54 DU 2 5 DU 3 s3;jþ1 0 s2;jþ1 2 0 ð25Þ 2 We can now calculate the shock speeds (three of them) as 2 s1;jþ1 3 2 ðF 1 ÞR ðF 1 ÞL ðU 1 ÞR ðU 1 ÞL 3 6 7 2 6 s 1 7 6 ðF 2 ÞR ðF 2 ÞL 7 4 2;jþ2 5 ¼ 6 ðU2 ÞR ðU2 ÞL 7 4 5 s3;jþ1 ðF 3 ÞR ðF 3 ÞL 2 ð26Þ ðU 3 ÞR ðU 3 ÞL This is an approximation based on ‘n’ waves for a hyperbolic system with ‘n’ equations (n ¼ 3 in the illustrated example considered above, as in the case of 1D Euler equations). Such simplification of the matrix (to a diagonal matrix) in the Rankine– Hugoniot condition leads to an elegant method which can be easily extended to any hyperbolic system of conservation laws, even for more involved systems such as the magneto-hydrodynamics equations. We term this method as MOVERS-n, representing n-wave approximation for the shock speed evaluation. 3.2. MOVERS-1 A simpler alternative to n-wave approximation to shock speed in MOVERS is a one-wave approximation. This changes the dissipation from a diagonal matrix form to a scalar form. One way of choosing a scalar dissipation, satisfying the Rankine–Hugoniot condition as derived before, is to choose one among all the n-waves representing the shock speed. A simple guide line for this choice is to find out which of the n-waves has the maximum possible information. It is still better if this information also reflects the expected behavior of characteristics converging at the shocks and diverging in expansion regions. Such a Central Rankine–Hugoniot Solver with one-wave approximation for the shock speed is named as MOVERS-1. To illustrate this idea, we consider the system of Euler equations in one dimension as follows. 3.2.1. MOVERS-1 for 1D Euler equations For the 1D Euler equations, the procedure in the previous approximation leads to three values of the shock speed as in the equation (26). For the n-wave approximation, each of the diagonal elements becomes the coefficient of dissipation for each of the equations in the system and this makes a a matrix. Instead, here we use a single wave approximation and hence seek a scalar dissipation coefficient. 776 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 To choose a single wave, we start with a criterion based on suitability of the numerical wave speed to represent the actual physical requirements for a shock wave and an expansion wave. The characteristics must converge on to the shock wave and must diverge at an expansion wave. Since the coefficient of numerical diffusion represents a wave speed, it is preferable for it to mimic these properties. We also consider the amount of information carried by each of the terms while choosing the single wave, as the best term will carry the maximum information. Numerical experiments for the 1D shock tube problem and the quasi-1D nozzle problem have shown that out of all the three components in (26), the third coefficient (for the energy equation) has the desirable properties. The wave speed for the mass conservation does not represent merging of characteristics at shocks while for the momentum equation, the second component of (26) can lead to problem of denominator becoming zero for equal left and right fluxes (for a steady discontinuity) even when a Riemann problem persists. The energy term has information of all the primitive variables and thus contains the maximum of information. The plot of a (which represents the wave speed) for quasi-1D nozzle flow shown in Fig. 2(a) changes sign at two points (at X  1:5 and X  2:1) where an expansion and a shock occur, respectively. It can be seen that wave speeds diverge at X ¼ 1:5 (with negative wave speed on the left and positive wave speed on the right of sonic point) and converge at X ¼ 2:1 (with positive wave speed on the left and negative wave speed on the right of sonic point), indicating that the choice of the third component (corresponding to energy equation) for the wave speed in the 1-wave approximation is the correct choice. The wave speed behavior for the unsteady flow in shock tube (Fig. 2(b)) also shows the same behavior, by a sign change from negative to positive for an expansion at ðX  4:5Þ and changing sign from positive to negative at the shock at ðX  þ4:5Þ. Therefore, the third component in (26) is chosen as the scalar coefficient of numerical dissipation for all the three equations. Thus, we can write 2 s1;jþ1 3 2 ðF 3 ÞR ðF 3 ÞL ðU 3 ÞR ðU 3 ÞL 3 7 6 2 6 s 1 7 6 ðF 3 ÞR ðF 3 ÞL 7 4 2;jþ2 5 ¼ 6 ðU3 ÞR ðU3 ÞL 7 5 4 s3;jþ1 ðF 3 ÞR ðF 3 ÞL 2 ð27Þ ðU 3 ÞR ðU 3 ÞL This effectively means choosing a scalar numerical dissipation coefficient based on Rankine–Hugoniot condition as a representative numerical wave speed sjþ1 ¼ 2 ðF 3 ÞR  ðF 3 ÞL ðU 3 ÞR  ðU 3 ÞL ð28Þ even for a system of equations. For system of equations representing ideal MHD flow, we again choose the energy term for the same reasons as described for Euler equations and for shallow water equations, we choose the momentum equation which has the maximum information. 3.3. MOVERS-L The two variants of the new scheme proposed here, namely MOVERS-n and MOVERS-1, require an entropy fix to avoid unphysical expansion shocks near the sonic points. Another variant MOVERS-L, in which the numerical dissipation is switched from that in the Local Lax–Friedrichs method to the numerical dissipation based on the Rankine–Hugoniot condition, based on a limiter, is presented here. The coefficient of numerical diffusion of this hybrid scheme is written as a ¼ aMOVERS þ /ðrÞ½aLLF  aMOVERS  ð29Þ where aLLF is the local maximum of the spectral radius of the flux Jacobian matrix and aMOVERS is the corrected wave speed obtained using MOVERS-n (or 1) (see the expressions (33)–(36) in the following subsection), using the R–H condition. Here, / 800 800 Wave Speed Vs Distance(Nozzle1d) Wave Speed Vs Distance(Sods Shocktube) 600 600 400 400 200 200 0 0 0 0.5 1 1.5 2 2.5 3 0 5 10 Fig. 2. Plot of wave speed ðaÞ vs. distance (X) for (a) quasi-1D converging diverging nozzle (steady flow); (b) for Sod’s shock tube problem. S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 777 is a standard limiter. This is to enforce the MOVERS to work only at the discontinuities while reverting to the dissipative (entropy condition satisfying) Local Lax–Friedrichs method elsewhere. We use the minmod limiter suggested by Yee [29] as defined below   /jþ1 r jþ1 ¼ minmodð1; r þj ; r jþ1 Þ 2 ð30Þ 2 The ratio of solution differences are defined as rþj ¼ U j  U j1 U jþ1  U j and r j ¼ U jþ1  U j U j  U j1 ð31Þ If the denominator in the above expression is small (when U j  U jþ1 ), then we modify r to prevent numerical overflow as rþj ¼ signðU jþ1  U j Þ U j  U j1 d if jU jþ1  U j j < d ð32Þ where d is a small parameter. The same approach is used for r j also. The limiter based Central Rankine–Hugoniot Solver, MOVERS-L, therefore, can capture steady discontinuities exactly and still avoid expansion shocks. For a second order accurate version of MOVERS-L, as the left and right states of a cell interface are obtained by reconstruction, the scheme switches from a second order accurate LLF method in smooth regions to second order accurate MOVERS in high gradient regions. Two kinds of limiters will be operating in this case – the first being a standard limiter used for second order reconstruction of conserved variables as seen in Eqs. (38) and (39) of Section 3.6. The other limiter (29) and (30) is specific to this scheme for switching over from MOVERS-n (or 1) to LLF method (here we use only minmod limiter for obtaining a convex combination). For practical applications, since second order accurate LLF method is used in smooth regions, the solution is not degraded by the limiter in these regions. By second order accurate LLF method is meant the LLF method with reconstruction for obtaining the left and right states of the conserved variable, as in Kurganov and Tadmor [13]. 3.4. Wave speed correction In both n-wave and 1-wave versions of the new algorithm, while calculating the wave speeds (expressions (26) and (28), respectively), if the denominator DU becomes small, the condition refers to an isentropic case and we assign the minimum in the characteristic speed spectrum ofhthe flux Jacobian as the wave speed. When there is a steady shock ðs ¼ 0Þ, since the disi crete update formula U jnþ1 ¼ U nj  DDxt F njþ1  F nj1 approximates the Euler equations in a consistent and stable manner, all the 2 2 above representative numerical wave speeds of (26) and (27) reduce to the speed of the shock, which is zero in the steady case. Thus, any version of the Central Rankine–Hugoniot Solver, MOVERS, will capture steady shocks exactly, without any numerical dissipation. Because of the form of the dissipation we have assumed, even when there is a steady contact discontinuity, this algorithm will capture it exactly (which is demonstrated in the section on results). The above strategy is employed only when the gradients are significant. Elsewhere, when the denominator DU is too small, the evaluated value of s can be nonphysically high. We use a wave speed correction algorithm to bring the wave speeds to physically admissible values. The wave speed correction algorithm is a min-max wave speed limiter with the absolute values of discrete wave speed being restricted within the spectrum of characteristic values (eigenvalues) of the hyperbolic system while still retaining the same direction of the computed uncorrected discrete wave speed sjþ1=2 ¼ ( F F j Ujþ1 Uj j if U j –U jþ1 kmin;jþ1=2 if U j ¼ U jþ1 jþ1 j ð33Þ sjþ1 ¼ signðsjþ1 Þkmax if jsjþ1 j P kmax ð34Þ sjþ1 ¼ signðsjþ1 Þkmin if jsjþ1 j 6 kmin ð35Þ 2 2 2 2 2 2 With these corrections, the steps to update the vector conservation laws involve computation of representative numerical wave speeds at all the points in the domain and correction of such wave speeds to make them lie within the bounds of spectral radius. Here, kmax and kmin refer to the local maxima and minima of the magnitudes of characteristic speeds defined by the eigenvalues of the flux Jacobian matrix. As an example, for the Euler equations in one dimension, we define kmax ¼ Maxðmaxðjuj; ju  aj; ju þ ajÞR ; maxðjuj; ju  aj; ju þ ajÞL Þ kmin ¼ Maxðminðjuj; ju  aj; ju þ ajÞR ; minðjuj; ju  aj; ju þ ajÞL Þ ð36Þ since the eigenvalues in this case are u  a, u and u þ a with u being the fluid velocity and a representing the sound speed. 3.5. Entropy fix for MOVERS-n and MOVERS-1 The numerical diffusion in both the n-wave variant and 1-wave variant becomes low at the expansive sonic points where wave speeds change the sign. To avoid any expansion shocks, an entropy fix is used, similar to Harten’s entropy correction [28], as given in the following expression 778 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 sj ¼ s2j þ d2 2d if jsj j < d where d ¼ jkmax and 0 < j < 1 ð37Þ As Local Lax–Friedrichs method is a good alternative to low dissipative schemes like Roe scheme or MOVERS in the regions of expansion fans with sonic points, it is also used as an entropy fix [6]. In the numerical experiments done in this work, it was found that reducing the numerical dissipation of the LLF method by half (i.e., j ¼ 0:5) works reasonably well, avoiding expansion shocks and non-smoothness. Since this entropy correction adds numerical diffusion everywhere, loss of accuracy in shocks can be expected. A suitable correction to this scheme can then be to retain sufficient dissipation in most of the domain (which includes the sonic points), except at the discontinuities. This is the motivation for the limiter based variant of the scheme MOVERS-L. 3.6. Second order accurate schemes The second order accurate schemes are constructed by assuming piecewise linear approximation for the conserved variables, as Uðx; tn Þ ¼ ½U nj þ ðU x Þnj ðx  xj Þ½xj1=2 ;xjþ1=2  ; xj1=2 ¼ xj  Dx 2 ð38Þ where U nj is the computed cell average, and ðU x Þnj are approximations to the exact derivatives, U x ðx; tn Þ, evaluated from the computed cell averages. The non-oscillatory behavior of the scheme depends on the appropriate choice of approximate derivatives, and we use the one parameter family of minmod limiter for this purpose   U j  U j1 U jþ1  U j1 U jþ1  U j ; ;b ; ðU x Þj ¼ minmod b Dx 2Dx Dx 16b62 ð39Þ The numerical values of fluxes and the primitive variables required for the update are computed at the cell interfaces from the above reconstruction procedure. Several benchmark problems of gas dynamics, shallow water flows and ideal magnetohydrodynamics are solved and results are presented in the next section. 4. Results and discussion 4.1. Results with MOVERS for Burgers equation MOVERS is first tested for the inviscid Burgers equation in one dimension, given by ou of ðuÞ 1 þ ¼ 0 where f ðuÞ ¼ u2 ot ox 2 ð40Þ The test case used models a shock wave and an expansion fan [4]. The results with the Local Lax–Friedrichs method are also presented for comparison, in Fig. 3. The shock wave is captured exactly with MOVERS and the expansion fan is captured well with Harten’s entropy fix. The results obtained with MOVERS for the 2D Burgers equation, given by ou of1 ðuÞ of2 ðuÞ 1 þ þ ¼ 0 with f 1 ðuÞ ¼ u2 and f 2 ðuÞ ¼ u ot ox oy 2 ð41Þ are shown in Fig. 4. This test case is taken from Spekreijse [30]. The shock is again captured exactly. 4.2. Results with MOVERS-n for Euler equations The Euler equations in one dimension are given by oU oF þ ¼0 ot ox ð42Þ 2 ð43Þ where 3 2 3 q qu 6 7 6 7 U ¼ 4 qu 5 and F ¼ 4 p þ qu2 5 qE pu þ quE Here U is the vector of conserved variables and F is the flux vector, with q representing the density, u the fluid velocity, p the pressure and E the total energy, defined by E ¼ qðcp1Þ þ 12 u2 . This system of equations is hyperbolic, with the eigenvalues given by k1 ¼ u  a, k2 ¼ u and k3 ¼ u þ a. The quasi-1D and multi-dimensional forms of the Euler equations are given in [3]. Numerical solution of inviscid compressible flows for the 1D case of steady contact discontinuity in a shock tube, the quasi-one-dimensional nozzle flow [31], the 2D cases of oblique shock reflection, supersonic flow over a ramp in a channel [32], S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 LLF EXACT 1.5 u 1 0.5 0.5 u 0 0 MOVERS EXACT 1.5 1 0.5 779 0 0 1 x 0.5 1 x Fig. 3. 1D Burgers equation solution: first order accurate results with MOVERS. 1 0.75 0.5 0.25 0 0 0.25 0.5 0.75 1 Fig. 4. 2D Burgers equation solution with first order accurate MOVERS using 32  32 points: contour plot of the conserved variable, u. 1 0.8 Entropy Fix No Entropy Fix 0.6 0.4 0.2 0 0 5 10 X Fig. 5. Steady contact discontinuity in a shock tube with first order accurate MOVERS-n, with and without entropy fix; density plot. supersonic flow over a forward-facing step in a channel [33], double Mach reflection [33], shear flow of Mach 3 over Mach 2 [34] and hypersonic flow across a half-cylinder [35] are obtained using MOVERS-n. The results consistently show high accuracy of the new solver in capturing flow features. It is worth noting that this is achieved without any Riemann solvers and complicated flux splittings, in a simple format of a central scheme. 780 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 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 POINTS 200 500 1.2 POINTS 50 100 3 0.2 0 0.5 1 1.5 2 2.5 3 Fig. 6. Nozzle flow: density plot with first order accurate MOVERS-n. 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0.2 0.4 0.6 0.8 1 0 0 0.2 0.4 0.6 0.8 1 Fig. 7. Horizontal shear flow (Mach 3 flow over Mach 2 flow), Mach contours with MOVERS-n (first order) on a 20  20 grid. 1 1 0.5 0.5 0 0 1 2 0 1 1 0.5 0.5 0 0 1 2 0 1 1 0.5 0.5 0 0 1 2 0 0 1 2 0 1 2 0 1 2 Fig. 8. Oblique shock reflection, density contours (0.9:0.1:2.7) with MOVERS-n on grids (a) 60  20, (b) 120  40 and (c) 240  80. The result obtained with MOVERS-n for the steady contact discontinuity in a shock tube is shown in Fig. 5. The contact discontinuity is captured exactly when the entropy fix is not used. The results for the standard quasi-1D nozzle flow are shown in Fig. 6. In the horizontal shear flow (2D) test case, the slip surface is captured exactly for a shear flow of Mach 3 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 781 over Mach 2 (Fig. 7), if the entropy fix is not used. The results for the two-dimensional oblique shock reflection problem (Fig. 8) on crude, medium and fine grids, show a consistent improvement in the resolution of the solutions, with grid refinement and from first order to second order accuracy. The trend persists for unsteady problems as well, as depicted in Figs. 9 and 10 for the supersonic flow over a step in a channel and the double Mach reflection problems and all the flow features are captured with high accuracy. Stagnation line pressure and Mach number plots for steady hypersonic flow (M = 6) over a halfcylinder in Fig. 11 highlight the improvement in solution over the parent LLF scheme. C p plots for two bench-mark cases of transonic flow over NACA0012 airfoil [36] are presented in Figs. 12–14, for M ¼ 0:85, a ¼ 1 and M ¼ 0:8, a ¼ 1:25 . MOVERSn without any entropy fix resolves the shocks very accurately in most of the cases, with at most one point in the shock. However, non-smoothness can be seen elsewhere in the plot. With entropy fix, the non-smoothness is removed (with j ¼ 0:25 in (37)) but the shocks get smeared. The results of the n-wave solver are compared with the results of Roe’s scheme in Figs. 15 and 16. The solution for supersonic flow over a ramp in a channel is shown in Fig. 17 using unstructured grids. The results presented demonstrate the efficiency of this version of the Central Rankine–Hugoniot Solver (MOVERS-n) in capturing the flow features accurately, clearly demonstrating its low diffusive nature. The steady shocks aligned with grid lines are captured exactly, without any numer- 1 0.5 0 0 0.5 1 1.5 2 2.5 3 0 0 0.5 1 1.5 2 2.5 3 1 0.5 Fig. 9. Density contours at time = 4 with MOVERS-n using structured mesh of (240  80) points for Mach 3 flow over a forward-facing step in channel. 1 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 1 0.5 0 0 Fig. 10. Density contours at time = 0.2 with MOVERS-n using structured mesh of (240  60) points for double Mach reflection due to Mach 10 shock moving through a 30° wedge. 782 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 60 8 LLF LLF 7 50 6 40 5 30 4 3 20 2 10 1 0 0 X X Fig. 11. Results along stagnation line for Mach 6 flow across half cylinder using structured mesh of 45  45 points with first order accurate MOVERS-n. 1 0.5 0 I order II order 0 0.2 0.4 0.6 0.8 1 Fig. 12. Pressure coefficient distribution on NACA0012 airfoil for Mach 0.85 flow at 1° angle of attack with MOVERS-n. 1 1 0.5 0.5 0 0 I order 0 0.2 0.4 0.6 0.8 I order 1 0 0.2 0.4 0.6 0.8 1 Fig. 13. Pressure coefficient distribution on NACA0012 airfoil for Mach 0.85 flow at 1° angle of attack with MOVERS-n: (a) without entropy fix; (b) with entropy fix. 783 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 1 1 0.5 0.5 0 0 I order 0 0.2 0.4 0.6 0.8 I order 1 0 0.2 0.4 0.6 0.8 1 Fig. 14. Pressure coefficient distribution on NACA0012 airfoil for Mach 0.8 flow at 1.25° angle of attack with MOVERS-n: (a) without entropy fix; (b) with entropy fix. 1 0.5 0 0 0.5 1 1.5 2 0 0.5 1 1.5 2 1 0.5 0 Fig. 15. Comparison of first order accurate solutions of MOVERS-n and Roe’s scheme for supersonic flow over a ramp in a channel. 2 2 0 0 0 0 Fig. 16. Comparison of first order accurate solutions of (a) Roe’s scheme with (b) MOVERS-n for Mach 6 flow across half cylinder. 784 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 17. Supersonic flow over a ramp in a channel: pressure contours obtained with MOVERS-n using an unstructured mesh of 7503 points. ical dissipation. The entropy fix is necessary to suppress unphysical behaviour and non-smoothness in the solution for this version of the algorithm. 4.3. Results with MOVERS-n for shallow water equations The one-dimensional form of shallow water equations in differential form can be written as U t þ FðUÞx ¼ S ð44Þ where U¼  h hu  ; F¼ " hu 2 hu þ 12 gh 2 # ; S¼  0 ghBx The system is hyperbolic with eigenvalues k1 ¼ u  given by F 0 ðUÞ ¼  0 1 u2 þ gh 2u   ð45Þ pffiffiffiffiffiffi pffiffiffiffiffiffi gh and k2 ¼ u þ gh for the one-dimensional flux Jacobian matrix, ð46Þ The non-homogeneous (source) term S is very common and quite often represents the water bed beneath. Numerical solution with MOVERS-n for some bench-mark problems are presented below. 785 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 1D Dam-break problem on a variable depth river-bed This is a test case [37] with variable river-bed and so the system of equations is not homogeneous, which means that a source term is present in the equation. For this test case, the river-bed is defined as BðxÞ ¼ ( 0:125 cos 10 p x  12 þ 1 0 if 0:4 6 x 6 0:6 otherwise and we have the following initial conditions uðxÞ ¼ 0 hðx; 0Þ ¼ ( 1  BðxÞ if 0 6 x 6 12 0:5  BðxÞ if 12 < x 6 1 The presence of discontinuity at x ¼ 0:5 is like a barrier separating two river heights. Reflection boundary conditions are taken into consideration assuming walls at x ¼ 0 and 1. Results presented in Fig. 18 show all the flow features being captured well. 4.3.1. 1D Steady state trans-critical flow with shock In this benchmark test case [38], we study the convergence towards steady flow over a hump. We consider the system with initial conditions uðx; 0Þ ¼ 0 hðx; 0Þ þ ZðxÞ ¼ H0 where H0 is the constant water level downstream provided by the boundary condition. We test our scheme to this steady state flow where the bathymetry is non-trivial and is given by ZðxÞ ¼ ( 0:2  0:05ðx  10Þ2 0 if 8 6 x 6 12; otherwise in a channel of length L ¼ 25 m and H ¼ 2 m, CFL ¼ 0:5. In this test case we impose an upstream boundary condition for the discharge as q ¼ 0:18 m2 =s and the downstream boundary condition for the water level is H0 ¼ 0:33 m. The final water level at time 200 s computed with mesh size Dx ¼ 0:125 m and residue for computation are compared with LLF method and displayed in Fig. 19. The numerical simulation of several rapidly varying two-dimensional flows are presented in the following sub-section. 4.3.2. 2D Partial dam-break problem This problem corresponds to a partial breach [39]. Consider a 2D partial dam-break problem with an unsymmetrical breach. The computational domain is defined in a 200 m long and 200 m wide channel. The non-symmetrical breach or sluice gates are 75-m wide, and the thickness is 10 m (see Fig. 20). A dam is assumed to fail instantaneously or the sluice gates are assumed to be opened instantly. The initial upstream water depth is h1 ¼ 10 m and downstream water depth is h0 ¼ 5 m. A grid of 40  40 cells is used for this problem. This test case is similar to the shock tube problem for the Euler equations. 1.8 1 LLF LLF 1.6 Exact Exact 0.9 1.4 1.2 0.8 1 0.8 0.7 0.6 0.6 0.4 0.2 0.5 0 0.4 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 Fig. 18. Dam-break flow with bottom topography on 100 grid points using first order accurate MOVERS-n at time 0.1 s. 1 786 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 0.45 2 10 LLF 0.4 0 10 0.35 0.3 10 LLF 0.25 EXACT 10 0.2 0.15 10 0.1 0.05 0 5 10 15 20 25 10 30 0 0.5 X 1 Iterations 1.5 2 4 x 10 Fig. 19. Transcritical steady flow with hydraulic jump on 120 points at t ¼ 200 s: first order accurate MOVERS-n compared with LLF. 30m 75m 95m 95m 200m Fig. 20. Problem sketch for partial dam breach. The numerical results displaying 3D views of the water surface elevation after dam failure and velocity contours at time t ¼ 7:2 s are presented in Fig. 21. At the instant the dam-breaks, water is released through the breach, forming a positive wave that propagates downstream and a negative wave that moves upstream. The results for steady oblique hydraulic jump compared with LLF scheme are depicted in Fig. 22. Fig. 23 presents the surface height variation for the oblique jump. We can see that all the flow features are resolved well. 4.4. Results with MOVERS-n for ideal magneto-hydrodynamic equations The one-dimensional ideal MHD equations are given by 3 qu1 2 2 1 2 7 6 qu1 þ p þ 2 B  Bx 7 6 7 6 7 6 qu1 u2  Bx By 7 6 7 6 q u u  B B 1 3 x z 7 o 6 7¼0 7þ 6 6 7 6 0 ot 6 Bx 7 ox 6 7 7 6 7 6 6 B 7 By u1  u2 Bx 7 6 6 y 7 7 6 7 6 7 6 B u  u B z 1 3 x 4 Bz 5 5 4   1 2 u1 E þ p þ 2 B  Bx ðu1 Bx þ u2 By þ u3 Bz Þ E q 3 6 qu1 7 7 6 7 6 6 qu2 7 7 6 6 o 6 qu3 7 7 2 2 ð47Þ S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 10 9 8 7 6 5 0 200 50 150 100 100 150 50 200 0 200 180 160 140 y(m) 120 100 80 60 40 20 0 0 50 100 150 200 x(m) Fig. 21. Solution of partial dam-break problem using first order accurate MOVERS-n at time 7.2 s with 40  40 points. LLF(1.0:0.1:2.1) 5 4 3 2 1 2 4 6 8 10 MOVERS-n(1.0:0.1:2.1) 5 4 3 2 1 2 4 6 8 10 Fig. 22. Solution of oblique jump problem using first order accurate MOVERS-n compared with LLF method. 787 788 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 2.5 2 1.5 1 0.5 6 15 4 10 2 5 0 0 Fig. 23. Solution of oblique jump problem using first order accurate MOVERS-n. Since the magnetic field in the flow direction is constant, both time and space derivatives of Bx are zero. The seven eigenvalues of the Jacobian matrix (see [40]) are k1 ¼ u  c f ; k2 ¼ u  cA ; k3 ¼ u  cs ; k5 ¼ u þ c f ; k6 ¼ u þ cA ; k7 ¼ u þ cs k4 ¼ u; ð48Þ Several standard problems in one and two dimensions are solved with the scheme MOVERS-n and results are presented here. 4.4.1. 1D Coplanar problem This 1D MHD problem with the initial Riemann conditions is given in [40]. The initial conditions ½q; u1 ; u2 ; p; By  on the leftðx 6 0Þ are [1, 0, 0, 1, 1] and on rightðx P 0Þ are [0.125, 0, 0, 0.1, 1]. Two hundred grid points are used in a domain of [50, 50]. The CFL number is chosen as 0.75 and the results given are at t ¼ 10 with Bx ¼ 0:75 and c ¼ 2. The results are shown for the domain [30:30] in Figs. 24 and 25 in comparison with LLF method and HLLC Riemann solver [41], well known for its accuracy. The results here are deliberately given for a very coarse grid for this problem to compare the schemes. MOVERS-n is seen to capture the discontinuities much better than LLF method while it improves even over the HLLC scheme, although being simpler. When the grids are made finer with 800 points, the results (shown in Fig. 26) are in close agreement with those of HLLC Riemann solver. 4.4.2. Riemann problem with seven discontinuities This is a one-dimensional problem which involves all the three components of velocity and magnetic field. It was suggested by Ryu et al. [42]. The spatial interval considered is (1, 1) on the x-axis. The initial conditions ½q; u1 ; u2 ; u3 ; p; By ; Bz  on the leftðx 6 0Þ are [1.08, 1.2, 0.01, 0.5, 0.95, 3.6, 2.0] and on rightðx P 0Þ are [1, 0, 0, 0, 1, 4, 2]. 1 1 HLLC HLLC 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 10 20 30 0 0 10 20 Fig. 24. Numerical simulation of ideal MHD problem (A) with first order accurate MOVERS-n and 200 points at time = 10. 30 789 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 1 1 HLLC HLLC 0.8 0.6 0.4 0.5 0.2 0 0 0 10 20 30 0 10 20 30 Fig. 25. Numerical simulation of ideal MHD problem (A) with first order accurate MOVERS-n and 200 points at time 10. Fig. 26. Numerical simulation of ideal MHD problem (A) with first order accurate MOVERS-n on 800 points at time 10. Fig. 27. Numerical simulation of ideal MHD problem (B) with first order accurate MOVERS-n at time 0.4 on 500 points. Fast shocks, rotational discontinuities and slow shocks propagate from each side of the contact discontinuity. The rotation of the magnetic field across the initial discontinuity generates two rotational discontinuities. All the above features are seen to be captured by the scheme using 500 points (Figs. 27 and 28) with the solution falling close to HLLC scheme but is much 790 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 28. Numerical simulation of ideal MHD problem (B) with first order accurate MOVERS-n at time 0.4 on 500 points. Fig. 29. Numerical simulation of ideal MHD problem (B) with first order accurate MOVERS-n at time 0.4 on 1500 points. better than LLF method (not shown). In this problem, due to the complex interaction of multiple waves, more amount of numerical diffusion is used for entropy fix. The MOVERS-n predictions merge with that of HLLC with fine grid of 1500 points (Fig. 29). 4.4.3. Orszag–Tang vortex problem This is a problem of evolution of a compressible vortex system [43,44] involving interaction between several shock waves traveling at different speeds. The computational domain is ½0; 2p  ½0; 2p with periodic boundary conditions in both x and y directions using a uniform grid of 288  288 points. The initial conditions are qðx; y; 0Þ ¼ c2 ; pðx; y; 0Þ ¼ c; v x ðx; y; 0Þ ¼ siny v y ðx; y; 0Þ ¼ sinx; Bx ðx; y; 0Þ ¼ siny; By ðx; y; 0Þ ¼ sin2x Value of c is taken as 5/3. The solutions for 1D ideal MHD equations shown in Figs. 24–28 are compared with Local Lax–Friedrichs Scheme (not shown) which is very dissipative and the HLLC Solver which is known to be very accurate. The solutions obtained with MOVERS-n are, as expected, comparable with the HLLC solutions although the former scheme is much simpler and free of Riemann solvers. The solutions for 2D-ideal MHD equations for the Orszag–Tang vortex system in Figs. 30–33 are very accurate and are comparable with reported results [44]. 4.5. Results with MOVERS-1 for Euler equations The solutions for 1D steady contact discontinuity problem, Sod’s shock tube test case, quasi-1D nozzle flow, supersonic flow over a step in a wind tunnel and hypersonic flow across half cylinder, obtained with MOVERS-1, are given in Figs. 34–37. S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 30. Numerical simulation of Orszag–Tang vortex system with first order accurate MOVERS-n at time = 1. Fig. 31. Numerical simulation of Orszag–Tang vortex system with first order accurate MOVERS-n at time = 1. Fig. 32. Numerical simulation of Orszag–Tang vortex system with first order accurate MOVERS-n at time = 3. 791 792 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 33. Numerical simulation of Orszag–Tang vortex system with first order accurate MOVERS-n at time = 3. a b Fig. 34. Density plot with first order accurate MOVERS-1 for: (a) steady 1D contact; (b) Sod’s 1D shock tube problem. a b 50 100 1.2 0.8 0.8 0.4 0.4 0 1 2 200 500 1.2 3 0 1 2 3 Fig. 35. Density plot with first order accurate MOVERS-1 for quasi-1D nozzle flow: (a) 50, 100 points; (b) 200, 500 points. The solutions are as good as those obtained with the n-wave variant or even better than some of them. As mentioned before, the choice of one of the three wave speeds (corresponding to the energy equation and which uses the maximum possible 793 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 36. Supersonic flow over forward-facing a step in channel: density contours with MOVERS-1 on 240  80 grid. 1.5 1.5 1.5 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.5 -1 -0.5 0 -1.5 -1 -0.5 0 -1.5 -1.5 -1 -0.5 0 Fig. 37. Hypersonic (M = 6) flow across a semi-cylinder: pressure contours with MOVERS-1 using unstructured mesh of 3897 points. information) is made based on reasons mentioned in Section 3.2. The shock solution of the quasi-1D nozzle flow consistently contains a maximum of only two points which is currently possible only by highly accurate Riemann solvers such as Roe scheme when only first order scheme is used. The solution for shear flow problem with 1-wave variant is almost identical to that of the n-wave variant (Fig. 7) and hence is not reproduced here. The scheme is used to solve benchmark problem of supersonic flow over a forward-facing step in a wind tunnel and the slip lines starting at the Mach stem from right above the step are clearly seen in Fig. 36. The computation with unstructured mesh is also found to be accurate as evident in the solution for hypersonic flow across a half cylinder presented in Fig. 37. It is found that the residue for the second order solution does not fall to very low values for this version of MOVERS. However, instead of calculating the wave speed at the interface directly with the immediate neighboring points as in Eq. (33), if the wave speeds are calculated at every cell center [18] as sj ¼ (F jþ1 F j1 U jþ1 U j1 if U j1 –U jþ1 kmin;j if U j1 ¼ U jþ1 ð49Þ sj ¼ signðsj Þkmax;j if jsj j P kmax;j ð50Þ sj ¼ signðsj Þkmin;j if jsj j 6 kmin;j ð51Þ 794 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 a b Fig. 38. Convergence history for MOVERS-1: (a) supersonic flow (M = 2) over a ramp in channel; (b) hypersonic flow (M = 6) across a semi-cylinder. ρ Fig. 39. Steady contact discontinuity in a shock tube with first order accurate MOVERS-L. a b ρ ρ Fig. 40. Quasi-1D nozzle flow with first order accurate MOVERS-L. S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 795 Fig. 41. Horizontal shear flow (Mach 3 flow over Mach 2 flow), Mach contours with first order accurate MOVERS-L on 40  20 grid. Fig. 42. Mach 2 flow over Ramp in channel: density contours with MOVERS-L. and the interface wave speed is then taken to be the maximum of left and right state values, the residual of the second order accurate scheme falls to much lower values. With such a change in computing the interface flux, the solutions agree closely with the parent methods MOVERS-(n,1) presented and hence only the convergence plots are given here in Fig. 38. 4.6. Results with MOVERS-L The results for MOVERS-L, a limiter based Central Rankine–Hugoniot Solver in which the dissipation switches from the one in the MOVERS-1(or n) to that in the LLF method, obtained with a minmod limiter, are presented in Figs. 39–43. Steady discontinuities are resolved accurately and all features of the flow are captured well, as can be seen from the results. Note that the high resolution is achieved here without the use of an entropy fix, based only on the flow variables. 5. MOVERS and upwind interpretation The results presented in the previous sections show that the Central Rankine–Hugoniot Solver developed in this work, MOVERS, is a very efficient scheme, capturing the steady discontinuities exactly when aligned with grid lines, and providing low diffusion solutions otherwise. The motivation with which this work started was to develop an accurate flow solver based on a simple central discretization scheme, taking the Local Lax–Friedrichs method as the building block. It is worth noting that MOVERS can also be interpreted as an upwind scheme. An upwind interpretation may give a different perspective for this new algorithm and may be useful in modifying the existing upwind codes in a simple way. In the following, an upwind interpretation for n-wave variant is presented. The cases of other variants can be treated in a similar fashion. Consider the case of a system of hyperbolic conservation laws in one dimension oU oF þ ¼0 ot ox ð52Þ Let us introduce a linearization of the flux vector as e þB F  AU ð53Þ e and B are constant matrices. We further simplify by assuming A e to be a diagonal matrix where A 2~ k1 60 6 e¼6 A 6 .. 4 . 0 0 ~k2 .. . 0  0 3 07 7 7 .. 7  . 5    ~kn  ð54Þ 796 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 Fig. 43. Transonic flow across NACA0012 airfoil: pressure coefficient distribution obtained with MOVERS-L using O-type mesh with 160  60 points. Then the system of conservation laws can be rewritten as e þ BÞ oU oU oð AU e oU ¼ 0 þ ¼ þA ot ox ot ox ð55Þ which is the approximate Riemann problem being solved at the cell interfaces, with the initial conditions as UL Uðx; t ¼ 0Þ ¼ ð56Þ UR Expanding, we obtain oU m ~ oU m þ km ¼ 0; ot ox m ¼ 1; . . . ; n ð57Þ The approximate Riemann problem consists of n linear convection equations. The n wave speeds need to be constants as per the linearization assumption made and this requirement can be satisfied by taking them as local averages (not necessarily arithmetic averages), across a cell interface. Now consider the flux Jacobian matrix for the original conservation laws, given by A¼ oF oU ð58Þ which has n eigenvalues, k1 ; . . . ; kn . If the physical situation presents a shock wave, then relevant wave speed will be the shock speed, defined by the Rankine–Hugoniot condition as DF ¼ sDU ð59Þ Since our aim is to enhance resolution of the shocks, we approximate the wave speeds by relevant speeds of the shock coming from the jump condition as ~km ¼ jsm;I j; m ¼ 1; . . . ; n ð60Þ leading to 2 js1;I j 0  6 0 js j  2;I 6 e ¼6 A .. 6 .. 4 . .  0 0 0 0 .. .    jsn;I j 3 7 7 7 7 5 ð61Þ Here, the subscript I denotes the cell interface. Using the notations DU ¼ U R  U L and DF  F R  F L , we can now use the approximate Flux Jacobian in the relation e DU DF ¼ A and the connection to MOVERS-n has now become apparent, if we use ð62Þ S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 sm;I ¼ DF m ; DU m m ¼ 1; . . . ; n 797 ð63Þ with the wave speed corrections introduced before. If the evaluated wave speeds are positive jsm;I j ¼ sm;I ; m ¼ 1; . . . ; n e¼A e and the cell interface flux becomes Then A FI ¼ 1 1e 1 1 ðF L þ F R Þ  AðU ðF L þ F R Þ  ðF R  F L Þ ¼ F L R  UL Þ ¼ 2 2 2 2 ð64Þ ð65Þ which is the upwind flux for positive wave speeds. If the evaluated wave speeds are negative jsm;I j ¼ sm;I ; m ¼ 1; . . . ; n ð66Þ Then e ¼ A e A ð67Þ and the cell interface flux becomes FI ¼ 1 1e 1 1 ðF L þ F R Þ þ AðU ðF L þ F R Þ þ ðF R  F L Þ ¼ F R R  UL Þ ¼ 2 2 2 2 ð68Þ which is again the upwind flux for negative wave speeds. Thus, MOVERS-n is also an upwind scheme, similar to the Flux Difference Splitting method of Roe, but the expressions are simpler and are not based on Roe averages. Then, this scheme can be interpreted as an upwind method with the wave speeds being evaluated based on the Rankine–Hugoniot jump condition, with wave speed corrections. The interpretation of other variants of the new algorithm can be done in a similar way, noting that MOVERS-L does not need an entropy fix and is a combination of LLF and either n-wave variant or 1-wave variant. It is also important to note that the upwind interpretation of MOVERS-1 is possible only if we take a very broad outlook, like treating the Local Lax–Friedrichs method also as an upwind scheme. To go beyond such interpretations, it is important to note that simplicity and accuracy in capturing flow features like discontinuities are worthy goals to strive for in developing new algorithms, and this is the motive behind the work presented in this work. With the above interpretation, the Central Rankine–Hugoniot Solvers developed in this work can be thought of as an advancement over both the existing central discretization methods (in capturing the steady discontinuities exactly) and upwind methods (being simpler and being at least equally accurate) and may also be thought of as uniting the two frameworks. 6. Conclusion An accurate numerical solver for hyperbolic partial differential equations is presented here. The new solver is in the framework of central schemes and is based on the direct discretization of Rankine–Hugoniot (jump) condition. The solver is simple and easily adoptable for any hyperbolic system. For hyperbolic vector conservation laws, three different variants of the new algorithm are proposed. MOVERS-n is based on the propagation of n-waves while MOVERS-1 is based on an approximation of a single wave, for the shock speed in the jump condition, to represent the coefficient of numerical dissipation. MOVERS-L is a hybrid limiter based solver smoothly varying between Local Lax–Friedrichs method and the method based on the jump condition, thus avoiding the need for an entropy fix. These different variants of the scheme are used to solve various bench-mark test cases of gas dynamics, shallow water flows and ideal magneto-hydrodynamics. The results show high accuracy in capturing flow features, with grid-aligned shock waves and contact discontinuities being resolved exactly. The variants of MOVERS can also be interpreted as upwind methods but are simpler, while retaining the accuracy in capturing discontinuities, and one of the variants avoids the use of any entropy fix due to the flow based adaptation of the numerical dissipation. References [1] C. Hirsch, Numerical Computation of Internal and External Flows: Fundamentals of Numerical Discretization, vol. 1, John Wiley & Sons, 1988. [2] C. Hirsch, Numerical Computation of Internal and External Flows: Computational Methods for Inviscid and Viscous Flows, vol. 2, John Wiley & Sons, 1990. [3] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, second ed., Springer, 1999. [4] Culbert B. Laney, Computational Gasdynamics, Cambridge University Press, 1998. [5] R.J. Leveque, Numerical Methods for Conservation Laws, Birkhäuser, 1990. [6] R.J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002. [7] P.L. Roe, Approximate Riemann solvers parameter vectors and difference schemes, J. Comput. Phys. 43 (1981) 357–372. [8] J.J. Quirk, A contribution to the great Riemann solver debate, Int. J. Numer. Methods Fluids 18 (1994) 555–574. [9] H. Nessyahu, E. Tadmor, Non-oscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys. 87 (1990) 408–463. [10] R.C. Swanson, E. Turkel, On central difference and upwind schemes, J. Comput. Phys. 101 (1992) 292–306. [11] R.C. Swanson, R. Radespiel, E. Turkel, On some numerical dissipation schemes, J. Comput. Phys. 147 (1998) 518–544. [12] Pierre Degond, Pierre-FranCok Peyrard, Giovanni Russo, Philippe Villedieu, Polynomial upwind schemes for hyperbolic systems, C.R. Acad. Sci. Paris, Serie I 328 (1999) 479–483. 798 S. Jaisankar, S.V. Raghurama Rao / Journal of Computational Physics 228 (2009) 770–798 [13] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, J. Comput. Phys. 160 (2000) 214–282. [14] P.D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Commun. Pure Appl. Math. 7 (1954) 159–193. [15] K.O. Friedrichs, Symmetric hyperbolic linear differential equations, Commun. Pure Appl. Math. 7 (1954) 345. [16] V.V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, J. Comput. Math. Phys., USSR 1 (1961) 267–279. [17] S. Jaisankar, S.V. Raghurama Rao, Diffusion regulation for Euler solvers, J. Comput. Phys. 221 (2007) 577–599. [18] S. Jaisankar, S.V. Raghurama Rao, A discrete Rankine–Hugoniot solver for hyperbolic conservation laws, in: Proceedings of the Seventh Asian Computational Fluid Dynamics Conference, Bangalore, India, November 2007. [19] P.D. Lax, B. Wendroff, Systems of conservation laws, Commun. Pure Appl. Math. 13 (1960) 217–237. [20] 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-81-1259, 1981. [21] Meng-Sing Liou, C.J. Steffen Jr., A new flux splitting scheme, J. Comput. Phys. 107 (1993) 23–39 (also NASA TM-104404, 1991). [22] Meng-Sing Liou, Ten years in the making – AUSM family, AIAA Paper No. AIAA-2001-2521, 2001 (also NASA/TM 2001-210977). [23] A. Jameson, Analysis and design of numerical schemes for gas dynamics. 2: Artificial diffusion and discrete shock structure, Int. J. Comput. Fluid Dynam. 5 (1995) 1–38. [24] 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, Int. J. Comput. Fluid Dynam. 4 (1995) 171–218. [25] S.V. Raghurama Rao, K. Balakrishna, An accurate shock capturing algorithm with a relaxation system (ASCARS) for hyperbolic conservation laws, AIAA Paper No. AIAA-2003-4115, 2003. [26] R. Courant, K.O. Friedrichs, Supersonic Flow and Shock Waves, Interscience Publishers, New York, 1948. p. 124. [27] J.C. Tannehil, D.A. Anderson, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, second ed., Taylor & Francis, 1997. [28] A. Harten, J.M. Hyman, Self-adjusting grid methods for one-dimensional hyperbolic conservation laws, J. Comput. Phys. 50 (1983) 235–269. [29] H.C. Yee, Construction of explicit and implicit symmetric TVD schemes and their applications, J. Comput. Phys. 68 (1987) 151–179. [30] S. Spekreijse, Multigrid solution of monotone second order discretizations of hyperbolic conservation laws, Math. Comput. 48 (179) (1987) 135–155. [31] John D. Anderson Jr., Computational Fluid Dynamics: The Basics with Applications, McGrawhill Inc., 1995. p. 356. [32] D. Levy, K.G. Powell, B. van Leer, Use of a rotated Riemann solver for two-dimensional Euler equations, J. Comput. Phys. 106 (1993) 201–214. [33] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115–173. [34] M. Manna, A three-dimensional high resolution upwind finite volume Euler solver, Technical Note 180, von Karman Institute for Fluid Dynamics, 1992. [35] Meng-Sing Liou, Mass flux schemes and connection to shock instability, J. Comput. Phys. 160 (2000) 623–648. [36] H. Viviand, Numerical solutions of two-dimensional reference test cases, in: Test Cases for Inviscid Flow Field Methods, AGARD 211, 1985. [37] Justin Hudson, Numerical techniques for shallow water equations, Numerical Analysis Report 2/99, University of Reading, 1999. [38] F. Jimenez, M. Hanif Chaudhry, Computation of supercritical free-surface flows, J. Hydraul. Eng. 144 (4) (1988). [39] R.J. Fennema, M.H. Chaudhry, Implicit methods for two-dimensional unsteady free surface flows, J. Hydraul. Res. 27 (3) (1989) 321–332. [40] M. Brio, C.C. Wu, An upwind differencing scheme for the equations of ideal magnetohydrodynamics, J. Comput. Phys. 75 (1988) 400–422. [41] K.F. Gurski, An HLLC-type approximate Riemann solver for ideal magnetohydrodynamics, SIAM J. Sci. Comput. 25 (6) (2004) 2165–2187. [42] D. Ryu, T. Jones, T. Frank, Numerical magnetohydrodynamics in astrophysics: algorithm and tests for multidimensional flow, Astrophys. J. 452 (1995) 785–796. [43] S.A. Orszag, C.M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, J. Fluid Mech. 90 (1979) 129. [44] J. Balbas, E. Tadmor, C.C. Wu, Non-oscillatory central schemes for one- and two-dimensional MHD equations I, J. Comput. Phys. 201 (20) (2004) 261– 285.