# A brief timeline of the significant research findings in Computational Fluid Dynamics (CFD) for flows with discontinuity (Shocks)

Shocks are familiar to us as waves that can pass through fluids at a finite speed. But when a shock is described as discontinuity of a flow, the idea becomes bit arcane. Discontinuity is nothing but an extremely thin region (a few nanometers) where the fluid properties change abruptly. This invisible wave is always associated with enormous powers to change flow properties of a certain volume in a blink of time. It can change a place of low temp. low pressure and zero speed into a region of high pressure and high temp. with a finite velocity. Shocks always pass at the speed of sound. This the speed at which the information also pass through the fluids. When an object is moving at a subsonic speed ( less than the speed of sound), the flow particles ahead of the object have the chance to know that the object is coming towards it and it change its position making the way for the object to pass easily. If enough time is not given to reach the information regarding the incoming object, then there arise discontinuities i.e. forms a shock wave around the object. This is an indispensable concept for the people who want to investigate and research (both analytically and experimentally) supersonic (higher than the speed of sound) flow. One who wants to study numerical techniques about compressible supersonic flow must have a robust and lucid idea about shock: how it propagate and why and when it is created. A chronology of the noteworthy research findings in the fields of CFD esp. shock capturing methods over the last seventy years is presented in this article. Exploring the earlier as well as the contemporary outcomes of the numerical schemes in an analogous manner helps to choose a suitable scheme for solving a certain fluid flow problem. And this timeline has been made in such a way that it represents various analogies among the schemes developed by the pioneers of CFD over the late twentieth century. Since capturing discontinuities i.e. the shocks in fluid flow are one of the most challenging and complex tasks in scientific computation, we concentrate on the schemes that can be used for capturing shock waves and other discontinuities. The effectiveness, as well as the drawbacks of the different numerical approaches, is illustrated in this paper from the viewpoint of their uses and applications.

Readers, who are not widely familiar with CFD and numerical analysis, are advised to take a look at the following reference before going any further.

De Souza, Althea. “How to–Understand Computational Fluid Dynamics Jargon.” (2005).

Introduction

In the mid-1940s, numerical approaches were first practiced as a tool for solving partial differential equations (PDEs) by John von Neumann [32, 33]. Over the past seventy years, scientific computation has appeared as a flexible tool to support theory and experiments, and numerical methods obtain a position at the heart of many of today’s advanced scientific computations for solving PDEs. A large portion of the theory of PDEs was developed in response to the models that originated in the physical sciences. The Laplace equation, the wave equation and the heat equation are ideal linear PDEs. Whereas governing equations such as the minimal surface equation, the Schrodinger equation, Euler equations and Navier–Stokes equations are canonical examples that are driving much of the study of nonlinear PDEs. In this paper, we provide comments and reviews on the development of numerical solver of the hyperbolic partial differential equations with an emphasis on shock capturing schemes.

Shock capturing numerical methods have seen revolutionary developments over the past 20 years. These are methods which deal with the numerical solutions of PDEs with discontinuous solutions. From the historical point of view, shock-capturing methods can be categorized into two general groups: (1) conventional methods i.e. low/first order schemes and (2) modern shock capturing methods i.e. high-resolution/higher-order schemes. Modern shock-capturing methods are generally upwind based in contrast to classical symmetric or central discretization. Upwind-type differencing schemes discretize hyperbolic partial differential equations by using differencing based on the direction determined by the sign of the characteristic speeds. Adaptive finite difference stencil is used in upwind schemes to find the direction of propagation of information in a flow field. On the other hand, symmetric or central schemes do not consider any information about the wave propagation in the discretization.

A stable calculation in presence of shock waves requires a certain amount of numerical dissipation, in order to avoid the formation of unphysical numerical oscillations. In the case of classical shock-capturing methods, numerical dissipation terms are usually linear and the same amount is uniformly applied at all grid points. Classical shock-capturing methods only exhibit accurate results in the case of a smooth and weak-shock solution, but when strong shock waves are present in the solution, non-linear instabilities and oscillations can arise across discontinuities. Modern shock-capturing methods have, however, a non-linear numerical dissipation, with an automatic feedback mechanism which regulates the extent of dissipation in any cell of the mesh, in accord to the gradients in the solution. These schemes have proven to be stable, non-oscillating and accurate even for problems containing strong shock waves.

Some of the well-known classical shock-capturing methods are as below which requires addition of artificial viscosity:

