Numerical Methods

Cerisse has implemented several numerical methods within the finite volume framework. While most of these methods are designed for high-speed flows, not all are limited to such cases. The method discussed here specifically addresses the Euler (convective) term of the equations. In compressible flows, high-frequency noise can accumulate when physical scales are not adequately resolved, especially near discontinuities. To mitigate this, many methods either directly filter or dissipate the noise, or use upwind-type stencils that indirectly introduce dissipation. It is important to note that no perfect numerical scheme exists; a compromise must be made between accuracy, speed, and stability. In the sections below, a brief overview of the rationale behind the numerical methods is provided. This description is not exhaustive, and the reader is encouraged to refer to the original papers for a more comprehensive explanation.

HLLC Riemann Solver

Cerisse employs the Harten-Lax-van Leer Contact (HLLC) solver, developed by Toro et al. (1994). The HLLC solver enhances approximate Riemann Solvers by incorporating the intermediate contact wave, a feature that is crucial for modelling reactive flows applications.

Scheme of the HLLC Riemann Solver with three waves propagating ar speed SL, SM and SR. Two acoustic waves and one contact. The waves separate four constant states ULU_L , ULāˆ—U_L^\ast, URāˆ—U_R^\ast, URU_R

The value of URPU_{RP} in the interface depends on the wave-speeds

