Riemann solvers for hyperbolic systems

Let $m\ge 1$, let $G\subseteq {\mathbb{R}}^m$ be open, and let $f:G\to {\mathbb{R}}^m$ be $C^3$ (Continuously differentiable maps, local inverses, and local diffeomorphisms; the $C^3$ convention is that of $C^k$ Euclidean maps and diffeomorphisms). The one-dimensional conservation system with state domain $G$ is $u_t+{\mathrm{\partial }}_{x}f(u)=0,u(t,x)\in G,$ where $u$ is the conserved state and $f$ its flux. Its flux Jacobian is $A(u)=Df(u)$. The system is strictly hyperbolic on an open set $N\subseteq G$ if, for every $u\in N$, $A(u)$ has $m$ distinct real eigenvalues, written in increasing order ${\lambda }_1(u)<\cdots <{\lambda }_m(u)$. An $i$th characteristic family on $N$ is a $C^2$ right-eigenvector field $r_i:N\to {\mathbb{R}}^m\setminus \{0\}$ with $A(u)r_i(u)={\lambda }_i(u)r_i(u)$. It is genuinely nonlinear on $N$ if $D{\lambda }_i(u)r_i(u)$ is nowhere zero there; after shrinking $N$ its sign is constant. It is linearly degenerate on $N$ if $D{\lambda }_i(u)r_i(u)=0$ throughout $N$. Rescaling $r_i$ by a nonzero constant changes the wave-strength parameter but not the family.

For $Q=(0,\mathrm{\infty })\times \mathbb{R}$, a locally bounded map $u:Q\to G$ is a distributional weak solution with initial data $u_0$ if ${\int }_Q(u\cdot {\phi }_t+f(u)\cdot {\phi }_x)dtdx=0$ for every vector test function $\phi \in C_c^{\mathrm{\infty }}(Q;{\mathbb{R}}^m)$ and $u(t,\cdot )\to u_0$ in $L_{\mathrm{l}\mathrm{o}\mathrm{c}}^1(\mathbb{R};{\mathbb{R}}^m)$ as $t\downarrow 0$. A Riemann initial datum is the step $u_0(x)=u_L$ for $x<0$ and $u_0(x)=u_R$ for $x>0$.

In this weak-solution clause, $u$ is required to be Lebesgue measurable and $f\circ u\in L_{\mathrm{l}\mathrm{o}\mathrm{c}}^1(Q;{\mathbb{R}}^m)$, so the test-function integral is defined. Local boundedness gives $u\in L_{\mathrm{l}\mathrm{o}\mathrm{c}}^1(Q;{\mathbb{R}}^m)$.

For a system as in Strictly hyperbolic conservation systems and characteristic families, a discontinuity with constant states $u_L,u_R\in G$ and speed $s\in \mathbb{R}$ is a Rankine–Hugoniot jump if $f(u_R)-f(u_L)=s(u_R-u_L)\mathrm{.}$ For a system strictly hyperbolic on $N\subseteq G$, states $u_L,u_R\in N$, and $1\le i\le m$, it is a strict Lax $i$-shock if it is a Rankine–Hugoniot jump and ${\lambda }_j(u_L)<s(j<i),{\lambda }_i(u_R)<s<{\lambda }_i(u_L),s<{\lambda }_j(u_R)(j>i)\mathrm{.}$ An empty family of inequalities is omitted. A linearly-degenerate contact has equal $i$-speeds on both sides, so it is not a strict Lax shock.

Let $f:G\to {\mathbb{R}}^m$ be $C^3$ and let $\bar{u}\in G$ be a state at which $A(\bar{u})=Df(\bar{u})$ has distinct real eigenvalues ${\lambda }_1^0<\cdots <{\lambda }_m^0$. For each $i$, there is a neighborhood $N_i$ of $\bar{u}$ and $C^2$ functions ${\lambda }_i:N_i\to \mathbb{R}$ and $r_i:N_i\to {\mathbb{R}}^m$ such that $A(u)r_i(u)={\lambda }_i(u)r_i(u)$, ${\lambda }_i(\bar{u})={\lambda }_i^0$, and the eigenvectors can be normalized by a fixed linear functional. After intersecting these finitely many neighborhoods, the eigenvalues remain real, simple, and ordered, and $r_1(u),\dots ,r_m(u)$ form a basis at each $u$.