• MacCormack method
• Lax–Wendroff method
• Beam–Warming method

Examples of modern shock-capturing schemes include

• Higher-order total variation diminishing (TVD) schemes first proposed by Harten,
• Flux-corrected transport scheme introduced by Boris and Book,
• Monotonic Upstream-centered Schemes for Conservation Laws (MUSCL) based on Godunov approach and introduced by van Leer,
• Various Essentially Non-Oscillatory schemes (ENO) proposed by Harten et al., and
• Piecewise Parabolic Method (PPM) proposed by Woodward and Colella.

Another important class of high-resolution schemes belongs to the approximate Riemann solvers including:

• Roe Solver
• HLLE Solver
• HLLC Solver
• Rotated-hybrid Riemann solvers

The timeline

In computational fluid dynamics (CFD), the upwind scheme is an eminent category of the numerical schemes for solving PDEs. Historically, the origin of upwind methods can be traced back to the work of Courant, Isaacson, and Rees who proposed the CIR method [1].

In 1959, Godunov is credited with introducing the first exact Riemann solver [2] for the Euler equations, by extending the previous CIR (Courant-Isaacson-Rees) method to non-linear systems of hyperbolic conservation laws. This category, known as Godunov-type methods, is very efficient for modern shock-capturing methods for the following reasons: the ability to automatically capture clean, sharp, and correctly positioned discontinuities such as shock waves and contact surfaces; their accuracy throughout the flow field and their intrinsic robustness.

In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This second-order finite difference method was introduced by Robert W. MacCormack in 1969[30]. The MacCormack method is elegant and easy to understand and program. The MacCormack method is a variation of the two-step Lax–Wendroff scheme that contains Predictor and Corrector Step but is much simpler in an application. The MacCormack method is well suited for nonlinear equations (Inviscid Burgers equation, Euler equations, etc.) The order of differencing can be reversed for the time step (i.e., forward/backward followed by backward/forward). For nonlinear equations, this procedure provides the best results. Unlike first-order upwind scheme, the MacCormack does not introduce diffusive errors in the solution. However, it is known to introduce dispersive errors (Gibbs phenomenon) in the region where the gradient is high.

In 1973, a conservative, first order, non-oscillatory shock-capturing scheme known as Flux-corrected transport (FCT) is developed by Boris, D. L. Book et al. for solving Euler equations and other hyperbolic equations. An FCT algorithm consists of two stages (a transport stage and a flux-corrected anti-diffusion stage) which are known as predictor-corrector structure [3]. The correction step removes the large dissipative error made in the predictor step, thus expose a solution with second- or third-order accuracy. During the second step the corrective fluxes are compared to the provisional solution values and limited where necessary, in order to ensure that no new extrema will arise, nor existing extrema grow. The comparison step makes the overall method nonlinear: the coefficients in the scheme depend on the solution itself, even when applied to a linear equation.

Later, in 1979, the MUSCL scheme (Monotonic Upstream-Centered Scheme for Conservation Laws) is introduced in the study of partial differential equations in a finite volume method by B. van Leer. MUSCL can provide highly accurate numerical solutions for a given system, even in cases where the solutions exhibit shocks, discontinuities, or large gradients. [4]. In this scheme B. van Leer constructed one kind of high-order total variation diminishing (TVD) method to obtain second-order spatial accuracy. This technique supersedes the idea of the piecewise constant approximation of Godunov’s scheme. Fluxes are calculated at each cell boundaries by obtaining slope limited and reconstructed left and right states. These fluxes can, in turn, be used as input to a Riemann solver. Alternatively, the fluxes can be used in Riemann-solver-free schemes, such as the Kurganov and Tadmor scheme which is outlined later. MUSCL methods are generally second-order accurate in smooth regions (although they can be formulated for higher orders) and provide good resolution, monotonic solutions around discontinuities. This scheme is straight-forward to implement and computationally efficient.

