Scalar Conservation Laws

Abstract

We present a theoretical aspect of conservation laws by using simplest scalar models with essential properties. We start by rewriting the general scalar conservation law as a quasilinear partial differential equation and solve it by method of characteristics. Here we come across with the notion of strong and weak solutions depending on the initial value of the problem. Taking into account a special initial data for the left and right side of a discontinuity point, we get the related Riemann problem. An illustration of this problem is provided by some examples. In the remaining part of the chapter, we extend this analysis to the gas dynamics given in the Euler system of equations in one dimension. The transformations of this system into the Lagrangian coordinates follow by applying a suitable change of coordinates which is one of the main issues of this section. We next introduce a first-order Godunov finite volume scheme for scalar conservation laws which leads us to write Godunov schemes in both Eulerian and Lagrangian coordinates in one dimension where, in particular, the Lagrangian scheme is reformulated as a finite volume method. Finally, we end up the chapter by providing a comparison of Eulerian and Lagrangian approaches.

1. Introduction

We present a general form of scalar conservation laws with further properties including some basic models and provide examples of computational methods for them. The equations described by

$${\partial }_{t}u+{\partial }_{x}f\left(u\right)=0\mathrm{,}t>0\mathrm{,}x\in \mathbb{R}$$E1

in one dimension are known as scalar conservation laws where $u=u\left(t,x\right)$ is the conserved quantity and $f=f\left(u\right)$ is the associated flux function depending on $t$ and $x$. Whenever an initial condition $u\left(0,x\right)=u_0\left(x\right)$ is attached to Eq. (1), the problem is called the Cauchy problem the solution of which is a content of this chapter. The outlook of chapter is as follows. We introduce basic concepts and provide particular examples of scalar conservation laws in the first part. The equation of gas dynamics in Eulerian coordinates in one dimension is the main issue of the second part. After providing further instruction for these equations, we provide a transformation of the Eulerian equations in the Lagrangian coordinates. In the final part, we give as an example of computational methods for conservation laws, the Godunov schemes for the Eulerian, and the Lagrangian coordinates, respectively.

1.1 Conservation laws: integral form and differential form

We start by investigating the relation of the equations in gas dynamics with conservation laws. We take into account the equation of conservation of mass in one dimension. The density and the velocity are assumed to be constant in the tube where $x$ is the distance and $\rho \left(t,x\right)$ is the density at the time $t$ and at the point $x$. Then if we integrate the density on $\left[x_1,x_2\right]$, we get total mass ${\int }_{x_1}^{x_2}\rho \left(t,x\right)\mathit{dx}$ at time $t$. Assigning the velocity by $u\left(t,x\right)$, then mass flux at becomes $\rho \left(t,x\right)u\left(t,x\right).$ It follows that the rate of change of the mass in $\left[x_1,x_2\right]$ is

$$\frac{d}{\mathit{dt}}{\int }_{x_1}^{x_2}\rho \left(t,x\right)\mathit{dx}=\rho \left(t,x_1\right)u\left(t,x_1\right)-\rho \left(t,x_2\right)u\left(t,x_2\right).$$E2

The last equation is called integral form of conservation law. Integrating this expression in time from $t_1$ to $t_2\mathrm{,}$ we get

$${\int }_{x_1}^{x_2}\rho \left(t_2,x\right)\mathit{dx}-{\int }_{x_1}^{x_2}\rho \left(t_1,x\right)\mathit{dx}={\int }_{t_1}^{t_2}\rho \left(t,x_1\right)u\left(t,x_1\right)\mathit{dt}-{\int }_{t_1}^{t_2}\rho \left(t,x_2\right)u\left(t,x_2\right)\mathit{dt}.$$E3

Using the fundamental theorem of calculus after reduction of Eq. (3), it follows that

$$\rho \left(t,x_2\right)u\left(t,x_2\right)-\rho \left(t,x_1\right)u\left(t,x_1\right)={\int }_{x_1}^{x_2}{\partial }_x\left(\rho \left(t,x\right)u\left(t,x\right)\right)\mathit{dx}.$$E4

As a result, we get

$${\int }_{t_1}^{t_2}{\int }_{x_1}^{x_2}\left\{{\partial }_t\rho \left(t,x\right)+{\partial }_x\rho (t,x\left)u\right(t,x)\right\}\mathit{dx}\mathit{dt}=0.$$E5

Here the end points of the integrations are arbitrary; that is, for any $\left[x_1,x_2\right]$ and $\left[t_1,t_2\right]$, the integrant must be zero. It follows that the conservation of mass yields

$${\partial }_t\rho +{\partial }_x\left(\mathit{u\rho }\right)=0\mathrm{,}$$E6

which is said to be the differential form of the conservation law.

1.2 A first-order quasilinear partial differential equations

A general solution to a quasilinear partial differential equation of the form

$$a\left(t,x,u\right){\partial }_{t}u+b\left(t,x,u\right){\partial }_{x}u=c\left(t,x,u\right)$$E7

where $a\mathrm{,}b\mathrm{,}c$ are non-zero and smooth on a given domain $D\in {\mathbb{R}}^3$ follows by the characteristic method where the characteristic curves are defined by

$$\frac{\mathit{dt}}{a\left(t,x,u\right)}=\frac{\mathit{dx}}{b\left(t,x,u\right)}=\frac{\mathit{du}}{c\left(t,x,u\right)}.$$E8

By applying a parametrization of $c$, the relation (8) is transformed to a system of ordinary differential equation (ODE):

$$\frac{\mathit{dt}}{\mathit{dc}}=a\left(t,x,u\right)\mathrm{,}\frac{\mathit{dx}}{\mathit{dc}}=b\left(t,x,u\right)\mathrm{,}\frac{\mathit{du}}{\mathit{dc}}=c\left(t,x,u\right).$$E9

In addition to these equations, if an initial condition $u_0=u\left(x_0\right)$ is also given, then by the existence theorem of ODE, there is a unique characteristic curve passing from each point $\left(t_0,x_0,u_0\right)$ leading to an integral surface which is the solution to Eq. (7).

Observe that the scalar conservation law (1) is a particular example of Eq. (7) if we assign $a\left(t,x,u\right)=1$, $b\left(t,x,u\right)u_x={\left(f\left(u\right)\right)}_x$, and $c\left(t,x,u\right)=0$. The conserved quantity can be observed by integrating equation (1) over $\left[x_0,x_1\right]$. Indeed

$$\begin{array}{cc}\frac{d}{\mathit{dt}}{\int }_{x_0}^{x_1}u\left(t,x\right)\mathit{dx}& ={\int }_{x_0}^{x_1}{\partial }_{t}u\left(t,x\right)\mathit{dx}=-{\int }_{x_0}^{x_1}f{\left(u\left(t,x\right)\right)}_x\mathit{dx}\\ & =f\left(u\left(t,x_1\right)\right)-f\left(u\left(t,x_0\right)\right)\\ & =\left[\text{inflow}\mathrm{at}\text{the}\text{point}{x}_1\right]-\left[\text{outflow}\mathrm{at}\text{the}\text{point}{x}_0\right].\end{array}$$E10

This means, the quantity $u\left(t,x\right)$ is conserved so that it depends on the difference of the flux functions between the points $x_0$ and $x_1.$

1.3 Strong (classical) solutions

We consider the initial value problem

$$\begin{array}{l}{\partial }_{t}u+{\partial }_x\left(f\left(u\right)\right)=0\mathrm{,}t>0\mathrm{,}x\in \mathbb{R}\\ u\left(0,x\right)=u_0\left(x\right)\mathrm{,}x\in \mathbb{R}\end{array}$$E11

where the initial data is assumed to be continuously differentiable, that is, $u_0\left(x\right)\in C^1\left(\mathbb{R}\right)$. Applying the chain rule to the relation (11), it follows that

$$\begin{array}{l}{\partial }_{t}u+f^{\prime }\left(u\right){\partial }_{x}u=0\mathrm{,}t>0\mathrm{,}x\in \mathbb{R}\mathrm{,}\\ u\left(0,x\right)=u_0\left(x\right)\mathrm{,}x\in \mathbb{R}\mathrm{,}\end{array}$$E12

where we define characteristic curves of Eq. (12) to be the solution of $\frac{d}{\mathit{dt}}x\left(t\right)=f^{\prime }\left(u\left(t,x\left(t\right)\right)\right)=f^{\prime }\left(u\right)$. Then a solution to the system (12) in a domain $\mathrm{\Omega }\in \mathbb{R}$ is said to be a strong (or classical) solution if it satisfies Eq. (11), and it is continuously differentiable on a domain $\mathrm{\Omega }\in \mathbb{R}.$ Let u be a strong solution and the initial data u₀ be differentiable. Observe that (12) is equivalent to a quasilinear form:

