Waves Speed Averaging Impact on Godunov type Schemes for Hyperbolic Equations with Discontinuous Coefficients: The linear scalar case

This paper deals with the waves speed averaging impact impact on Godunov type schemes for linear scalar hyperbolic equations with discontinuous coefficients. In many numerical schemes of Godunov type used in fluid dynamics, electromagnetic, electro-hydrodynamic problems and so on, usually a Riemann problem needs to be solved to estimate fluxes. The exact solution is generally not possible to obtain, but good approximations are provided in many situations like Roe and HLLC Riemann solvers in fluids. However all these solvers assume that the acoustic waves speed are continuous by considering some averaging. This could unfortunately lead to a wrong solution as we will show in this paper for the linear scalar case. Providing a Riemann solver in the general case of non-linear hyperbolic systems with discontinuous waves speed is a very hard task, therefore in this paper and as a first step, we focus on the linear and scalar case. In a previous work we proposed for such problems a Riemann solution that takes into account the discontinuities of the waves speed, we provided a numerical argument to show the validity of the solution. In this paper, first a new argument using regularization technique is provided to reinforce the validity of the proposed solution. Then, the corresponding Godunov scheme is derived and the effect of waves speed averaging is clearly demonstrated with a clear connection to the distribution product phenomenon.


Introduction
In many numerical methods such as finite volume, Discontinuous Galerkin (DG), Discontinuous finite volume, and so on, estimation of numerical fluxes at cells (subcells) faces is the most important part of the numerical scheme. The accuracy of the method depends on the accuracy of the flux estimation. For the convective fluxes, generally a Riemann problem is considered and then an approximation Riemann solver is used. This leads to a stable upwind numerical schemes. This approach was first proposed by Godunov [1], consequently we refer to such methods by Godunov type methods. Depending of the problem to solve, many Riemann solver approximations where developed. Among the most popular in computational fluid dynamic we can cite the Roe solver [2,3], the HLL Riemann solver [7] and the HLLC solver [6]. For the Roe solver, the Jacobian matrix is averaged in such a way that hyperbolicity, consistency with the exact Jacobian and conservation across discontinuities still fulfilled. This solver has been modified [4], [5] to overcome the shortcoming for low-density flows. The HLL solver [7] solves the original nonlinear flux to take nonlinearity into account. It has a major drawback however due to the space averaging process, the contact discontinuities, shear waves and material interfaces are not captured. To remedy this problem, the HLLC solver [6] was proposed by adding the missing wave to the structure. However all these methods assume that the waves speed are continuous across the left and right states of the Riemann problem (through the cell interfaces of the mesh) by applying diverse averaging process. This is not true in general; typical situations are recirculation for turbulent flows and transitions from subsonic to supersonic for transonic regimes. The impact of this averaging on the obtained numerical methods has never been addressed. To handle this important issue, we focus in this paper, as a first step, on the scalar and linear case. A Riemann solver of scalar hyperbolic linear equation with discontinuous coefficient is developed and published in a proceeding [8] providing a numerical argument of its validity, this was based on a first idea developed in [11]. This solver takes into account the discontinuities of waves speed. To reinforce the validity of the proposed solver and in absence of a rigorous mathematical proof, another strong argument based on a regularization strategy is proposed in this work. Then, and to assess the impact of the averaging process, a Godunov type scheme is derived using the proposed Riemann solver and compared to the case of averaged waves speed. The test will demonstrate as well that the critical situation is the one where we face a product of distributions, which is not defined in the classical theory of distributions, that could explain the non-validity of the waves speed averaging in such cases.