In 1981, Roe introduced a Flux Differencing Scheme which is largely used, even now, because of its accuracy, quality and mathematical clarity [5]. Concerning the FDS schemes, they are based on the difference between the decomposition of fluxes, constructed on an approximated solution of the local Riemann problem between two adjacent states. The basic premise of this problem is that changes in a flow can be transmitted only through entropy waves and acoustic waves, and only at some given speeds, which represent the eigenvalues of the governing non-linear equation system. Unfortunately, Roe’s scheme admits rarefaction shocks that do not satisfy the entropy condition and sometimes gives rise to spurious solutions, such as carbuncle phenomena, odd–even decoupling. In upwind differences, we must, in the first place, establish which way the wind blows, more precisely, we must determine in which direction each of a variety of signals moves through the computational grid. For this purpose, a physical model of the interaction between computational cells is needed. One such model is the Riemann approach that is discussed neighboring cells interact through discrete, finite-amplitude waves. The nature, propagation speed, and amplitude of these waves are found by solving, exactly or approximately, Riemann’s initial-value problem for the discontinuity at the cell interface. The numerical technique of distinguishing between the influence of the forward- and the backward-moving waves is called flux-difference splitting, and the best and older examples of this kind of methods are the methods of Roe [5] and Osher [6]. The characteristics of approximate Riemann solvers, that were investigated by Osher and Roe, are used as a ‘building block’ for the construction of modern shock-capturing schemes.

Flux Vector Splitting (FVS) scheme first drew the attention of the researchers by the pioneering work of Steger and Warming [7]. Later, Van Leer’s flux vector splitting scheme [13] with the implicit relaxation algorithms became famous for its efficiency, simplicity and ability to capture the sharp shock. The Van Leer flux splitting method utilizes a different premise from the Roe’s technique. Since fluid flow contains both positive and negative wave speeds (so that eigenvalue information can pass both upstream and downstream), the flux can be split into two components that are the positive and the negative parts according to the sign of the eigenvalues of the coefficient matrix. Each was discretized using relatively upwind stencils to maintain stability and accuracy. There are many possible ways to split the flux term by basing the splitting on the eigenvalue structure or some similar representation of the flow:

• Riemann approach which was discussed before.
• Boltzmann approach

In Boltzmann approach, the interaction of neighboring cells is accomplished through the mixing of pseudo-particles that move in and out of each cell according to a given velocity distribution. In contrast, the FVS has proved to be a simple and useful technique for arriving at upwind differencing and is preeminently suited for use in implicit schemes. These approaches are considered as the first level of upwind schemes because these approaches achieve upwinding by splitting the flux vector. The identification of the upwind directions is done with less effort than in the Godunov-type methods. Unfortunately, the simplicity of these splittings come at the price of reduced accuracy due to numerical diffusion. Concerning the Steger-Warming and Van Leer schemes, it has been shown that they cannot resolve intermediate characteristic fields.

Harten–Lax–Van Leer (HLL) scheme is a direct approximation of the numerical flux to compute the Godunov flux [8]. HLL scheme is very efficient, robust, has an entropy satisfaction property, resolves isolated shock exactly and preserves positivity. Harten, Lax, and van Leer showed the way to construct a simple approximate Riemann solution which contains only one intermediate step. This construction assumes theoretical bounds on the smallest and largest signal velocities in the exact Riemann solution.

In the same year, Ami Harten introduced the concept of Total variation diminishing (TVD) [9]. TVD characteristics are especially attractive to capture sharper shock by avoiding spurious or non-physical oscillations when the variation of field variable is discontinuous. TVD scheme enables sharper shock predictions on coarse grids saving computation time by eliminating the necessity of fine grids and as the scheme preserves monotonicity there are no spurious oscillations in the solution thus resolves the issue of false shock predictions in other numerical schemes.

Essentially Non-Oscillatory (ENO), another scheme constructed by Harten et. al. in 1987[10], which yield uniform high order accuracy in smooth regions and resolve discontinuities in the derivatives sharply. An essentially non-oscillatory piecewise polynomial solution is constructed by averaging of approximate solution over each cell. Unlike standard finite difference methods this procedure uses an adaptive stencil of grid points and, consequently, the resulting schemes are highly nonlinear. In 1994, Liu, Osher, and Chan introduced a new version of ENO for a third order finite volume version which is called WENO. Both ENO and WENO schemes use the idea of adaptive stencils to automatically achieve high order accuracy and non-oscillatory property near discontinuities.

In 1988 Davis proposed a number of algorithms for obtaining the priori bounds that are applied on the smallest and largest signal velocities i.e. wave-speeds in the exact Riemann solution in HLL scheme [11]. In [11], he also used his approximate Riemann solutions with the approach of van Leer et al. to construct a second-order schemes.