$${\partial }_{t}u+\lambda \left(u\right){\partial }_{x}u=0\mathrm{,}$$E13

with $\lambda \left(u\right)=f^{\prime }\left(u\right)$. Applying the method of characteristics to Eq. (13), the partial differential equation is transformed to a system of ordinary differential equations. We consider the characteristic curve passing through the point $\left(0,x_0\right):$

$$\begin{array}{l}{\partial }_{t}x=\lambda \left(u\left(t,x\left(t\right)\right)\right)\\ x\left(0\right)=x_0.\end{array}$$E14

Along this characteristic curve,

$${\partial }_{t}u\left(t,x\left(t\right)\right)={\partial }_{t}u\left(t,x\left(t\right)\right)+{\partial }_{t}x{\partial }_{x}u\left(t,x\left(t\right)\right)={\partial }_{t}u+\lambda \left(u\right){\partial }_{x}u=0$$E15

is satisfied, that is, $u$ is constant. Hence, the characteristic curves are straight lines satisfying

$$x=x_0+\lambda \left(u_0\left(x_0\right)\right)t=0.$$E16

Hence we can define smooth solutions by $u\left(t,x\right)=u_0\left(x_0\right)$. If the slope of the characteristics is $m_{\mathit{char}}=\frac{1}{\lambda \left(u_0\left(x_i\right)\right)}\mathrm{,}$ then depending on the behavior of $\lambda$, the solution takes different forms. If $\lambda \left(u_0\left(x\right)\right)$ is increasing, then the slopes of the characteristics are decreasing. As a result, the characteristics do not intersect, and thus solution can be defined for all $t$ which is greater than zero. On the other hand, if $\lambda \left(u_0\left(x\right)\right)$ is decreasing, then the slopes of the characteristics will be increasing which implies that the characteristics intersect at some point. But at the intersection point, solution cannot take both values $u_0\left(x_1\right)$ and $u_0\left(x_2\right)$. Therefore, we cannot define the strong solution for all $t>0$.

1.4 Linear advection equation

The basic example of the scalar conservation law is the linear advection equation. It can be obtained by setting $a\left(t,x,u\right)=1$, $b\left(t,x,u\right)=\lambda$, and $c\left(t,x,u\right)=0$ in Eq. (7). The flux function takes the form $f\left(u\right)=\mathit{\lambda u}$ where $\lambda$ is a constant. Then the following quasilinear partial differential equation

$${\partial }_{t}u+\lambda {\partial }_{x}u=0$$E17

is a linear advection equation. Similar to Eqs. (11) and (12), an initial value problem for linear advection equation is described by

$$\begin{array}{l}{\partial }_{t}u+{\partial }_{x}f\left(u\right)=0\mathrm{,}-\infty <x<\infty \mathrm{,}t\ge 0\mathrm{,}\\ u\left(0,x\right)=u_0\left(x\right)=f\left(x_0\right)\mathrm{,}-\infty <x<\infty .\end{array}$$E18

Applying the method of characteristics, it follows that $\frac{\mathit{dt}}{1}=\frac{\mathit{dx}}{\lambda }=\frac{\mathit{du}}{0}$ or equivalently

$$u=c_1\mathrm{,}\frac{\mathit{dx}}{\mathit{dt}}=\lambda =c_1\mathrm{,}x=c_{1}t+c_2\mathrm{,}$$E19

where $c_1$ and $c_2$ are constant and $x-\mathit{\lambda t}=c_2.$ As a conclusion, the solution is

$$u\left(t,x\right)=u_0\left(x-\mathit{\lambda t}\right)\mathrm{,}t\ge 0.$$E20

Here $\lambda$ is the $\text{wave speed}$, and the characteristic lines $x-\mathit{\lambda t}=c_2$ are wavefronts which are constants.

1.5 Burgers’ equation

Burgers’ equation is the simplest nonlinear partial differential equation and is the one of the most common models used in the scalar conservation laws and fluid dynamics. The classical Burgers’ equation is described by

$${\partial }_{t}u+u{\partial }_{x}u=\nu {\partial }_{\mathit{xx}}u\mathrm{,}$$E21

where $\nu {\partial }_{\mathit{xx}}u$ is the viscosity term. Equation (21) can be considered as a combination of nonlinear wave motion and linear diffusion term so that it is balance between time evolution, nonlinearity, and diffusion. The term $u{\partial }_{x}u$ is a convection term that may have an effect to wave breaking, and the term $\nu {\partial }_{\mathit{xx}}u$ is a diffusion term that may cause to efface the wave breaking and to flatten discontinuities, and thus we expect to achieve a smooth solution. We try to find a traveling wave solution of Eq. (21) of the form

$$u\left(t,x\right)=g\left(\xi \right)=g\left(x-\lambda \mathit{t}\right)\mathrm{,}\text{with}\xi =x-\mathit{\lambda t}\mathrm{,}$$E22

where $g$ and $\lambda$ are to be determined. Applying the chain rule, we get

$${\partial }_{t}u=-\lambda g^{\prime }\left(\xi \right)\mathrm{,}{\partial }_{x}u=g^{\prime }\left(\xi \right)\mathrm{,}{\partial }_{\mathit{xx}}u=g″\left(\xi \right).$$E23

Plugging these terms in Eq. (21), we get

$$-\lambda g^{\prime }\left(\xi \right)+g\left(\xi \right)g^{\prime }\left(\xi \right)-\nu g″\left(\xi \right)=0.$$E24

Taking integration with respect to $\xi$ gives

$$-\lambda \mathit{g}+\frac{1}{2}{g}^2-\nu g^{\prime }=C\mathrm{,}C:\mathit{\text{constant}}.$$E25

Rewriting Eq. (25) by

$$\left(g-g_1\right)\left(g-g_2\right)=g^2-2\lambda \mathit{g}-2C=2\mathit{\nu }\mathit{dg}/\mathit{d\xi }\mathrm{,}$$E26

it follows that $g_{1\mathrm{,}2}=\lambda \pm \sqrt{{\lambda }^2+2C}$. Supposing that $g_1\mathrm{,}{g}_2$ are real implies $g_1>g_2$. Using separation of variable and then integrating equation (26), we get

$$g\left(\xi \right)=\frac{g_1+g_2{e}^{\left(\frac{g_1-g_2}{2\nu }\right)\xi }}{1+e^{\left(\frac{g_1-g_2}{2\nu }\right)\xi }}=\frac{g_1+g_2}{2}-\frac{g_1-g_2}{2}tanh\left(\frac{g_1-g_2}{4\nu }\xi \right)$$E27

As a result the explicit form of traveling wave solution of Eq. (21) becomes

$$u\left(t,x\right)=\lambda -\frac{g_1-g_2}{2}tanh\left(\frac{1}{4\nu }\left(g_1-g_2\right)\left(x-\mathit{\lambda t}\right)\right)$$E28

where $\lambda =\frac{g_1+g_2}{2}$ is the wave speed. We can observe that $\underset{\xi \to -\infty }{lim}g\left(\xi \right)=g_1$ and $\underset{\xi \to \infty }{lim}g\left(\xi \right)=g_2$ with $g^{\prime }\left(\xi \right)<0$ for all $\xi$. This means the solution $g\left(\xi \right)$ decreases monotonically with $\xi$ from the value $g_1$ to $g_2$. At $\xi =0$, $u=\frac{g_1+g_2}{2}=\lambda \mathrm{,}$ that is the wave form $g\left(\xi \right)$ travels from left to right with speed $\lambda$ equal to the average value of its asymptotic values. The solution resembles to a shock form as it connects the asymptotic states $g_1$ and $g_2$. Without the viscosity term, the solutions to Burgers equation allow shock forms which finally break. The diffusion term prevents incrementally deformation of the wave and its breaking by withstanding the nonlinearity. As a conclusion, there exists a balance between nonlinear advection term and the linear diffusion term. The wave form is notably affected by the diffusion coefficient $\nu$. If $\nu$ is smaller, then the transition layer between two asymptotic values of solution is sharper. In the limit $\nu \to 0$, the solutions converge to the step shock wave solutions to the inviscid Burgers’ equation.

Remark. If the initial data is smooth and very small, then the $u_{\mathit{xx}}$ term is negligible compared to other terms before the beginning of wave breaking. As the wave breaking starts, the $u_{\mathit{xx}}$ term raises faster than $u_x$ term. After a while, the term $u_{\mathit{xx}}$ becomes comparable to the other terms so that it keeps the solution smooth, giving rise to avoid breakdown solutions.