Riemann Solution for Hyperbolic Equation with Discontinuous Coefficient
Consider the following scalar linear hyperbolic equation with discontinuous coefficient, The initial condition ϕ 0 (x) and the coefficient a(x) are only bounded functions and can be discontinuous. From the theoretical point of view, the well-posedness of this type of problems is studied in [9]. It is shown that the more critical case is when the solution ϕ and the coefficient a(.) are discontinuous at the same location which leads to a product of distributions (for instance if a(.) is some Heaviside function and a Dirac function resulting from the derivative of ϕ). This product is not defined in the classical space of distributions which is not an algebra. The wellposdeness of the problem is then studied in a more appropriate space of generalized functions introduced by J.F Colombeau , known as well as the Colombeau algebra. For more details we refer to [13,14,9]. We recalled this critical case to make later the connection between this theoretical problem of product of distributions and its impact on the numerical aspect. Now, let's define the Riemann problem associated with problem (1) In this equation, the acoustic waves speed a(.) is discontinuous which is not taken into account in the existing Riemann solvers where acoustic waves speed are assumed to be continuous in the vicinity of the origin and some averaging is applied to ensure this assumption as we mentioned before. In the proceedings [8], a Riemann solver (a solution to problem (2)) is proposed based on the following observations of different situations. Case 1: a L > 0 and a R > 0 we have propagation of the discontinuity (of initial condition) to the right and we do not need to consider what happening within the fan defined by the two acoustic waves, because they will catch up if a L > β R and if a L < β R an expansion will appear. Case 2: a L < 0 and a R < 0 similar the previous case with a propagation of the discontinuity to the left. Case 3: a L < 0 and a R > 0 we have propagation of the discontinuity to the left and the right simultaneously, and we need to determine what happened within the fan defined by the two acoustics waves. We assume that a constant state appears and its expression will be given below. Case 4: a L > 0 and a R < 0 in this case we have opposite acoustic waves speed and then the discontinuity will remain blocked, which means there is no propagation. Based on the above observations, the Riemann solution of problem (2) is given by Where the expression of the constant λ is given by

On the proof
In this section we will first recall the numerical argument presented in [8] with more details and then we propose another strong argument that supports the validity of the proposed Riemann solver. This argument is based first on a regularization strategy that allows computing the exact solution using characteristics technique and finally obtain the exact solution of the original problem as an asymptotic limit of the regularized solution.

Numerical argument
To prove the validity of solution (3) and formula (4) at least numerically, the Riemann problem (2) is solved using a centered second order finite volume scheme stabilized by a first order artificial viscosity, which is equivalent in this case to a finite difference scheme. Let's set a n i = a(t n , x i ) ϕ n i = ϕ(t n , x i ) Then the explicit numerical scheme is give by Where ε is the diffusion coefficient that will be set to 0.1 for the whole simulation.

regularization approach argument
Here we propose to do a regularization of the initial condition ϕ 0 and the coefficient a(.), in such a way we can compute the exact regularized solution by the characteristics technique, and then the exact solution of the initial problem is obtained by the asymptotic limit.
Let's replace ϕ 0 and a(.) by: Where Where β ε = (a R − a L )/2ε 2 The regularization is achieved by connecting the two states by a quadratic functions. Now, consider the same problem with the regularized functions: The exact solution is then obtained using the characteristics technique. First we need to solve the equation Where y(.) ε is a curve passing by the point x at time t The solution of (7) is then given by Equation (8) is solved using the separable variable technique and by distinguish- ε, +∞[. The obtained solution is given by Where Notice that γ 1 = −γ 2 and γ 3 = −γ 4 , which implies that expressions of y ε (.) in the intervals [−ε, −ε + 1 σ 1 [ and [−ε + 1 σ 1 , < 0[ are identical and the same applies to intervals [0, ε − 1 σ 2 [ and [ε − 1 σ 2 , ε[ Therefore the expression of y ε (.) could be written in a more simple form as: The exact solution as we mentioned above, is then given by ϕ ε (x, t) = ϕ ε 0 (y(0 : x, t). To pass to the limit in ε to obtain the solution of the initial problem is not easy to perform analytically, this is done numerically instead. The regularized exact solution is implemented and the solution is plotted for different values of ε and compared to the proposed Riemann solver. Note that the implementation is not straightforward, indeed for a given point (x, t) the characteristics passing by this point could jump from and interval to another before crossing the x axes. Therefore we need to pay attention to apply the appropriate formula in (11). Figure (5) shows the regularized exact solution for different values of ε. We can see that it is converging to the proposed Riemann solver as ε goes to zero. This demonstrates clearly the validity of the proposed solution. In the following an answer to the relevant question concerning the impact of waves speed averaging on numerical schemes that use Riemann solvers is given.