The HLLE (Harten, Lax, van Leer and Einfeldt) solver, proposed by Einfeldt in the paper [12], is an approximate solution to the Riemann problem which is only based on the integral form of the conservation laws. The stability and robustness of the HLLE solver are closely linked with the signal velocities and a single central average state. However, this scheme is very efficient, robust, has an entropy satisfaction property, resolves isolated shock exactly and preserves positivity.

Later, several attempts had been made by E. F. Toro to understand and solve the phenomenon regarding contact discontinuity and region inside the two opposite acoustic waves.

In 1989, Toro modeled the constant-covolume equation of state in the solution of the Riemann problem [14]. He used a single algebraic (non-linear) equation to find the pressure between the acoustic waves. In the same year, he developed another algorithm where the intercell fluxes are defined by a weighted average [15] to obtain higher accuracy without solving ‘generalized’ Riemann problems or adding an anti-diffusive term. An advantage of the present method is its simplicity. It also has the potential for efficiency, because it is well suited to the use of approximations in the solution of the associated Riemann problem.

In 1991, He used another weighted average flux numerical method to solve the time-dependent Euler equations by using Riemann solver locally [16]. Here, the present linearized Riemann solver is tied with exact Riemann solver in an adaptive fashion for severe flow regimes.

In [17], three topics on modern shock capturing methods for the time-dependent Euler equations were addressed by Toro, mainly: a Weighted Average Flux Method; then an efficient, robust, non-iterative exact Riemann solver was presented; and an improved version of the Harten-Lax-van Leer Riemann solver. Also, a very simple linearized Riemann solver is presented together with a Riemann solver adaptation procedure.

In 1993, in order to combine the advantages FDS and FVS, Liou and Steffen (1993) developed a hybrid flux splitting method that is able to recognize contact waves [18]. This scheme termed as AUSM (advection upstream splitting method) has received much attention by the computational community. It is developed as a numerical inviscid flux function for solving a general system of conservation equations. It is based on the upwind concept and was motivated to provide an alternative approach to other upwind methods, such as the Godunov method, FDS method by Roe, FDS scheme by Solomon and Osher, FVS method by Van Leer, FVS by Steger and Warming. The AUSM first recognizes that the inviscid flux consists of two physically distinct parts, i.e., convective and pressure fluxes. The former is associated with the flow (advection) speed, while the latter with the acoustic speed; or respectively classified as the linear and nonlinear fields. Currently, the convective and pressure fluxes are formulated using the eigenvalues of the flux Jacobian matrices. This scheme later substantially improved in [19, 20] to yield a more accurate and robust version.

The HLLC (Harten-Lax-van Leer-Contact) solver was introduced by Toro in 1994 [21]. It restores the missing rarefaction wave by some estimates, like linearization, these can be simple but also more advanced exists like using the Roe average velocity for the middle wave speed. In 1994 Quirk proposed in his paper [22] that HLLC are quite robust and efficient but somewhat more diffusive. In addition, Quirk proposed a strategy to use combined fluxes so that a dissipative approach can be used in the shock regions.

The first WENO scheme is constructed in 1994 by Liu, Osher and Chan [23] for a third order finite volume version. Instead of choosing the “smoothest” stencil to pick one interpolating polynomial for the ENO reconstruction, they use a convex combination of all candidates to attain the essentially non-oscillatory property, while additionally gaining one order of improvement in accuracy. The resulting weighted ENO schemes are based on cell averages and a TVD Runge-Kutta time discretization. In 1996, third and fifth order finite difference WENO schemes in multi-space dimensions are generated by Jiang and Shu [24]. The key idea in WENO schemes is a linear combination of lower order fluxes based on local characteristic decompositions to find a higher order approximation and to evade spurious oscillatory. The key idea lies at the approximation level, where a nonlinear adaptive procedure is used to automatically choose the locally smoothest stencil, hence avoiding crossing discontinuities in the interpolation procedure as much as possible. The main plus point of this scheme is its proficiency to accomplish arbitrarily high order formal accuracy in smooth regions while maintaining stable, non-oscillatory and sharp discontinuity transitions.

