MUSCL scheme

Last updated

In the study of partial differential equations, the MUSCL scheme is a finite volume method that can provide highly accurate numerical solutions for a given system, even in cases where the solutions exhibit shocks, discontinuities, or large gradients. MUSCL stands for Monotonic Upstream-centered Scheme for Conservation Laws (van Leer, 1979), and the term was introduced in a seminal paper by Bram van Leer (van Leer, 1979). In this paper he constructed the first high-order, total variation diminishing (TVD) scheme where he obtained second order spatial accuracy.

Contents

The idea is to replace the piecewise constant approximation of Godunov's scheme by reconstructed states, derived from cell-averaged states obtained from the previous time-step. For each cell, slope limited, reconstructed left and right states are obtained and used to calculate fluxes at the cell boundaries (edges). These fluxes can, in turn, be used as input to a Riemann solver , following which the solutions are averaged and used to advance the solution in time. Alternatively, the fluxes can be used in Riemann-solver-free schemes, which are basically Rusanov-like schemes.

Linear reconstruction

1D advective equation
u
t
+
u
x
=
0
{\displaystyle u_{t}+u_{x}=0}
, with step wave propagating to the right. Shows the analytical solution along with a simulation based upon a first order upwind spatial discretization scheme. StepFirstOrdUpwind.png
1D advective equation , with step wave propagating to the right. Shows the analytical solution along with a simulation based upon a first order upwind spatial discretization scheme.

We will consider the fundamentals of the MUSCL scheme by considering the following simple first-order, scalar, 1D system, which is assumed to have a wave propagating in the positive direction,

Where represents a state variable and represents a flux variable.

The basic scheme of Godunov uses piecewise constant approximations for each cell, and results in a first-order upwind discretisation of the above problem with cell centres indexed as . A semi-discrete scheme can be defined as follows,

This basic scheme is not able to handle shocks or sharp discontinuities as they tend to become smeared. An example of this effect is shown in the diagram opposite, which illustrates a 1D advective equation with a step wave propagating to the right. The simulation was carried out with a mesh of 200 cells and used a 4th order Runge–Kutta time integrator (RK4).

To provide higher resolution of discontinuities, Godunov's scheme can be extended to use piecewise linear approximations of each cell, which results in a central difference scheme that is second-order accurate in space. The piecewise linear approximations are obtained from

Thus, evaluating fluxes at the cell edges we get the following semi-discrete scheme

1D advective equation
u
t
+
u
x
=
0
{\displaystyle u_{t}+u_{x}=0}
, with step wave propagating to the right. Shows the analytical solution along with a simulation based upon a second order, central difference spatial discretization scheme. StepSecOrdCentralDifference.png
1D advective equation , with step wave propagating to the right. Shows the analytical solution along with a simulation based upon a second order, central difference spatial discretization scheme.

where and are the piecewise approximate values of cell edge variables, i.e.,

Although the above second-order scheme provides greater accuracy for smooth solutions, it is not a total variation diminishing (TVD) scheme and introduces spurious oscillations into the solution where discontinuities or shocks are present. An example of this effect is shown in the diagram opposite, which illustrates a 1D advective equation , with a step wave propagating to the right. This loss of accuracy is to be expected due to Godunov's theorem. The simulation was carried out with a mesh of 200 cells and used RK4 for time integration.

An example of MUSCL type left and right state linear-extrapolation. LinExtrap.jpg
An example of MUSCL type left and right state linear-extrapolation.

MUSCL based numerical schemes extend the idea of using a linear piecewise approximation to each cell by using slope limited left and right extrapolated states. This results in the following high resolution, TVD discretisation scheme,

Which, alternatively, can be written in the more succinct form,

The numerical fluxes correspond to a nonlinear combination of first and second-order approximations to the continuous flux function.

The symbols and represent scheme dependent functions (of the limited extrapolated cell edge variables), i.e.,

where, using downwind slopes:

and