Waves speed averaging effect
In this section we will demonstrate the effect of waves speed averaging, a process widely used in numericals schemes of hyperbolic equations, by comparing results of a Godunov scheme that uses the proposed Riemann solver with the one using the averaging process.

Godunov scheme
Let's derive the Godunov numerical schemes for the linear scalar hyperbolic equations with discontinuous coefficients using the proposed Riemann solver and the one for the case of waves speed averaged. Let h and ∆t be the space and time steps, then set x i−1/2 = (i − 1/2)h and t n = n∆t, and The Godunov scheme approximation is then given by where ϕ n+1,L i−1/2 and ϕ n+1,R i−1/2 are obtained by the classical projection process and are given in the following for both sachems.

Non-averaging waves speed case
The terms above are given by the following formulas using the propose Riemann solver and depending on the sign of a n i−1 and a n i Case a n i−1 < 0 and a n i > 0 with Case a n i−1 < 0 and a n i < 0 Case a n i−1 > 0 and a n i < 0 Case a n i−1 > 0 and a n i > 0

Averaging waves speed case
For the case where the waves speed are averaged, the Godunov scheme is obtained in the same way with formulas reduced to only two cases as a result of the averaging process, namely we have; Case of α n i− 1 2 = 0.5(a n i + a n i−1 ) > 0 Case of α n i− 1 2 = 0.5(a n i + a n i−1 ) < 0

Numerical tests
To demonstrate the effect of waves speed averaging we first consider the coefficient a(.) to be a two states function and set to Which means that we only considering the case of of the appearing the of intermediate constant state λ in the proposed Riemann solver. As for the initial condition, first we consider a sinusoidal function, in this case we avoid to have a product of distribution of type Heaviside times a Dirac function (not defined in the classical space of distributions) which was the main motivation to the use of the Colombeau algebra for the mathematical analyse of the problem. Results of the test are depicted in Figures (10,11). As we can see both schemes give the same results, this means that in this case averaging waves speed allows to well capture the intermediate state with a correct value. Now to analyse the case of the discontinuity of the coefficient a(.) and the initial condition at the same location that leads to a product of distribution, we consider the same function a(.) as before, and we add and subtracted a value to the sinusoidal function to create discontinuity at the origin as showed in Figure (12). The result of the propagation is depicted in Figures (13, 14), which show a sensitive difference between the two solutions. Figures (15, 16, 17) confirm the result for a discontinuous initial polynomial function. as we can see the scheme obtained by waves speed averaging is capable to predict the plateau (the intermediate state) but with a completely wrong value. In another hand this agreed with the theory of generalized function algebra that tries to handle the situation of product of distributions, while the averaging process in the numerical scheme is avoiding this situation.

Conclusions
The paper dealt with the impact of waves speed averaging process commonly used in Godunov type numerical schemes. First a Riemann solver for scalar hyperbolic equations with discontinuous coefficients, published first in a conference paper, is presented. After recalling the numerical argument to demonstrate the validity of the solver provided in the conference paper, a new argument based on a regularization procedure and asymptomatic analysis is introduced to reinforce the validity of the proposed solution in absence of a rigorous mathematical proof. The mathematical analysis of the PDE model is briefly discussed pointing out the problem of product of distributions that could results from the discontinuity of the coefficient and the solution at the same location that require the use of Colombeau algebra for analyses. To demonstrate the impact of waves speed averaging, a Gudunov scheme using the proposed Riemann solver is derived and results are compared to the same scheme using waves speed averaging. The tests show clearly a significant discrepancy of the solutions in the case of discontinuous coefficient with discontinuous solution at the same location. This shows first that the averaging process could lead to a wrong solution and second it makes a connection to the phenomenon of product of distributions. Note that this work is a first step toward a construction of a Riemann solver for non-linear hyperbolic equations and systems with discontinuous coefficients and its application for instance for inviscid fluxes estimation in the Navier-Stokes and electromagnetic equations discretization.