In 2000, a Riemann free solver, second-order, high-resolution scheme is presented, known as Kurganov and Tadmor (KT) central scheme, which uses MUSCL reconstruction [25]. A precursor to the Kurganov and Tadmor (KT) central scheme is the Nessyahu and Tadmor (NT) central scheme [26]. KT Central Scheme is a fully discrete method that is straight forward to implement and can be used on scalar and vector problems and can be viewed as a modification of the Lax-Friedrichs scheme. The algorithm is based upon central differences with comparable performance to Riemann type solvers when used to obtain solutions for PDE’s describing systems that exhibit high-gradient phenomena. The KT scheme extends the NT scheme and has a smaller amount of numerical viscosity than the original NT scheme. It also has the added advantage that it can be implemented as either a fully discrete or semi-discrete scheme.

Kim et al. (2003) developed a carbuncle-free FDS scheme by employing the HLL-type splitting that introduces a multi-dimensional dissipation term to overcome shock-instability without tuning parameters while maintaining the accuracy of Roe’s scheme [27]. For the hydrodynamic case, latest research results showed the possibility to avoid the iterations to calculate the exact solution for the Euler equations [28].

Rotated-hybrid Riemann solvers were introduced by Nishikawa and Kitamura [29] in order to overcome the carbuncle problems of the Roe solver and the excessive diffusion of the HLLE solver at the same time. A robust and accurate Riemann solver was developed by combining the Roe solver and the HLLE/Rusanov solvers: in two orthogonal directions, the two Riemann solvers can be combined into a single Roe-type solver (the Roe solver with modified wave speeds). In particular, the one that is derived from the Roe and HLLE solvers, called Rotated-RHLL solver, is extremely robust (carbuncle-free for all possible test cases on both structured and unstructured grids) and accurate (as accurate as the Roe solver for the boundary layer calculation).

Concluding Remarks

The high–resolution extensions of classical methods usually acquire all the features of the first–order formulations. These extensions also adopt the aspects that impair the solution. Classical methods and their extensions are not satisfactory for accurate resolution in the regions of the steep gradient. Non-physical effects (such as soil the solution by smudging or spurious oscillations) are often produced in these schemes. Since linear methods cannot offer a non-oscillatory solution higher than first order, a number of techniques have been developed that largely overcome this problem over the last seven decades. To avoid spurious or non-physical oscillations where shocks are present, schemes that exhibit a Total Variation Diminishing (TVD) characteristic are especially attractive. Two high-resolution schemes that are proving to be particularly effective for shock capturing and avoiding unnecessary oscillations are MUSCL and the WENO method. Mainly, two type of approaches have been employed to capture discontinuities in these high-order methods: artificial viscosity and solution limiting. The  former  involves  user  specified  parameters,  and  latter  often  causes iterative  convergence  to  stall. Another desirable feature of any shock-capturing approach is accuracy-preservation away from discontinuities.  Although MUSCL methods provide good resolution and monotonic solutions around discontinuities, for problems comprising both shocks and complex smooth solution structure, WENO schemes can provide higher accuracy than second-order schemes along with good resolution around discontinuities.

Though the formulation of novel numerical procedures has significant impact on modern developments of the solver of solving fluid and heat-related problems, still it requires effort to find an approach which is robust, parameter-free, accuracy preserving, and convergent for steady flow problems. In recent decades, the parallel development of numerical methods is carried on with full attention. The rise of nonlinear stochastic PDEs with a wide range of applications, from financial models and data assimilation in atmospheric sciences to material sciences and biological models will, therefore, assimilate stochastic aspects and the numerical methods will be multifaceted intrinsically.

References