The function is a limiter function that limits the slope of the piecewise approximations to ensure the solution is TVD, thereby avoiding the spurious oscillations that would otherwise occur around discontinuities or shocks - see Flux limiter section. The limiter is equal to zero when and is equal to unity when . Thus, the accuracy of a TVD discretization degrades to first order at local extrema, but tends to second order over smooth parts of the domain.

The algorithm is straight forward to implement. Once a suitable scheme for has been chosen, such as the Kurganov and Tadmor scheme (see below), the solution can proceed using standard numerical integration techniques.

Kurganov and Tadmor central scheme

A precursor to the Kurganov and Tadmor (KT) central scheme, (Kurganov and Tadmor, 2000), is the Nessyahu and Tadmor (NT) a staggered central scheme, (Nessyahu and Tadmor, 1990). It is a Riemann-solver-free, second-order, high-resolution scheme that uses MUSCL reconstruction. It 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 Rusanov flux (also called the local Lax-Friedrichs flux) supplemented with high order reconstructions. 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. Here we consider the semi-discrete scheme.

The calculation is shown below:

1D advective equation
u
t
+
u
x
=
0
{\displaystyle u_{t}+u_{x}=0}
, with step wave propagating to the right. Shows the analytical solution along with a simulation based upon the Kurganov and Tadmor central scheme with SuperBee limiter. StepKTsuperbee.png
1D advective equation , with step wave propagating to the right. Shows the analytical solution along with a simulation based upon the Kurganov and Tadmor central scheme with SuperBee limiter.

Where the local propagation speed, , is the maximum absolute value of the eigenvalue of the Jacobian of over cells given by

and represents the spectral radius of

Beyond these CFL related speeds, no characteristic information is required.

