Modified Semi-Lagrangian Godunov-Type Method Without Numerical Viscosity for Shocks

1. Introduction

Compressible flows simulation presents a significant challenge due to the existence of discontinuities in the fluid flow, including shocks, contact waves, and a wide range of continuous flow scales. To address these challenges, various high-resolution schemes have been developed to accurately resolve smooth flow scales while simultaneously capturing discontinuities. Notable examples include the artificial viscosity method [1,2], the total variation diminishing (TVD) method [3], the essentially non-oscillatory (ENO) scheme [4], and the weighted essentially non-oscillatory (WENO) scheme [5].

Harten et al. developed the ENO scheme [4,6,7,8], the first successful high-order method to allow the spatial discretization of hyperbolic conservation laws.This property is very useful for the numerical simulation of hyperbolic conservation laws because high-order numerical methods often cause false vibrations, especially near shocks or other interruptions. The spatial discretization of finite volumes of ENO has been studied [4], where it has been shown that it is uniformly accurate in the capturing of the location of any discontinuity in the flow. Subsequently, Shu and Osher [9,10] developed a finite difference ENO scheme. The main idea of the ENO scheme is to select the interpolation point stencil that suppresses the oscillations, i.e., to choose the best solution varying stencil, to approximate the cell boundary flow, and to avoid large false oscillations caused by interpolating data between discontinuities.

WENO schemes have been studied [5,11,12] to deal with potential numerical instability when selecting ENO stencils. The WENO method applies a convex combination of all ENO stencils, which obtains higher precision in smooth regions than the ENO method while maintaining the ENO property at discontinuities.

The mentioned above methods belong to the class of Euler-type methods. These methods [1,2,3,4,5,6,7,8,9,10,11,12,13] exhibit numerical viscosity in the presence of shocks. Additionally, there exist methods [14,15,16,17] that incorporate exact Riemann solutions while still being classified as Euler-type methods, which also exhibit numerical viscosity for shocks. Godunov et al. [14] proposed mitigating this issue by using moving grids; however, this approach significantly complicates the problem, particularly in multi-dimensional cases.

For highly rarefied gases, in particular, the direct simulation Monte Carlo (DSMC) method is applied, in the numerical scheme of which the fact of instantaneous change in flow parameters on shock waves is used [18].

In [19], we proposed a Lagrangian Godunov-type method that eliminates the numerical viscosity for shocks and provided a detailed description of the wall conditions. The method utilizes a fixed homogeneous grid. This paper presents a modified method that involves altering the formula for the density calculation in rarefaction waves. It was found that this formula allows obtaining good results for all types of test problems solved in [19]. While in work [19], two different formulas were used to obtain successful results: the linear law of the distribution of flow parameters in a rarefaction wave for the Shu and Osher shock tube problem [20] and the nonlinear law for the remaining problems of consideration. In addition, the following innovation is proposed here, which consists of considering an additional condition on the direction of the shock wave at the stage of calculating the acoustic substep. Moreover, this innovation, introduced into the scheme of the method, does not affect the results of the test problems already considered in [19], but manifests itself only in the numerical solution of the two last test cases of the interaction of two Riemann problems, presented in this work.

2. Materials and Methods

In the case of the one-dimensional scalar hyperbolic conservation law [13], the equations are given as follows:

$$\begin{gathered} {\displaystyle \frac{\partial \rho }{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \rho }{\partial \mathrm{x}}}+\rho {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \mathrm{p}}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \epsilon }{\partial \mathrm{x}}}+{\displaystyle \frac{\mathrm{p}}{\rho }}{\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}=0. \end{gathered}$$

(1)

Here, u is the convective velocity of flow, $\rho$ is the density, p is the pressure, $\epsilon$ is the internal specific energy, x is the direction, and t is the time.

The pressure is determined using the equation of state

$$\mathrm{p}=(\kappa -1)\rho \epsilon .$$

(2)

where $\kappa$ is the ratio of specific heat coefficients at constant pressure and volume. For air, $\kappa =1.4$.

Godunov et al. [14] proposed an iterative exact Riemann solver, which was implemented in a finite-difference scheme. Later, other researchers adopted a control volume method instead, while we employed the Lagrangian approach.

In this paper, as well as in a previous paper [19], the following method is proposed: The convection and acoustic stages are considered separately and are classified as distinct physical processes because the local acoustic velocity (i.e., sound speed) can vary significantly from the flow convection velocity. Therefore, we solve the following system during the acoustic stage:

$$\begin{gathered} {\displaystyle \frac{\partial \rho }{\partial \mathrm{t}}}+\rho {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{t}}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial \mathrm{p}}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{t}}}+{\displaystyle \frac{\mathrm{p}}{\rho }}{\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}=0. \end{gathered}$$

(3)

and the Equation (2). The system (3) is derived from (1) by discarding the convective terms and can be expressed in a matrix form as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \rho }{\partial \mathrm{t}}}\\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{t}}}\\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{t}}}\end{array}\right\}=\left[\mathrm{A}\right]\left\{\begin{array}{c}{\displaystyle \frac{\partial \rho }{\partial \mathrm{x}}}\\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}\\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{x}}}\end{array}\right\}$$

(4)

where $\left[\mathrm{A}\right]$ is the Jacobian matrix:

$$\left[\mathrm{A}\right]=\left[\begin{array}{ccc}0& \rho & 0\\ {\displaystyle \frac{(\kappa -1)\epsilon }{\rho }}& 0& \kappa -1\\ 0& (\kappa -1)\epsilon & 0\end{array}\right].$$

(5)

The matrix $\left[\mathrm{A}\right]$ can be expressed as

$$\left[\mathrm{A}\right]=\left[\mathrm{R}\right]\left[\mathrm{\Lambda }\right]\left[\mathrm{L}\right],$$

(6)

where $[\mathrm{\Lambda }]$ is the diagonal matrix of eigenvalues of the matrix $\left[\mathrm{A}\right]$.

$$\left[\mathrm{\Lambda }\right]=\left[\begin{array}{ccc}-\mathrm{c}& 0& 0\\ 0& 0& 0\\ 0& 0& \mathrm{c}\end{array}\right].$$

(7)

Here $\mathrm{c}$ is the sound velocity, which is determined as follows:

$$\mathrm{с}=\sqrt{\kappa {\displaystyle \frac{\mathrm{p}}{\rho }}},$$

(8)

A corresponding left-eigenvector matrix $[\mathrm{L}]$, determined as

$$\left[\mathrm{L}\right]=\left[\begin{array}{ccc}-{\displaystyle \frac{\mathrm{c}}{\kappa \rho }}& 1& -{\displaystyle \frac{\kappa -1}{\mathrm{c}}}\\ 1& 0& -{\displaystyle \frac{\rho \kappa }{{\mathrm{c}}^2}}\\ {\displaystyle \frac{\mathrm{c}}{\kappa \rho }}& 1& {\displaystyle \frac{\kappa -1}{\mathrm{c}}}\end{array}\right]$$

(9)

and a corresponding right eigenvectors matrix $[\mathrm{R}]$, defined as

