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 ULā , ULāā, URāā, URā
The value of URPā in the interface depends on the wave-speeds
The corresponding flux FRPā is just FRPā=F(URPā).
The Average State
Assuming the sonic waves speeds, SLā and SRā are known, we need SMā to estimate the intermediate average states Uā. 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ā
and the pressure
plāā=prāā=pā
where q=u the velocity normal to the discontinuity. The region between the sonic waves has constant pressure pā and normal velocity uā. To calculate the value of the contact wave speed, the Euler equations across the Riemann fan should be solved, resulting in
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ā)
The density of the inytermediate left state is:
Ļlāā=ĻlāSLāāSMāSLāāqlāā
and solutions to the intermediate-left state (in conserved variables)
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 l or L to r and R, 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ā)
The Sonic Wave Speed Estimates
To compute the intermediate states the sonic wave speeds (SLā,SRā) are needed. Following Batten et al. (1997), the wave speeds can be obtained from
SLā=min(qlāāclā,q~āāc~)
and
SRā=min(qrā+crā,q~ā+c~)
where q~ā=u~nxā+v~nyā+w~nzā and the average state is defined by
u~=(1+rĻā)(ulā+rĻāurā)ā
H~=(1+rĻā)(Hlā+rĻāHrā)ā
and
H~=CpāT~+21ā(u~2+v~2+w~2)
where rĻā is the ratio of densities
rĻā=Ļrā/Ļ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
Under construction
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:
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:
for simplicity fā²ā”dĻ/dx, and all derivatives evaluated at the cell face i+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/2, xi+2ā=3Īx/2 , etc. For example at xiā2ā=ā5Īx/2, we get
We look for a vector of weights αā=(α0ā,α1ā,...α5ā)such that that αāā ΦT=fdesiredā provides desired approximation for interpolation Ļi+1/2ā or gradients at the interface.
For example to obtain the first derivative, all terms should cancel except terms involving Īxfā² that should be 1, therefore the desired function should be fdesiredā=(0,1,0,0,0)
The system to solve wold be:
αāAT=fdesiredā
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:
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)
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ā
To achieve high order, an interpolation is built from node values. A fourth order interpolation, for example
The derivatives are obtained by conventional cell-centred central derivatives at iā1, i , i+1, etc. For example, using second order, the derivative would be a combination of derivatives of the form:
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ā 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ā ) :
S^Ļā=V1āā«SĻādV
otherwise the convergence is limited to second order.
Convergence of 4thcentral finite difference approach using MMS solving the 3D Navier-Stokes
at three different Reynolds numbers for a perfect gas with γ=1.4 and Pr=0.7.
The relevant python scripts are in tools/numerics
Convergence of 4thcentral finite volume approach using MMS solving the 3D Navier-Stokes
at three different Reynolds numbers for a perfect gas with γ=1.4 and Pr=0.7.
The relevant python scripts are in tools/numerics
The integration of the FV approach is done in symbolic form. This can be extremely costly in full Navier-Stokes with current symbolic python packages, Sympy, despite the simple analytic functions. To aid this, the density was kept constant to evaluate the MMS in the finite volume case
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:
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
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)
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
where ε is a small offset to ensure the denominator is never zero, while PTVDā=ā£Ļi+1āāĻiāā£+ā£ĻiāāĻiā1āā£and PJSTā=Ļi+1āā2Ļiā+Ļ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 :
Ļ=ĻĻā+ĻPāĻĻ2ā+ĻP2āā
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ā=kāāwkāFi+2(k)ā
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,...,r compute the polynomial reconstruction on each stencil:
F^(k) using third order or similar. For each stencil, compute a smoothness indicatorβkāā (same as in WENO-JS):
In practice, the indicators are precomputed using finite difference formulas. For example:
βkā=m=1ā2ācmā(ĪmF^(k))2
These expressions are designed to measure the oscillation or variation of F^(k)(x), and are minimized when the function is smooth over stencil k.
TENO introduces a normalized sensorĻ (like in WENO-Z):
Ļ=ā£Ī²0āāβrāā£
Then define the exponential cutoff function for each stencil:
Ļkā=exp(āĪ»(βkā/Ļ+ε)qā)
Turn off stencils with large smoothness indicators (discontinuities), and use optimal linear weights to combine the rest:
γkā={dkā,0,Ā āifĀ Ļkā>Ī“otherwiseā
Normalise the weights
Ļkā=ājāγjāγkāā
And finally use the same functional form as WENO
Fi+1/2ā=kāāwkāFi+2(k)ā
Under construction
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.
ātāEtāā+āxjāā(Ļe+Ļk+p)ujāā=0
where Etā=Ļe+Ļ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
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+2ĪtāRHS(Un)
Un+1=Un+ĪtRHS(Un+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.
[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
[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
[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
[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