• [1] Courant, Richard; Isaacson, E; Rees, M. (1952). “On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Differences”. Comm. Pure Appl. Math. 5: 243..255. doi:10.1002/cpa.3160050303.
• [2] Godunov, S. K. (1959), “A difference scheme for numerical computation of discontinuous solution of hyperbolic equation”, Math. Sbornik 47: 271–306
• [3] Jay P. Boris and David L. Book, “Flux-corrected transport, I: SHASTA, a fluid transport algorithm that works”, J. Comput. Phys. 11, pp. 38 (1973).
• [4] van Leer, B. (1979), Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov’s Method, J. Com. Phys., 32, 101–136.)
• [5] Roe, P. L. (1981). Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of computational physics, 43(2), 357-372.
• [6] Osher, “Numerical solution of singular perturbation problems and hyperbolic systems of conservation laws”, North-Holland Mathematical Studies 47 (1981), 179
• [7] L. Steger and R. F. Warming, ‘Flux vector splitting of the inviscid gas dynamics equations with  application to finite-difference methods’, J. Comput. Phys., 40, 263-293 (1981)
• [8] Harten A, Lax PD, van Leer B (1983) On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review 25:35-61
• [9] Harten, Ami (1983), “High resolution schemes for hyperbolic conservation laws”, J. Comput. Phys. 49: 357–393, doi:10.1006/jcph.1997.57131
• [10] Harten, A., Engquist, B., Osher, S., & Chakravarthy, S. R. (1987). Uniformly high order accurate essentially non-oscillatory schemes, III. In Upwind and High-Resolution Schemes (pp. 218-290). Springer Berlin Heidelberg.
• [11] Davis SF (1988) Simplified second-order Godunov-type methods. SIAM J Sci and Star Comput 9:445
• [12] Einfeldt, B. (1988), “On Godunov-type methods for gas dynamics”, SIAM J. Numer. Anal. 25 (2): 294–318, doi:10.1137/0725021
• [13] K. Anderson, J. L. Thomas and B. van Leer (1986), “Comparison of finite volume flux vector splittings for the Euler equations”, AIAA Journal,24, 1453–1460.
• [14] Toro EF (1989a) A fast Riemann solver with constant covolume applied to the random choice method. Int J Numer Methods in Fluids 9:1145.
• [15] Toro EF (1989b) A weighted average flux method for hyperbolic conservation laws. Proc Royal Soc London A 423:401
• [16] Toro EF (1991) A linearized Riemann solver for the Euler equations of gas dynamics. Proc Roy Soc London A 434:683.
• [17] Toro EF (1992b) The weighted average flux method applied to the Euler equations. Phil Trans Roy Soc London A 341:499
• [18] Liou, M.-S. and Steffen, C., “A New Flux Splitting Scheme,” J. Comput. Phys., Vol. 107, 23-39, 1993.
• [19] Liou, M.-S., “A Sequel to AUSM: AUSM+” J. Comput. Phys., Vol. 129, 364-382, 1996.
• [20] Wada, Y. and Liou, M.-S., “An Accurate and Robust Flux Splitting Scheme for Shock and Contact Discontinuities,” SIAM J. Scientific Computing, Vol. 18, 633-657, 1997. 3
• [21] Toro, E. F.; Spruce, M.; Speares, W. (1994), “Restoration of the contact surface in the HLL-Riemann solver”, Shock Waves 4: 25–34, doi:10.1007/BF01414629
• [22] Quirk, J. J. (1994), “A contribution to the great Riemann solver debate”
• [23] Liu, Xu-Dong, Stanley Osher, and Tony Chan. “Weighted essentially non-oscillatory schemes.” Journal of computational physics 115.1 (1994): 200-212.
• [24] Jiang, Guang-Shan, and Chi-Wang Shu. “Efficient Implementation of Weighted ENO Schemes.” Journal of Computational Physics 126.1 (1996): 202-228.
• [25] Kurganov, Alexander and Eitan Tadmor (2000), New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations, J. Comp. Phys.,160, 241–282.
• [26] Nessyahu, Haim, and Eitan Tadmor. “Non-oscillatory central differencing for hyperbolic conservation laws.” Journal of computational physics 87.2 (1990): 408-463.
• [27] Kim, Sung-soo, et al. “Cures for the shock instability: development of a shock-stable Roe scheme.” Journal of Computational Physics 185.2 (2003): 342-374.
• [28] Wu, Y.Y.; Cheung, K.F. (2008), “Explicit solution to the exact Riemann problem and application in nonlinear shallow-water equations”, Int. J. Numer. Meth. Fluids 57 (11): 1649–1668, doi:10.1002/fld.1696
• [29] Nishikawa, H.; Kitamura, K. (2008), “Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers”, J. Comput. Phys. 227 (4): 2560–2581, doi:10.1016/j.jcp.2007.11.003
• [30] MacCormack, R. W., The Effect of viscosity in hypervelocity impact cratering, AIAA Paper, 69-354 (1969).
• [32] von Neumann, J.; Richtmyer, R. D. (March 1950). “A Method for the Numerical Calculation of Hydrodynamic Shocks”. Journal of Applied Physics 21 (3): 232–237. doi:10.1063/1.1699639.
• [33] Charney, J. G.; Fjörtoft, R.; Neumann, J. (1950). “Numerical Integration of the Barotropic Vorticity Equation”. Tellus 2 (4): 237–254. doi:10.1111/j.2153-3490.1950.tb00336.x.