1.6 Inviscid Burgers’ equation

Whenever $\nu =0$, Eq. (21) is called the inviscid Burgers’ equation. This equation can be obtained by substituting $f\left(u\right)=u^2/2$ in the scalar conservation law (1), that is

$${\partial }_{t}u+{\partial }_x\left(u^2/2\right)={\partial }_{t}u+u{\partial }_{x}u=0.$$E29

Observe that $f\left(u\right)$ is a nonlinear function of $u$; thus, the inviscid Burgers’ equation is a nonlinear equation. Equation (29) is now equivalent to Eq. (17) with $\lambda =u$. We know the solution of Eq. (17); so, plugging $\lambda =u$ into the relation (20) implies that the solution of Eq. (29) is

$$u\left(t,x\right)=f\left(x-\mathit{ut}\right)=u_0\left(x-\mathit{ut}\right).$$E30

Recall that the characteristic speed $\lambda$ is constant for linear advection equation; that is, the characteristic curves become parallel for Eq. (17). In contrast, for the inviscid Burgers’ equation (29), the characteristic speed $\lambda =u$ depends on $u$. As a result the characteristic lines are not parallel. If we apply the implicit function theorem to Eq. (29), the solution can be written as a function of $t$ and $x$ as $u_0$ is differentiable. More particularly, differentiating Eq. (30) with respect to $t$, we get

$${\partial }_{t}u=-u_0^{\prime }\left(u_{t}t+u\right)\Rightarrow {\partial }_{t}u=-\frac{u_0^{\prime }u}{1+u_0^{\prime }t}\mathrm{;}$$E31

and differentiating equation (30) with respect to $x$, we get

$${\partial }_{x}u=u_0^{\prime }\left(1-u_{x}t\right)\Rightarrow {\partial }_{x}u=\frac{u_0^{\prime }}{1+u_0^{\prime }t}.$$E32

Thus, substituting Eqs. (31) and (32) in (29), we can recover the inviscid Burgers’ equation. Consequently, the relations (31) and (32) imply that the solutions of Eq. (1) and particularly of Eq. (29) depend on the initial value $u_0$. It can be observed that whenever $u_0^{\prime }\left(x\right)>0$, then by Eq. (32), ${\partial }_{x}u$ decreases in time because $1+u_{0^{\prime }}t>0$ for $t>0$. In other words, the profile of the wave flattens as time increases. On the other hand, whenever $u_0^{\prime }\left(x\right)<0$, then ${\partial }_{x}u$ increases in time as $1+u_{0}t<0.$ Hence $u_x$ in Eq. (32) tends to $\infty$ as $1+u_0^{\prime }t$ approaches to zero. As a result, wave profile become sharp after some time. For further details on the Burgers’ equations, we refer the reader to [12, 13, 22] and the references therein.

1.7 Shock waves

Let the constants u_L and u_R are given with a linear function, $\phi \left(t\right)=\mathit{\lambda t}$. Then

$$u\left(t,x\right)=\left\{\begin{array}{cc}{u}_R\text{if}& x>\mathit{\lambda t}\mathrm{,}\\ u_L\text{if}& x<\mathit{\lambda t}\mathrm{,}\end{array}\right)$$E33

is a simple example of discontinuous solution of the conservation law (11). If $u_L\ne u_R$, the relation (33) is called a shock wave connecting $u_L$ to $u_R$ with shock speed $\lambda$. As an example, if we take into account the characteristics of the inviscid Burgers’ equations which are of the form $\frac{\mathit{dx}}{\mathit{dt}}=u\left(t,x\right)\mathrm{,}$ it follows that

$$x\left(t\right)=u_0\left(x_0\right)t+x_0$$E34

where $u_0\left(x\right)=u\left(0,x\right)$ and $x_0=x\left(0\right)$; thus, the characteristics are straight lines. Depending on the behavior of these characteristics, we have two cases. If $u_L>u_R\mathrm{,}$ characteristics intersect, the solution will have an infinite slope, and the wave will break; as a result a shock is obtained. This is illustrated in Figure 1. On the other hand, if $u_R>u_L$, the characteristics do not intersect, and hence a region without characteristic will appear which is physically unacceptable. This is shown in Figure 2. We get rid of this by introducing the rarefaction waves.

Figure 1.

For the initial value $u_L>u_R$, characteristics, and shock wave.

Figure 2.

For the initial value $u_R>u_L$, characteristics and rarefaction waves.

1.8 Rarefaction waves

A rarefaction wave is a strong solution which is a union of characteristic lines. A rarefaction fan is a collection of rarefaction waves. These waves are constant on the characteristic line $x-x_0=\mathit{\alpha t}$. Here $\alpha \in \left(f^{\prime }\left(u_L\right),f^{\prime }\left(u_R\right)\right)$ where $u_L$ and $u_R$ are the values of $u$ at the edge of the rarefaction wave fan. If moreover $f^{\prime }$ is invertible, then the solution $u=u\left(t,x\right)$ satisfies

$$u\left(x,t\right)={\left(f^{\prime }\right)}^{-1}\left(\frac{x-x_0}{t}\right).$$E35

If, for instance, f is convex, then the rarefaction waves are increasing. If we consider again the inviscid Burgers’ equation with the initial values, then the region without characteristics in Figure 2 will be covered by rarefaction solution which is described by