URP={ULifĀ SL>0ULāˆ—ifĀ SL≤0<SMURāˆ—ifĀ SM≤0≤SRURifĀ SR<0U_{RP} = \begin{cases} U_L & \text{if } S_L > 0 \\ U_L^* & \text{if } S_L \leq 0 < S_M \\ U_R^* & \text{if } S_M \leq 0 \leq S_R \\ U_R & \text{if } S_R < 0 \end{cases}

The corresponding flux FRPF_{RP} is just FRP=F(URP)F_{RP}=F(U_{RP}).

The Average State

Assuming the sonic waves speeds, SLS_L and SRS_R are known, we need SMS_M to estimate the intermediate average states Uāˆ—U^\ast. The normal velocity and the pressure do not change across a contact discontinuity (mechanical equilbrium) and therefore the normal velocity is the contact wave speed

SM=ulāˆ—=urāˆ—=uāˆ—S_M=u_l^{\ast}=u_r^{\ast}=u^{\ast}

and the pressure

plāˆ—=prāˆ—=pāˆ—p_l^{\ast}=p_r^{\ast}=p^{\ast}

where q=uq = u the velocity normal to the discontinuity. The region between the sonic waves has constant pressure pāˆ—p^{\ast} and normal velocity uāˆ—u^{\ast}. To calculate the value of the contact wave speed, the Euler equations across the Riemann fan should be solved, resulting in

SM=ρrqr(SRāˆ’qr)āˆ’Ļlql(SLāˆ’ql)+plāˆ’prρr(SRāˆ’qr)āˆ’Ļl(SLāˆ’ql)S_M = \frac{\rho_r q_r(S_R - q_r)-\rho_l q_l(S_L - q_l) + p_l - p_r}{\rho_r (S_R - q_r)-\rho_l (S_L - q_l)}

From the contact wave speed, the Rankine-Hugoniot conditions are applied to each acoustic wave to find the average state. In the left wave, the jump relations are:

FLāˆ—āˆ’FL=SL(ULāˆ—āˆ’UL)F_L^\ast - F_L = S_L (U_L^\ast - U_L)

The density of the inytermediate left state is:

ρlāˆ—=ρlSLāˆ’qlSLāˆ’SM\rho_l^\ast = \rho_l \frac{S_L - q_l}{S_L - S_M }

and solutions to the intermediate-left state (in conserved variables)

(ρu)lāˆ—=(SLāˆ’ql)ρul+(pāˆ—āˆ’pl)SLāˆ’SM(\rho u)_l^{\ast} = \frac{(S_L-q_l)\rho u_l +(p^{\ast}- p_l)}{S_L-S_M}

and

elāˆ—=(SLāˆ’ql)ρElāˆ’plql+pāˆ—SMSLāˆ’SMe_l^{\ast} = \frac{(S_L-q_l)\rho E_l -p_lq_l+p^{\ast}S_M}{S_L-S_M}

The procedure to compute the intermediate-right state solution is analogous, but applying the Rankine-Hugoniot condition in the right sonic wave, that is interchanging the subscripts ll or LL to rr and RR, respectively. in previous equations. From the two intermediate states the flux may be obtained and replaced in the numerical scheme. The first-order scheme is

Fi+1/2=F(URP)F_{i+1/2}= F(U_{RP})

The Sonic Wave Speed Estimates

To compute the intermediate states the sonic wave speeds (SL,SRS_L,S_R) are needed. Following Batten et al. (1997), the wave speeds can be obtained from

SL=min⁔(qlāˆ’cl,q~āˆ’c~)S_L = \min(q_l - c_l, \tilde{q} - \tilde{c})

and

SR=min⁔(qr+cr,q~+c~)S_R = \min(q_r + c_r, \tilde{q} + \tilde{c})

where q~=u~nx+v~ny+w~nz\tilde{q}=\tilde{u} n_x+\tilde{v} n_y +\tilde{w} n_z and the average state is defined by

u~=(ul+rρur)(1+rρ)\tilde{u}=\frac{(u_l+r_{\rho} u_r)}{(1+r_{\rho})}
H~=(Hl+rρHr)(1+rρ)\tilde{H}=\frac{(H_l+r_{\rho} H_r)}{(1+r_{\rho})}

and

H~=CpT~+12(u~2+v~2+w~2)\tilde{H}=C_p \tilde{T} + \frac{1}{2} \left( \tilde{u}^2+\tilde{v}^2+\tilde{w}^2 \right)

where rρr_{\rho} is the ratio of densities

rρ=ρr/ρlr_{\rho}=\sqrt{\rho_r/\rho_l}

High order extension

The above scheme is formally first-order, a high order extension can be build by using a Total Variation Diminishing (TVD) reconstruction

Rusanov Scheme

The Rusanov flux is a simple upwind flux that requires a single wave-speed estimate. In the current implementation, it is very compact and can be used to perform quick tests. The numerical flux function is typically given by:

Fi+1/2=12(Fi+1+Fi)āˆ’12(∣λi+1∣+∣λi∣)(Ui+1āˆ’Ui)F_{i+1/2}= \frac{1}{2} \left( F_{i+1} + F_{i} \right) - \frac{1}{2} \left( \left| \lambda_{i+1} \right| + \left| \lambda_i \right| \right) (U_{i+1} - U_i)

where Ī»i\lambda_i ​ is the local characteristic speed at the i-th cell (often taken as the maximum eigenvalue of the Jacobian of the flux function).

Central differences

The simplest approach is to create a central interpolation based on values of the fluxes evaluated at cell points

Fi+1/2=L(..Fiāˆ’1,Fi,Fi+1,Fi+2..)=āˆ‘k=āˆ’HHFiāˆ’kαk+HF_{i+1/2}= \mathcal{L} \left ( .. F_{i-1}, F_{i}, F_{i+1}, F_{i+2} .. \right) = \sum_{k= -H}^{H} F_{i-k} \alpha_{k+H}

where H would be the number of cells in the stencil, which woud correspond to the order. For example, for a 4th order estimate, H would be 2 and the formulation would be:

Fi+1/2=α0Fiāˆ’1+α1Fi+α2Fi+1+α3Fi+2F_{i+1/2}= \alpha_0 F_{i-1} + \alpha_1 F_{i} + \alpha_2 F_{i+1} + \alpha_3 F_{i+2}

where Fi≔F(Ui)F_i \equiv F(U_i). The interpolation coefficents can be obtained from Taylor series around the cell face i+1/2i+1/2, see next sections.

Coefficents for interpolation and first derivatives (finite differences)

To obtain the coefficients, we evaluate the function using a series expansion around cell face i+1/2i+1/2, which will correspond to x=0x=0

Ļ•(x)ā‰ˆĻ•i+1/2+xf′+x22!f′′+x33!f′′′+x44!f′′′′+x55!f′′′′′+O(x6)\phi(x) \approx \phi_{i+1/2} + x f' + \frac{x^2}{2!} f'' + \frac{x^3}{3!} f''' + \frac{x^4}{4!} f'''' + \frac{x^5}{5!} f''''' + \mathcal{O}(x^6)

for simplicity f′≔dĻ•/dxf' \equiv {d \phi}/{dx}, and all derivatives evaluated at the cell face i+1/2i+1/2. For symmetric interpolation of order 2, we need 1 point at the left and another one at the right (for a total of 2) . For order 4 , 4 points are required, etc. Evaluating the expression at node points: xi+1=Ī”x/2x_{i+1}= \Delta x/2, xi+2=3Ī”x/2x_{i+2}= 3 \Delta x/2 , etc. For example at xiāˆ’2=āˆ’5Ī”x/2x_{i-2}= -5 \Delta x/2, we get

Ļ•iāˆ’2ā‰ˆĻ•i+1/2āˆ’5Ī”x2f′+(5/2Ī”x)22!fā€²ā€²āˆ’(5/2Ī”x)33!f′′′+(5/2Ī”x)44!fā€²ā€²ā€²ā€²āˆ’(5/2Ī”x)55!f′′′′′+O(Ī”x6)\phi_{i-2} \approx \phi_{i+1/2} - \frac{5 \Delta x}{2} f' + \frac{(5/2 \Delta x) ^2}{2!} f'' - \frac{(5/2 \Delta x)^3}{3!} f''' + \frac{(5/2 \Delta x)^4}{4!} f'''' - \frac{(5/2 \Delta x)^5}{5!} f''''' + \mathcal{O}(\Delta x^6)

Rearranging coefficents and remove truncation error for clarity, we obtain the six neighbouring points

Ļ•iāˆ’2=Ļ•i+1/2āˆ’52Ī”xf′+258(Ī”x)2fā€²ā€²āˆ’12548(Ī”x)3f′′′+625384(Ī”x)4fā€²ā€²ā€²ā€²āˆ’31253840(Ī”x)5f′′′′′\phi_{i-2} = \phi_{i+1/2} - \frac{5}{2} \Delta x f' + \frac{25}{8} (\Delta x)^2 f'' - \frac{125 }{48} (\Delta x)^3 f''' + \frac{625}{384} (\Delta x)^4f'''' - \frac{3125}{3840} (\Delta x)^5 f'''''
Ļ•iāˆ’1=Ļ•i+1/2āˆ’32Ī”xf′+98(Ī”x)2fā€²ā€²āˆ’2748(Ī”x)3f′′′+81384(Ī”x)4fā€²ā€²ā€²ā€²āˆ’2433840(Ī”x)5f′′′′′\phi_{i-1} = \phi_{i+1/2} - \frac{3}{2} \Delta x f' + \frac{9}{8} (\Delta x)^2 f'' - \frac{27}{48} (\Delta x)^3 f''' + \frac{81}{384} (\Delta x)^4f'''' - \frac{243}{3840} (\Delta x)^5 f'''''
Ļ•i=Ļ•i+1/2āˆ’12Ī”xf′+18(Ī”x)2fā€²ā€²āˆ’148(Ī”x)3f′′′+1384(Ī”x)4fā€²ā€²ā€²ā€²āˆ’13840(Ī”x)5f′′′′′\phi_{i} = \phi_{i+1/2} - \frac{1}{2} \Delta x f' + \frac{1}{8} (\Delta x)^2 f'' - \frac{1}{48} (\Delta x)^3 f''' + \frac{1}{384} (\Delta x)^4f'''' - \frac{1}{3840} (\Delta x)^5 f'''''
Ļ•i+1=Ļ•i+1/2+12Ī”xf′+18(Ī”x)2f′′+148(Ī”x)3f′′′+1384(Ī”x)4f′′′′+13840(Ī”x)5f′′′′′\phi_{i+1} = \phi_{i+1/2} + \frac{1}{2} \Delta x f' + \frac{1}{8} (\Delta x)^2 f'' + \frac{1}{48} (\Delta x)^3 f''' + \frac{1}{384} (\Delta x)^4f'''' + \frac{1}{3840} (\Delta x)^5 f'''''
Ļ•i+2=Ļ•i+1/2+32Ī”xf′+98(Ī”x)2f′′+2748(Ī”x)3f′′′+81384(Ī”x)4f′′′′+2433840(Ī”x)5f′′′′′\phi_{i+2} = \phi_{i+1/2} + \frac{3}{2} \Delta x f' + \frac{9}{8} (\Delta x)^2 f'' + \frac{27}{48} (\Delta x)^3 f''' + \frac{81}{384} (\Delta x)^4f'''' + \frac{243}{3840} (\Delta x)^5 f'''''
Ļ•i+3=Ļ•i+1/2+52Ī”xf′+258(Ī”x)2f′′+12548(Ī”x)3f′′′+625384(Ī”x)4f′′′′+31253840(Ī”x)5f′′′′′\phi_{i+3} = \phi_{i+1/2} + \frac{5}{2} \Delta x f' + \frac{25}{8} (\Delta x)^2 f'' + \frac{125}{48} (\Delta x)^3 f''' + \frac{625}{384} (\Delta x)^4f'''' + \frac{3125}{3840} (\Delta x)^5 f'''''

This can be arranged in a matrix system

Af=ΦA f = \Phi
A=(1āˆ’5/225/8āˆ’125/8625/384āˆ’3125/38401āˆ’3/29/8āˆ’27/4881/384āˆ’243/38401āˆ’1/21/8āˆ’1/481/384āˆ’1/384011/21/81/481/3841/384013/29/827/4881/384243/384015/225/8125/48625/3843125/3840)A = \left( \begin{array}{cccccc} 1 & -5/2 & 25/8 & -125/8 & 625/384 & -3125/3840 \\ 1 & -3/2 & 9/8 & -27/48 & 81/384 & -243/3840 \\ 1 & -1/2 & 1/8 & -1/48 & 1/384 & -1/3840 \\ 1 & 1/2 & 1/8 & 1/48 & 1/384 & 1/3840 \\ 1 & 3/2 & 9/8 & 27/48 & 81/384 & 243/3840 \\ 1 & 5/2 & 25/8 & 125/48 & 625/384 & 3125/3840 \end{array} \right)

with the array of unknowns

f=(Ļ•i+1/2Ī”xf′(Ī”x)2f′′(Ī”x)3f′′′(Ī”x)4f′′′′(Ī”x)5f′′′′′)f= \left( \begin{array}{c} \phi_{i+1/2} \\ \Delta x f' \\ (\Delta x)^2 f'' \\ (\Delta x)^3 f''' \\ (\Delta x)^4 f'''' \\ (\Delta x)^5 f''''' \\ \end{array} \right)

and the "known" information at nodes

Φ=(Ļ•iāˆ’2Ļ•iāˆ’1Ļ•iĻ•i+1Ļ•i+2Ļ•i+3)\Phi = \left( \begin{array}{c} \phi_{i-2} \\ \phi_{i-1} \\ \phi_{i} \\ \phi_{i+1} \\ \phi_{i+2} \\ \phi_{i+3} \\ \end{array} \right)

We look for a vector of weights α‾=(α0,α1,...α5)\underline{\alpha} = (\alpha_0, \alpha_1, ... \alpha_5 ) such that that α‾⋅ΦT=fdesired\underline{\alpha} \cdot \Phi^T = f_{desired} provides desired approximation for interpolation Ļ•i+1/2\phi_{i+1/2} or gradients at the interface.

For example to obtain the first derivative, all terms should cancel except terms involving Ī”xf′\Delta x f' that should be 1, therefore the desired function should be fdesired=(0,1,0,0,0)f_{desired} =( 0, 1, 0, 0, 0 )

The system to solve wold be:

α‾AT=fdesired\underline{\alpha} A^T = f_{desired}

which can be solved directly (check script solve.py in tools/numerics) and the solution is the required coefficents expressed as fractions. For example, the fourth and sixth order derivatives evaluated at cell faces are:

dĻ•dx∣i+1/2=1Ī”x(āˆ’124Ļ•i+2+98Ļ•i+1āˆ’98Ļ•i+124Ļ•iāˆ’1)+O(Ī”x)4\left.{\frac{d\phi}{dx}} \right|_{i+1/2} = \frac{1}{\Delta x } \left( \frac{-1}{24}\phi_{i+2} + \frac{9}{8} \phi_{i+1} - \frac{9}{8} \phi_{i} + \frac{1}{24} \phi_{i-1}\right) + \mathcal{O}(\Delta x)^4
dĻ•dx∣i+1/2=1Ī”x(3640Ļ•i+3āˆ’25384Ļ•i+2+7564Ļ•i+1āˆ’7564Ļ•i.+25384Ļ•iāˆ’1āˆ’3640Ļ•iāˆ’2)+O(Ī”x)6\left.{\frac{d\phi}{dx}} \right|_{i+1/2} = \frac{1}{\Delta x } \left( \frac{3}{640} \phi_{i+3} - \frac{25}{384}\phi_{i+2} + \frac{75}{64} \phi_{i+1} - \frac{75}{64} \phi_{i} . + \frac{25}{384} \phi_{i-1} - \frac{3}{640} \phi_{i-2} \right) + \mathcal{O}(\Delta x)^6

To obtain symmetrical interpolations, set fdesired=(1,0,0,0,0,0)f_{desired} = (1,0,0,0,0,0) and the expressions are

Ļ•i+1/2=āˆ’116Ļ•i+2+916Ļ•i+1+916Ļ•iāˆ’116Ļ•iāˆ’1+O(Ī”x)4\phi_{i+1/2} = - \frac{1}{16}\phi_{i+2} + \frac{9}{16} \phi_{i+1} + \frac{9}{16} \phi_{i} - \frac{1}{16} \phi_{i-1} + \mathcal{O}(\Delta x)^4
Ļ•i+1/2=3256Ļ•i+3āˆ’25256Ļ•i+2+75128Ļ•i+1+75128Ļ•iāˆ’25256Ļ•iāˆ’1+3256Ļ•iāˆ’2+O(Ī”x)6\phi_{i+1/2} = \frac{3}{256} \phi_{i+3} - \frac{25}{256}\phi_{i+2} + \frac{75}{128} \phi_{i+1} + \frac{75}{128} \phi_{i} - \frac{25}{256} \phi_{i-1} + \frac{3}{256} \phi_{i-2} + \mathcal{O}(\Delta x)^6

Coefficents for interpolation and first derivatives (finite volume)

The above coefficients are based on a finite difference, where

ϕi+1=ϕ(xi+1)\phi_{i+1} = \phi(x_{i+1})

However, Cerisse uses the finite volume approach, where the solution obtained

Ļ•^i=1Ī”x∫xiāˆ’1/2xi+1/2Ļ•(x)dx\hat{\phi}_{i} = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} \phi(x) dx