The above flux calculation is most frequently called Lax-Friedrichs flux (though it's worth mentioning that such flux expression does not appear in Lax, 1954 but rather on Rusanov, 1961).

An example of the effectiveness of using a high resolution scheme is shown in the diagram opposite, which illustrates the 1D advective equation , with a step wave propagating to the right. The simulation was carried out on a mesh of 200 cells, using the Kurganov and Tadmor central scheme with Superbee limiter and used RK-4 for time integration. This simulation result contrasts extremely well against the above first-order upwind and second-order central difference results shown above. This scheme also provides good results when applied to sets of equations - see results below for this scheme applied to the Euler equations. However, care has to be taken in choosing an appropriate limiter because, for example, the Superbee limiter can cause unrealistic sharpening for some smooth waves.

The scheme can readily include diffusion terms, if they are present. For example, if the above 1D scalar problem is extended to include a diffusion term, we get

for which Kurganov and Tadmor propose the following central difference approximation,

Where,

Full details of the algorithm (full and semi-discrete versions) and its derivation can be found in the original paper (Kurganov and Tadmor, 2000), along with a number of 1D and 2D examples. Additional information is also available in the earlier related paper by Nessyahu and Tadmor (1990).

Note: This scheme was originally presented by Kurganov and Tadmor as a 2nd order scheme based upon linear extrapolation. A later paper (Kurganov and Levy, 2000) demonstrates that it can also form the basis of a third order scheme. A 1D advective example and an Euler equation example of their scheme, using parabolic reconstruction (3rd order), are shown in the parabolic reconstruction and Euler equation sections below.

Piecewise parabolic reconstruction

An example of MUSCL type state parabolic-reconstruction. ParabolicExtrap.jpg
An example of MUSCL type state parabolic-reconstruction.

It is possible to extend the idea of linear-extrapolation to higher order reconstruction, and an example is shown in the diagram opposite. However, for this case the left and right states are estimated by interpolation of a second-order, upwind biased, difference equation. This results in a parabolic reconstruction scheme that is third-order accurate in space.

We follow the approach of Kermani (Kermani, et al., 2003), and present a third-order upwind biased scheme, where the symbols and again represent scheme dependent functions (of the limited reconstructed cell edge variables). But for this case they are based upon parabolically reconstructed states, i.e.,

and

1D advective equation
u
t
+
u
x
=
0
{\displaystyle u_{t}+u_{x}=0}
, with step wave propagating to the right. Shows the analytical solution along with a simulation based upon the Kurganov and Tadmor Central Scheme with parabolic reconstruction and van Albada limiter. StepParabolicKTalbada.png
1D advective equation , with step wave propagating to the right. Shows the analytical solution along with a simulation based upon the Kurganov and Tadmor Central Scheme with parabolic reconstruction and van Albada limiter.

Where = 1/3 and,

and the limiter function , is the same as above.

Parabolic reconstruction is straight forward to implement and can be used with the Kurganov and Tadmor scheme in lieu of the linear extrapolation shown above. This has the effect of raising the spatial solution of the KT scheme to 3rd order. It performs well when solving the Euler equations, see below. This increase in spatial order has certain advantages over 2nd order schemes for smooth solutions, however, for shocks it is more dissipative - compare diagram opposite with above solution obtained using the KT algorithm with linear extrapolation and Superbee limiter. This simulation was carried out on a mesh of 200 cells using the same KT algorithm but with parabolic reconstruction. Time integration was by RK-4, and the alternative form of van Albada limiter, , was used to avoid spurious oscillations.

Example: 1D Euler equations

For simplicity we consider the 1D case without heat transfer and without body force. Therefore, in conservation vector form, the general Euler equations reduce to

where

and where is a vector of states and is a vector of fluxes.

The equations above represent conservation of mass, momentum, and energy. There are thus three equations and four unknowns, (density) (fluid velocity), (pressure) and (total energy). The total energy is given by,

where represents specific internal energy.

In order to close the system an equation of state is required. One that suits our purpose is

where is equal to the ratio of specific heats for the fluid.

We can now proceed, as shown above in the simple 1D example, by obtaining the left and right extrapolated states for each state variable. Thus, for density we obtain

where

Similarly, for momentum , and total energy . Velocity , is calculated from momentum, and pressure , is calculated from the equation of state.

Having obtained the limited extrapolated states, we then proceed to construct the edge fluxes using these values. With the edge fluxes known, we can now construct the semi-discrete scheme, i.e.,

The solution can now proceed by integration using standard numerical techniques.

The above illustrates the basic idea of the MUSCL scheme. However, for a practical solution to the Euler equations, a suitable scheme (such as the above KT scheme), also has to be chosen in order to define the function .

High resolution simulation of Euler equations based on G A Sod's 'Shock Tube' problem. Shows the analytical solutions along with simulated (2nd order) solutions based upon the Kuganov and Tadmor Central Scheme with Linear Extrapolation and Ospre limiter. SodProbKTospre.png
High resolution simulation of Euler equations based on G A Sod's 'Shock Tube' problem. Shows the analytical solutions along with simulated (2nd order) solutions based upon the Kuganov and Tadmor Central Scheme with Linear Extrapolation and Ospre limiter.

The diagram opposite shows a 2nd order solution to G A Sod's shock tube problem (Sod, 1978) using the above high resolution Kurganov and Tadmor Central Scheme (KT) with Linear Extrapolation and Ospre limiter. This illustrates clearly demonstrates the effectiveness of the MUSCL approach to solving the Euler equations. The simulation was carried out on a mesh of 200 cells using Matlab code (Wesseling, 2001), adapted to use the KT algorithm and Ospre limiter. Time integration was performed by a 4th order SHK (equivalent performance to RK-4) integrator. The following initial conditions (SI units) were used:

High resolution simulation of Euler equations based on G A Sod's 'Shock Tube' problem - SI units. Shows the analytical solutions along with simulated (3rd order) solutions based upon the Kurganov and Tadmor Central Scheme with parabolic reconstruction and van Albada limiter. SodProbParabolicKTalbada.png
High resolution simulation of Euler equations based on G A Sod's 'Shock Tube' problem - SI units. Shows the analytical solutions along with simulated (3rd order) solutions based upon the Kurganov and Tadmor Central Scheme with parabolic reconstruction and van Albada limiter.

The diagram opposite shows a 3rd order solution to G A Sod's shock tube problem (Sod, 1978) using the above high resolution Kurganov and Tadmor Central Scheme (KT) but with parabolic reconstruction and van Albada limiter. This again illustrates the effectiveness of the MUSCL approach to solving the Euler equations. The simulation was carried out on a mesh of 200 cells using Matlab code (Wesseling, 2001), adapted to use the KT algorithm with Parabolic Extrapolation and van Albada limiter. The alternative form of van Albada limiter, , was used to avoid spurious oscillations. Time integration was performed by a 4th order SHK integrator. The same initial conditions were used.

Various other high resolution schemes have been developed that solve the Euler equations with good accuracy. Examples of such schemes are,

More information on these and other methods can be found in the references below. An open source implementation of the Kurganov and Tadmor central scheme can be found in the external links below.

See also

Related Research Articles

<span class="mw-page-title-main">Laplace's equation</span> Second-order partial differential equation

In mathematics and physics, Laplace's equation is a second-order partial differential equation named after Pierre-Simon Laplace, who first studied its properties. This is often written as

<span class="mw-page-title-main">Stress–energy tensor</span> Tensor describing energy momentum density in spacetime

The stress–energy tensor, sometimes called the stress–energy–momentum tensor or the energy–momentum tensor, is a tensor physical quantity that describes the density and flux of energy and momentum in spacetime, generalizing the stress tensor of Newtonian physics. It is an attribute of matter, radiation, and non-gravitational force fields. This density and flux of energy and momentum are the sources of the gravitational field in the Einstein field equations of general relativity, just as mass density is the source of such a field in Newtonian gravity.

Linear elasticity is a mathematical model of how solid objects deform and become internally stressed due to prescribed loading conditions. It is a simplification of the more general nonlinear theory of elasticity and a branch of continuum mechanics.

In the calculus of variations, a field of mathematical analysis, the functional derivative relates a change in a functional to a change in a function on which the functional depends.

<span class="mw-page-title-main">Euler equations (fluid dynamics)</span> Set of quasilinear hyperbolic equations governing adiabatic and inviscid flow

In fluid dynamics, the Euler equations are a set of quasilinear partial differential equations governing adiabatic and inviscid flow. They are named after Leonhard Euler. In particular, they correspond to the Navier–Stokes equations with zero viscosity and zero thermal conductivity.

In theoretical physics, the Batalin–Vilkovisky (BV) formalism was developed as a method for determining the ghost structure for Lagrangian gauge theories, such as gravity and supergravity, whose corresponding Hamiltonian formulation has constraints not related to a Lie algebra. The BV formalism, based on an action that contains both fields and "antifields", can be thought of as a vast generalization of the original BRST formalism for pure Yang–Mills theory to an arbitrary Lagrangian gauge theory. Other names for the Batalin–Vilkovisky formalism are field-antifield formalism, Lagrangian BRST formalism, or BV–BRST formalism. It should not be confused with the Batalin–Fradkin–Vilkovisky (BFV) formalism, which is the Hamiltonian counterpart.

<span class="mw-page-title-main">Large eddy simulation</span> Mathematical model for turbulence

Large eddy simulation (LES) is a mathematical model for turbulence used in computational fluid dynamics. It was initially proposed in 1963 by Joseph Smagorinsky to simulate atmospheric air currents, and first explored by Deardorff (1970). LES is currently applied in a wide variety of engineering applications, including combustion, acoustics, and simulations of the atmospheric boundary layer.

In physics and fluid mechanics, a Blasius boundary layer describes the steady two-dimensional laminar boundary layer that forms on a semi-infinite plate which is held parallel to a constant unidirectional flow. Falkner and Skan later generalized Blasius' solution to wedge flow, i.e. flows in which the plate is not parallel to the flow.

In numerical methods, total variation diminishing (TVD) is a property of certain discretization schemes used to solve hyperbolic partial differential equations. The most notable application of this method is in computational fluid dynamics. The concept of TVD was introduced by Ami Harten.

The Newman–Penrose (NP) formalism is a set of notation developed by Ezra T. Newman and Roger Penrose for general relativity (GR). Their notation is an effort to treat general relativity in terms of spinor notation, which introduces complex forms of the usual variables used in GR. The NP formalism is itself a special case of the tetrad formalism, where the tensors of the theory are projected onto a complete vector basis at each point in spacetime. Usually this vector basis is chosen to reflect some symmetry of the spacetime, leading to simplified expressions for physical observables. In the case of the NP formalism, the vector basis chosen is a null tetrad: a set of four null vectors—two real, and a complex-conjugate pair. The two real members often asymptotically point radially inward and radially outward, and the formalism is well adapted to treatment of the propagation of radiation in curved spacetime. The Weyl scalars, derived from the Weyl tensor, are often used. In particular, it can be shown that one of these scalars— in the appropriate frame—encodes the outgoing gravitational radiation of an asymptotically flat system.

In mathematics, the cylindrical harmonics are a set of linearly independent functions that are solutions to Laplace's differential equation, , expressed in cylindrical coordinates, ρ (radial coordinate), φ (polar angle), and z (height). Each function Vn(k) is the product of three terms, each depending on one coordinate alone. The ρ-dependent term is given by Bessel functions (which occasionally are also called cylindrical harmonics).

The Cauchy momentum equation is a vector partial differential equation put forth by Cauchy that describes the non-relativistic momentum transport in any continuum.

The hybrid difference scheme is a method used in the numerical solution for convection–diffusion problems. It was introduced by Spalding (1970). It is a combination of central difference scheme and upwind difference scheme as it exploits the favorable properties of both of these schemes.

An axial fan is a type of fan that causes gas to flow through it in an axial direction, parallel to the shaft about which the blades rotate. The flow is axial at entry and exit. The fan is designed to produce a pressure difference, and hence force, to cause a flow through the fan. Factors which determine the performance of the fan include the number and shape of the blades. Fans have many applications including in wind tunnels and cooling towers. Design parameters include power, flow rate, pressure rise and efficiency.

Lagrangian field theory is a formalism in classical field theory. It is the field-theoretic analogue of Lagrangian mechanics. Lagrangian mechanics is used to analyze the motion of a system of discrete particles each with a finite number of degrees of freedom. Lagrangian field theory applies to continua and fields, which have an infinite number of degrees of freedom.

<span class="mw-page-title-main">Central differencing scheme</span>

In applied mathematics, the central differencing scheme is a finite difference method that optimizes the approximation for the differential operator in the central node of the considered patch and provides numerical solutions to differential equations. It is one of the schemes used to solve the integrated convection–diffusion equation and to calculate the transported property Φ at the e and w faces, where e and w are short for east and west. The method's advantages are that it is easy to understand and implement, at least for simple material relations; and that its convergence rate is faster than some other finite differencing methods, such as forward and backward differencing. The right side of the convection-diffusion equation, which basically highlights the diffusion terms, can be represented using central difference approximation. To simplify the solution and analysis, linear interpolation can be used logically to compute the cell face values for the left side of this equation, which is nothing but the convective terms. Therefore, cell face values of property for a uniform grid can be written as:

Unsteady flows are characterized as flows in which the properties of the fluid are time dependent. It gets reflected in the governing equations as the time derivative of the properties are absent. For Studying Finite-volume method for unsteady flow there is some governing equations >

The upwind differencing scheme is a method used in numerical methods in computational fluid dynamics for convection–diffusion problems. This scheme is specific for Peclet number greater than 2 or less than −2

In fluid mechanics, fluid flow through porous media is the manner in which fluids behave when flowing through a porous medium, for example sponge or wood, or when filtering water using sand or another porous material. As commonly observed, some fluid flows through the media while some mass of the fluid is stored in the pores present in the media.

In astrophysics, the Chandrasekhar virial equations are a hierarchy of moment equations of the Euler equations, developed by the Indian American astrophysicist Subrahmanyan Chandrasekhar, and the physicist Enrico Fermi and Norman R. Lebovitz.

References

    Further reading