$$u\left(t,x\right)=\left\{\begin{array}{lll}{u}_L& \text{if}& x/t\le f^{\prime }\left(u_L\right)\mathrm{,}\\ {\left(f^{\prime }\right)}^{-1}\left(x/t\right)& \text{if}& f^{\prime }\left(u_L\right)\le x/t\le f^{\prime }\left(u_R\right)\mathrm{,}\\ u_R& \text{if}& f^{\prime }\left(u_R\right)\le x/t.\end{array}\right.$$E36

An illustration of rarefaction waves and rarefaction fan in Eq. (36) is given in Figure 3.

Figure 3.

Rarefaction fan.

Remark. Whenever characteristics intersect, we may have multiple valued solution or no solution; but we have no more classical (strong) solution. To get rid of this situation, we introduce a more wide-ranging notion of solution, the weak solution, in the next part. By this arrangement, we may have non-differentiable and even discontinuous solutions.

1.9 Weak solution

Weak solutions occur whenever there is no smooth (classical) solution. These solutions may not be differentiable or even not continuous. Considering $\varphi :\mathbb{R}\times {\mathbb{R}}^{+}\to \mathbb{R}$ as a smooth test function with a compact support and multiplying the scalar conservation law (1) by this test function $\varphi$, it follows after integration by parts that

$$\begin{array}{l}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\varphi {\partial }_{t}u+\varphi {\partial }_{x}f\left(u\right)\mathit{dx}\mathit{dt}\\ ={\int }_{-\infty }^{\infty }{\left.\mathit{\varphi u}\right|}_0^{\infty }\mathit{dx}-{\int }_0^{\infty }{\int }_{-\infty }^{\infty }u{\partial }_t\varphi \mathit{dx}\mathit{dt}+{\int }_0^{\infty }{\left.\mathit{\varphi f}\left(u\right)\right|}_{-\infty }^{\infty }\mathit{dt}-{\int }_0^{\infty }{\int }_{-\infty }^{\infty }f\left(u\right){\partial }_x\varphi \mathit{dx}\mathit{dt}\\ {\left.=-{\int }_0^{\infty }{\int }_{-\infty }^{\infty }u{\partial }_t\varphi \mathit{dx}\mathit{dt}-{\int }_0^{\infty }{\int }_{-\infty }^{\infty }f\left(u\right){\partial }_x\varphi \mathit{dx}\mathit{dt}-{\int }_{-\infty }^{\infty }\mathit{u\varphi }\right|}_{t=0}\mathit{dx}.\end{array}$$E37

Putting the initial condition $u_0\left(x\right)=u\left(0,x\right)$ to the above relation, it follows that

$${\int }_0^{\infty }{\int }_{-\infty }^{\infty }u{\varphi }_t+f\left(u\right){\varphi }_x\mathit{dx}\mathit{dt}+{\int }_{-\infty }^{\infty }u\left(0,x\right)\varphi \left(x\right)\mathit{dx}=0.$$E38

Observe that there are no more derivatives of $u$ and $f$ which may lead less smoothness. In other words, the smoothness requirement is reduced for finding a solution. Thus, the function $u\left(t,x\right)$ is said to be the weak solution of the initial value problem (11) if the relation (38) satisfied for all test function $\varphi .$ Here it is significant to note that $u$ needs not be smooth or continuous to satisfy Eq. (38). Consequently, by weak solutions, we extend the solutions so that discontinuous solutions may also be covered. However, in general weak solutions are not unique. We can also notice that strong solutions are also weak solutions and a weak solution which is continuous and piecewise differentiable is also strong solution.

1.10 Riemann problem

The Riemann problem is a Cauchy problem with a particular initial value which consists a conservation law together with piecewise constant data having a single discontinuity. We consider the Riemann problem for a convex flux described by

$$\begin{array}{l}{\partial }_{t}u+{\partial }_x\left(f\left(u\right)\right)=0\mathrm{,}x\in \mathbb{R}\mathrm{,}t\in {\mathbb{R}}_{+}\mathrm{,}\\ u\left(0,x\right)=\left\{\begin{array}{cc}{u}_L\text{if}& x<0\mathrm{,}\\ u_R\text{if}& x>0.\end{array}\right.\end{array}$$E39

The solution is a set of shock and rarefaction waves depending on the relation between $u_L$ and $u_R.$ There are two cases to investigate:

Case 1: $\left(u_L>u_R\right)$ A shock is obtained because the left-hand side wave moves faster than the right-hand side one. Thus the solution

$$u\left(t,x\right)=\left\{\begin{array}{cc}{u}_L\text{if}& x/t<\lambda \mathrm{,}\\ u_R\text{if}& x/t>\lambda \mathrm{,}\end{array}\right.$$E40

is a shock wave satisfying the shock speed $\lambda =\frac{f\left(u_R\right)-f\left(u_L\right)}{u_R-u_L}.$

Case 2: ($u_L<u_R$) The solution given in Case 1 is also a solution for this case. In addition, we have rarefaction solutions of the form (36) illustrated by Figure 3.

1.11 Rankine-Hugoniot jump condition

A jump discontinuity along the characteristic line is controlled by the Rankine-Hugoniot jump condition. Integrating the scalar conservation law (1) in $\left[x_1,x_2\right]$, it follows that

$$\frac{d}{\mathit{dt}}{\int }_{x_1}^{x_2}{\left.u\left(t,x\right)\mathit{dx}+f\left(u\right)\right|}_{x_1}^{x_2}=0.$$E41

Suppose that there is a discontinuity at the point $x=\xi \left(t\right)\in \left(x_1,x_2\right)$ where $u$ and $u^{\prime }$ are continuous on the $\left[x_1,\xi \left(t\right)\right)$ and $\left(\xi \left(t\right),x_2\right]$, respectively. Suppose also that whenever $x_1\to \xi {\left(t\right)}^{-}$ and $x_2\to \xi {\left(t\right)}^{+}$, their limits exist. Next, Eq. (41) can be rewritten as

$$\frac{d}{\mathit{dt}}{\int }_{x_1}^{\xi \left(t\right)}u\left(t,x\right)\mathit{dx}+\frac{d}{\mathit{dt}}{\int }_{\xi \left(t\right)}^{x_2}u\left(t,x\right)\mathit{dx}=-\left(f\left(t,x_2\right)-f\left(t,x_1\right)\right).$$E42

By the fundamental theorem of calculus, the relations (41) and (42) yield

$$u\left({\xi }^{-},x\right){\xi }^{\prime }\left(t\right)-u\left({\xi }^{+},x\right){\xi }^{\prime }\left(t\right)+\frac{d}{\mathit{dt}}{\int }_{x_1}^{\xi \left(t\right)}{u}_t\left(t,x\right)\mathit{dx}+\frac{d}{\mathit{dt}}{\int }_{\xi \left(t\right)}^{x_2}{u}_t\left(t,x\right)\mathit{dx}.$$E43

Taking the limit whenever $x_1\to \xi {\left(t\right)}^{-}$ and $x_2\to \xi {\left(t\right)}^{+}$, it follows that

$${\xi }^{\prime }\left(t\right)\left(x_2-x_1\right)=f\left(x_2\right)-f\left(x_1\right)\Rightarrow \lambda ={\xi }^{\prime }\left(t\right)=\frac{f\left(x_2\right)-f\left(x_1\right)}{x_2-x_1}.$$E44

The relation (44) is said to be the Rankine-Hugoniot jump condition. Geometrical meaning of the Rankine-Hugoniot jump condition is that the shock speed is the slope of the secant line through the points $\left(u_L,f\left(u_L\right)\right)$ and $\left(u_R,f\left(u_R\right)\right)$ on the graph of $f$.

1.12 Entropy functions

Entropy and entropy flux are defined for attaining physically meaningful solutions. If $u$ is the smooth solution of the conservation law (1), then the relation

$${\partial }_{t}G\left(u\right)+{\partial }_{x}F\left(u\right)=0$$E45

is satisfied for continuously differentiable functions $G$ and $F$ where the pair $\left(G,F\right)$ is called as entropy pair so that $G$ is entropy and $F$ is entropy flux. If in addition $u$ is smooth, then Eq. (45) becomes

$$G^{\prime }\left(u\right){\partial }_{t}u+F^{\prime }\left(u\right){\partial }_{x}u=0$$E46

which looks like to the scalar conservation law (1). Indeed, if we multiply Eq. (1) by $G^{\prime }\left(u\right)$, it follows that

$$G^{\prime }\left(u\right){\partial }_{t}u+G^{\prime }\left(u\right)f^{\prime }\left(u\right){\partial }_{x}u=0.$$E47

It follows that Eqs. (46) and (47) are equivalent with $F^{\prime }\left(u\right)=G^{\prime }\left(u\right)f^{\prime }\left(u\right).$ Here the function $u\left(t,x\right)$ is said to be the entropy solution of Eq. (1) if

$${\partial }_{t}G\left(u\right)+{\partial }_{x}F\left(u\right)\le 0$$

holds for all convex entropy pairs $\left(G\left(u\right),F\left(u\right)\right)$.

1.13 Entropy condition

Weak solutions to conservation laws may contain discontinuities as a result of a discontinuity in the initial data or of characteristics that cross each other or because of the jump conditions which are satisfied across the discontinuities. Although the Rankine-Hugoniot jump condition is satisfied, the uniqueness of the solution may always not be guaranteed. In order to eliminate the nonphysical solutions among the weak solutions, we need an additional condition, so-called entropy condition. It is described by the following: A discontinuity propagating with the characteristic speed $\lambda$ given by the Rankine-Hugoniot jump condition satisfies the entropy condition if holds.

$$f^{\prime }\left(u_L\right)>\lambda >f^{\prime }\left(u_R\right)$$E48

Example 1.1. The weak solutions to conservation laws need not be unique. If we write the inviscid Burgers’ equation in quasilinear form and multiply by $2u$, we obtain $2u{\partial }_{t}u+2u^2{\partial }_{x}u=0$. In conservative form it becomes

$${\partial }_t\left(u^2\right)+{\partial }_x\left(\frac{2}{3}{u}^3\right)=0\mathrm{,}\text{with}f\left(u^2\right)=\frac{2}{3}{\left(u^2\right)}^{3/2}.$$E49

The inviscid Burgers’ equation and Eq. (49) have exactly the same smooth solutions. But their weak solutions are different. A shock traveling speed for the inviscid Burgers’ equation is ${\lambda }_1=\left(u_L+u_R\right)/2$; however for Eq. (49), we have ${\lambda }_2=(\frac{2}{3}\left(\frac{u_L^3-u_R^3}{u_L^2-u_R^2}\right)$. That is ${\lambda }_1\ne {\lambda }_2$ whenever $u_L\ne u_R$, and thus these two equations have different weak solutions.

Example 1.2. We first consider the initial value problem for $u_L>u_R$ given by

$${\partial }_{t}u+{\partial }_x\left(u^2/2\right)=0\mathrm{,}{u}_0=\left\{\begin{array}{cc}1& \text{if}x\le 0\mathrm{,}\\ 0& \text{if}x>0.\end{array}\right.$$E50

Applying the method of characteristics for $t>0$, it follows that

$$\frac{\mathit{du}}{\mathit{dt}}=0\mathrm{,}\frac{\mathit{dx}}{\mathit{dt}}=\left\{\begin{array}{cc}1& \text{if}x\le 0\mathrm{,}\\ 0& \text{if}x>0.\end{array}\right.$$E51

Next if we integrate Eq. (51) with respect to $t$, we get the characteristic curves

$$x=\left\{\begin{array}{cc}t-c& \text{if}x\le 0\mathrm{,}\\ b& \text{if}x>0\mathrm{,}\end{array}\right.$$E52

where $c>0$ and $b$ are constants. Due to the discontinuity at the point $x=0$, there is no strong (classical) solution. The speed of propagation is $\lambda =\frac{u_L+u_R}{2}=0.5.$ Moreover, the weak solution for $t\le \lambda =0.5$ becomes

$$u\left(t,x\right)=\left\{\begin{array}{cc}1& \text{if}\frac{x}{t}\le 0.5\\ 0& \text{if}\frac{x}{t}>0.5\end{array}\right.\mathrm{,}$$E53

which satisfies both the jump condition and the entropy condition as $u_L=1>u_R=0$. The characteristic curves can be observed in Figure 4.

Figure 4.

For initial value $u_L>u_R$, the characteristic solutions.

Example 1.3. We now interchange the roles of $u_L$ and $u_R$ of the Example 1.2 so that $u_L<u_R$ to get an initial value problem:

$${\partial }_{t}u+{\partial }_x\left(u^2/2\right)=0\mathrm{,}{u}_0=\left\{\begin{array}{cc}0& \text{if}x\le 0\mathrm{,}\\ 1& \text{if}x>0.\end{array}\right.$$E54

By the method of characteristics, we obtain a solution

$$u_1\left(t,x\right)=\left\{\begin{array}{cc}0& \text{if}\frac{x}{t}\le 1\\ 1& \text{if}\frac{x}{t}>1\end{array}\right.$$E55

which is a classical (strong) solution on both sides of the characteristic line $\frac{x}{t}=1$. Since it satisfies the Rankine-Hugoniot jump condition along the discontinuity curve, it is a weak solution. However, the entropy condition is not satisfied. It yields an empty region between the characteristic lines shown in Figure 4. In order to cover this empty state, we consider another solution described by

$$u_2\left(t,x\right)=\left\{\begin{array}{lll}0& \text{if}& x\le 0\mathrm{,}\\ \frac{x}{t}& \text{if}& 0\le \frac{x}{t}\le 1\mathrm{,}\\ 1& \text{if}& \frac{x}{t}\ge 1\end{array}\right.$$E56

which satisfies both jump and entropy conditions. Here we can observe the rarefaction fan arising on the interval $0\le \frac{x}{t}\le 1$. An illustration of this solution is supplied in Figure 5.

Figure 5.

For initial value $u_L<u_R$, characteristic solutions $u_1\left(t,x\right)$ and $u_2\left(x,t\right)$ with rarefaction fan.

2. The gas dynamic equations in one dimension

The equation of fluid dynamics can be represented in Eulerian and Lagrangian forms. Eulerian coordinates are related to the coordinates of a fixed observer. On the other hand, Lagrangian coordinates are in usual related to the local flow velocity. That is, due to the velocity taking different values in different parts of the fluid, the change of coordinates is different from one point to another one.

2.1 Eulerian coordinates

The equations of gas dynamics in Eulerian coordinates can be written in the following conservative forms:

$$\left\{\begin{array}{l}{\partial }_t\left(\rho \right)+{\partial }_x\left(\mathit{\rho u}\right)=0\mathrm{,}\\ {\partial }_t\left(\mathit{\rho u}\right)+{\partial }_x\left(\rho u^2+p\right)=0\mathrm{,}\\ {\partial }_t\left(\mathit{\rho e}\right)+{\partial }_x\left(\left(\mathit{\rho e}+p\right)u\right)=0\end{array}\right.$$E57

where we ignored the heat conduction. If we denote

$$U=\left(\begin{array}{c}\rho \\ \mathit{\rho u}\\ \mathit{\rho e}\end{array}\right)\mathrm{,}F\left(U\right)=\left(\begin{array}{c}\mathit{\rho e}\\ \rho u^2+p\\ \mathit{\rho eu}+\mathit{pu}\end{array}\right)\mathrm{,}$$E58

then Eq. (57) can be written by

$${\partial }_{t}U+{\partial }_{x}F\left(U\right)=0$$E59

where $\rho$ is density, $p$ is pressure, $u$ is velocity, and $e$ is the specific internal energy.

2.2 Hyperbolicity of the Euler system

If we do not neglect the heat conduction, then the $U$ and $F$ terms in Eq. (59) become

$$U=\left(\begin{array}{c}\rho \\ \mathit{\rho u}\\ E\end{array}\right)\text{and}F=\left(\begin{array}{c}\mathit{\rho u}\\ \rho u^2+p\\ \left(E+p\right)u\end{array}\right)\mathrm{,}$$E60

where $E$ is total energy such that $E=\frac{1}{2}\rho u^2+\mathit{\rho e}\mathrm{,}$ $e=\frac{p}{\left(\delta -1\right)\rho }$, and for perfect gases $\delta =c_p/c_v$ is the ratio of specific heats. Rewriting Eq. (59) in quasilinear form, we get

$${\partial }_{t}U+A\left(U\right){\partial }_{x}U=0\mathrm{,}$$E61

where $A\left(U\right)=\frac{\partial F}{\partial U}$ is the Jacobian matrix. The eigenvalues of $A\left(U\right)$ then are ${\lambda }_1=u\mathrm{,}{\lambda }_2=u-a\mathrm{,}{\lambda }_3=u+a$ where $a$ is the sound speed given by $a=\sqrt{\frac{\mathit{\delta p}}{\rho }}$. Moreover the corresponding eigenvectors are

$$E^{\left(1\right)}=\left(\begin{array}{c}1\\ u\\ \frac{1}{2}{u}^2\end{array}\right)\mathrm{,}{E}^{\left(2\right)}=\left(\begin{array}{c}1\\ u-a\\ H-\mathit{ua}\end{array}\right)\mathrm{,}{E}^{\left(3\right)}=\left(\begin{array}{c}1\\ u+a\\ H+\mathit{ua}\end{array}\right)$$E62

which are real, and the eigenvectors are linearly independent implying that the Euler equations for perfect gases are hyperbolic.

2.3 Rankine-Hugoniot conditions for the Euler system

Using the results in the previous part, the Rankine-Hugoniot jump conditions for the Euler system will be of the form

$$\begin{array}{c}s\left({\rho }_1-{\rho }_2\right)=m_2-m_1\mathrm{,}\\ s\left(m_2-m_1\right)=\frac{m_2^2}{{\rho }_2}+p_2-\frac{m_1^2}{{\rho }_1}-p_1\mathrm{,}\\ s\left({\rho }_2{E}_2-{\rho }_1{E}_1\right)=m_2{H}_2-H_1{m}_1\mathrm{,}\end{array}$$E63

where the indices $1$ and $2$ refer to the left and right of the shock, respectively, and $s$ denotes the wave speed.

2.4 Riemann problem for the Euler system

The Riemann problem for the one-dimensional Euler equation (57) is represented by

$$\begin{array}{l}{\partial }_{t}U+{\partial }_x\left(F\left(U\right)\right)=0\mathrm{,}x\in \mathbb{R}\mathrm{,}t>0\mathrm{,}\\ U\left(0,x\right)=U_0\left(x\right)=\left\{\begin{array}{cc}{U}_L\text{if}& x<0\mathrm{,}\\ U_R\text{if}& x>0.\end{array}\right.\end{array}$$E64

The reader is addressed to the references [18, 24] for further details.

2.5 Lagrangian coordinates

We aim to transform the equations of gas dynamics (57) given in the Eulerian coordinates into the Lagrangian coordinates for one-dimensional case. We start denoting by $\mathbf{u}=\mathbf{u}\left(t,x\right)$ the velocity field of the fluid flow and consider the differential system

$$\frac{\mathit{dx}}{\mathit{dt}}=\mathbf{u}\left(t,x\right).$$E65

We set the following change of coordinates from Euler coordinates to Lagrange coordinates for space and time as $\left(t,x\right)\to \left(t^{\prime },\xi \right)$ where $\xi =\left({\xi }_1,{\xi }_2,{\xi }_3\right)\in {\mathbb{R}}^3$ so that

$$t^{\prime }=t\mathrm{,}\frac{\partial x\left(t^{\prime },\xi \right)}{\partial t^{\prime }}=\mathbf{u}\left(t^{\prime },x\left(t^{\prime },\xi \right)\right)\mathrm{,}{x}_0=x\left(0,\xi \right)=\xi .$$E66

It follows that $\left(t^{\prime },\xi \right)=\left(t,\left({\xi }_1,{\xi }_2,{\xi }_3\right)\right)$ are the Lagrangian coordinates associated with the velocity field $u$. We set

$$J\left(t,\xi \right)=det\left(\frac{\partial x_i}{\partial {\xi }_j}\left(t,\xi \right)\right)\mathrm{,}$$E67

which gives

$$\frac{\partial J}{\partial t}\left(t,\xi \right)=J\left(t,\xi \right)\left(\mathbf{div}\mathbf{u}\right)\left(t,x\left(t,\xi \right)\right)\mathrm{,}\text{where}\mathrm{,}\mathbf{div}\mathbf{u}=\sum _{j=1}^3\frac{\partial u_j}{\partial x_j}.$$E68

It follows by some algebraic manipulations that the gas dynamic equations become

$$\left\{\begin{array}{l}{\partial }_t\left(\mathit{\rho J}\right)=0\mathrm{,}\left(\text{Conservation}\text{of}\text{mass}\right)\mathrm{,}\\ {\partial }_t\left(\mathit{\rho uJ}\right)+{\partial }_{\xi }\left(p\right)=0\mathrm{,}\left(\text{Conservation}\text{of}\text{momentum}\right)\mathrm{,}\\ {\partial }_t\left(\mathit{\rho eJ}\right)+{\partial }_{\xi }\left(\mathit{pu}\right)=0\mathrm{,}\left(\text{Conservation}\text{of}\text{energy}\right).\end{array}\right.$$E69

In order to derive a more convenient form of the system (69), we derive firstly the equation of conservation of mass:

$$\mathit{\rho J}={\rho }_0=\rho \left(0,\xi \right)$$E70

where ${\rho }_0\left(\xi \right)=\rho \left(0,\xi \right).$ Assuming that $\rho >0$, we introduce the specific volume $\tau =1/\rho$, and by using Eq. (68) we get

$$J={\rho }_0\tau \mathrm{,}\text{and}{\partial }_{t}J=J{\partial }_{x}u={\partial }_{\xi }u$$E71

which yields

$${\rho }_0{\partial }_t\tau -{\partial }_{\xi }u=0.$$E72

Hence the second and third equations of Eq. (69) become

$$\begin{array}{l}{\rho }_0{\partial }_{t}u+{\partial }_{\xi }p=0\mathrm{,}\left(\text{Conservation}\text{of}\text{momentum}\right)\mathrm{,}\\ {\rho }_0{\partial }_{t}e+{\partial }_{\xi }\left(\mathit{pu}\right)=0\mathrm{,}\left(\text{Conservation}\text{of}\text{energy}\right).\end{array}$$E73

Moreover, we define a mass variable $m$ by

$$m\left(\xi \right)={\int }_0^{\xi }\rho \left(0,y\right)\mathit{dy}\mathrm{,}\text{or}\text{equivalently}\mathrm{,}\mathit{dm}=\rho \left(0,\xi \right)\mathit{d\xi }={\rho }_0\mathit{d\xi }.$$E74

Finally, using Eqs. (69) and (73), the Euler system (57) can be written in Lagrangian coordinates with the mass variable in the form

$$\left\{\begin{array}{l}{\partial }_t\tau -{\partial }_{m}u=0\mathrm{,}\\ {\partial }_{t}u+{\partial }_{m}p=0\mathrm{,}\\ {\partial }_{t}e+{\partial }_m\left(\mathit{pu}\right)=0\mathrm{,}\end{array}\right.$$E75

where $p=p\left(\tau ,\xi \right)=p\left(\tau ,e-u^2/2\right)$. If we set $V=\left(\begin{array}{c}\tau \\ u\\ e\end{array}\right)\mathrm{,}F\left(V\right)=\left(\begin{array}{c}-u\\ p\\ \mathit{pu}\end{array}\right)$ with $\tau >0\mathrm{,}u\in \mathbb{R}\mathrm{,}e-u^2/2>0\mathrm{,}$ we obtain a scalar conservation law of the form

$${\partial }_{t}V+{\partial }_{m}F\left(V\right)=0$$E76

which is strictly hyperbolic. This can be verified by checking the Jacobian matrix of the flux calculated with respect to the variables $\left(\tau ,u,e\right)$

$$\left(\begin{array}{ccc}0& -1& 0\\ p_{\tau }& -{\mathit{up}}_{\epsilon }& p_{\epsilon }\\ {\mathit{up}}_{\tau }& p-u^2{p}_{\epsilon }& {\mathit{up}}_{\epsilon }\end{array}\right)$$E77

with $e=\epsilon +\frac{1}{2}{u}^2$. The eigenvalues are ${\sigma }_1=-\sqrt{p_{\tau }-{\mathit{pp}}_{\epsilon }}<{\sigma }_2=0<{\sigma }_3=\sqrt{p_{\tau }-{\mathit{pp}}_{\epsilon }}$ so that they are all distinct, and thus the system is strictly hyperbolic.

In fact there are different versions of the gas dynamics in Lagrangian coordinates. In this part we followed the approaches stated in [9, 10, 12]. For further details we cite these works with references therein.

2.6 Rankine-Hugoniot conditions for the Lagrangian system

Similarly as in the Euler system, the Rankine-Hugoniot jump conditions for the Lagrangian system (79) are of the form

$$\begin{array}{l}\sigma \left({\tau }_1-{\tau }_0\right)=-\left(u_1-u_0\right)\mathrm{,}\\ \sigma \left(u_1-u_0\right)=p_1-p_0\mathrm{,}\\ \sigma \left(e_1-e_0\right)=p_1{u}_1-p_0{u}_0\mathrm{,}\end{array}$$E78

where $\sigma$ denotes the speed of propagation of the discontinuity with respect to the mass variable.

Remark. The Eulerian and Lagrangian Rankine-Hugoniot relations are equivalent. Moreover, Eulerian entropy relations are equivalent to all Lagrangian entropy relations (see [9] for further detail).

Example 2.1. For simplicity of notation, we take $\left(t,x\right)$ as the Lagrangian coordinates. Then the system of equations

$$\left\{\begin{array}{l}{\partial }_t\tau -{\partial }_{x}u=0\mathrm{,}\\ {\partial }_{t}u+{\partial }_{x}p\left(\tau \right)=0\mathrm{,}\end{array}\right.$$E79

is a one-dimensional isentropic gas dynamics in Lagrangian coordinates which is also known as $p$-system. It is the simplest nontrivial example of a nonlinear system of conservation laws. Here $\tau$ is the specific volume, $u$ is the velocity, and the pressure $p=p\left(\tau \right)$ is given as a function of $\tau$ by

$$p\left(\tau \right)={\mathit{\kappa \tau }}^{-\gamma }\mathrm{,}\gamma >0\mathrm{,}\kappa =\frac{{\left(\gamma -1\right)}^2}{4\gamma }.$$E80

The system (79) is equivalent to

$${\partial }_{t}V+{\partial }_{x}f\left(V\right)=0\mathrm{,}\text{with}V=\left(\begin{array}{c}\tau \\ u\end{array}\right)\mathrm{,}f\left(V\right)=\left(\begin{array}{c}-u\\ p\left(\tau \right)\end{array}\right)\mathrm{,}$$E81

where $\tau >0$ and $\left(\tau ,u\right)\in {\mathbb{R}}^2.$ If we assume that $p^{\prime }\left(\tau \right)<0$, it follows that the Jacobian matrix of $f$

$$J\left(f\right)=\left(\begin{array}{cc}0& -1\\ p^{\prime }\left(\tau \right)& 0\end{array}\right)$$E82

has two real distinct eigenvalues ${\sigma }_1=-\sqrt{(-p^{\prime }\left(\tau \right)}<{\sigma }_2=\sqrt{-p^{\prime }\left(\tau \right)}$. In other words, the system (81) is strictly hyperbolic. On the other hand, for the case $p^{\prime }\left(\tau \right)>0$, it becomes elliptic. Moreover, one can verify that the solutions of the $p$-system (79) and the Euler system (57) are equivalent.

3. Godunov schemes

The Godunov scheme deals with solving the Riemann problem forward in time for each grid cell and then taking the mean value over these cells. The Riemann problem is solved per mesh point at each time step iteratively. If there are no strong shock discontinuities, this process may cost much and will not be effective. To get rid of such a situation, we establish approximate Riemann solvers that are easier to implement and also low cost to use. Eulerian and Lagrangian Godunov schemes are current Godunov scheme in literature. Both have advantages and disadvantages depending on the structure of the problem. A brief comparison of the method for these two approaches is presented in the last part of the chapter. In this work we will not go further in numerical examples and details of these methods; instead, we aim to present a general form of Godunov schemes for gas dynamics in Eulerian and Lagrangian coordinate. Before introducing these, we present a first-order Godunov scheme for scalar conservation laws.

3.1 First-order Godunov scheme

Consider the scalar conservation law (1). Godunov scheme is a numerical scheme which takes advantage of analytical solutions of the Riemann problem for the conservation law (1). The numerical flux functions are evaluated at the spatial steps $x_{j-1/2}$ and $x_{j+1/2}$ by handling the solutions of the Riemann problem. On each grid cell ${\mathcal{I}}_i=\left[x_{j-1/2},x_{j+1/2}\right]\mathrm{,}$ we have a piecewise constant function. The Riemann problem for (1) for the left and right sides of ${\mathcal{I}}_i$ are described by

$$u_L\left(x\right)=\left\{\begin{array}{cc}{u}_{j-1}^n& \mathrm{;}x<0\mathrm{,}\\ u_j^n& \mathrm{;}x>0\mathrm{,}\end{array}\right.u_R\left(x\right)=\left\{\begin{array}{cc}{u}_j^n& \mathrm{;}x<0\mathrm{,}\\ u_{j+1}^n& \mathrm{;}x>0\mathrm{,}\end{array}\right.$$E83

respectively. These two solutions to the Riemann problem will be the numerical solution $\tilde{u}\left(t,x\right)$. Once establishing the solution over the mesh $\left[t^n,t^{n+1}\right]$, we approximate the solution at the next time step $t^{n+1}$ by the average value

$$U_j^{n+1}=\frac{1}{\mathrm{\Delta }x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}\tilde{u}\left(x,t^{n+1}\right)\mathit{dx}.$$E84

Proceeding this process, we define the solution $\tilde{u}\left(x,t^{n+1}\right)$ iteratively. Then $U_j^{n+1}$ can be calculated by using the integral form of the conservation law (1) in the following way: We integrate (1) for $u\left(t,x\right)$ over each grid cell $\left[t_n,t_{n+1}\right]\times I_j:$

$$\begin{array}{ll}{\int }_{x_{j-1/2}}^{x_{j+1/2}}{{\displaystyle \tilde{u}}}^n\left(x,t^{n+1}\right)\mathit{dx}& -{\int }_{x_{j-1/2}}^{x_{j+1/2}}{{\displaystyle \tilde{u}}}^n\left(x,t^n\right)\mathit{dx}\\ & ={\int }_{t^n}^{t^{n+1}}f\left({{\displaystyle \tilde{u}}}_{j-1/2}^n\right)\mathit{dt}-{\int }_{t^n}^{t^{n+1}}f\left({{\displaystyle \tilde{u}}}_{j+1/2}^n\right)\mathit{dt}.\end{array}$$E85

Dividing both parts by $\mathrm{\Delta }x$ and using the fact that $\tilde{u}\left(x,t^n\right)=u_j^n$ is constant at the end points $x_{j-1/2}$ and $x_{j+1/2}$, we get

$$u_j^{n+1}=u_j^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\left(f\left({{\displaystyle \tilde{u}}}_{j-1/2}^n\right)-f\left({{\displaystyle \tilde{u}}}_{j+1/2}^n\right)\right).$$E86

Thus, Godunov method is a conservative numerical scheme. It can be restated in an alternative form. Assigning the constant value of $u_j^n$ at the points $x_{j-1/2}$ and $x_{j+1/2}$ by $u^{\ast }(U_{j-1}^n,U_j^n)$ and $u^{\ast }(U_j^n,U_{j+1}^n)$, respectively, the numerical flux functions become

$$\begin{array}{l}f\left({\tilde{u}}_{j-1/2}^n\right)=f\left(u^{\ast }\left(U_{j-1}^n,U_j^n\right)\right)=\mathcal{F}\left(U_{j-1}^n,U_j^n\right)\mathrm{,}\\ f\left({\tilde{u}}_{j+1/2}^n\right)=f\left(u^{\ast }\left(U_j^n,U_{j+1}^n\right)\right)=\mathcal{F}\left(U_j^n,U_{j+1}^n\right).\end{array}$$E87

Therefore, a first-order Godunov method takes the form

$$U_j^{n+1}=U_j^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\left(\mathcal{F}\left(U_j^n,U_{j+1}^n\right)-\mathcal{F}\left(U_{j-1}^n,U_j^n\right)\right).$$E88

Here the constant value of ${\tilde{u}}^n$ depends on the initial data. In other words, the Godunov method considers the Riemann problem as constant in each grid interval ${\mathcal{I}}_i$. It follows that, at the subsequent time stage, the exact solutions of the problem are picked as the numerical fluxes at the grid boundary.

The Godunov method is consistent with the exact solution of the Riemann problem for the conservation law (1). If we suppose that $u_j^n=u_j^{n+1}=\overline{u}$, then ${\tilde{u}}_{j+1/2}^n=\overline{u}$ and $\mathcal{F}\left(\overline{u},\overline{u}\right)=f\left(\overline{u}\right)$. For the stability, CFL condition requires that

$$\underset{x\in \mathcal{R}\mathrm{,}t>0}{sup}\mid f^{\prime }\left(u\left(t,x\right)\right)\mid \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\le 1$$E89

for each $u_j^n$. Next, if assigning $u^{\ast }$ as the intermediate value over the grid ${\mathcal{I}}_i$ in the Riemann solution, it implies that

$$u^{\ast }\left(u_L,u_R\right)=\left\{\begin{array}{cc}{u}_L\mathrm{,}& \lambda >0\mathrm{,}\\ u_R\mathrm{,}& \lambda <0\mathrm{,}\end{array}\right.$$E90

where $\lambda$ is the wave propagation speed. Hence the numerical flux for Godunov’s method can be generalized by

$$f\left(u_L,u_R\right)=\left\{\begin{array}{cc}\underset{u_L\le u\le u_R}{min}f\left(u\right)\mathrm{,}& \text{if}{u}_L\le u_R\mathrm{,}\\ \underset{u_L\ge u\ge u_R}{max}f\left(u\right)\mathrm{,}& \text{if}{u}_R<u_L.\end{array}\right.$$E91

For numerical illustration of Godunov schemes, we cite the articles [14, 20, 27].

3.2 Godunov method in Eulerian coordinates

We consider Eq. (59) with (60). The eigenvalues of $F^{\prime }\left(U\right)$ are ${\sigma }_1=u-c<{\sigma }_2=u<{\sigma }_3=u+c$. Then the Riemann problem at the point $x_{i+1/2}$ between the states $U_i$ and $U_{i+1}$ which is solved by the Godunov scheme can be written by

$$\left\{\begin{array}{l}{\rho }_i^{n+1}={\rho }_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}_i}\left({\left(\mathit{\rho u}\right)}_{i+1/2}^n-{\left(\mathit{\rho u}\right)}_{i-1/2}^n\right)\\ {\left(\mathit{\rho u}\right)}_i^{n+1}={\left(\mathit{\rho u}\right)}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}_i}\left({\left(\rho u^2+p\right)}_{i+1/2}^n-{\left(\rho u^2+p\right)}_{i-1/2}^n\right).\\ {\left(\mathit{\rho e}\right)}_i^{n+1}={\left(\mathit{\rho e}\right)}_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{x}_i}\left({\left(\left(\mathit{\rho e}+p\right)u\right)}_{i+1/2}^n-{\left(\left(\mathit{\rho e}+p\right)u\right)}_{i-1/2}^n\right)\end{array}\right.$$E92

3.3 Godunov method in Lagrangian coordinates

Consider the initial condition for a quantity $v$ given by the mean value

$$v_i^0=\frac{1}{\mathrm{\Delta }{\xi }_i}{\int }_{{\xi }_{i-1/2}}^{{\xi }_{i+1/2}}v\left(\xi ,0\right)\mathit{d\xi }.$$E93

The eigenvalues satisfy ${\sigma }_1<{\sigma }_2=0<{\sigma }_3.$ Setting $u_{i+1/2}$ and $p_{i+1/2}$ as the values of $u$ and $p$ at the contact discontinuity between $V_i^n$ and $V_{i+1}^n$, it follows that

$$F\left(\left(w_R\left(0;V_i^n,V_{i+1}^n\right)\right)\right)={\left(-u_{i+j/2}^n,p_{i+j/2}^n,{\left(\mathit{pu}\right)}_{i+j/2}^n\right)}^T.$$E94

Then Godunov scheme for the Lagrangian coordinates takes the form

$$\left\{\begin{array}{l}{\tau }_i^{n+1}={\tau }_i^n+\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{m}_i}\left(u_{i+1/2}^n-u_{i-1/2}^n\right)\\ u_i^{n+1}=u_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{m}_i}\left(p_{i+1/2}^n-p_{i-1/2}^n\right)\\ e_i^{n+1}=e_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{m}_i}\left({\left(pu\right)}_{i+1/2}^n-{\left(pu\right)}_{i-1/2}^n\right)\end{array}\right.$$E95

where

$$\mathrm{\Delta }{m}_i={\rho }_i^0\mathrm{\Delta }{\xi }_i\mathrm{,}{p}_i^n=p\left({\tau }_i^n,{\epsilon }_i^n\right)\mathrm{,}{\epsilon }_i^n=e_i^n-\frac{{\left(u_i^n\right)}^2}{2}.$$E96

If we now consider the moving coordinates, Godunov scheme can also be derived equivalently by the following. Setting $x_{i+1/2}={\xi }_{i+1/2}$ with the approximation of $u=\mathit{dx}/\mathit{dt}$, it follows that the Eulerian coordinate $x_{i+1/2}$ of the interface ${\xi }_{i+1/2}$ at $t_n$ is upgraded with respect to

$$x_{i+1/2}^{n+1}=x_{i+1/2}^n+\mathrm{\Delta }{\mathit{tu}}_{i+1/2}^n.$$E97

Next we deduce

$${\rho }_i^n\left(x_{i+1/2}^n-x_{i-1/2}^n\right)=\mathrm{\Delta }{m}_i$$E98

by a simple induction process. Hence the Lagrangian Godunov schemes become

$$\left\{\begin{array}{c}\mathrm{\Delta }{m}_i={\rho }_i^0\left(x_{i+1/2}^0-x_{i-1/2}^0\right)\\ x_{i+1/2}^{n+1}=x_{i+1/2}^n+\mathrm{\Delta }{\mathit{tu}}_{i+1/2}^n\end{array}\right.$$E99

with

$$\left\{\begin{array}{l}{\rho }_i^{n+1}={\left(x_{i+1/2}^{n+1}-x_{i-1/2}^{n+1}\right)}^{-1}\mathrm{\Delta }{m}_i\\ u_i^{n+1}=u_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{m}_i}\left(p_{i+1/2}^n-p_{i-1/2}^n\right).\\ e_i^{n+1}=e_i^n-\frac{\mathrm{\Delta }t}{\mathrm{\Delta }{m}_i}\left({\left(pu\right)}_{i+1/2}^n-{\left(pu\right)}_{i-1/2}^n\right)\end{array}\right)$$E100

Notice that the Lagrangian Godunov schemes can be reformulated as a finite volume method. Equation (100) can be written in conservative form:

$${\partial }_t\left(\mathit{\phi J}\right)+{\partial }_{\xi }f=0.$$E101

If we integrate these equations on $\left({\xi }_{i-1/2},{\xi }_{i-1/2}\right)$ it follows that

$$\frac{d}{\mathit{dt}}{\int }_{x_{i-1/2}}^{x_{i+1/2}}\mathit{\phi d\xi }+\left(f_{i+1/2}-f_{i-1/2}\right)=0.$$E102

Here we omit the dependency of $f\mathrm{,}\phi$ and $x$ on $t$. Moreover, if we suppose that $\phi$ is constant in each cell $\left({\xi }_{i-1/2},{\xi }_{i-1/2}\right)$, it follows by an explicit one-step method that is

$$\mathrm{\Delta }{x}_i^{n+1}{\phi }_i^{n+1}=\mathrm{\Delta }{x}_i^n{\phi }_i^n-\mathrm{\Delta }t\left(f_{i+1/2}^n-f_{i-1/2}^n\right).$$E103

Moreover, if $\left(\rho ,u,e\right)$ are constant in each cell with $v=u$, we get the Godunov scheme:

$$\left\{\begin{array}{l}\mathrm{\Delta }{x}_i^n{\rho }_i^n=\mathrm{\Delta }{m}_i\\ \mathrm{\Delta }{m}_i^{n+1}{u}_i^{n+1}=\mathrm{\Delta }{m}_i^n{u}_i^n-\mathrm{\Delta }t\left(p_{i+1/2}^n-p_{i-1/2}^n\right)\\ \mathrm{\Delta }{m}_i^{n+1}{e}_i^{n+1}=\mathrm{\Delta }{m}_i^n{e}_i^n-\mathrm{\Delta }t\left({\left(pu\right)}_{i+1/2}^n-{\left(pu\right)}_{i-1/2}^n\right)\end{array}\right.$$E104

provided $\left(u_{i+1/2}^n,p_{i+1/2}^n\right)$ are determined by the solution of the Riemann problem, which is the desired result.

3.4 Comparison of Eulerian and Lagrangian schemes

In the literature there are two types of Godunov schemes: the Eulerian and Lagrangian. To compare one with the other, both have advantages and disadvantages. These are briefly listed in the following:

3.4.1 Eulerian approach

It is more nature; that is the properties of a flow field are described as functions of the coordinates which are in the natural physical space and time. The flow is determined by examining the behavior of the functions. Eulerian coordinates correspond to the coordinates of a fixed observer. This approach is ease of implementation and computation. The computational grids derived from the geometry constraints are generated in advance. The computational cells are fixed in space, and the fluid particles move across the cell interfaces. Since the Eulerian schemes consider the implementation at the nodes of a fixed grid, this may lead to spurious oscillations for the problems like diffusion-dominated transport equations. By adding artificial diffusion, one can get rid of these oscillations; however the nature of the problem may differ from the original one. Besides, refining the grids may also lead to remove numerical oscillations, but this process may augment the computation cost. Besides, while refining the grids, it may cause restriction of the size of time step which is limited by CFL condition. This restriction does not occur in Lagrangian case.

3.4.2 Lagrangian approach

It is based on the notion of mass coordinate denoted by $m\left(\xi \right)$. An important feature of the mass coordinate is that two segments have the same length if the mass contained in these segments is the same. This leads to face with a disadvantage; that is, at each iteration time step, the problem has to be converted from the natural coordinate system to the mass coordinate system. Once the solution at the next step is known, it has to be remapped into the natural coordinate system. As a result, this process raises the cost of the computation. Lagrangian coordinates are associated to the local flow velocity. In other words, as the velocity has different values in different parts of the fluid, then the change of coordinates is different from one point to another one in Lagrangian coordinates. Thus Lagrangian coordinates are equivalent to the Eulerian coordinates at another time. Lagrangian description states the motions and properties of the given fluid particles as they travel to different locations. Hence the computational grid points are precisely fluid particles. Since the particle paths in steady flow coincide with the streamlines, no fluid particles will cross the streamlines. Hence, there is no convective flux across cell boundaries, and the numerical diffusion is minimized. As a result, Godunov method in a Lagrangian grid is easier to handle. Moreover, in the case of higher schemes, the subsonic character of the flow makes the transformation much easier than in Eulerian schemes. Lagrangian schemes consider the implementation in a grid that moves with the flow which is an advantage for the problems like the transport equations since the advective and diffusion terms can separately be examined.

Apart from the two main approaches, there is another method which is a combination of both, so-called Eulerian-Lagrangian methods. It combines the advantages and eliminates disadvantages of both approaches to get a more efficient method. For further details we address the reader to the reference in the next part.

Notes

We have tried to present only the theoretical aspects of scalar conservation laws with some basic models and provide some examples of computational methods for the scalar models. There are plenty of contributors to the subject; however, we just cite some important of these and the references therein. Scalar conservation laws are thoroughly studied in particular in [12]; for a more general introduction including systems, see [13, 15, 18, 19, 22] and the references therein. There are some important works related to the concept of entropy provided by [7, 15, 16]. A more precise study of the shock and rarefaction waves can be found in [23]. A simple analysis for inviscid Burgers’ equation is done by [21]. The readers who are deeply interested in systems of conservation laws and the Riemann problem should see [8, 13, 15, 22, 24]. A well-ordered work of the propagation and the interaction of nonlinear waves are provided by [26]. We refer the reader to the papers [1, 17] for the theory of hyperbolic conservation laws on spacetime geometries and finite volume analysis with different aspects. A widely introductory material for finite difference and finite volume schemes to scalar conservation laws can be found in [18]. In this chapter we have studied the one-dimensional gas dynamics on the Eulerian and Lagrangian coordinates. For the detail on the Lagrangian conservation laws, we refer [10]; moreover for both Eulerian and Lagrangian conservation laws, we cite [11]. The proof of the equivalency of the Euler and Lagrangian equations for weak solutions is given in [25]. There are several numerical works for Lagrangian approach; some of the basic works on Lagrangian schemes are given in [2, 3, 4, 5, 6]. We refer the reader to the book [7] for a detailed analysis of the mathematical standpoint of compressible flows. Moreover Godunov-type schemes are precisely analyzed in [14, 27]; whereas, Lagrangian Godunov schemes can be found in [2, 12, 20]. As a last word, we must cite [9] as a recent and more general book consisting of scalar and system approaches of both Eulerian and Lagrangian conservation laws with theoretical and numerical parts which can be a basic source for the curious readers.