In a finite difference, the interpolated value at i+1/2i+1/2 is (assuming fourth order)

Ļ•i+1/2=α0Ļ•iāˆ’1+α1Ļ•i+α2Ļ•i+1+α3Ļ•i+2\phi_{i+1/2}= \alpha_0 \phi_{i-1} + \alpha_1 \phi_{i} + \alpha_2 \phi_{i+1} + \alpha_3 \phi_{i+2}

while in a finite volume, the interpolated value would be

Ļ•i+1/2=α^0Ļ•^iāˆ’1+α^1Ļ•^i+α^2Ļ•^i+1+α^3Ļ•^i+2\phi_{i+1/2}= \hat{\alpha}_0 \hat{\phi}_{i-1} + \hat{\alpha}_1 \hat{\phi}_{i} + \hat{\alpha}_2 \hat{\phi}_{i+1} + \hat{\alpha}_3 \hat{\phi}_{i+2}

In common second order approaches, both coefficients αk=α^k\alpha_k = \hat{\alpha}_k are the same α=α^\alpha=\hat{\alpha}. However, this changes for high order approximations.

To build the matrix of coefficients, we re-use the Taylor expansion of the previous section, expanding from the face

Ļ•(x)ā‰ˆĻ•i+1/2+xf′+x22!f′′+x33!f′′′+x44!f′′′′+x55!f′′′′′+O(x6)\phi(x) \approx \phi_{i+1/2} + x f' + \frac{x^2}{2!} f'' + \frac{x^3}{3!} f''' + \frac{x^4}{4!} f'''' + \frac{x^5}{5!} f''''' + \mathcal{O}(x^6)