$$\left[\mathrm{R}\right]=\left[\begin{array}{ccc}-{\displaystyle \frac{\rho }{2\mathrm{c}}}& {\displaystyle \frac{\kappa -1}{\kappa }}& {\displaystyle \frac{\rho }{2\mathrm{c}}}\\ {\displaystyle \frac{1}{2}}& 0& {\displaystyle \frac{1}{2}}\\ -{\displaystyle \frac{\mathrm{c}}{2\kappa }}& -{\displaystyle \frac{(\kappa -1)\epsilon }{\kappa \rho }}& {\displaystyle \frac{\mathrm{c}}{2\kappa }}\end{array}\right].$$

(10)

Since

$$\left[\mathrm{L}\right]\left[\mathrm{R}\right]=\left[\mathrm{E}\right],$$

(11)

where $\left[\mathrm{E}\right]$ is a diagonal identity matrix, multiplying Equation (4) from the left by $\left[\mathrm{L}\right]$ and applying (11) results in the following

$$\left[\mathrm{L}\right]\left\{\begin{array}{c}{\displaystyle \frac{\partial \rho }{\partial \mathrm{t}}}\\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{t}}}\\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{t}}}\end{array}\right\}=\left[\mathrm{\Lambda }\right]\left[\mathrm{L}\right]\left\{\begin{array}{c}{\displaystyle \frac{\partial \rho }{\partial \mathrm{x}}}\\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}\\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{x}}}\end{array}\right\}.$$

(12)

If $\left[\mathrm{L}\right]=[\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}]$, from (12) we can get the characteristic form

$$\begin{gathered} {\displaystyle \frac{\partial {\mathrm{w}}_1}{\partial \mathrm{t}}}-\mathrm{c}{\displaystyle \frac{\partial {\mathrm{w}}_1}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial {\mathrm{w}}_2}{\partial \mathrm{t}}}=0, \\ {\displaystyle \frac{\partial {\mathrm{w}}_3}{\partial \mathrm{t}}}+\mathrm{c}{\displaystyle \frac{\partial {\mathrm{w}}_3}{\partial \mathrm{x}}}=0, \end{gathered}$$

(13)

where

$$\begin{gathered} {\mathrm{w}}_1=-{\displaystyle \frac{\mathrm{c}*}{\kappa \rho *}}\rho +\mathrm{u}-{\displaystyle \frac{\kappa -1}{\mathrm{c}*}}\epsilon , \\ {\mathrm{w}}_2=\rho -{\displaystyle \frac{\kappa \rho *}{\mathrm{c}{*}^2}}\epsilon , \\ {\mathrm{w}}_3={\displaystyle \frac{\mathrm{c}*}{\kappa \rho *}}\rho +\mathrm{u}+{\displaystyle \frac{\kappa -1}{\mathrm{c}*}}\epsilon . \end{gathered}$$

(14)

Here, an asterisk represents constant values. Equation (13) describes the motion of two waves with velocities $-\mathrm{c}$ and $\mathrm{c}$. The waves transfer the values ${\mathrm{w}}_1$ and ${\mathrm{w}}_3$ as given in (14), respectively.

The proposed Lagrangian approach can be used to solve the system of Equation (13). However, as in the previous paper [19], we use the Godunov exact solver to solve the full nonlinear systems (1) and (2). We then subtract the convective velocity from the solution to obtain the solution of Equation (3). For more details, see below.

We accept the time step ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}$ at the acoustic stage in accordance with the Courant criterion

$${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={\mathrm{k}}_{\mathrm{c}}{\displaystyle \frac{\mathrm{h}}{{\mathrm{c}}_{\mathrm{s}}}},$$

(15)

where ${\mathrm{c}}_{\mathrm{s}}$ is the shock acoustic velocity and ${\mathrm{k}}_{\mathrm{c}}$ is integer. For example, if ${\mathrm{k}}_{\mathrm{c}}=3$, then the shock spreads through three grid cells during the acoustic stage.

During the convection stage, only the system with convective terms is solved:

$$\begin{gathered} {\displaystyle \frac{\partial \rho }{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \rho }{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \mathrm{u}}{\partial \mathrm{x}}}=0, \\ {\displaystyle \frac{\partial \epsilon }{\partial \mathrm{t}}}+\mathrm{u}{\displaystyle \frac{\partial \epsilon }{\partial \mathrm{x}}}=0. \end{gathered}$$

(16)

as well as Equation (13).

The time step during the convection stage is taken to be similar to that of (15):

$${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={\mathrm{k}}_{\mathrm{u}}{\displaystyle \frac{\mathrm{h}}{{\mathrm{u}}_{\mathrm{s}}}},$$

(17)

where ${\mathrm{u}}_{\mathrm{s}}$ is the shock convective velocity and ${\mathrm{k}}_{\mathrm{u}}$ is integer. For example, if ${\mathrm{k}}_{\mathrm{u}}=2$, then the shock spreads through two grid cells at the convection stage. If ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}\ne {}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}$, then the marching time steps ${}_{\mathrm{\Delta }}\mathrm{t}$ is chosen so that

$$\begin{gathered} {}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={\mathrm{n}}_{\mathrm{c}}{}_{\mathrm{\Delta }}\mathrm{t}, \\ {}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={\mathrm{n}}_{\mathrm{u}}{}_{\mathrm{\Delta }}\mathrm{t}, \end{gathered}$$

(18)

where ${\mathrm{n}}_{\mathrm{c}}$ and ${\mathrm{n}}_{\mathrm{u}}$ are integers.

2.1. A Methodological Scheme During the Acoustics Stage

During the acoustics stage, the system of Equation (3) is solved. Although the wave equation system (13) cannot be derived for a nonlinear problem, the fundamental properties of the linear system (13) are retained in the nonlinear system (3). The topologies of the nonlinear and linear systems are identical. This phenomenon corresponds to either a shock wave or a rarefaction wave [14]. Unfortunately, the variables and for the nonlinear system (3) are unknown. However, the exact solution of the Riemann problem for the nonlinear systems (1) and (2) can be determined using the Godunov method [14]. This solution is represented using the large variables: ${\mathrm{R}}_{\mathrm{L}}$ and ${\mathrm{R}}_{\mathrm{R}}$ for density, ${\mathrm{U}}_{\mathrm{L}}={\mathrm{U}}_{\mathrm{R}}$ for speed, ${\mathrm{E}}_{\mathrm{L}}$ and ${\mathrm{E}}_{\mathrm{R}}$ for internal specific energy, ${\mathrm{P}}_{\mathrm{L}}={\mathrm{P}}_{\mathrm{R}}$ for pressure, and ${\mathrm{D}}_{\mathrm{L}}$ and ${\mathrm{D}}_{\mathrm{R}}$ for the propagation frequency of shock or rarefaction waves. The subscripts $\mathrm{L}$ and $\mathrm{R}$ denote the left and right cells adjacent to the edge between them.