Fix $\gamma >1$. In one space dimension the conserved state is $U=(\rho ,m,E)\in {\mathbb{R}}^3$, where $\rho$ is density, $m$ momentum, and $E$ total energy per unit volume. The admissible-state set is the open set ${\mathcal{G}}_{\gamma }=\{(\rho ,m,E):\rho >0,E>m^2\mathrm{/}(2\rho )\}\mathrm{.}$ The ideal-gas pressure on this set is $p(U)=(\gamma -1)(E-\frac{m^2}{2\rho })>0,$ the velocity is $u=m\mathrm{/}\rho$, and the flux is $F(U)=(m,\frac{m^2}{\rho }+p,\frac{m}{\rho }(E+p))\mathrm{.}$ The sound speed is $c=\sqrt{\gamma p\mathrm{/}\rho }$. These definitions apply only on ${\mathcal{G}}_{\gamma }$; a state with nonpositive density or nonpositive internal energy is inadmissible. The set is open because $\rho$ and $E-m^2\mathrm{/}(2\rho )$ are continuous on $\rho >0$ and the inequalities are strict; all displayed flux functions are smooth on this domain.

For states $u_L,u_R$ in the open state domain $G$ of Strictly hyperbolic conservation systems and characteristic families, a Roe matrix is a map $\tilde{A}:G\times G\to {\mathbb{R}}^{m\times m}$ such that, for all $u,u_L,u_R\in G$, $\tilde{A}(u,u)=Df(u),\tilde{A}(u_L,u_R)(u_R-u_L)=f(u_R)-f(u_L),$ and $\tilde{A}(u_L,u_R)$ is diagonalizable over $\mathbb{R}$ with real eigenvalues. The map $\tilde{A}$ is also continuous at each diagonal pair $(u,u)$. If $\tilde{A}=R\mathrm{\Lambda }{R}^{-1}$, define $\mathrm{\mid }\tilde{A}\mathrm{\mid }=R\mathrm{\mid }\mathrm{\Lambda }\mathrm{\mid }{R}^{-1}$ by replacing each diagonal eigenvalue with its absolute value. The Roe flux associated with this matrix is ${\widehat{f}}_{\mathrm{R}\mathrm{o}\mathrm{e}}(u_L,u_R)=\frac{1}{2}(f(u_L)+f(u_R))-\frac{1}{2}\mathrm{\mid }\tilde{A}(u_L,u_R)\mathrm{\mid }(u_R-u_L)\mathrm{.}$ The jump identity is an algebraic conservation property. This definition does not assert entropy admissibility of the approximate wave pattern.