Integrating

I=āˆ«Ļ•(x)dx=Ļ•i+1/2x+x22f′+x33!f′′+x44!f′′′+x55!f′′′′+x66!f′′′′′I = \int \phi(x) dx = \phi_{i+1/2} x + \frac{x^2}{2} f' + \frac{x^3}{3!} f'' + \frac{x^4}{4!} f''' + \frac{x^5}{5!} f'''' + \frac{x^6}{6!} f'''''

The finite volume representations (in 1D) are just differences of the above, for example to obtain the value at iāˆ’2i-2:

Ļ•^iāˆ’2=1Ī”x[I(xiāˆ’3/2)āˆ’I(xiāˆ’5/2)]=Ļ•i+1/2+12[(3/2)2āˆ’(5/2)2]Ī”xfā€²āˆ’13![(3/2)3āˆ’(5/2)3](Ī”x)2f′′+14![(3/2)4āˆ’(5/2)4](Ī”x)3fā€²ā€²ā€²āˆ’15![(3/2)5āˆ’(5/2)5](Ī”x)4f′′′′+16![(3/2)6āˆ’(5/2)6](Ī”x)5f′′′′′\hat{\phi}_{i-2} = \frac{1}{\Delta x } \left[ I(x_{i-3/2}) - I(x_{i-5/2}) \right] =\phi_{i+1/2} + \frac{1}{2}\left[ (3/2)^2 - (5/2)^2 \right] \Delta x f' - \frac{1}{3!} \left[ (3/2)^3 - (5/2)^3 \right] (\Delta x)^2 f'' \\ + \frac{1}{4!} \left[ (3/2)^4 - (5/2)^4 \right] (\Delta x)^3 f''' - \frac{1}{5!} \left[ (3/2)^5 - (5/2)^5 \right] (\Delta x)^4 f'''' + \frac{1}{6!} \left[ (3/2)^6 - (5/2)^6 \right](\Delta x)^5 f'''''

Rearranging the coefficents we obtain the six neighbouring finite volume points estimated from derivatives and values at the interface.

Ļ•^iāˆ’2=Ļ•i+1/2āˆ’52(Ī”x)f′+196(Ī”x)2fā€²ā€²āˆ’6524(Ī”x)3f′′′+211120(Ī”x)4fā€²ā€²ā€²ā€²āˆ’665720(Ī”x)5f′′′′′\hat{\phi}_{i-2} = \phi_{i+1/2} - \frac{5}{2} (\Delta x) f' + \frac{19}{6} (\Delta x)^2 f'' - \frac{65}{24} (\Delta x)^3 f''' + \frac{211}{120} (\Delta x)^4 f'''' - \frac{665}{720} (\Delta x)^5 f'''''
Ļ•^iāˆ’1=Ļ•i+1/2āˆ’32(Ī”x)f′+76(Ī”x)2fā€²ā€²āˆ’1524(Ī”x)3f′′′+31120(Ī”x)4fā€²ā€²ā€²ā€²āˆ’63720(Ī”x)5f′′′′′\hat{\phi}_{i-1} = \phi_{i+1/2} - \frac{3}{2} (\Delta x) f' + \frac{7}{6} (\Delta x)^2 f'' - \frac{15}{24} (\Delta x)^3 f''' + \frac{31}{120} (\Delta x)^4 f'''' - \frac{63}{720} (\Delta x)^5 f'''''
Ļ•^i=Ļ•i+1/2āˆ’12Ī”xf′+16(Ī”x)2fā€²ā€²āˆ’124(Ī”x)3f′′′+1120(Ī”x)4fā€²ā€²ā€²ā€²āˆ’1720(Ī”x)5f′′′′′\hat{\phi}_{i} = \phi_{i+1/2} - \frac{1}{2}\Delta x f' + \frac{1}{6} (\Delta x)^2 f'' - \frac{1}{24} (\Delta x)^3 f''' + \frac{1}{120} (\Delta x)^4 f'''' - \frac{1}{720} (\Delta x)^5 f'''''
Ļ•^i+1=Ļ•i+1/2+12Ī”xf′+16(Ī”x)2f′′+124(Ī”x)3f′′′+1120(Ī”x)4f′′′′+1720(Ī”x)5f′′′′′\hat{\phi}_{i+1} = \phi_{i+1/2} + \frac{1}{2}\Delta x f' + \frac{1}{6} (\Delta x)^2 f'' + \frac{1}{24} (\Delta x)^3 f''' + \frac{1}{120} (\Delta x)^4 f'''' + \frac{1}{720} (\Delta x)^5 f'''''
Ļ•^i+2=Ļ•i+1/2+32(Ī”x)f′+76(Ī”x)2f′′+1524(Ī”x)3f′′′+31120(Ī”x)4f′′′′+63720(Ī”x)5f′′′′′\hat{\phi}_{i+2} = \phi_{i+1/2} + \frac{3}{2} (\Delta x) f' + \frac{7}{6} (\Delta x)^2 f'' + \frac{15}{24} (\Delta x)^3 f''' + \frac{31}{120} (\Delta x)^4 f'''' + \frac{63}{720} (\Delta x)^5 f'''''
Ļ•^i+3=Ļ•i+1/2+52(Ī”x)f′+196(Ī”x)2f′′+6524(Ī”x)3f′′′+211120(Ī”x)4f′′′′+665720(Ī”x)5f′′′′′\hat{\phi}_{i+3} = \phi_{i+1/2} + \frac{5}{2} (\Delta x) f' + \frac{19}{6} (\Delta x)^2 f'' + \frac{65}{24} (\Delta x)^3 f''' + \frac{211}{120} (\Delta x)^4 f'''' + \frac{665}{720} (\Delta x)^5 f'''''

Following the same process as with finite difference we obtain the matrix coefficents (note that is very similar to the finite differenced one, but numbers in columns above 3 change)

A=(1āˆ’5/219/6āˆ’65/24211/120āˆ’665/7201āˆ’3/27/6āˆ’15/2431/120āˆ’63/7201āˆ’1/21/6āˆ’1/241/120āˆ’1/72011/21/61/241/1201/72013/27/615/2431/12063/72015/219/665/24211/120665/720)A = \left( \begin{array}{cccccc} 1 & -5/2 & 19/6 & -65/24 & 211/120 & -665/720 \\ 1 & -3/2 & 7/6 & -15/24 & 31/120 & -63/720 \\ 1 & -1/2 & 1/6 & -1/24 & 1/120 & -1/720 \\ 1 & 1/2 & 1/6 & 1/24 & 1/120 & 1/720 \\ 1 & 3/2 & 7/6 & 15/24 & 31/120 & 63/720 \\ 1 & 5/2 & 19/6 & 65/24 & 211/120 & 665/720 \end{array} \right)