In the proposed method, ${\mathrm{w}}_1$ and ${\mathrm{w}}_3$ (14) are replaced by ${\overline{\mathrm{w}}}_1{=\text{{}\mathrm{P}}_{\mathrm{L}}{,\mathrm{U}}_{\mathrm{L}},{\mathrm{R}}_{\mathrm{L}}{,\mathrm{E}}_{\mathrm{L}}\text{}}$ and ${\overline{\mathrm{w}}}_3{=\text{{}\mathrm{P}}_{\mathrm{R}}{,\mathrm{U}}_{\mathrm{R}},{\mathrm{R}}_{\mathrm{R}}{,\mathrm{E}}_{\mathrm{R}}\text{}}$, which are each transported by the local acoustic velocities ${\mathrm{C}}_{\mathrm{L}}={\mathrm{D}}_{\mathrm{L}}-{\mathrm{U}}_{\mathrm{L}}$ and ${\mathrm{C}}_{\mathrm{R}}={\mathrm{D}}_{\mathrm{R}}-{\mathrm{U}}_{\mathrm{R}}$.

First, we must solve the Riemann problem using the Godunov method [14] for each pair of cells ($\mathrm{i}$ and $\mathrm{i}+1$) with the corresponding data ($\left\{{\mathrm{p}}_{\mathrm{i}},{\mathrm{u}}_{\mathrm{i}},{\rho }_{\mathrm{i}},{\epsilon }_{\mathrm{i}}\right\}$ and $\left\{{\mathrm{p}}_{\mathrm{i}+1},{\mathrm{u}}_{\mathrm{i}+1},{\rho }_{\mathrm{i}+1},{\epsilon }_{\mathrm{i}+1}\right\}$), respectively. Therefore, the values of the large variables $\text{{}\mathrm{P},\mathrm{U},{\mathrm{R}}_{\mathrm{L}},{\mathrm{E}}_{\mathrm{L}},{\mathrm{D}}_{\mathrm{L}},{\mathrm{R}}_{\mathrm{R}},{\mathrm{E}}_{\mathrm{R}},{\mathrm{D}}_{\mathrm{R}}\text{}}$ are obtained and assigned to the grid cells as follows:

$${{\mathrm{P}}_{\mathrm{L}}}_{\mathrm{i}+1}^{\mathrm{k}-1}=\mathrm{P},{{\mathrm{U}}_{\mathrm{L}}}_{\mathrm{i}+1}^{\mathrm{k}-1}=\mathrm{U},{{\mathrm{R}}_{\mathrm{L}}}_{\mathrm{i}+1}^{\mathrm{k}-1}={\mathrm{R}}_{\mathrm{L}},{{\mathrm{E}}_{\mathrm{L}}}_{\mathrm{i}+1}^{\mathrm{k}-1}={\mathrm{E}}_{\mathrm{L}},{{\mathrm{C}}_{\mathrm{L}}}_{\mathrm{i}+1}^{\mathrm{k}-1}={\mathrm{D}}_{\mathrm{L}}-\mathrm{U},$$

(19)

which causes the vector ${\overline{\mathrm{w}}}_1$ of the wave to spread to the left, along the negative direction of the x-axis; and:

which constitute the vector ${\overline{\mathrm{w}}}_1$ of the wave spreading along the negative direction of the x-axis (to the left); and:

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}=\mathrm{P},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}=\mathrm{U},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}={\mathrm{R}}_{\mathrm{R}},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}={\mathrm{E}}_{\mathrm{R}},{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}={\mathrm{D}}_{\mathrm{R}}-\mathrm{U},$$

(20)

which causes the vector ${\overline{\mathrm{w}}}_3$ of the wave to spread to the right, along the positive direction of the x-axis. In the case of a rarefaction wave, the variable D represents the wave velocity of the slower characteristic. For instance, for the right rarefaction wave, D_R is determined as

$${\mathrm{D}}_{\mathrm{R}}=\mathrm{U}+{\mathrm{c}}_{\mathrm{i}+1}-{\displaystyle \frac{(\kappa -1)}{2}}({\mathrm{u}}_{\mathrm{i}+1}-\mathrm{U}),$$

(21)

where ${\mathrm{c}}_{\mathrm{i}+1}$ is the sound velocity for $\mathrm{i}+1$ grid cell, which is determined according to the formula (8). For the left rarefaction wave ${\mathrm{D}}_{\mathrm{L}}$ is given as follows:

$${\mathrm{D}}_{\mathrm{L}}=\mathrm{U}-{\mathrm{c}}_{\mathrm{i}}-{\displaystyle \frac{(\kappa -1)}{2}}({\mathrm{u}}_{\mathrm{i}}-\mathrm{U}).$$

(22)

Examine a wave spreading in the positive x-direction (toward the right). The grid cells are examined in pairs. The coordinates of the cells after transfer during the acoustic stage for the right wave ${\overline{\mathrm{w}}}_3$ are determined as follows:

$$\begin{gathered} {\mathrm{x}}_1={\mathrm{x}}_{\mathrm{i}}+{{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}}_{\mathrm{\Delta }}{\mathrm{t}}_{\mathrm{c}}, \\ {\mathrm{x}}_2={\mathrm{x}}_{\mathrm{i}+1}+{{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}}_{\mathrm{\Delta }}{\mathrm{t}}_{\mathrm{c}}. \end{gathered}$$

(23)

The solution to the problem is expressed in the following form:

(1)

If the conditions

