Sand box Approximation Schemes

From CFD-Wiki

This content is retrieved from Wikipedia and will be modified

Please excuse any mistakes "Michail 13:19, 7 November 2011 (MST)"

S.K. Godunov's theorem

In numerical analysis and computational fluid dynamics, Godunov's theorem — also known as Godunov's order barrier theorem — is a mathematical theorem important in the development of the theory of high resolution schemes for the numerical solution of partial differential equations.

The theorem states that:

Linear numerical schemes for solving partial differential equations (PDE's), having the property of not generating new extrema (monotone scheme), can be at most first-order accurate.

Professor Sergei K. Godunov originally proved the theorem as a Ph.D. student at Moscow State University. It is his most influential work in the area of applied and numerical mathematics and has had a major impact on science and engineering, particularly in the development of methodologies used in computational fluid dynamics (CFD) and other computational fields. One of his major contributions was to prove the theorem (Godunov, 1954; Godunov, 1959), that bears his name.

The theorem

We generally follow Wesseling (2001).

Aside

Assume a continuum problem described by a PDE is to be computed using a numerical scheme based upon a uniform computational grid and a one-step, constant step-size, M grid point, integration algorithm, either implicit or explicit. Then if $x_{j} = j\,\Delta x \$ and $t^{n} = n\,\Delta t \$ , such a scheme can be described by

$\sum\limits_m^{M} {\beta _m } \varphi _{j + m}^{n + 1} = \sum\limits_m^{M} {\alpha _m \varphi _{j + m}^n }. \quad \quad ( 1)$

It is assumed that $\beta _m \$ determines $\varphi _j^{n + 1} \$ uniquely. Now, since the above equation represents a linear relationship between $\varphi _j^{n } \$ and $\varphi _j^{n + 1} \$ we can perform a linear transformation to obtain the following equivalent form,

$\varphi _j^{n + 1} = \sum\limits_m^{M} {\gamma _m \varphi _{j + m}^n }. \quad \quad ( 2)$

Theorem 1: Monotonicity preserving

The above scheme of equation (2) is monotonicity preserving if and only if

$\gamma _m \ge 0,\quad \forall m . \quad \quad ( 3)$

Proof - Godunov (1959)

Case 1: (sufficient condition)

Assume (3) applies and that $\varphi _j^n \$ is monotonically increasing with $j \$ .

Then, because $\varphi _j^n \le \varphi _{j + 1}^n \le \cdots \le \varphi _{j + m}^n$ it therefore follows that $\varphi _j^{n + 1} \le \varphi _{j + 1}^{n + 1} \le \cdots \le \varphi _{j + m}^{n + 1} \$ because

$\varphi _j^{n + 1} - \varphi _{j - 1}^{n + 1} = \sum\limits_m^{M} {\gamma _m \left( {\varphi _{j + m}^n - \varphi _{j + m - 1}^n } \right)} \ge 0 . \quad \quad ( 4)$

This means that monotonicity is preserved for this case.

Case 2: (necessary condition)

For the same monotonically increasing $\varphi_j^n \quad$ , assume that $\gamma _p^{} < 0 \$ for some $p \$ and choose

$\varphi _i^n = 0, \quad i < k;\quad \varphi _i^n = 1, \quad i \ge k . \quad \quad ( 5)$

Then from equation (2) we get

$\varphi _j^{n + 1} - \varphi _{j-1}^{n+1} = \sum\limits_m^M {\gamma _m } \left( {\varphi _{j + m}^{n} - \varphi _{j + m - 1}^{n} } \right) = \left\{ {\begin{array}{*{20}c} {0,} & {\left[ {j + m \ne k} \right]} \\ {\gamma _m ,} & {\left[ {j + m = k} \right]} \\ \end{array}} \right . \quad \quad ( 6)$

Now choose $j=k-p \$ , to give

$\varphi _{k-p}^{n + 1} - \varphi _{k-p-1}^{n + 1} = {\gamma _p \left( {\varphi _{k}^n - \varphi _{k - 1}^n } \right)} < 0 , \quad \quad ( 7)$

which implies that $\varphi _j^{n + 1} \$ is NOT increasing, and we have a contradiction. Thus, monotonicity is NOT preserved for $\gamma _p < 0 \$ , which completes the proof.

Theorem 2: Godunov’s Order Barrier Theorem

Linear one-step second-order accurate numerical schemes for the convection equation

${{\partial \varphi } \over {\partial t}} + c{ { \partial \varphi } \over {\partial x}} = 0 , \quad t > 0, \quad x \in \mathbb{R} \quad \quad (10)$

cannot be monotonicity preserving unless

$\sigma = \left| c \right|{{\Delta t} \over {\Delta x}} \in \mathbb{ N} , \quad \quad (11)$

where $\sigma \$ is the signed Courant–Friedrichs–Lewy condition (CFL) number.

Proof - Godunov (1959)

Assume a numerical scheme of the form described by equation (2) and choose

$\varphi \left( {0,x} \right) = \left( {{x \over {\Delta x}} - {1 \over 2}} \right)^2 - {1 \over 4}, \quad \varphi _j^0 = \left( {j - {1 \over 2}} \right)^2 - {1 \over 4} . \quad \quad (12)$

The exact solution is

$\varphi \left( {t,x} \right) = \left( {{{x - ct} \over {\Delta x}} - {1 \over 2}} \right)^2 - {1 \over 4} . \quad \quad (13)$

If we assume the scheme to be at least second-order accurate, it should produce the following solution exactly

$\varphi _j^1 = \left( {j - \sigma - {1 \over 2}} \right)^2 - {1 \over 4}, \quad \varphi _j^0 = \left( {j - {1 \over 2}} \right)^2 - {1 \over 4}. \quad \quad (14)$

Substituting into equation (2) gives:

$\left( {j - \sigma - {1 \over 2}} \right)^2 - {1 \over 4} = \sum\limits_m^{M} {\gamma _m \left\{ {\left( {j + m - {1 \over 2}} \right)^2 - {1 \over 4}} \right\}}. \quad \quad (15)$

Suppose that the scheme IS monotonicity preserving, then according to the theorem 1 above, $\gamma _m \ge 0 \$ .

Now, it is clear from equation (15) that

$\left( {j - \sigma - {1 \over 2}} \right)^2 - {1 \over 4} \ge 0, \quad \forall j . \quad \quad (16)$

Assume $\sigma > 0, \quad \sigma \notin \mathbb{ N} \$ and choose $j \$ such that $j > \sigma > \left( j - 1 \right) \$ . This implies that $\left( {j - \sigma } \right) > 0 \$ and $\left( {j - \sigma - 1} \right) < 0 \$ .

It therefore follows that,

$\left( {j - \sigma - {1 \over 2}} \right)^2 - {1 \over 4} = \left( j - \sigma \right) \left(j - \sigma - 1 \right) < 0, \quad \quad (17)$

which contradicts equation (16) and completes the proof.

The exceptional situation whereby $\sigma = \left| c \right|{{\Delta t} \over {\Delta x}} \in \mathbb{N} \$ is only of theoretical interest, since this cannot be realised with variable coefficients. Also, integer CFL numbers greater than unity would not be feasible for practical problems.

References

Godunov, Sergei K. (1954), Ph.D. Dissertation: Different Methods for Shock Waves, Moscow State University.
Godunov, Sergei K. (1959), A Difference Scheme for Numerical Solution of Discontinuous Solution of Hydrodynamic Equations, Math. Sbornik, 47, 271-306, translated US Joint Publ. Res. Service, JPRS 7226, 1969.
Wesseling, Pieter (2001), Principles of Computational Fluid Dynamics, Springer-Verlag.

Total variation diminishing

In numerical methods, total variation diminishing (TVD) is a property of certain discretization schemes used to solve hyperbolic partial differential equations. The concept of TVD was introduced by Ami Harten.^[1]

In systems described by partial differential equations, such as the following hyperbolic advection equation,

$\frac{\part u}{\part t} + a\frac{\part u}{\part x} = 0,$

the total variation (TV) is given by,

$TV = \int \left| \frac{\part u}{\part x} \right| dx ,$

and the total variation for the discrete case is,

$TV = \sum_j \left| u_{j+1} - u_j \right| .$

A numerical method is said to be total variation diminishing (TVD) if,

$TV \left( u^{n+1}\right) \leq TV \left( u^{n}\right) .$

A system is said to be monotonicity preserving if the following properties are maintained as a function of t:

No new local extrema can be created within the solution spatial domain,
The value of a local minimum is non-decreasing, and the value of a local maximum is non-increasing.

Template:Harvnb proved the following properties for a numerical scheme,

A monotone scheme is TVD, and

A TVD scheme is monotonicity preserving.

Monotone schemes are attractive for solving engineering and scientific problems because they do not provide non-physical solutions.

Godunov's theorem proves that only first order linear schemes preserve monotonicity and are therefore TVD. Higher order linear schemes, although more accurate for smooth solutions, are not TVD and tend to introduce spurious oscillations (wiggles) where discontinuities or shocks arise. To overcome these drawbacks, various high-resolution, non-linear techniques have been developed, often using flux/slope limiters.

References

↑ Harten, Ami (1983), "High resolution schemes for hyperbolic conservation laws", J. Comput. Phys. 49: 357–393, Template:Doi

Flux limiters

Flux limiters are used in high resolution schemes — numerical schemes used to solve problems in science and engineering, particularly fluid dynamics, described by partial differential equations (PDE's). They are used in high resolution schemes, such as the MUSCL scheme, to avoid the spurious oscillations (wiggles) that would otherwise occur with high order spatial discretization schemes due to shocks, discontinuities or sharp changes in the solution domain. Use of flux limiters, together with an appropriate high resolution scheme, make the solutions total variation diminishing (TVD).

Note that flux limiters are also referred to as slope limiters because they both have the same mathematical form, and both have the effect of limiting the solution gradient near shocks or discontinuities. In general, the term flux limiter is used when the limiter acts on system fluxes, and slope limiter is used when the limiter acts on system states.

How they work

The main idea behind the construction of flux limiter schemes is to limit the spatial derivatives to realistic values - for scientific and engineering problems this usually means physically realisable values. They are used in high resolution schemes for solving problems described by PDEs and only come into operation when sharp wave fronts are present. For smoothly changing waves, the flux limiters do not operate and the spatial derivatives can be represented by higher order approximations without introducing non-real oscillations. Consider the 1D semi-discrete scheme below,

$\frac{d u_i}{d t} + \frac{1}{\Delta x_i} \left[ F \left( u_{i + \frac{1}{2}} \right) - F \left( u_{i - \frac{1}{2}} \right) \right] =0,$

where, $F \left( u_{i + \frac{1}{2}} \right) \$ and $F \left( u_{i - \frac{1}{2}} \right) \$ represent edge fluxes for the ith cell. If these edge fluxes can be represented by low and high resolution schemes, then a flux limiter can switch between these schemes depending upon the gradients close to the particular cell, as follows,

$F \left( u_{i + \frac{1}{2}} \right) = f^{low}_{i + \frac{1}{2}} - \phi\left( r_i \right) \left( f^{low}_{i + \frac{1}{2}} - f^{high}_{i + \frac{1}{2}} \right)$ ,

$F \left( u_{i - \frac{1}{2}} \right) = f^{low}_{i - \frac{1}{2}} - \phi\left( r_{i-1} \right) \left( f^{low}_{i - \frac{1}{2}} - f^{high}_{i - \frac{1}{2}} \right)$ ,

where

$f^{low} = \$ low resolution flux,

$f^{high} = \$ high resolution flux,

$\phi\ (r) = \$ flux limiter function,

and $r\$ represents the ratio of successive gradients on the solution mesh, i.e.,

$r_{i} = \frac{u_{i} - u_{i-1}}{u_{i+1} - u_{i}}$ .

The limiter function is constrained to be greater than or equal to zero, i.e., $r \ge 0$ . Therefore, when the limiter is equal to zero (sharp gradient, opposite slopes or zero gradient), the flux is represented by a low resolution scheme. Similarly, when the limiter is equal to 1 (smooth solution), it is represented by a high resolution scheme. The various limiters have differing switching characteristics and are selected according to the particular problem and solution scheme. No particular limiter has been found to work well for all problems, and a particular choice is usually made on a trial and error basis.

Limiter functions

The following are common forms of flux/slope limiter function, $\phi\ (r)$ :

CHARM [not 2nd order TVD] (Zhou, 1995)

$\phi_{cm}(r)=\left\{ \begin{array}{ll} \frac{r\left(3r+1\right)}{\left(r+1\right)^{2}}, \quad r>0, \quad\lim_{r\rightarrow\infty}\phi_{cm}(r)=3 \\ 0 \quad \quad\, , \quad r\le 0 \end{array}\right.$

HCUS [not 2nd order TVD] (Waterson & Deconinck, 1995)

$\phi_{hc}(r) = \frac{ 1.5 \left(r+\left| r \right| \right)}{ \left(r+2 \right)} ; \quad \lim_{r \rightarrow \infty}\phi_{hc}(r) = 3$ .

HQUICK [not 2nd order TVD] (Waterson & Deconinck, 1995)

$\phi_{hq}(r) = \frac{2 \left(r + \left|r \right| \right)}{ \left(r+3 \right)} ; \quad \lim_{r \rightarrow \infty}\phi_{hq}(r) = 4$ .

Koren (Koren, 1993)

$\phi_{kn}(r) = \max \left[ 0, \min \left(2 r, \left(2 + r \right)/3, 2 \right) \right]; \quad \lim_{r \rightarrow \infty}\phi_{kn}(r) = 2$ .

minmod - symmetric (Roe, 1986)

$\phi_{mm} (r) = \max \left[ 0 , \min \left( 1 , r \right) \right] ; \quad \lim_{r \rightarrow \infty}\phi_{mm}(r) = 1$ .

monotonized central (MC) - symmetric (van Leer, 1977)

$\phi_{mc} (r) = \max \left[ 0 , \min \left( 2 r, 0.5 (1+r), 2 \right) \right] ; \quad \lim_{r \rightarrow \infty}\phi_{mc}(r) = 2$ .

Osher (Chatkravathy and Osher, 1983)

$\phi_{os} (r) = \max \left[ 0 , \min \left( r, \beta \right) \right], \quad \left(1 \leq \beta \leq 2 \right) ; \quad \lim_{r \rightarrow \infty}\phi_{os} (r) = \beta$ .

ospre - symmetric (Waterson & Deconinck, 1995)

$\phi_{op} (r) = \frac{1.5 \left(r^2 + r \right) }{\left(r^2 + r +1 \right)} ; \quad \lim_{r \rightarrow \infty}\phi_{op} (r) = 1.5$ .

smart [not 2nd order TVD] (Gaskell & Lau, 1988)

$\phi_{sm}(r) = \max \left[ 0, \min \left(2 r, \left(0.25 + 0.75 r \right), 4 \right) \right] ; \quad \lim_{r \rightarrow \infty}\phi_{sm}(r) = 4$ .

superbee - symmetric (Roe, 1986)

$\phi_{sb} (r) = \max \left[ 0, \min \left( 2 r , 1 \right), \min \left( r, 2 \right) \right] ; \quad \lim_{r \rightarrow \infty}\phi_{sb} (r) = 2$ .

Sweby - symmetric (Sweby, 1984)

$\phi_{sw} (r) = \max \left[ 0 , \min \left( \beta r, 1 \right), \min \left( r, \beta \right) \right], \quad \left(1 \leq \beta \leq 2 \right) ; \quad \lim_{r \rightarrow \infty}\phi_{sw} (r) = \beta$ .

UMIST (Lien & Leschziner, 1994)

$\phi_{um}(r) = \max \left[ 0, \min \left(2 r, \left(0.25 + 0.75 r \right), \left(0.75 + 0.25 r \right), 2 \right) \right] ; \quad \lim_{r \rightarrow \infty}\phi_{um}(r) = 2$ .

van Albada 1 - symmetric (van Albada, et al., 1982)

$\phi_{va1} (r) = \frac{r^2 + r}{r^2 + 1 } ; \quad \lim_{r \rightarrow \infty}\phi_{va1} (r) = 1$ .

van Albada 2 Alternative form [not 2nd order TVD] used on high spatial order schemes (Kermani, 2003)

$\phi_{va2} (r) = \frac{2 r}{r^2 + 1} ; \quad \lim_{r \rightarrow \infty}\phi_{va2} (r) = 0$ .

van Leer - symmetric (van Leer, 1974)

$\phi_{vl} (r) = \frac{r + \left| r \right| }{1 + \left| r \right| } ; \quad \lim_{r \rightarrow \infty}\phi_{vl} (r) = 2$ .

All the above limiters indicated as being symmetric, exhibit the following symmetry property,

$\frac{ \phi \left( r \right)}{r} = \phi \left( \frac{1}{r} \right)$ .

This is a desirable property as it ensures that the limiting actions for forward and backward gradients operate in the same way.

File:LimiterRegion.png

Admissible limiter region for second-order TVD schemes.

Unless indicated to the contrary, the above limiter functions are second order TVD. This means that they are designed such that they pass through a certain region of the solution, known as the TVD region, in order to guarantee stability of the scheme. Second-order, TVD limiters satisfy at least the following criteria:

$r \le \phi(r) \le 2r, \left( 0 \le r \le 1 \right) \$ ,
$1 \le \phi(r) \le r, \left( 1 \le r \le 2 \right) \$ ,
$1 \le \phi(r) \le 2, \left( r > 2 \right) \$ ,
$\phi(1) = 1 \$ ,

The admissible limiter region for second-order TVD schemes is shown in the Sweby Diagram opposite (Sweby, 1984), and plots showing limiter functions overlaid onto the TVD region are shown below. In this image, plots for the Osher and Sweby limiters have been generated using $\beta = 1.5$ .

Generalised minmod limiter

An additional limiter that has an interesting form is the van-Leer's one-parameter family of minmod limiters (van Leer, 1979; Harten and Osher, 1987; Kurganov and Tadmor, 2000). It is defined as follows

$\phi_{mg}(u,\theta)=\textrm{minmod}\left(\theta\frac{u_{i}-u_{i-1}}{\Delta x},\;\frac{u_{i+1}-u_{i-1}}{2\Delta x},\;\theta\frac{u_{i+1}-u_{i}}{\Delta x}\right),\quad\theta\in\left[1,2\right],$

where the multivariable minmod limiter is defined as

$\textrm{minmod}\left(z_{1},z_{2},\cdots\right):=\left\{ \begin{array}{cc} \min_{j}\quad & \textrm{if}\quad z_{j}>0\quad \forall j\\ \max_{j\quad} & \textrm{if}\quad z_{j}<0\quad \forall j\\ 0 & \textrm{otherwise}\end{array}\right. .$

Note: $\phi_{mg} \$ is most dissipative for $\theta=1, \$ when it reduces to $\phi_{mm}, \$ and is least dissipative for $\theta=2 \$ .

References

Chakravarthy, S.R. & S. Osher (1983), "High resolution applications of the Osher upwind scheme for the Euler equations", Proc. AIAA 6th Computational Fluid Dynamics Conference, AIAA Paper 83-1943, at 363–373
Gaskell, P.H. & A.K.C. Lau (1988), "Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm", Int. J. Num. Meth. Fluids 8 (6): 617–641, DOI:10.1002/fld.1650080602
Harten, A. & S. Osher (1987), "Uniformly high-order accurate nonoscillatory schemes. I", SIAM J. Numer. Anal. 24 (2): 279–309, DOI:10.1137/0724022
Hirsch, C. (1990), Numerical Computation of Internal and External Flows. Volume 2: Computational Methods for Inviscid and Viscous Flows, Wiley
Kermani, M.J.; A.G. Gerber & J.M. Stockie (2003), "Thermodynamically Based Moisture Prediction Using Roe’s Scheme", 4th Conference of Iranian AeroSpace Society, Amir Kabir University of Technology, Tehran, Iran, January 27–29
Koren, B. (1993), "A robust upwind discretisation method for advection, diffusion and source terms", in Vreugdenhil, C.B. & B. Koren, Numerical Methods for Advection-Diffusion Problems, Braunschweig: Vieweg, ISBN 3528076453, at 117
Kurganov, A. & E. Tadmor (2000), Solution of Two-Dimensional Riemann problems for Gas Dynamics without Riemann Problem Solvers, Report by Dept. of Mathematics, Univ. Michigan Available on-line at: CiteSeer.
Lien, F.S. & M.A. Leschziner (1994), "Upstream monotonic interpolation for scalar transport with application to complex turbulent flows", Int. J. Num. Meth. Fluids 19 (6): 527–548, DOI:10.1002/fld.1650190606
Leonard, B.P.; M.A. Leschziner & J. McGuirk (1978), "The QUICK algorithm: a uniformly 3rd-order finite-difference method for highly convective flows", Proc. 1st Conf. on Numerical Methods in Laminar & Turbulent Flow, Swansea, at 807
Roe, P.L. (1986), "Characteristic-based schemes for the Euler equations", Ann. Rev. Fluid Mech. 18: 337–365, DOI:10.1146/annurev.fl.18.010186.002005
Sweby, P.K. (1984), "High resolution schemes using flux-limiters for hyperbolic conservation laws", SIAM J. Num. Anal. 21 (5): 995–1011, DOI:10.1137/0721062
Van Albada, G.D.; B. Van Leer & W.W. Roberts (1982), "A comparative study of computational methods in cosmic gas dynamics", Astron. Astrophysics 108: 76–84
Van Leer, B. (1974), "Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme", J. Comp. Phys. 14 (4): 361–370, DOI:10.1016/0021-9991(74)90019-9
Van Leer, B. (1977), "Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow", J. Comp. Phys. 23 (3): 263–275, DOI:10.1016/0021-9991(77)90094-8
Van Leer, B. (1979), "Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method", J. Comp. Phys. 32: 101–136, DOI:10.1016/0021-9991(79)90145-1
Waterson, N.P. & H. Deconinck (1995), A unified approach to the design and application of bounded higher-order convection schemes
Zhou, G. (1995), Numerical simulations of physical discontinuities in single and multi-fluid flows for arbitrary Mach numbers, Goteborg, Sweden: Chalmers Univ. of Tech.

Sand box Approximation Schemes

From CFD-Wiki

Contents

S.K. Godunov's theorem

The theorem

See also

References

Further reading

Total variation diminishing

See also

References

Further reading

Flux limiters

How they work

Limiter functions

Generalised minmod limiter

See also

References

Further reading

Views

My wiki

wiki navigation

wiki search

wiki toolbox