and we obtain the coefficients solving the same system

α‾AT=fdesired\underline{\alpha} A^T = f_{desired}

The derivatives in the face are

dĻ•dx∣i+1/2=1Ī”x(āˆ’112Ļ•^i+2+54Ļ•^i+1āˆ’54Ļ•^i+112Ļ•^iāˆ’1)+O(Ī”x)4\left.{\frac{d\phi}{dx}} \right|_{i+1/2} = \frac{1}{\Delta x } \left( \frac{-1}{12}\hat{\phi}_{i+2} + \frac{5}{4} \hat{\phi}_{i+1} - \frac{5}{4} \hat{\phi}_{i} + \frac{1}{12} \hat{\phi}_{i-1}\right) + \mathcal{O}(\Delta x)^4
dĻ•dx∣i+1/2=1Ī”x(190Ļ•^i+3āˆ’536Ļ•^i+2+4936Ļ•^i+1āˆ’4936Ļ•^i+536Ļ•^iāˆ’1āˆ’190Ļ•^iāˆ’2)+O(Ī”x)6\left.{\frac{d\phi}{dx}} \right|_{i+1/2} = \frac{1}{\Delta x } \left( \frac{1}{90} \hat{\phi}_{i+3} - \frac{5}{36}\hat{\phi}_{i+2} + \frac{49}{36} \hat{\phi}_{i+1} - \frac{49}{36} \hat{\phi}_{i} + \frac{5}{36} \hat{\phi}_{i-1} - \frac{1}{90} \hat{\phi}_{i-2} \right) + \mathcal{O}(\Delta x)^6

Interpolation schemes are similarly found by setting fdesired=(1,0,0,0,0,0)f_{desired} = (1,0,0,0,0,0) and the result is:

Ļ•i+1/2=āˆ’112Ļ•^i+2+712Ļ•^i+1+712Ļ•^iāˆ’112Ļ•^iāˆ’1+O(Ī”x)4\phi_{i+1/2} = - \frac{1}{12}\hat{\phi}_{i+2} + \frac{7}{12} \hat{\phi}_{i+1} + \frac{7}{12}\hat{\phi}_{i} - \frac{1}{12} \hat{\phi}_{i-1} + \mathcal{O}(\Delta x)^4
Ļ•i+1/2=160Ļ•^i+3āˆ’215Ļ•^i+2+3760Ļ•^i+1+3760Ļ•^iāˆ’215Ļ•^iāˆ’1+160Ļ•^iāˆ’2+O(Ī”x)6\phi_{i+1/2} = \frac{1}{60} \hat{\phi}_{i+3} - \frac{2}{15} \hat{\phi}_{i+2} + \frac{37}{60} \hat{\phi}_{i+1} + \frac{37}{60} \hat{\phi}_{i} - \frac{2}{15} \hat{\phi}_{i-1} + \frac{1}{60} \hat{\phi}_{i-2} + \mathcal{O}(\Delta x)^6

These formulas are used to evaluate derivatives in the cell faces for viscous terms in the direction normal to the faces (for example in the x-direction).

Cross-derivatives

In the viscous terms, cross-derivative terms appear, i.e. the derivative in the y direction has to be compute in the x face

āˆ‚Ļ•āˆ‚y∣i+1/2\left.{\frac{\partial \phi }{\partial y}} \right|_{i+1/2}

To achieve high order, an interpolation is built from node values. A fourth order interpolation, for example

āˆ‚Ļ•āˆ‚y∣i+1/2=α0āˆ‚Ļ•āˆ‚y∣iāˆ’1+α1āˆ‚Ļ•āˆ‚y∣i+α2āˆ‚Ļ•āˆ‚y∣i+2+α3āˆ‚Ļ•āˆ‚y∣i+2\left.{\frac{\partial \phi}{\partial y}} \right|_{i+1/2} = \alpha_0 \left.{\frac{\partial \phi}{\partial y}} \right|_{i-1} + \alpha_1 \left.{\frac{\partial \phi}{\partial y}} \right|_{i} + \alpha_2 \left.{\frac{\partial \phi}{\partial y}} \right|_{i+2} +\alpha_3 \left.{\frac{\partial \phi}{\partial y}} \right|_{i+2}

The derivatives are obtained by conventional cell-centred central derivatives at iāˆ’1i-1, ii , i+1i+1, etc. For example, using second order, the derivative would be a combination of derivatives of the form:

āˆ‚Ļ•āˆ‚y∣iāˆ’1=Ļ•iāˆ’1,j+1āˆ’Ļ•iāˆ’1,jāˆ’1Ī”y\left.{\frac{\partial \phi}{\partial y}} \right|_{i-1} = \frac{ \phi_{i-1,j+1} - \phi_{i-1,j-1} }{\Delta y}

Method of Manufactured Solutions

Using the Method of Manufactured Solutions (MMS), it is possible to calculate the convergence of a numerical method implementation. By introducing a known analytical solution, Ļ•e\phi_e and calculating the corresponding source term, the theoretical convergence can be evaluated within both finite difference and finite volume frameworks. If the solution is sufficiently smooth, convergence is observed in the finite volume approach even when a quadrature method is not used for flux integration. This is consistent with Motheau and Wakefield (2021) observations. The order of convergence remains robust across a wide range of Reynolds numbers, from 0.01 to 100 (see Figures below). As a note, in the finite volume context the analytical solution and the source term, both have to be integrated (as the finite volume solution converges to Ļ•^→ϕ^e\hat{\phi} \rightarrow \hat{\phi}_e ) :

S^ρ=1V∫SρdV\hat{S}_\rho = \frac{1}{V} \int S_\rho dV

otherwise the convergence is limited to second order.

Convergence of 4th central finite difference approach using MMS solving the 3D Navier-Stokes at three different Reynolds numbers for a perfect gas with γ=1.4\gamma=1.4 and Pr=0.7. The relevant python scripts are in tools/numerics
Convergence of 4th central finite volume approach using MMS solving the 3D Navier-Stokes at three different Reynolds numbers for a perfect gas with γ=1.4\gamma=1.4 and Pr=0.7. The relevant python scripts are in tools/numerics

Skew-symmetric

Numerical errors associated with discretisation can be categorised into truncation and aliasing errors (Kravchenko and Moin, 1997; Lilly, 1965). The concept of numerical order alone is insufficient to fully characterize performance. Key properties such as dissipation, dispersion, and conservation are strongly influenced by the discretization scheme used for the convective term. To illustrate this, consider a one-dimensional scalar equation and three possible formulations for the nonlinear, hyperbolic term:

āˆ‚Uāˆ‚t+āˆ‚HUāˆ‚x=0\frac{\partial U}{\partial t} + \frac{\partial H U}{\partial x} = 0
Hdiv=āˆ‚UVāˆ‚xH_{div} =\frac{\partial UV }{\partial x}
Hconv=Uāˆ‚Vāˆ‚x+Vāˆ‚Uāˆ‚xH_{conv} = U \frac{\partial V }{\partial x} + V \frac{\partial U }{\partial x}
Hskew=12(Hconv+Hdiv)=12āˆ‚UVāˆ‚x+12(Uāˆ‚Vāˆ‚x+Vāˆ‚Uāˆ‚x)H_{skew} = \frac{1}{2} \left( H_{conv}+ H_{div} \right) = \frac{1}{2} \frac{\partial UV }{\partial x} +\frac{1}{2} \left( U \frac{\partial V }{\partial x} + V \frac{\partial U }{\partial x} \right)

Although the three forms above are equivalent at the continuous level, their discretizations differ significantly in terms of intrinsic properties and performance.

The skew-symmetric form, when used with centered schemes, has been demonstrated to conserve quadratic quantities of interest—such as kinetic energy in the incompressible limit . This conservation is attributed to the reduction of aliasing errors. Furthermore, a Fourier analysis of the three forms reveals that the skew-symmetric formulation possesses superior built-in de-aliasing characteristics (Blaisdell et al. 1996) .

The method implemented here is the approach of Ducros et al. (2000) , which capitalises on the built in de-aliasing property of the skew-symmetric operator of centred schemes, while ensuring local conservation by employing the flux-based formulation:

The flux can be derived in a convective an pressure term

F=UV+Fp=Fadv+FpF = UV +F_p = F^{adv} + F_p

A second order scheme can be then

Fi+1/2adv,skew=14(Ui+Ui+1)(Vi+Vi+1)F^{adv,skew}_{i+1/2} = \frac{1}{4} \left( U{i} + U_{i+1} \right) \left( V_{i} + V_{i+1} \right)

Despite the built-in de-aliasing properties of skew-schemes, they remain part of the family of centered schemes. As such, they can still face stability challenges. Additional strategies to mitigate oscillations and handle shock waves is needed.

Artificial dissipation

A way to stabilise the mechanism is through two terms following Jameson et al. (1981)

āˆ‚Fāˆ‚xā‰ˆāˆ‚Fāˆ‚x∣num+āˆ‚2α2Uāˆ‚x2+āˆ‚4α4Uāˆ‚x4\frac{\partial F}{\partial x} \approx \left .\frac{\partial F}{\partial x} \right |_{num} + \frac{\partial^2 \alpha_2 U}{\partial x^2} + \frac{\partial^4 \alpha_4 U}{\partial x^4}

The second-derivative term is used to capture discontinuities (hereafter shock term) and the fourth-derivative is to control high-frequency noise (damping term). The shock term acts near discontinuities and the damping in smooth parts of the flow. Using the same conservative form as before, the flux is modified by

Fi+1/2adv=Fi+1/2skew+α2āˆ‚Uāˆ‚x∣i+1/2+ α4āˆ‚3Uāˆ‚x3∣i+1/2F^{adv}_{i+1/2} = F^{skew}_{i+1/2} + \left . \alpha_2 \frac{\partial U}{\partial x} \right |_{i+1/2} +\left . \ \alpha_4 \frac{\partial^3 U}{\partial x^3} \right |_{i+1/2}

The terms can be rewritten using differences, the shock term is

α2āˆ‚Uāˆ‚x∣i+1/2ā‰ˆĪ±2Ī”Ui+1/2Ī”x\left . \alpha_2 \frac{\partial U}{\partial x} \right |_{i+1/2} \approx \alpha_2 \frac{\Delta U_{i+1/2} }{\Delta x}

where Ī”Ui+1/2=Ui+1āˆ’Ui\Delta U_{i+1/2} = U_{i+1}-U_i. The flux modification is then (grouping constants)

Fi+1/2shock=ϵi+1/2(2)(Ui+1āˆ’Ui)F^{shock}_{i+1/2} = \epsilon_{i+1/2}^{(2)} \left( U_{i+1} - U_i \right)

The damping term is similarly

α4āˆ‚3Uāˆ‚x3∣i+1/2ā‰ˆĪ±4Ī”xāˆ‚2Ī”Uāˆ‚x2∣i+1/2\left . \alpha_4 \frac{\partial^3 U}{\partial x^3} \right |_{i+1/2} \approx \frac{\alpha_4}{\Delta x} \left . \frac{\partial^2 \Delta U}{\partial x^2} \right |_{i+1/2}

Using central differences for āˆ‚2Ī”U/āˆ‚x2{\partial^2 \Delta U}/{\partial x^2} as

āˆ‚2Ī”Uāˆ‚x2=Ī”Ui+3/2āˆ’2Ī”Ui+1/2+Ī”Uiāˆ’1/2Ī”x2+O(Ī”x2)\frac{\partial^2 \Delta U}{\partial x^2} = \frac{\Delta U_{i+3/2} - 2 \Delta U_{i+1/2} + \Delta U_{i-1/2} }{\Delta x^2} + \mathcal{O}(\Delta x^2)

and replacing the difference, we get

Fi+1/2damp=ϵi+1/2(4)(Ui+2āˆ’3Ui+1+3Uiāˆ’Uiāˆ’1)F^{damp}_{i+1/2} = \epsilon_{i+1/2}^{(4)} \left( U_{i+2} - 3 U_{i+1} + 3 U_i - U_{i-1} \right)

For high-order schemes, the damping term accuracy can be increased by using a high-order central scheme:

āˆ‚2Ī”Uāˆ‚x2=āˆ’1/12Ī”Ui+5/2+4/3Ī”Ui+3/2āˆ’5/2Ī”Ui+1/2+4/3Ī”Uiāˆ’1/2āˆ’1/12Ī”Uiāˆ’3/2Ī”x2+O(Ī”x4)\frac{\partial^2 \Delta U}{\partial x^2} = \frac{ -1 /12 \Delta U_{i+5/2} + 4/3 \Delta U_{i+3/2} - 5/2 \Delta U_{i+1/2} + 4/3 \Delta U_{i-1/2} -1/12 \Delta U_{i-3/2} }{\Delta x^2} + \mathcal{O}(\Delta x^4)

the flux is similarly written as

Fi+1/2damp=ϵi+1/2(4)12(āˆ’Ui+3+17Ui+2āˆ’46Ui+1+46Uiāˆ’17Uiāˆ’1+Uiāˆ’2)F^{damp}_{i+1/2} = \frac{ \epsilon_{i+1/2}^{(4)} }{12} \left( -U_{i+3} + 17 U_{i+2} - 46 U_{i+1} + 46 U_i - 17 U_{i-1} + U_{i-2} \right)