For a two-state problem for a flux $f$ on the open state domain $G$ of Strictly hyperbolic conservation systems and characteristic families, choose explicit speeds $S_L<S_R$. A signal speed means the propagation speed of one wave in the particular Riemann solution being approximated. If the HLL model is used as an enclosing approximation to that fan, its admissibility hypothesis is that every such speed lies in $[S_L,S_R]$; this definition supplies no general wave-speed bound. Set $f_L=f(u_L)$ and $f_R=f(u_R)$. The HLL flux is ${\widehat{f}}_{\mathrm{H}\mathrm{L}\mathrm{L}}=\{\begin{array}{ll}{f}_L,& 0\le S_L,\\ {\displaystyle \frac{S_R{f}_L-S_L{f}_R+S_L{S}_R(u_R-u_L)}{S_R-S_L}},& S_L<0<S_R,\\ f_R,& S_R\le 0.\end{array}$ In the middle case the model has three constant regions separated by the two waves $S_L,S_R$ and one intermediate state. The wave-speed enclosure is a separate hypothesis whenever an approximation claim uses it.

For the scalar case of Strictly hyperbolic conservation systems and characteristic families, use its definition of a weak solution on $Q=(0,\mathrm{\infty })\times \mathbb{R}$. An entropy pair on $G$ is a function $\eta \in C^2(G)$ that is convex on each interval component of $G$, together with $q\in C^1(G)$ satisfying $q^{\mathrm{\prime }}(u)={\eta }^{\mathrm{\prime }}(u)f^{\mathrm{\prime }}(u)$. A weak solution $u$ is entropy-admissible for this pair when $\eta \circ u,q\circ u\in L_{\mathrm{l}\mathrm{o}\mathrm{c}}^1(Q)$ and ${\mathrm{\partial }}_t\eta (u)+{\mathrm{\partial }}_{x}q(u)\le 0$ in distributions; equivalently, ${\int }_Q(\eta (u){\phi }_t+q(u){\phi }_x)dtdx\ge 0$ for every nonnegative $\phi \in C_c^{\mathrm{\infty }}(Q)$. For a single jump of speed $s$, the entropy inequality requires $q(u_R)-q(u_L)-s(\eta (u_R)-\eta (u_L))\le 0$.

On a neighborhood where an eigenvector field $r_i$ is $C^2$, the local $i$-rarefaction curve through $u_0$ is the unique local solution $\frac{d{R}_i}{d\sigma }(\sigma ;u_0)=r_i(R_i(\sigma ;u_0)),R_i(0;u_0)=u_0\mathrm{.}$ The parameter $\sigma$ is the signed integral-curve parameter for the chosen normalization of $r_i$. Existence on a common small interval for nearby $u_0$ follows from Nearby initial values share one Picard–Lindelöf time interval and one state cylinder; the $C^2$ flow on a compact local tube follows from A $C^2$ vector field has a $C^2$ flow on a common compact tube. No global integral curve is asserted.

For a conservation system as in Strictly hyperbolic conservation systems and characteristic families, constant states $u_L,u_R\in G$, and a constant speed $s$, define $u(t,x)=u_L$ for $x<st$ and $u(t,x)=u_R$ for $x>st$, and let $u_0(x)=u_L$ for $x<0$ and $u_0(x)=u_R$ for $x>0$ (values on the jump lines are immaterial). Then $u$ is a distributional weak solution with initial datum $u_0$ in the sense of Strictly hyperbolic conservation systems and characteristic families if and only if $f(u_R)-f(u_L)=s(u_R-u_L)\mathrm{.}$

Let $G\subseteq {\mathbb{R}}^m$ be open, let $f:G\to {\mathbb{R}}^m$ be $C^3$, and let $\bar{u}\in G$ be a state where $A(\bar{u})=Df(\bar{u})$ has $m$ distinct real eigenvalues ${\lambda }_1(\bar{u})<\cdots <{\lambda }_m(\bar{u})$. Fix $i\in \{1,\dots ,m\}$ and its $i$th eigenpair $({\lambda }_i(\bar{u}),r_i(\bar{u}))$. There are a neighborhood $N_0$ of $\bar{u}$, $\eta >0$, and $C^1$ functions $v_i(\sigma ,u_0),s_i(\sigma ,u_0)$ for $u_0\in N_0$, $\mathrm{\mid }\sigma \mathrm{\mid }<\eta$, such that, with ${\mathrm{\ell }}_i(u_0)$ the dual left eigenvector normalized by ${\mathrm{\ell }}_i(u_0)\cdot r_i(u_0)=1$, ${\mathrm{\ell }}_i(u_0)\cdot v_i(\sigma ,u_0)=1,S_i(\sigma ;u_0)=u_0+\sigma v_i(\sigma ,u_0),$ and for $\sigma \mathrm{\ne }0$ the pair $(u_0,S_i(\sigma ;u_0))$ satisfies Rankine–Hugoniot with speed $s_i(\sigma ,u_0)$. Also $S_i(0;u_0)=u_0$ and ${\mathrm{\partial }}_{\sigma }{S}_i(0;u_0)=r_i(u_0)$. If the $i$th field is genuinely nonlinear on $N_0$ and $d_i(u):=D{\lambda }_i(u)r_i(u)$ has constant nonzero sign, then all sufficiently small nonzero strengths with $\sigma d_i(u_0)<0$ give strict Lax $i$-shocks.

For $U=(\rho ,m,E)\in {\mathcal{G}}_{\gamma }$, put $u=m\mathrm{/}\rho$, $H=(E+p)\mathrm{/}\rho$, and $c=\sqrt{\gamma p\mathrm{/}\rho }$. The flux Jacobian has the three simple real eigenvalues and right eigenvectors ${\lambda }_{-}=u-c,r_{-}=(1,u-c,H-uc{)}^T,$ ${\lambda }_0=u,r_0=(1,u,u^2\mathrm{/}2{)}^T,$ ${\lambda }_{+}=u+c,r_{+}=(1,u+c,H+uc{)}^T\mathrm{.}$ The middle field is linearly degenerate; the acoustic fields are genuinely nonlinear, with $D{\lambda }_{-}{r}_{-}=-(\gamma +1)c\mathrm{/}(2\rho )$ and $D{\lambda }_{+}{r}_{+}=(\gamma +1)c\mathrm{/}(2\rho )$.

Fix $\gamma >1$ and admissible left/right primitive states $({\rho }_K,u_K,p_K)$, $K\in \{L,R\}$, with ${\rho }_K,p_K>0$ and $c_K=\sqrt{\gamma p_K\mathrm{/}{\rho }_K}$. Set $A_K=\frac{2}{(\gamma +1){\rho }_K},B_K=\frac{\gamma -1}{\gamma +1}{p}_K\mathrm{.}$ For a trial star pressure $p>0$, the wave pressure response from side $K$ is $f_K(p)=\{\begin{array}{ll}{\displaystyle \frac{2c_K}{\gamma -1}}[{\left({\displaystyle \frac{p}{p_K}}\right)}^{(\gamma -1)\mathrm{/}(2\gamma )}-1],& 0<p\le p_K,\\ (p-p_K)\sqrt{{\displaystyle \frac{A_K}{p+B_K}}},& p\ge p_K\mathrm{.}\end{array}$ The two formulas agree at $p=p_K$. The star-pressure residual is $\mathrm{\Phi }(p)=f_L(p)+f_R(p)+u_R-u_L,p>0.$ This is a scalar root equation; it does not itself define a Riemann solution.

For any supplied Roe matrix, the Roe flux is consistent: ${\widehat{f}}_{\mathrm{R}\mathrm{o}\mathrm{e}}(u,u)=f(u)$. For a finite periodic array define $U_j^{new}=U_j-\nu ({\widehat{f}}_{j+1\mathrm{/}2}-{\widehat{f}}_{j-1\mathrm{/}2}),$ where ${\widehat{f}}_{j+1\mathrm{/}2}={\widehat{f}}_{\mathrm{R}\mathrm{o}\mathrm{e}}(U_j,U_{j+1})$ and indices are modulo $n$, with $n\ge 1$. For any scalar $\nu$, the update preserves the sum of cell states. The Roe jump property is the defining matrix identity in Roe linearization and Roe numerical flux and is not an entropy assertion.

When $S_L<0<S_R$, the HLL intermediate state $u_{\ast }=\frac{S_R{u}_R-S_L{u}_L-(f_R-f_L)}{S_R-S_L}$ together with the assigned model flux $f_{\ast }=f_L+S_L(u_{\ast }-u_L)$ satisfies the two model balances $f_{\ast }-f_L=S_L(u_{\ast }-u_L),f_R-f_{\ast }=S_R(u_R-u_{\ast })\mathrm{.}$ This $f_{\ast }$ is the middle formula in HLL two-wave approximate flux. In general $f_{\ast }\mathrm{\ne }f(u_{\ast })$, even when $u_{\ast }\in G$; the model is not asserted to be a weak solution of the original conservation law. For equal states the HLL flux is consistent. A finite-volume update using the same HLL flux at each periodic interface is conservative: for any finite periodic array with adjacent pairs in the flux domain, define $U_j^{new}=U_j-\nu ({\widehat{f}}_{\mathrm{H}\mathrm{L}\mathrm{L}}(U_j,U_{j+1})-{\widehat{f}}_{\mathrm{H}\mathrm{L}\mathrm{L}}(U_{j-1},U_j))\mathrm{.}$ Here $n\ge 1$. The result holds for any scalar $\nu$ and assumes the supplied speeds; it does not prove that they enclose a nonlinear Riemann fan.

For a side state $K=({\rho }_K,u_K,p_K)$ and any trial pressure $p>0$, the rarefaction/shock branches defined in Ideal-gas rarefaction and shock pressure response give the scalar velocity response $f_K(p)$ along the corresponding acoustic wave. For a left 1-wave the star velocity is $u_{\ast }=u_L-f_L(p)$; for a right 3-wave it is $u_{\ast }=u_R+f_R(p)$. The star density on a rarefaction branch is ${\rho }_{\ast }(p)={\rho }_K(p\mathrm{/}{p}_K{)}^{1\mathrm{/}\gamma },$ and on a compressive shock branch, with $q=p\mathrm{/}{p}_K>1$ and $\beta =(\gamma -1)\mathrm{/}(\gamma +1)$, it is ${\rho }_{\ast }(p)={\rho }_K\frac{q+\beta }{\beta q+1}\mathrm{.}$ The shock branches satisfy the strict acoustic Lax inequalities. When the two pressure responses match at a positive $p$, their star states are joined by a linearly-degenerate contact at speed $u_{\ast }$.