$$\left|{{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right|<{\mathrm{e}}_1,\left|{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right|<{\mathrm{e}}_1$$

(24)

are met, then for all grid cells satisfying condition ${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$, the solution at the next time step k is trivial and can be expressed as

$$\begin{gathered} {{\mathrm{P}}_{\mathrm{j}}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}+{{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}),{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}+{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}), \\ {{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}+{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}),{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}+{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}). \end{gathered}$$

(25)

(2)

If the condition

$$\left({{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)<-{\mathrm{e}}_1,({\mathrm{x}}_2-{\mathrm{x}}_1)>{\mathrm{e}}_2$$

(26)

are satisfied, a weak shock occurs, which does not overtake the solution from the preceding cell. Therefore, the solution is given as follows:

$$\begin{gathered} {{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{E}}_{\mathrm{R}\mathrm{i}}^{\mathrm{k}-1},\mathrm{if}{\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_1+{\mathrm{x}}_2)/2, \\ {{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},\mathrm{if}({\mathrm{x}}_1+{\mathrm{x}}_2)/2<{\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2; \end{gathered}$$

(27)

(3)

If

$$\left({{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{C}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)\ge {\mathrm{e}}_1$$

(28)

then

(3.1) if

$$\left({{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)<{\mathrm{e}}_1,$$

then we have the shock.

(3.1 a) If

$$\left({{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}+{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)\ge 0,$$

then

- if

$$({\mathrm{x}}_1-\mathrm{h}/2)\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_1+\mathrm{h}/2)$$

then

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{E}}_{\mathrm{R}\mathrm{i}}^{\mathrm{k}-1},$$

(29)

- if

$$({\mathrm{x}}_1+\mathrm{h}/2)<{\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$$

then

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}.$$

(30)

(3.1 b) If

$$\left({{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}+{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)<0,$$

then

- if

$$({\mathrm{x}}_2-\mathrm{h}/2)\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_2+\mathrm{h}/2)$$

then

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1};$$

(31)

- if

$${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}<({\mathrm{x}}_2-\mathrm{h}/2)$$

then

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1},{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{E}}_{\mathrm{R}\mathrm{i}}^{\mathrm{k}-1}.$$

(32)

The modification of the method proposed in this article, compared to work [19], consists of introducing the conditions “(3.1 a)” and “(3.1 b)”. Instead of checking these conditions, in [19] formulas (29), (30) were calculated, and formulas (31), (32) were not.

(3.2) If

$$\left({{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)\ge {\mathrm{e}}_1,$$

then a rarefaction wave occurs, and the solution for ${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$ is determined from the Riemann invariants of the hyperbolic conservation law as

$${{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}+{\displaystyle \frac{{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{{\mathrm{U}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1),$$

(33)

$${{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\left(\sqrt{{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}}+{\displaystyle \frac{\sqrt{{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}}-\sqrt{{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1)\right)}^2,$$

(34)

$${{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}={\left({\left({{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)}^{{\displaystyle \frac{\kappa -1}{2}}}+{\displaystyle \frac{{\left({{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}+1}^{\mathrm{k}-1}\right)}^{\frac{\kappa -1}{2}}-{\left({{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}-1}\right)}^{\frac{\kappa -1}{2}}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1)\right)}^{{\displaystyle \frac{2}{\kappa -1}}},$$

(35)

$${{\mathrm{P}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}=(\kappa -1){{\mathrm{R}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}{{\mathrm{E}}_{\mathrm{R}}}_{\mathrm{j}}^{\mathrm{k}}.$$

(36)

It will be shown below that formula (35) using instead of the formula from [19] allows us to obtain a universal scheme of the method, with good results being obtained for all test problems considered in this work. In this case, parasitic sawtooth oscillations do not arise in the case of using formula (35) when solving the Shu and Osher problem, as well as test examples about the interaction of two Riemann problems.

Under conditions (1)–(3), if the solution has already been assigned to the grid cell at the given substep, it is replaced with a new solution obtained using the exact Riemann solver. In this instance, the solution already assigned will be denoted with the subscript R, while the next solution applied to this cell will be denoted with the subscript L The resulting new solution comes from the result of the Riemann problem with subscript R. Furthermore, when assigning a solution to the grid cell, the condition is checked to ensure that the spreading of these solutions with a local acoustic velocity does not overwrite the solution of the forward shock, if such a shock wave exists.

For the wave spreading in the negative x-direction (toward the left), and with transfer values ${\overline{\mathrm{w}}}_1{=\text{{}\mathrm{P}}_{\mathrm{L}}{,\mathrm{U}}_{\mathrm{L}}{,\mathrm{R}}_{\mathrm{L}}{,\mathrm{E}}_{\mathrm{L}}\text{}}$, the conditions and expressions are derived in a similar manner.

After simulating wave spreading in both the left and right directions, each grid cell contains sets of transfer values, ${\overline{\mathrm{w}}}_1{{=\text{{}\mathrm{P}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{U}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{R}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{E}}_{\mathrm{L}}}^{\mathrm{k}}\text{}}$ and ${\overline{\mathrm{w}}}_3{{=\text{{}\mathrm{P}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{U}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{R}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{E}}_{\mathrm{R}}}^{\mathrm{k}}\text{}}$, respectively, along with the solution from the previous time substep ${\text{{}\mathrm{p}}^{\mathrm{k}-1}{,\mathrm{u}}^{\mathrm{k}-1},{\rho }^{\mathrm{k}-1},{\epsilon }^{\mathrm{k}-1}\text{}}$ (see Figure 1). The Riemann problem for these values is solved three times using the Godunov exact solver. The first instance of solving the problem uses ${\overline{\mathrm{w}}}_3{{=\text{{}\mathrm{P}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{U}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{R}}_{\mathrm{R}}}^{\mathrm{k}}{{,\mathrm{E}}_{\mathrm{R}}}^{\mathrm{k}}\text{}}$ as the left-state data and ${\text{{}\mathrm{p}}^{\mathrm{k}-1}{,\mathrm{u}}^{\mathrm{k}-1},{\rho }^{\mathrm{k}-1},{\epsilon }^{\mathrm{k}-1}\text{}}$ as the right-state data (see Figure 2), the resulting right solution is saved for the later calculations. The second instance of solving the problem uses ${\text{{}\mathrm{p}}^{\mathrm{k}-1}{,\mathrm{u}}^{\mathrm{k}-1},{\rho }^{\mathrm{k}-1},{\epsilon }^{\mathrm{k}-1}\text{}}$ as the left-state data and ${\overline{\mathrm{w}}}_1{{=\text{{}\mathrm{P}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{U}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{R}}_{\mathrm{L}}}^{\mathrm{k}}{{,\mathrm{E}}_{\mathrm{L}}}^{\mathrm{k}}\text{}}$ as the right-state data (see Figure 3). The resulting left solution is then assigned as the right-state data for the third problem (see Figure 4). The previously saved solution is assigned to the data of the third problem on the left (see Figure 4). Finally, in the third instance, the problem is solved using the previously mentioned data (see Figure 4), resulting in the values of the large variables ${\text{{}\mathrm{P}}^{\mathrm{k}*}{,\mathrm{U}}^{\mathrm{k}*}{{,\mathrm{R}}_{\mathrm{L}}}^{\mathrm{k}*}{{,\mathrm{E}}_{\mathrm{L}}}_{\mathrm{k}*}{,\mathrm{R}}_{\mathrm{R}}\mathrm{k}*{{,\mathrm{E}}_{\mathrm{R}}}^{\mathrm{k}*}\text{}}$ once again.

Since the pressure and convective velocity are identical on both the left and right sides, their values after the acoustic stage are given as

$${\mathrm{p}}_{\mathrm{с}\mathrm{i}}^{\mathrm{k}}={\mathrm{P}}_{\mathrm{i}}^{\mathrm{k}*},{\mathrm{u}}_{\mathrm{с}\mathrm{i}}^{\mathrm{k}}={\mathrm{U}}_{\mathrm{i}}^{\mathrm{k}*}.$$

(37)

The density and energy values are determined as follows:

$$\begin{gathered} {\rho }_{\text{с}\mathrm{i}}^{\mathrm{k}}{{=\mathrm{R}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}*},{\epsilon }_{\text{с}\mathrm{i}}^{\mathrm{k}}{{=\mathrm{E}}_{\mathrm{R}}}_{\mathrm{i}}^{\mathrm{k}*},\mathrm{if}{\mathrm{U}}_{\mathrm{i}}^{\mathrm{k}*}\ge 0, \\ {\rho }_{\text{с}\mathrm{i}}^{\mathrm{k}}{{=\mathrm{R}}_{\mathrm{L}}}_{\mathrm{i}}^{\mathrm{k}*},{\epsilon }_{\text{с}\mathrm{i}}^{\mathrm{k}}{{=\mathrm{E}}_{\mathrm{L}}}_{\mathrm{i}}^{\mathrm{k}*},\mathrm{if}{\mathrm{U}}_{\mathrm{i}}^{\mathrm{k}*}<0. \end{gathered}$$

(38)

The solution at the acoustics stage is ${\text{{}\mathrm{p}}_{\mathrm{c}}^{\mathrm{k}}{,\mathrm{u}}_{\mathrm{c}}^{\mathrm{k}},{\rho }_{\mathrm{c}}^{\mathrm{k}},{\epsilon }_{\mathrm{c}}^{\mathrm{k}}\text{}}$.

2.2. A Methodological Scheme During the Convection Stage

During the convection stage, the system of Equations (2) and (16) is solved. In the proposed method, the value ${\text{{}\mathrm{p}}_{\mathrm{c}}^{\mathrm{k}}{,\mathrm{u}}_{\mathrm{c}}^{\mathrm{k}},{\rho }_{\mathrm{c}}^{\mathrm{k}},{\epsilon }_{\mathrm{c}}^{\mathrm{k}}{\text{}}=\text{{}\mathrm{p}}_{\mathrm{u}}^{\mathrm{k}-1}{,\mathrm{u}}_{\mathrm{u}}^{\mathrm{k}-1},{\rho }_{\mathrm{u}}^{\mathrm{k}-1},{\epsilon }_{\mathrm{u}}^{\mathrm{k}-1}\text{}}$ is transported with the local convective velocity ${\mathrm{u}}_{\mathrm{c}}^{\mathrm{k}}$. In this approach, the results obtained from the acoustic stage serve as the initial input for the convection stage. Moreover, the grid cells were examined in pairs at the convection stage. The coordinates of the cells after transfer during the convection stage were determined as follows:

$$\begin{gathered} {\mathrm{x}}_1={\mathrm{x}}_{\mathrm{i}}+{{\mathrm{u}}_{\mathrm{u}\mathrm{i}}^{\mathrm{k}-1}}_{\mathrm{\Delta }}{\mathrm{t}}_{\mathrm{u}}, \\ {\mathrm{x}}_2={\mathrm{x}}_{\mathrm{i}+1}+{{\mathrm{u}}_{\mathrm{u}\mathrm{i}+1}^{\mathrm{k}-1}}_{\mathrm{\Delta }}{\mathrm{t}}_{\mathrm{u}}. \end{gathered}$$

(39)

The variables here refer to the convection stage and further we omit subscript u.

The solution during the convection stage is expressed in the following form:

(1)

If the conditions

$$\left|{\mathrm{u}}_{\mathrm{i}+1}-{\mathrm{u}}_{\mathrm{i}}\right|<{\mathrm{e}}_1,\left|{\rho }_{\mathrm{i}+1}-{\rho }_{\mathrm{i}}\right|<{\mathrm{e}}_1,\left|{\epsilon }_{\mathrm{i}+1}-{\epsilon }_{\mathrm{i}}\right|<{\mathrm{e}}_1$$

(40)

are met, then for all grid cells satisfying condition ${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$, the solution at the next time moment $\mathrm{k}$ is trivial and can be expressed as:

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1}+{\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1}),{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1}+{\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1}),{\rho }_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({\rho }_{\mathrm{i}}^{\mathrm{k}-1}+{\rho }_{\mathrm{i}+1}^{\mathrm{k}-1}),{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\displaystyle \frac{1}{2}}({\epsilon }_{\mathrm{i}}^{\mathrm{k}-1}+{\epsilon }_{\mathrm{i}+1}^{\mathrm{k}-1}).$$

(41)

When applying solution (41) to the grid cells, the separate condition is verified to ensure that its spreading with the convection velocity does not overwrite the solution of the forward shock, if such a shock wave exists. At the start of the convection stage, the positions of the shock boundaries at the next time step are determined.

(2)

If the conditions

$$({\mathrm{u}}_{\mathrm{i}+1}-{\mathrm{u}}_{\mathrm{i}})<{\mathrm{e}}_1$$

(42)

are satisfied, then we have the shock.

(2.1) If

$$({\mathrm{x}}_2-{\mathrm{x}}_1)>{\mathrm{e}}_2,$$

(43)

then the solution is defined as follows:

(2.1 a) if

$${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_1+{\mathrm{x}}_2)/2$$

then

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}}^{\mathrm{k}-1},$$

(44)

(2.1 b) if

$$({\mathrm{x}}_1+{\mathrm{x}}_2)/2<{\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$$

then

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}+1}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}+1}^{\mathrm{k}-1},$$

(45)

(2.2) if

$$({\mathrm{x}}_2-{\mathrm{x}}_1)\le {\mathrm{e}}_2,$$

(46)

then

(2.2 a) if

$$\left|{\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1}\right|<{\mathrm{e}}_1$$

then the discontinuity boundary is first determined as

$$\mathrm{x}={\mathrm{x}}_{\mathrm{i}}+\mathrm{h}/2+{({\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1}+{\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1})}_{\mathrm{\Delta }}{\mathrm{t}}_{\mathrm{u}},$$

(47)

and the solution is applied to cell $\mathrm{j}$ to the left of the discontinuity boundary

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}}^{\mathrm{k}-1},$$

(48)

and in cell $\mathrm{j}+1$ to the right of the discontinuity boundary, the following values are applied:

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}+1}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}+1}^{\mathrm{k}-1},$$

(49)

(2.2 b) if

$$\left|{\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1}\right|\ge {\mathrm{e}}_1$$

then the solution is expressed as follows:

- if

$${\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1}>{\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1}\mathrm{and}({\mathrm{x}}_1-\mathrm{h}/2)\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_1+\mathrm{h}/2)$$

then

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}}^{\mathrm{k}-1},$$

(50)

- if

$${\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}-1}<{\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1}\mathrm{and}({\mathrm{x}}_2-\mathrm{h}/2)\le {\mathrm{x}}_{\mathrm{j}}\le ({\mathrm{x}}_2+\mathrm{h}/2)$$

then

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{p}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1},{\rho }_{\mathrm{j}}^{\mathrm{k}}={\rho }_{\mathrm{i}+1}^{\mathrm{k}-1},{\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\epsilon }_{\mathrm{i}+1}^{\mathrm{k}-1}.$$

(51)

(3)

If

$$({\mathrm{u}}_{\mathrm{i}+1}-{\mathrm{u}}_{\mathrm{i}})\ge {\mathrm{e}}_1,$$

(52)

then the rarefaction wave occurs, and the solution for ${\mathrm{x}}_1\le {\mathrm{x}}_{\mathrm{j}}\le {\mathrm{x}}_2$ is determined from the Riemann invariants of the hyperbolic conservation law as

$${\mathrm{u}}_{\mathrm{j}}^{\mathrm{k}}={\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1}+{\displaystyle \frac{{\mathrm{u}}_{\mathrm{i}+1}^{\mathrm{k}-1}-{\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}-1}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1),$$

(53)

$${\epsilon }_{\mathrm{j}}^{\mathrm{k}}={\left(\sqrt{{\epsilon }_{\mathrm{i}}^{\mathrm{k}-1}}+{\displaystyle \frac{\sqrt{{\epsilon }_{\mathrm{i}+1}^{\mathrm{k}-1}}-\sqrt{{\epsilon }_{\mathrm{i}}^{\mathrm{k}-1}}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1)\right)}^2,$$

(54)

$${\rho }_{\mathrm{j}}^{\mathrm{k}}={\left({\left({\rho }_{\mathrm{i}}^{\mathrm{k}-1}\right)}^{{\displaystyle \frac{\kappa -1}{2}}}+{\displaystyle \frac{{\left({\rho }_{\mathrm{i}+1}^{\mathrm{k}-1}\right)}^{\frac{\kappa -1}{2}}-{\left({\rho }_{\mathrm{i}}^{\mathrm{k}-1}\right)}^{\frac{\kappa -1}{2}}}{{\mathrm{x}}_2-{\mathrm{x}}_1}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_1)\right)}^{{\displaystyle \frac{2}{\kappa -1}}},$$

(54)

$${\mathrm{p}}_{\mathrm{j}}^{\mathrm{k}}=(\kappa -1){\rho }_{\mathrm{j}}^{\mathrm{k}}{\epsilon }_{\mathrm{j}}^{\mathrm{k}}.$$

(56)

It will be shown below that formula (55) using instead of the formula from [19] allows us to obtain a universal scheme of the method, with good results being obtained for all test problems considered in this work. In this case, parasitic sawtooth oscillations do not arise in the case of using formula (55) when solving the Shu and Osher problem, as well as test examples about the interaction of two Riemann problems.

Under conditions (1)–(3), if the solution has already been applied to the grid cell at the given substep, it is replaced with a new one when the pressure of the new solution is greater than the pressure of the solution already in the grid cell.

2.3. Wall Boundary Conditions

The main concept behind the choice of the following boundary conditions is to study the waves reflection mechanism from the wall. At the acoustics stage, when conditions (1)–(3) are met, the wave ${\overline{\mathrm{w}}}_3{=\text{{}\mathrm{P}}_{\mathrm{R}}{,\mathrm{U}}_{\mathrm{R}},{\mathrm{R}}_{\mathrm{R}}{,\mathrm{E}}_{\mathrm{R}}\text{}}$ spreads to the right with velocity ${\mathrm{C}}_{\mathrm{R}}$. Upon reaching the wall, it is reflected and transforms into a wave ${\overline{\mathrm{w}}}_1{=\text{{}\mathrm{P}}_{\mathrm{R}},-{\mathrm{U}}_{\mathrm{R}},{\mathrm{R}}_{\mathrm{R}}{,\mathrm{E}}_{\mathrm{R}}\text{}}$, which then spreads to the left of the wall with velocity $-{\mathrm{C}}_{\mathrm{R}}$. In this scenario, when the wave ${\overline{\mathrm{w}}}_1$ is distributed across the grid cells, the dam-break problem is solved using the Godunov method. In this algorithm, the variables on the left correspond to the data already applied to the cell, while the variables on the right represent ${\overline{\mathrm{w}}}_1{=\text{{}\mathrm{P}}_{\mathrm{R}},-{\mathrm{U}}_{\mathrm{R}},{\mathrm{R}}_{\mathrm{R}}{,\mathrm{E}}_{\mathrm{R}}\text{}}$. As a result, new values of the large variables ${\text{{}\mathrm{P}}^{*}{,\mathrm{U}}^{*},{\mathrm{R}}_{\mathrm{L}}^{*},{\mathrm{E}}_{\mathrm{L}}^{*},{\mathrm{D}}_{\mathrm{L}}^{*},{\mathrm{R}}_{\mathrm{R}}^{*},{\mathrm{E}}_{\mathrm{R}}^{*},{\mathrm{D}}_{\mathrm{R}}^{*}\text{}}$ are obtained and applied to the grid cells as

$${\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{P}}_{\mathrm{i}}^{*},{\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{U}}_{\mathrm{i}}^{*},{\rho }_{\mathrm{i}}^{\mathrm{k}}={\mathrm{R}}_{\mathrm{L}\mathrm{i}}^{*},{\epsilon }_{\mathrm{i}}^{\mathrm{k}}={\mathrm{E}}_{\mathrm{L}\mathrm{i}}^{*},{\mathrm{c}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{D}}_{\mathrm{L}\mathrm{i}}^{*}-{\mathrm{U}}_{\mathrm{i}}^{*}.$$

(57)

For the wave ${\overline{\mathrm{w}}}_1{=\text{{}\mathrm{P}}_{\mathrm{L}}{,\mathrm{U}}_{\mathrm{L}},{\mathrm{R}}_{\mathrm{L}}{,\mathrm{E}}_{\mathrm{L}}\text{}}$, which transforms into the wave ${\overline{\mathrm{w}}}_3{=\text{{}\mathrm{P}}_{\mathrm{L}}{,-\mathrm{U}}_{\mathrm{L}},{\mathrm{R}}_{\mathrm{L}}{,\mathrm{E}}_{\mathrm{L}}\text{}}$ after encountering a wall

$${\mathrm{p}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{P}}_{\mathrm{i}}^{*},{\mathrm{u}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{U}}_{\mathrm{i}}^{*},{\rho }_{\mathrm{i}}^{\mathrm{k}}={\mathrm{R}}_{\mathrm{R}\mathrm{i}}^{*},{\epsilon }_{\mathrm{i}}^{\mathrm{k}}={\mathrm{E}}_{\mathrm{R}\mathrm{i}}^{*},{\mathrm{c}}_{\mathrm{i}}^{\mathrm{k}}={\mathrm{D}}_{\mathrm{R}\mathrm{i}}^{*}-{\mathrm{U}}_{\mathrm{i}}^{*}.$$

(58)

During the convection stage, a wave reverses direction upon encountering the wall. Under conditions (1)–(3) of this stage, the wave $\overline{\mathrm{w}}=\left\{{\mathrm{p}}^{\mathrm{k}},{\mathrm{u}}^{\mathrm{k}},{\rho }^{\mathrm{k}},{\epsilon }^{\mathrm{k}}\right\}$ spreads to the right with velocity ${\mathrm{u}}^{\mathrm{k}}$. Upon reaching the wall, it is reflected and transforms into a wave $\overline{\mathrm{w}}=\left\{{\mathrm{p}}^{\mathrm{k}},-{\mathrm{u}}^{\mathrm{k}},{\rho }^{\mathrm{k}},{\epsilon }^{\mathrm{k}}\right\}$ that moves to the left with velocity $-{\mathrm{u}}^{\mathrm{k}}$. In this scenario, when the wave $\overline{\mathrm{w}}$ is distributed across the grid cells, the dam-break problem is solved using the Godunov method. In this algorithm, the variables on the left correspond to the data already applied to the cell, while the variables on the right represent $\overline{\mathrm{w}}=\left\{{\mathrm{p}}^{\mathrm{k}},-{\mathrm{u}}^{\mathrm{k}},{\rho }^{\mathrm{k}},{\epsilon }^{\mathrm{k}}\right\}$. As a result, new values of the large variables ${\text{{}\mathrm{P}}^{*}{,\mathrm{U}}^{*},{\mathrm{R}}_{\mathrm{L}}^{*},{\mathrm{E}}_{\mathrm{L}}^{*},{\mathrm{D}}_{\mathrm{L}}^{*},{\mathrm{R}}_{\mathrm{R}}^{*},{\mathrm{E}}_{\mathrm{R}}^{*},{\mathrm{D}}_{\mathrm{R}}^{*}\text{}}$ are obtained and applied to the grid cells as shown in (57). If the wave $\overline{\mathrm{w}}=\left\{{\mathrm{p}}^{\mathrm{k}},{\mathrm{u}}^{\mathrm{k}},{\rho }^{\mathrm{k}},{\epsilon }^{\mathrm{k}}\right\}$ spreads to the left, a similar line of reasoning yields the result expressed in formula (58).

3. Results

The first example is the standard shock-tube problem of Sod [21], with the following initial conditions:

$$\begin{gathered} {\mathrm{p}}_{\mathrm{j}}^0=1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=1,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<0, \\ {\mathrm{p}}_{\mathrm{j}}^0=0.1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=0.125,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 0. \end{gathered}$$

(59)

The modeling region was defined as $\mathrm{x}\in \left[-4.5,5.5\right]$. The grid consisted of 100 cells, with a grid spacing of $\mathrm{h}=0.1$. The time step for the acoustics stage was determined based on the Courant criterion (15), with ${\mathrm{k}}_{\mathrm{c}}=1$, meaning the shock spreads through a single grid cell during this stage. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=2$, indicating that the shock spreads through two grid cells during this stage. The marching time step was set to ${}_{\mathrm{\Delta }}\mathrm{t}=0.02020929$. The acoustics stage was executed over six time steps, with ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}=6{}_{\mathrm{\Delta }}\mathrm{t}$, followed by the convection stage, which was carried out after ten time steps, with ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}=10{}_{\mathrm{\Delta }}\mathrm{t}$.

The results after 110 marching time steps are presented in Figure 5, Figure 6 and Figure 7, where they are compared with the exact solution and the data obtained using the previous author's method [19]. Figure 5, Figure 6 and Figure 7 show the agreement between the results for the previous and modified methods and demonstrate that the numerical viscosity is not present in the shock. However, due to the accumulation of round-off errors, the shock is delayed by one grid cell at the observed time.

The second example is the Riemann problem of Lax [22], with the following initial conditions:

$$\begin{gathered} {\mathrm{p}}_{\mathrm{j}}^0=3.52773,{\mathrm{u}}_{\mathrm{j}}^0=0.69888,{\rho }_{\mathrm{j}}^0=0.445,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<0, \\ {\mathrm{p}}_{\mathrm{j}}^0=0.571,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=0.5,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 0. \end{gathered}$$

(60)

The modeling region was defined as $\mathrm{x}\in \left[-8,6\right]$. The grid consisted of 140 cells, with a grid spacing of h = 0.1. The time step for the acoustics stage was determined based on the Courant criterion (15), with ${\mathrm{k}}_{\mathrm{c}}=2$, meaning the shock spreads through two grid cells during this stage. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=3$, indicating that the shock spreads through three grid cells during this stage. The marching time step was set to ${}_{\mathrm{\Delta }}\mathrm{t}=0.2$. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results are presented in Figure 8, Figure 9 and Figure 10, where they are compared with the exact solution and the data obtained using the previous author's method [19] at the time $\mathrm{t}=2.0$. Figure 8, Figure 9 and Figure 10 show the agreement between the results for the previous and modified methods and demonstrate that the numerical viscosity is not present in the shock. However, due to the accumulation of round-off errors, the shock is delayed by one grid cell at the observed time.

The third example is the Shu–Osher shock-tube problem.

The Shu–Osher problem [20] simulates the propagation of a normal shock front within a one-dimensional inviscid flow, incorporating artificial density fluctuations. The simulation is initialized with the following conditions:

$$\begin{gathered} {\mathrm{p}}_{\mathrm{j}}^0=10.3333,{\mathrm{u}}_{\mathrm{j}}^0=2.629369,{\rho }_{\mathrm{j}}^0=3.857143,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<1/8, \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{p}}_{\mathrm{j}}^0=1.0,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=1+0.2\sin (16\mathrm{\pi }\mathrm{x}),\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 1/8. \end{gathered}$$

(61)

The modeling region was defined as $\mathrm{x}\in \left[0,1\right]$. Three types of grids were utilized: a coarse grid with 192 points, a fine grid with 384 points, and a finer grid with 1,000 points. The time step for the acoustics stage was determined based on the Courant criterion (15), with${\mathrm{k}}_{\mathrm{c}}=2$, meaning the shock spreads through two grid cells during this stage—for all grids. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=5$, indicating that the shock spreads through five grid cells during this stage—for all grids. The marching time step was set to ${}_{\mathrm{\Delta }}\mathrm{t}=0.010270978683$ for the coarse grid, ${}_{\mathrm{\Delta }}\mathrm{t}=0.005135489341$ for the fine grid, and ${}_{\mathrm{\Delta }}\mathrm{t}=0.001972027907$ for the finer grid. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results are presented in Figure 11, Figure 12 and Figure 13, where they are compared with the "exact" solution [20] at t = 0.178. The results indicate that the fine and finer grids provide greater accuracy compared to the coarse grid. The proposed method has greater accuracy than the previous one [19], as shown in Figure 12. However, due to the accumulation of round-off errors, the shock is delayed by a few grid cells for the coarse grid results at the observed time.

The fourth example is the double expansion wave problem [23], with initial conditions specified as follows:

$$\begin{gathered} \hspace{1em}\hspace{1em}{\mathrm{p}}_{\mathrm{j}}^0={10}^5,{\mathrm{u}}_{\mathrm{j}}^0=-100,{\rho }_{\mathrm{j}}^{\mathrm{k}}=1.2,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<0, \\ {\mathrm{p}}_{\mathrm{j}}^0={10}^5,{\mathrm{u}}_{\mathrm{j}}^0=100,{\rho }_{\mathrm{j}}^0=1.2,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 0. \end{gathered}$$

(62)

The modeling region was defined as $\mathrm{x}\in \left[-0.5,0.5\right]$. Two types of grids were utilized: a coarse grid with 100 points, and a finer grid with 1,000 points; the grid step was $\mathrm{h}=0.01$ and $\mathrm{h}=0.001$, respectively. The time step for the acoustics stage was determined based on the Courant criterion (15), with${\mathrm{k}}_{\mathrm{c}}=7$, meaning the shock spreads through seven grid cells during this stage—for all grids. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=2$, indicating that the shock spreads through two grid cells during this stage—for all grids. The marching time step was set to ${}_{\mathrm{\Delta }}\mathrm{t}=0.0002$ for the coarse grid, and ${}_{\mathrm{\Delta }}\mathrm{t}=0.00002$ for the finer grid. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results after 4 marching time steps for the coarse grid and after 40 marching time steps for the finer grid (t=0.0008) are presented in Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 , where they are compared with the exact solution. Figure 14, Figure 15 and Figure 16 demonstrate a strong correlation between the exact and numerical solutions for the double expansion wave problem, while Figure 17, Figure 18 and Figure 19 exhibit an even higher degree of agreement.

The fifth example presents a one-dimensional problem defined by specific initial conditions

$${\mathrm{p}}_{\mathrm{j}}^0=101325.0,{\mathrm{u}}_{\mathrm{j}}^0=600.0,{\rho }_{\mathrm{j}}^0=1.25.$$

(63)

In this case, at time $\mathrm{t}=0$, a thin wall was lowered at the channel section $\mathrm{x}=0$.

The modeling region was defined as $\mathrm{x}\in \left[-5,5\right]$. The grid consisted of 100 cells, with a grid spacing of $\mathrm{h}=0.1$. The time step for the acoustics stage was determined based on the Courant criterion (15), with${\mathrm{k}}_{\mathrm{c}}=2$, meaning the shock spreads through two grid cells during this stage. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=5$, indicating that the shock spreads through five grid cells during this stage. The marching time step was set to ${}_{\mathrm{\Delta }}\mathrm{t}=0.00079$. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results, compared with the exact solution at $\mathrm{t}=0.00474$, are presented in Figure 20, Figure 21 and Figure 22. Figure 20, Figure 21 and Figure 22 show the same results for the previous [19] and modified methods. In this case, the shock does not induce any delay, nor is there any numerical viscosity associated with it.

Furthermore, we consider two (sixth and seventh) test examples of the interaction of two Riemann problems with different initial conditions. The initial conditions for the sixth example for the simulation are

$$\begin{gathered} {\mathrm{p}}_{\mathrm{j}}^0=10,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^{\mathrm{k}}=10,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<0, \\ {\mathrm{p}}_{\mathrm{j}}^0=1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^{\mathrm{k}}=1,\mathrm{if}0\le {\mathrm{x}}_{\mathrm{j}}<0.4, \\ {\mathrm{p}}_{\mathrm{j}}^0=0.1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=0.125,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 0.4. \end{gathered}$$

(64)

This is a more complicated problem proposed by Sod [21] (second and third lines of the formula), to which another discontinuity decay was added (using the first line of the formula (64)).

The modeling region was defined as $\mathrm{x}\in \left[-2;2\right]$. The grid consisted of 400 cells, with a grid spacing of $\mathrm{h}=0.01$. The time step for the acoustics stage was determined based on the Courant criterion (15), with ${\mathrm{k}}_{\mathrm{c}}=2$, meaning the shock spreads through two grid cells during this stage. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=2$, indicating that the shock spreads through two grid cells during this stage. The marching time step was set to _Δ$\mathrm{t}=0.021$. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results are presented in Figure 23, Figure 24, Figure 25, Figure 26, Figure 27 and Figure 28, comparing the proposed method with the previous approach [19] and the data obtained by Harten [3] using the ULT1C scheme at $\mathrm{t}=0.63$. It is important to note that more advanced numerical schemes, such as WENO [24,25], exhibit even higher numerical viscosity for shocks. Consequently, the ULT1C scheme was selected for comparison. Figure 23, Figure 24, Figure 25, Figure 26, Figure 27 and Figure 28 illustrate the absence of numerical viscosity for the shock, in contrast to Harten's method. The proposed scheme allows us to obtain more universal and accurate results than the previous method [19], and it has no sawtooth oscillations (dashed blue line).

The seventh example is the interaction of two Riemann problems with the initial conditions given by (65).

$$\begin{gathered} {\mathrm{p}}_{\mathrm{j}}^0=100,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^{\mathrm{k}}=5,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}<0,, \\ {\mathrm{p}}_{\mathrm{j}}^0=1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^{\mathrm{k}}=1,\mathrm{if}0\le {\mathrm{x}}_{\mathrm{j}}<0.4, \\ {\mathrm{p}}_{\mathrm{j}}^0=0.1,{\mathrm{u}}_{\mathrm{j}}^0=0,{\rho }_{\mathrm{j}}^0=0.125,\mathrm{if}{\mathrm{x}}_{\mathrm{j}}\ge 0.4. \end{gathered}$$

(65)

Also, this is a more complicated problem proposed by Sod [21] (second and third lines of the formula), to which another discontinuity decay (stronger than in (64)) was added (using the first line of the formula (65)).

The modeling region was defined as $\mathrm{x}\in \left[-2;2\right]$. The grid consisted of 400 cells, with a grid spacing of $\mathrm{h}=0.01$. The time step for the acoustics stage was determined based on the Courant criterion (15), with ${\mathrm{k}}_{\mathrm{c}}=3$, meaning the shock spreads through three grid cells during this stage. The time step for the convection stage was determined based on equation (17), with ${\mathrm{k}}_{\mathrm{u}}=11$, indicating that the shock spreads through eleven grid cells during this stage. The marching time step was set to _Δ$\mathrm{t}=0.0236$. The acoustics and convection stages were executed during each marching time step, i.e., ${}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{c}}={}_{\mathrm{\Delta }}\mathrm{t}_{\mathrm{u}}={}_{\mathrm{\Delta }}\mathrm{t}$.

The results are presented in Figure 23, Figure 24, Figure 25, Figure 26, Figure 27 and Figure 28, comparing the proposed method with the previous approach [19] and the data obtained by Harten [3] using the ULT1C scheme at $\mathrm{t}=0.2124$. The proposed scheme allows to obtain the more universal and accurate results than the previous method [19] and it have no sawtooth oscillations (dashed blue line). While the TVD scheme shows strong numerical viscosity of the solution at the shock wave.

Figure 29. Density plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 29. Density plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 30. Magnified view of density plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 30. Magnified view of density plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 31. Velocity plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 31. Velocity plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 32. Magnified view of velocity plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 32. Magnified view of velocity plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 33. Pressure plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 33. Pressure plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 34. Magnified view of pressure plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

Figure 34. Magnified view of pressure plot for the two Riemann problems interaction (65): the ULT1C scheme of Harten (red dash dotted line), the previous method [19] (black solid line), and the present method (blue dashed line).

4. Discussion

The numerical solutions for the aforementioned well-known test cases were compared with the exact solutions, data from Harten’s method, and the results of the author's previous approach. The obtained results lead to the following conclusions:

(1) The primary advantage of the proposed method over the well-known TVD, ENO, and WENO methods is the elimination of the numerical viscosity (diffusion) for shocks;

(2) Over time, shocks or rarefaction waves may spread slightly slower or faster than in the exact solution. This discrepancy arises from rounding the wave front positions to the grid cell accuracy, a consequence of using a fixed homogeneous grid.

(3) The proposed scheme allows us to obtain more universal and accurate results than the previous method [19], and it has no sawtooth oscillations.

The accuracy of the proposed method for handling flow discontinuities is determined by the iterative Godunov exact solver. In the calculations, an iteration termination criterion of less than 10^-5 is used to ensure convergence by measuring the pressure difference between consecutive iterations.

The directions for further research are connected with one-dimensional, two-dimensional Sedov problem modeling and two-dimensional Riemann problem simulations.