The parameters ϵi+1/2(2)\epsilon_{i+1/2}^{(2)} and ϵi+1/2(4)\epsilon_{i+1/2}^{(4)} control the second and fourth order dissipation.

ϵi+1/2(2)=k(2)∣λi+1/2∣ψi+1/2\epsilon_{i+1/2}^{(2)}= k^{(2)} | \lambda_{i+1/2} | \psi_{i+1/2}

where ψ\psi is the shock/discontinuty detector (with value of 1 close to jumps) and Ī»\lambda is the eigenvalue. The sensor based on a variable Ļ•\phi is :

ψi=2āˆ£Ļ•i+1āˆ’2Ļ•i+Ļ•iāˆ’1∣PJST+PTVD+ε\psi_{i} = 2 \frac{|\phi_{i+1} - 2\phi_i + \phi_{i-1} |}{P_{JST} + P_{TVD} + \varepsilon}

where ε\varepsilon is a small offset to ensure the denominator is never zero, while PTVD=āˆ£Ļ•i+1āˆ’Ļ•i∣+āˆ£Ļ•iāˆ’Ļ•iāˆ’1∣P_{TVD} = |\phi_{i+1} -\phi_{i} | + |\phi_{i} -\phi_{i-1} |and PJST=Ļ•i+1āˆ’2Ļ•i+Ļ•iāˆ’1P_{JST} = \phi_{i+1} - 2\phi_i + \phi_{i-1}

The original formulation works with a sensor on pressure. Cerisse implements the improved approach of Bouheraoua (2014) by including an additional density sensor and coupling the two as :

ψ=ψρ2+ψP2ψρ+ψP\psi = \frac{\psi_\rho^2 + \psi_P^2}{\psi_\rho + \psi_P}

WENO and TENO

Weighted Essentially Non-Oscillatory (WENO) methods, introduced by Liu et al. (1994), employ a nonlinear adaptive procedure to automatically select the locally smoothest stencil. This approach aims to avoid using stencils that cross discontinuities when interpolating the interface flux. The WENO family encompasses various variations, which can be further classified. Despite these differences, all WENO methods share a common feature: the interface flux is expressed as a linear combination of fluxes derived from the stencils.

Fi+1/2=āˆ‘kwkFi+2(k)F_{i+1/2} = \sum_k w_k F^{(k)}_{i+2}

Reconstruction in characteristic variables improves performance, as the post-shock oscillations are reduced. WENO is known to be excessively dissipative in smooth parts of the flow

TENO

Like WENO, TENO uses multiple lower-order candidate stencils to build a high-order approximation. But instead of blending them with nonlinear weights, TENO uses a cutoff function to selectively activate only smooth stencils.

Given a set of candidate stencils, k=0,...,rk = 0, ... ,r compute the polynomial reconstruction on each stencil: F^(k)\hat{F}^{(k)} using third order or similar. For each stencil, compute a smoothness indicator βk\beta_k​ (same as in WENO-JS):

βk=āˆ‘l=1r∫xiāˆ’1xi+1Ī”x2lāˆ’1(dldxlF^(k)(x))2dx\beta_k = \sum_{l=1}^{r} \int_{x_{i-1}}^{x_{i+1}} \Delta x^{2l-1} \left( \frac{d^l}{dx^l} \hat{F}^{(k)}(x) \right)^2 dx

In practice, the indicators are precomputed using finite difference formulas. For example:

βk=āˆ‘m=12cm(Ī”mF^(k))2\beta_k = \sum_{m=1}^2 c_m \left( \Delta^m \hat{F}^{(k)} \right)^2

These expressions are designed to measure the oscillation or variation of F^(k)(x)\hat{F}^{(k)}(x), and are minimized when the function is smooth over stencil kk.

TENO introduces a normalized sensor Ļ„\tau (like in WENO-Z):

Ļ„=∣β0āˆ’Ī²r∣\tau = | \beta_0 - \beta_r |

Then define the exponential cutoff function for each stencil:

χk=exp⁔(āˆ’(βk/Ļ„+ε)qĪ»)\chi_k = \exp\left( - \frac{ \left( {\beta_k}/{\tau + \varepsilon} \right)^q }{\lambda} \right)

Turn off stencils with large smoothness indicators (discontinuities), and use optimal linear weights to combine the rest:

γk={dk,if χk>Ī“0,Ā otherwise\gamma_k = \begin{cases} d_k, & \text{if } \chi_k > \delta \\ 0, Ā  & \text{otherwise} \end{cases}

Normalise the weights

ωk=γkāˆ‘jγj\omega_k = \frac{\gamma_k}{\sum_j \gamma_j}

And finally use the same functional form as WENO

Fi+1/2=āˆ‘kwkFi+2(k)F_{i+1/2} = \sum_k w_k F^{(k)}_{i+2}

KEEP

The central KEEP (non-dissipative and physically-consistent kinetic energy and entropy preserving) schemes for compressible flows (Kuya and Kawai 2020) are based on splitting the energy equation.

āˆ‚Etāˆ‚t+āˆ‚(ρe+ρk+p)ujāˆ‚xj=0\frac{\partial E_t }{\partial t} + \frac{\partial (\rho e + \rho k + p) u_j}{\partial x_j} = 0

where Et=ρe+ρkE_t = \rho e + \rho k is the toral energy (internal plus kinetic energy). In the inviscid limit, from the momentum equation is possible to derive the kinetic energy equation

āˆ‚Ļkāˆ‚t+āˆ‚Ļujkāˆ‚xj+uāˆ‚pāˆ‚xj=0\frac{\partial \rho k }{\partial t} + \frac{\partial \rho u_j k }{\partial x_j} + u \frac{\partial p }{\partial x_j} = 0

which implies

āˆ‚Ļeāˆ‚t+āˆ‚Ļujeāˆ‚xj+pāˆ‚ujāˆ‚xj=0\frac{\partial \rho e }{\partial t} + \frac{\partial \rho u_j e }{\partial x_j} + p \frac{\partial u_j }{\partial x_j} = 0

using the fundamental equation of thermodynamics in differential form

de=Tdsāˆ’pρ2dρd e = T d s - \frac{p}{\rho^2} d \rho

it can be shown (see original paper) that if mass is conserved and internal energy follows the above expression, then entropy is conserved

āˆ‚Ļsāˆ‚t+āˆ‚Ļujsāˆ‚xj=0\frac{\partial \rho s }{\partial t} + \frac{\partial \rho u_j s }{\partial x_j} = 0

The KEEP scheme still need a shock capturing term.

Split terms

Given a function f=abf = a b, there are two forms to split the derivate

Divergence form of the derivative

āˆ‚fāˆ‚x=āˆ‚abāˆ‚x\frac{\partial f}{\partial x} = \frac{\partial a b}{\partial x}

Quadratic split

āˆ‚fāˆ‚x=Ā 12(aāˆ‚bāˆ‚x+bāˆ‚aāˆ‚x+āˆ‚abāˆ‚x)\frac{\partial f}{\partial x} =\ \frac{1}{2}\left( a \frac{\partial b }{\partial x} + b \frac{\partial a }{\partial x} + \frac{\partial a b}{\partial x} \right)

These formulations are analytically equivalent. With a function f=abcf = a b c, there are three forms to split the derivative

Divergence form

āˆ‚fāˆ‚x=āˆ‚abcāˆ‚x\frac{\partial f}{\partial x} = \frac{\partial a b c}{\partial x}

The quadratic form

āˆ‚fāˆ‚x=Ā 12(aāˆ‚bcāˆ‚x+bcāˆ‚aāˆ‚x+āˆ‚abcāˆ‚x)\frac{\partial f}{\partial x} =\ \frac{1}{2}\left( a \frac{\partial b c}{\partial x} + b c \frac{\partial a }{\partial x} + \frac{\partial abc }{\partial x} \right)

The cubic form

āˆ‚fāˆ‚x=Ā 14(aāˆ‚bcāˆ‚x+bāˆ‚acāˆ‚x+cāˆ‚abāˆ‚x+abāˆ‚cāˆ‚x+acāˆ‚bāˆ‚x+bcāˆ‚aāˆ‚x+Ā āˆ‚abcāˆ‚x)\frac{\partial f}{\partial x} =\ \frac{1}{4}\left( a \frac{\partial b c}{\partial x} + b \frac{\partial a c }{\partial x} + c \frac{\partial a b}{\partial x} + a b \frac{\partial c }{\partial x} + a c \frac{\partial b }{\partial x} + b c \frac{\partial a }{\partial x} +\ \frac{\partial abc }{\partial x} \right)

Time Marching

Expict Runge-Kutta RK2

A Runge-Kutta 2-stage, 2nd-order method is a time integration scheme that uses two evaluations (stages) of the right-hand side function (RHS) per timestep and achieves second-order accuracy in time.

Un+1/2=Un+Δt2RHS(Un)U^{n+1/2} = U^n + \frac{\Delta t}{2} RHS(U^n)
Un+1=Un+ΔtRHS(Un+1/2)U^{n+1} = U^n + \Delta t RHS(U^{n+1/2})

Strong Stability Preserving Runge-Kutta

The Strong Stability Preserving Runge-Kutta with n-stage and m-order, SSPRK(n,m), scheme is a time integration method designed to preserve the strong stability properties (e.g., total variation diminishing, monotonicity) of certain spatial discretizations when applied to hyperbolic PDEs

Some spatial discretizations (like TVD schemes) are non-oscillatory and stable under forward Euler time stepping with a small enough timestep. SSP Runge-Kutta schemes extend this stability to higher-order time integrators by writing the method as a convex combination of forward Euler steps.

U(0)=Un,U(i)=āˆ‘j=0iāˆ’1(αi,jU(j)+Ī”t βi,jRHS(U(j))),i=1,2,…,s,Un+1=U(s).\begin{aligned} U^{(0)} &= U^n, \\ U^{(i)} &= \sum_{j=0}^{i-1} \left( \alpha_{i,j} U^{(j)} + \Delta t \, \beta_{i,j} RHS(U^{(j)}) \right), \quad i = 1, 2, \dots, s, \\ U^{n+1} &= U^{(s)}. \end{aligned}

To ensure strong stability preservation, the method must satisfy:

  • αi,j≄0\alpha_{i,j} \geq 0, βi,j≄0\beta_{i,j} \geq 0

  • āˆ‘j=0iāˆ’1αi,j=1\sum_{j=0}^{i-1} \alpha_{i,j} = 1 (convex combination)

  • Each stage U(i)U^{(i)} is a convex combination of forward Euler steps

References

[1] Morinishi, Y. (1995). Conservative properties of finite difference schemes for incompressible flow. Center for Turbulence Research Annual Research Briefs

[2] Blaisdell, G., Spyropoulos, E., and Qin, J. (1996). The effect of the formulation of nonlinear terms on aliasing errors in spectral methods. Applied Numerical Mathematics, 21(3):207–219

[3] Ducros, F., Laporte, F., SoulĆØres, T., Guinot, V., Moinat, P., and Caruelle, B. (2000). High order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows. Journal of Computational Physics, 161(1):114–139

[4] Kravchenko, A. and Moin, P. (1997). On the effect of numerical errors in large eddy simulations of turbulent flows. Journal of Computational physics, 131(2):310–322

[5] Lilly, D. K. (1965). On the computational stability of numerical solutions of time-dependent non-linear geophysical fluid dynamics problems. Monthly Weather Review, 93(1):11–25

[6] Liu, X.-D., Osher, S., and Chan, T. (1994). Weighted essentially non-oscillatory schemes. Journal of Computational physics, 115(1):200–212

[7] Yuichi Kuya, Soshi Kawai, (2020) A stable and non-dissipative kinetic energy and entropy preserving (KEEP) scheme for non-conforming block boundaries on Cartesian grids, Computers & Fluids, Volume 200, 104427

[8] Fu, L., Hu, X. Y., and Adams, N. A. (2017). Targeted ENO schemes with tailored resolution property for hyperbolic conservation laws. Journal of Computational Physics, 349:97–121

[9] Toro, E. F., Spruce, M., and Speares, W. (1994). Restoration of the contact surface in the HLL-Riemann solver. Shock waves, 4(1):25–34

[10] Bouheraoua, L. (2014). Simulation aux grandes Ʃchelles et modƩlisation de la combustion supersonique. PhD thesis, Rouen, INSA.

[11] Shu, C-W. (1988). Total-Variation-Diminishing Time Discretizations. SIAM Journal of Scientific and Statistical Computing, 9(6):1073-1084

[12] Rusanov, V. V. E. (1962). The calculation of the interaction of non-stationary shock waves and obstacles. USSR Computational Mathematics and Mathematical Physics, 1(2), 304-320.

[13] Motheau E, Wakefield J. (2021). On the numerical accuracy in finite-volume methods to accurately capture turbulence in compressible flows. Int J Numer Meth Fluids, 93, 3020–3033.

[14] Jameson, A, Schmidt, W, Turkel, E (1981) Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes AIAA 1981-1259

[15] Batten, P, Clarke, N, Lambert, C and Causon D M (1997) On the Choice of Wavespeeds for the HLLC Riemann Solver. SIAM Journal on Scientific Computing 1997 18:6, 1553-1570

Last updated