\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 32, pp. 1--18.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University - San Marcos.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2015/32\hfil Existence of solutions to Riemann problems] {Existence of solutions to the Riemann problem for a model of two-phase flows} \author[M. D. Thanh, D. H. Cuong \hfil EJDE-2015/32\hfilneg] {Mai Duc Thanh, Dao Huy Cuong} \address{Mai Duc Thanh \newline Department of Mathematics, International University, Vietnam National University - HCMC, Quarter 6, Linh Trung Ward, Thu Duc District, Ho Chi Minh City, Vietnam} \email{mdthanh@hcmiu.edu.vn} \address{Dao Huy Cuong \newline Nguyen Huu Cau High School, 07 Nguyen Anh Thu, Trung Chanh Ward, Hoc Mon District, Ho Chi Minh City, Vietnam. \newline Department of Mathematics and Computer Science, University of Science, Vietnam National University - HCMC, 227 Nguyen Van Cu str., District 5, Ho Chi Minh City, Vietnam} \email{cuongnhc82@gmail.com} \thanks{Submitted November 19, 2014. Published February 5, 2015.} \subjclass[2000]{35L65, 35L67, 76T10, 76N10} \keywords{Two-phase flow; nonconservative; source term; jump relation; \hfill\break\indent shock; Riemann problem} \begin{abstract} We study the existence of solutions of the Riemann problem for a model of two-phase flows. The model has the form of a nonconservative hyperbolic system of balance laws. Based on a phase decomposition approach, we obtain all the wave curves. By developing an analytic method, we can establish a system of nonlinear algebraic equations for each solution of the Riemann problem. The system is under-determined and can be parameterized by the volume fraction in one phase. Therefore, an argument relying on the Implicit-Function Theorem leads us to the existence of solutions of the Riemann problem for the model for sufficiently large initial data. Furthermore, the structure of the Riemann solutions obtained by this method can also be obtained. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \allowdisplaybreaks \section{Introduction} \label{IN} In this article we consider the existence and the structure of solutions of the following Riemann problem for the model of two-phase flows, \begin{equation} \begin{gathered} \partial_t(\alpha\rho) + \partial_x(\alpha\rho u) = 0,\\ \partial_t(\alpha\rho u) + \partial_x(\alpha(\rho u^2 + p))= p \partial_x \alpha,\\ \partial_t(\beta\theta) + \partial_x(\beta\theta v) = 0,\\ \partial_t(\beta\theta v) + \partial_x(\beta(\theta v^2 + q)) = -p \partial_x \alpha,\\ \partial_t\theta + \partial_x(\theta v) =0,\quad x\in\mathbb{R}, t>0,\\ \end{gathered}\label{1.1} \end{equation} where $\alpha,\rho, u, p$ stand for the volume fraction, density, velocity, and pressure in the first phase of the flow, called \emph{Phase I}; and $\beta,\theta, v, q$ stand for the volume fraction, density, velocity, and pressure in the second phase of the flow, called \emph{Phase II}. The volume fractions satisfy \begin{equation} \alpha+\beta=1.\label{1.2} \end{equation} System \eqref{1.1} is obtained from the full model of two-phase flows, see \cite{BaerNunziato, BzilMenikoffSonKapilaSteward}, by assuming that the flow is isentropic in both phases. The first and the second equations of \eqref{1.1} describe the balance of mass and momentum in Phase I, while the third and the four equations of \eqref{1.1} describe the balance of mass and momentum in Phase; the fifth equation of \eqref{1.1}, called the \emph{compaction dynamics equation}, represents the evolution of the volume fractions. Observe that the third and the fifth equations of \eqref{1.1} yield \begin{equation} \partial_t\beta + v\partial_x\beta=0. \label{1.3} \end{equation} Note that the first version of the compaction dynamics equation \eqref{1.3} is used in \cite{ BaerNunziato}. The current compaction dynamics equation ( the last equation of \eqref{1.1}) was proposed in \cite{BzilMenikoffSonKapilaSteward}, and is suitable for our analysis as it is \emph{conservative}. System \eqref{1.1} has the form of nonconservative systems of balance laws, where the weak solutions will be understood in the sense of \emph{nonconservative products}, see \cite{DalMasoLeFlochMurat}. It has been known that the system \eqref{1.1} is hyperbolic, but not strictly hyperbolic, since the characteristic fields may coincide. Moreover, the fact that the dimension of the unknown function $U$ is large (five components) makes it hard to solve the Riemann problem for large initial data. However, the system \eqref{1.1} possesses a very interesting property: four characteristic fields involve quantities only in one phase. This allows the waves associated with these characteristic fields to change only in one phase and remain constants in the other. In an earlier work, \cite{ThanhJMAA}, a phase decomposition approach for studying the Riemann problem for \eqref{1.1} was proposed. Then, based on a geometrical approach where the positions of related various curves in a phase plane can be determined and compared, the author established several results on the existence of solutions of the Riemann problem for \eqref{1.1}. In this work, we will rely on an analytical method, instead, to establish existence results of Riemann solutions of \eqref{1.1}. For this aim, we will show that each Riemann solution, whose structure can be determined using the phase decomposition method, corresponds to a set of nonlinear algebraic equations. Then, we will show that these equations can be parameterized as functions of the volume fraction in one phase, says, Phase I. By applying an Implicit Function Theorem, we can obtain the existence of Riemann solutions. It is worth to note that the locality in applying Implicit Function Theorem here is constraint to only the gas volume fractions, while all other quantities of the Riemann data (density, pressure, temperature, velocity) can be taken relatively large. Moreover, we note that since each solution of the Riemann problem is corresponding to a system of nonlinear algebraic equations, computational methods can be developed for computing the exact solutions. Consequently, this work can be useful for developing numerical methods such as the Godunov method to approximate the initial-value problem for the model under study. Hyperbolic models in nonconservative forms have attracted the attention of many authors. The earlier works concerning nonconservative systems were carried out in \cite{IsaacsonTemple92, LeFloch88,LeFloch89, MarchesinPaes-Leme}. The Riemann problem for the model of a fluid in a nozzle with discontinuous cross-section was considered in \cite{LeFlochThanh03.2} for the isentropic case, and in \cite{Thanh09} for the non-isentropic case. The Riemann problem for for shallow water equations with discontinuous topography were solved in \cite{LeFlochThanh07, LeFlochThanhJCP}. The Riemann problem for a general system in nonconservative form was studied by \cite{GoatinLeFloch}. The Riemann problem for a model of two-phase flows was studied in \cite{ SchwendemanetalJCP, ThanhJMAA}. Two-fluid models of two-phase flows were studied in \cite{KeyfitzSandersSever, ThanhNARWA12}. Numerical approximations for two-phase flows were considered in \cite{ CoqueletalM2AN2009,Bouchutbook, KarniDuenas,ThanhKroenerNam, ThanhKroenerChalons,ThanhCNSNS14, ThanhKJM}. See also the references therein. The organization of this article is as follows. In Section 2 we recall basic properties of the system \eqref{1.1}, and the jump relations by using the phase decomposition approach. In Section 3 we establish the existence results of solutions of the Riemann problem. \section{Preliminaries} \subsection{Characteristic fields} Throughout, we assume for simplicity that the fluid in each phase is isentropic and ideal, where the equation of state is $$ p=p(\rho)=\kappa\rho^\gamma, \quad q=q(\theta)=l\theta^\delta,\quad \kappa,l>0,\; \gamma,\delta>1. $$ The system \eqref{1.1} can be re-written as \begin{equation} \begin{gathered} \partial_t\rho + u\partial_x\rho+\rho\partial_xu +\frac{\rho(u-v)}{\alpha}\partial_x\alpha = 0, \\ \partial_tu+h'(\rho)\partial_x\rho+u\partial_xu=0,\\ \partial_t\theta+v\partial_x\theta+\theta\partial_xv=0,\\ \partial_tv+k'(\theta)\partial_x\theta+v\partial_xv +\frac{p-q}{\beta\theta}\partial_x\alpha=0,\\ \partial_t\alpha + v\partial_x\alpha =0, \quad x\in\mathbb{R}, \;t>0,\\ \end{gathered} \label{2.1} \end{equation} where $$ h'(\rho)=\frac{p'(\rho)}{\rho}, \quad k'(\theta)=\frac{q'(\theta)}{\theta}. $$ Thus, if we choose the unknown vector function to be of the form $$ U=(\rho,u,\theta,v,\alpha), $$ we can re-write the system \eqref{1.1} as a system of balance laws in nonconservative form as \begin{equation} U_t+A(U)U_x=0, \label{2.2} \end{equation} where $$ A(U)=\begin{pmatrix} u&\rho&0 &0&\frac{\rho(u-v)}{\alpha}\\ h'(\rho)&u&0&0&0\\ 0&0&v&\theta&0\\ 0&0&k'(\theta)&v&\frac{p-q}{\beta\theta}\\ 0&0&0&0&v\\ \end{pmatrix}. $$ The characteristic equation of the matrix $A(U)$ is $$ (v-\lambda)((u-\lambda)^2-p')((v-\lambda)^2-q')=0, $$ which admits five roots as \begin{equation} \begin{gathered} \lambda_1(\rho,u)=u-\sqrt{p'(\rho)},\quad \lambda_2(\rho,u)=u+\sqrt{p'(\rho)},\\ \lambda_3(\theta,v)=v-\sqrt{q'(\theta)},\quad \lambda_4(\theta,v)=v+\sqrt{q'(\theta)},\quad \lambda_5(v)=v. \end{gathered}\label{2.3} \end{equation} The corresponding right eigenvectors can be chosen as \begin{equation} \begin{gathered} r_1(\rho,u)=\mu\begin{pmatrix}-\rho\\ \sqrt{p'(\rho)}\\ 0\\ 0\\ 0 \end{pmatrix}, \quad r_2(\rho,u)=\mu \begin{pmatrix}\rho\\ \sqrt{p'(\rho)}\\ 0\\ 0\\ 0 \end{pmatrix},\\ r_3(\theta,v)=\nu \begin{pmatrix}0\\ 0\\ -\theta\\ \sqrt{q'(\theta)}\\ 0 \end{pmatrix}, \quad r_4(\theta,v)=\nu \begin{pmatrix}0\\ 0\\ \theta\\ \sqrt{q'(\theta)}\\ 0\\ \end{pmatrix},\\ r_5(U)=\begin{pmatrix} -(u-v)^2\rho\beta q(\theta)\\ (u-v)p'(\rho)q'(\theta)\beta\\ (q(\theta)-p(\rho))((u-v)^2-p'(\rho))\alpha\\ 0\\ ((u-v)^2-p'(\rho))\alpha\beta q'(\theta)\\ \end{pmatrix},\\ \end{gathered}\label{2.4} \end{equation} where $$ \mu=\frac{2\sqrt{p'(\rho)}}{p''(\rho)\rho + 2p'(\rho)},\quad \nu=\frac{2\sqrt{q'(\theta)}}{q''(\theta)\theta + 2q'(\theta)}. $$ It is not difficult to check that the eigenvectors $r_i, i=1,2,3,4,5$ are linearly independent. Thus, the system is hyperbolic. Furthermore, it holds that $$ \lambda_3<\lambda_5<\lambda_4. $$ It is interesting that the eigenvalues $\lambda_5$ may coincide with either $\lambda_1$ or $\lambda_2$ on a certain hyper-surface of the phase domain, called the \emph{sonic surface} or \emph{resonant surface}. We call the \emph{supersonic region} to be the one in which \begin{equation} |u-v|>c:=\sqrt{p'(\rho)}, \label{2.5} \end{equation} the \emph{subsonic region} is the one in which $|u-v|\sqrt{p'(\rho)}\},\\ G_2=\{(\rho,u)| \quad |u-v_0| < \sqrt{p'(\rho)}\},\\ G_3=\{(\rho,u)| \quad u-v_0< -\sqrt{p'(\rho)}\},\\ \mathcal{C}_\pm = \{(\rho,u)| \quad u-v_0=\pm \sqrt{p'(\rho)}\},\\ \mathcal{C}=\mathcal{C}_+\cup\mathcal{C}_-. \end{gathered}\label{2.6} \end{equation} \begin{figure}[htb] \begin{center} \includegraphics[width=0.5\textwidth]{fig1} % rhou-plan.png \end{center} \caption{Projection of the phase domain hyper-plane $v\equiv v_0$}\label{rhouplan} \end{figure} Then, the states $U$ such that $(\rho,u)\in G_1, G_3$, belong to the supersonic region; the states $U$ such that $(\rho,u)\in G_2$, belong to the subsonic region; and the states $U$ such that $(\rho,u)\mathcal{C}_\pm$ belong to the sonic surface. On the other hand, it is not difficult to verify that \begin{equation} \begin{gathered} D\lambda_i(U)\cdot r_i(U)=1,\quad i=1,2,3,4,\\ D\lambda_5(U)\cdot r_5(U)=0, \end{gathered}\label{2.7} \end{equation} so that the first, second, third, fourth characteristic fields $(\lambda_i(U),r_i(U))$, $i=1,2,3,4$, are genuinely nonlinear, while the fifth characteristic field $(\lambda_5(U),r_5(U))$ is linearly degenerate. \subsection{Rarefaction waves} Rarefaction waves of the system \eqref{2.2}, and therefore of \eqref{1.1}, are the continuous piecewise-smooth self-similar solutions of \eqref{1.1} associated with nonlinear characteristic fields, which have the form $$ U(x,t)=V(\xi),\quad \xi=\frac{x}{t},\quad t>0, x\in\mathbb{R}. $$ Substituting this into \eqref{2.2}, we can see that rarefaction waves are solutions of the following initial-value problem for ordinary differential equations \begin{equation} \begin{gathered} \frac{d V(\xi)}{d \xi}= r_i(V(\xi)),\quad \xi\ge \lambda_i(U_0),\quad i=1,2,3,4,\\ V(\lambda_i(U_0))=U_0. \end{gathered} \label{2.8} \end{equation} Thus, the integral curve of the first characteristic field is given by \begin{equation} \begin{gathered} \frac{d \rho(\xi)}{d \xi}= \frac{-2\sqrt{p'(\rho)}}{ p''(\rho)\rho + 2p'(\rho)}\rho(\xi)<0,\\ \frac{d u(\xi)}{d \xi}= \frac{2\sqrt{p'(\rho)}}{p''(\rho)\rho + 2p'(\rho)}\sqrt{p'(\xi)}>0,\\ \frac{d \theta(\xi)}{d \xi}= \frac{d v(\xi)}{d \xi} = \frac{d \alpha(\xi)}{d \xi}= 0. \end{gathered} \label{2.9} \end{equation} This implies that \emph{$\theta, v,\alpha$ are constant through $1$-rarefaction waves}, $\rho$ is strictly decreasing with respect to $\xi$, and $u$ is strictly increasing with respect to $\xi$. Moreover, since $\rho$ is strictly monotone though 1-rarefaction waves, we can use $\rho$ as a parameter of the integral curve \begin{equation} \frac{d u}{d\rho}= \frac{-\sqrt{p'(\rho)}}{\rho}. \label{2.10} \end{equation} The integral curve \eqref{2.10} determines the \emph{forward} curve of 1-rarefaction wave $\mathcal{R}_1(U_0)$ consisting of all right-hand states that can be connected to the left-hand state $U_0$ using $1$-rarefaction waves \begin{equation} \mathcal{R}_1(U_0):\quad u=\omega_1((\rho_0,u_0);\rho) := u_0-\int_{\rho_0}^{\rho} \frac{\sqrt{p'(y)}}{ y} dy,\quad \rho\le \rho_0, \label{2.11} \end{equation} where $\rho\le \rho_0$ follows from the condition that the characteristic speed must be increasing through a rarefaction fan. Similarly, $\theta, v,\alpha$ \emph{are constant through $2$-rarefaction waves}. The \emph{backward } curve of 2-rarefaction wave $\mathcal{R}_2(U_0)$ consisting of all left-hand states that can be connected to the right-hand state $U_0$ using $2$-rarefaction waves is given by \begin{equation} \mathcal{R}_2(U_0):\quad u=\omega_2((\rho_0,u_0);\rho) :=u_0+\int_{\rho_0}^{\rho}\frac{\sqrt{p'(y)}}{ y} dy, \quad\rho\le \rho_0. \label{2.12} \end{equation} In the same way, $\rho, u,\alpha$ \emph{are constant through $3$- and $4$-rarefaction waves}. The \emph{forward} curve of 3-rarefaction wave $\mathcal{R}_3(U_0)$ consisting of all right-hand states that can be connected to the left-hand state $U_0$ using $3$-rarefaction waves is given by \begin{equation} \mathcal{R}_3(U_0):\quad v=\omega_3((\theta_0,v_0);\theta) :=v_0-\int_{\theta_0}^{\theta} \frac{\sqrt{q'(y)}}{y} dy, \quad\theta\le \theta_0. \label{2.13} \end{equation} The \emph{backward } curve of 4-rarefaction wave $\mathcal{R}_4(U_0)$ consisting of all left-hand states that can be connected to the right-hand state $U_0$ using $4$-rarefaction waves is given by \begin{equation} \mathcal{R}_4(U_0):\quad v=\omega_4((\theta_0,v_0);\theta) :=v_0+\int_{\theta_0}^{\theta}\frac{\sqrt{q'(y)}}{y} dy, \quad \theta\le \theta_0. \label{2.14} \end{equation} \subsection{Jump relations for shock waves} A discontinuity (shock or contact wave) of \eqref{1.1} is a weak solution (in the sense of nonconservative products) and is of the usual form \begin{equation} U(x,t)=\begin{cases} U_-, &\text{for } x<\sigma t,\\ U_+, &\text{for } x>\sigma t, \end{cases}\label{2.15} \end{equation} for some constant states $U_\pm$ and a constant shock speed $\sigma$. This discontinuity satisfies the generalized Rankine-Hugoniot relations for a given family of Lipschitz paths. The generalized Rankine-Hugoniot relations corresponding to any family of Lipschitz path for the conservative equations in \eqref{1.1} must coincide with the usual canonical ones. In particular, it holds that \begin{equation} \begin{gathered} -\sigma[\beta\theta] + [\beta\theta v] = 0, \\ -\sigma[\theta] + [\theta v] =0, \end{gathered} \label{2.16} \end{equation} where $\sigma$ is the shock speed, $[A]=A_+-A_-$, and $A_\pm$ denote the values on the right and left of the jump on the quantity $A$. This yields \begin{equation} \begin{gathered} \theta(v-\sigma) =M=\text{constant},\\ M[\beta] = 0. \end{gathered}\label{2.17} \end{equation} The second equation of \eqref{2.17} implies that either $M=0$ or $[\beta]=0$. Since $\theta> 0$, one obtains the following conclusion: across any discontinuity \eqref{2.15} of \eqref{1.1} \begin{equation} \text{ either $[\beta]=0$, or $v=\sigma$, a constant}. \label{2.18} \end{equation} It is derived from \eqref{2.18} that if $[\beta]=0$, then the volume fractions remain constant across the discontinuity. The system \eqref{1.1} is therefore reduced to the two independent sets of isentropic gas dynamics equations in both phases \begin{equation} \begin{gathered} \partial_t\rho + \partial_x(\rho u) = 0,\\ \partial_t(\rho u) + \partial_x( \rho u^2 + p)= 0,\\ \partial_t\theta + \partial_x(\theta v) = 0, \\ \partial_t(\theta v) + \partial_x(\theta v^2 + q) = 0,\quad x\in\mathbb{R}, t>0.\\ \end{gathered} \label{2.19} \end{equation} This implies that $\theta, v,\alpha$ \emph{are constant through $1$- and $2$-shock waves, while $\rho, u,\alpha$ are constant through $3$- and $4$-shock waves}. Given a left-hand state $U_0$, let us denote by $\mathcal{S}_i(U_0), i=1,3$ the \emph{forward} shock curves consisting of all right-hand states $U$ that can be connected to the left-hand state $U_0$ by an $i$-Lax shock, $i=1,3$, and by $\mathcal{S}_j(U_0), j=2,4$ the \emph{backward} shock curves consisting of all left-hand states $U$ that can be connected to the right-hand state $U_0$ by a $j$-Lax shock, $j=2,4$. These curves are given by: \begin{equation} \begin{gathered} \mathcal{S}_{1}(U_0) : \quad u =\omega_1((\rho_0,u_0);\rho) := u_0 - \Big(\frac{(p-p_0)(\rho-\rho_0)}{\rho_0\rho}\Big)^{1/2}, \quad\rho > \rho_0,\\ \mathcal{S}_{2}(U_0) : \quad u=\omega_2((\rho_0,u_0);\rho) := u_0 + \Big(\frac{(p-p_0)(\rho-\rho_0)}{\rho_0\rho}\Big)^{1/2}, \quad\rho > \rho_0,\\ \mathcal{S}_3(U_0) : \quad v =\omega_3((\theta_0,v_0);\theta) :=v_0 - \Big(\frac{(q-q_0)(\theta-\theta_0)}{\theta_0\theta}\Big)^{1/2}, \quad\theta > \theta_0,\\ \mathcal{S}_4(U_0) : \quad v =\omega_4((\theta_0,v_0);\theta) :=v_0 + \Big(\frac{(q-q_0)(\theta-\theta_0)}{\theta_0\theta}\Big)^{1/2}, \quad\theta > \theta_0. \end{gathered} \label{2.20} \end{equation} From \eqref{2.11}--\eqref{2.14} and \eqref{2.20}, we can now define the \emph{forward wave curves } issuing from $U_0$ as \begin{equation} \begin{gathered} \mathcal{W}_1(U_0) =\mathcal{R}_1(U_0)\cup\mathcal{S}_{1}(U_0),\\ \mathcal{W}_3(U_0) =\mathcal{R}_3(U_0)\cup\mathcal{S}_{3}(U_0), \end{gathered} \label{2.21} \end{equation} and the \emph{backward wave curves} issuing from $U_0$ by \begin{equation} \begin{gathered} \mathcal{W}_2(U_0) =\mathcal{R}_2(U_0)\cup\mathcal{S}_{2}(U_0),\\ \mathcal{W}_4(U_0) =\mathcal{R}_4(U_0)\cup\mathcal{S}_{4}(U_0). \end{gathered} \label{2.22} \end{equation} These curves are parameterized in such a way that the velocity is given as a function of the density in each phase, under the form $u=\omega_i(U_0;\rho), \rho>0, i=1,2,3,4$. It is not difficult to check that $\omega_1,\omega_3$ are strictly decreasing; and that $\omega_2,\omega_4$ are strictly increasing. Summarizing the above argument, we get the following result. \begin{lemma}\label{lem31} Through an $i$-wave (shock or rarefaction), $i=1,2$, the quantities $\theta, v,\alpha$ are constant. Through a $j$-wave (shock or rarefaction), $j=3,4$, the quantities $\rho, u,\alpha$ are constant. The wave curves $\mathcal{W}_i(U_0), i=1,2,3,4$, associated with the genuinely nonlinear characteristic fields issuing from a given state $U_0$ are given by \eqref{2.21} and \eqref{2.22}. \end{lemma} \subsection{Jump relations for contact waves} Contact discontinuities correspond to the second of \eqref{2.18} where $[\beta]\ne 0$. Any contact discontinuity associated with the fifth characteristic field can be characterized as follows. \begin{theorem}[{\cite[Theorem 3.3]{ThanhJMAA}}] \label{theo33} Let $U$ be a contact discontinuity of the form \eqref{2.15} associated with the linearly degenerate characteristic field $(\lambda_5,r_5)$; that is, $[\beta]\ne 0$ and $U_\pm$ belong to the same trajectory of the integral field of the 5fth characteristic field. Then, $U$ is a weak solution of \eqref{1.1} in the sense of nonconservative products. Moreover, this contact discontinuity $U$ satisfies the jump relations in the usual form \begin{equation} \begin{gathered} v_{\pm} = \sigma,\quad [\alpha\rho(u-v)] = 0,\\ [(u-v)^2+2h]=0,\quad [mu + \alpha p+ \beta q] = 0, \end{gathered} \label{2.23} \end{equation} where $m$ is the constant $m=\alpha\rho(u-v)$. \end{theorem} Thus, whenever a state on one side of a contact discontinuity is fixed, the state on the other side $U$ that can be connected with $U_0$ by a contact discontinuity must satisfy the equations \begin{equation} \begin{gathered} \alpha\rho(u-v) =\alpha_{g0}\rho_0(u_0-v_{}):=m ,\\ (u-v)^2+2h=(u_0-v_0)^2+2h_{g0}, \end{gathered} \label{2.24} \end{equation} and \begin{equation} q =\frac{\beta_0q_0-[mu + \alpha p]}{\beta}. \label{2.25} \end{equation} So, we can define the \emph{fifth wave curve} $\mathcal{W}_5(U_0) $ to be the curve of contact waves issuing from $U_0$. It is the set of states $U$ that can be connected to $U_0$ by a contact discontinuity, where $U$ can be defined by the equations \eqref{2.24} and \eqref{2.25}. \section{Existence and structure of Riemann solutions} In this section, we will show that the Riemann problem for \eqref{1.1} possesses a solution made up a finite number of waves (shocks, rarefaction waves, and contacts) separated by constant states, only. The case where the solution containing coinciding waves is much more complicated and it will be the topic for future developments. \subsection*{Notation} In this section, we use the following notation: \begin{itemize} \item[(i)] $W_k(U_i,U_j)$ ($S_k(U_i,U_j), R_k(U_i,U_j)$) denotes a $k$-wave ($k$-shock, or $k$-rar\-efaction wave, respectively) connecting the left-hand state $U_i$ to the right-hand state $U_j$; \item[(ii)] $W_m(U_i,U_j)\oplus W_n(U_j,U_k)$ indicates that there is an $m$-wave from the left-hand state $U_i$ to the right-hand state $U_j$, followed by an $n$-wave from the left-hand state $U_j$ to the right-hand state $U_k$; \item[(iii)] $ U_\pm=(\rho_{\pm},u_\pm)$ and $ V_\pm=(\theta_{\pm},v_0)$ denote the left- and right-hand states of the contact in Phase I and Phase II of a Riemann solution under consideration, respectively. \end{itemize} \subsection{Phase decomposition} We can see from the above that \begin{itemize} \item Through $i$-waves, $i=1,2$, the quantities in the phase II: $V=(\theta,v)$ and $\alpha,\beta$ remain constant; \item Through $j$-waves, $j=3,4$, the quantities in the phase I: $U=(\rho,u)$ and $\alpha,\beta$ remain constant. \end{itemize} So, in any Riemann solution: \begin{itemize} \item[(i)] Quantities in the phase I involve only in the 1st, 2nd and 5fth characteristic fields; \item[(ii)] Quantities in the phase II involve only in the 3rd, 4th and 5fth characteristic fields. \end{itemize} Only 5-contacts involve quantities in both phases. So, they can serve as a ``bridge'' connecting the two phases. We can consider the wave structure in each phase separately. \subsection*{Components of Phase I in Riemann solutions} Quantities in Phase I involve only waves from the 1st, 2nd, and and 5fth fields. However, although $$ \lambda_1(U)<\lambda_2(U), $$ $\lambda_5(U)$ can be in any order with $\lambda_1(U)$ and $\lambda_2(U)$. In our phase decomposition, the contact waves play a key role for constructing Riemann solutions. Following \cite{LeFlochThanh03.2, Thanh09, IsaacsonTemple92, IsaacsonTemple95}, we also impose the following admissibility condition: \begin{itemize} \item[(AC)] The contact wave must remain in the same subsonic or supersonic region. \end{itemize} This criterion means that the left-hand and right-hand states of the admissible contact will remain in the same subsonic or supersonic region. Therefore, our construction of Riemann solutions will rely on the location of the contact in the supersonic or subsonic region as follows. Let $U_-$ be the state on the left, and $U_+$ be the state on the right of the contact. \begin{itemize} \item[(A)] The contact belongs to the supersonic region $\lambda_5<\lambda_1$; \item[(B)] The contact belongs to the subsonic region $\lambda_1<\lambda_5<\lambda_2$; \item[(C)] The contact belongs to the supersonic region $\lambda_5>\lambda_2$. \end{itemize} In the next subsection, we build three classes of Riemann solutions, having different configurations, for these three cases. \subsection*{Components of Phase II in Riemann solutions} In Phase II, one has $\lambda_3<\lambda_5<\lambda_4$. The Riemann solutions in Phase II always have the form $$ W_3(V_L,V-)\oplus W_5(V_-,V_+)\oplus W_4(V_+,V_R), $$ see Figure \ref{fig41}. \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig2a} % Sol1xtII.png \quad \includegraphics[width=0.45\textwidth]{fig2b} % Sol1II.png \end{center} \caption{Riemann solutions in the phase II} \label{fig41} \end{figure} So, the quantities in Phase II involve only four constant states $V_L, V_\pm, V_R$. Since $$ V_-=(\theta_{-},v_0)\in\mathcal{W}_3(V_L),\quad V_+=(\theta_{+},v_0)\in\mathcal{W}_4(V_R). $$ Therefore, \begin{equation} v_0=v_\pm=\omega_3(V_L;\theta_{-}) =\omega_4(V_R;\theta_{+}), \label{3.1} \end{equation} where $\omega_3$ and $\omega_4$ are defined by \eqref{2.13}, \eqref{2.14}, and \eqref{2.20}. \subsection{Construction of Solutions} \subsection*{Case A: The contact belongs to the supersonic region $\lambda_1>\lambda_5$} The solution contains a 5-contact on the left of both waves in nonlinear families, see Figure \ref{fig42} (left). As before, $U_\pm, V_\pm$ denote the states on the left and right of the contact in Phase I and Phase II, respectively. One has $$ U_-=U_{L}. $$ Let $\{U_1\}=\mathcal{W}_1(U_+)\cap\mathcal{W}_2(U_R)$. \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig3a} % Sol1xtI.png \quad \includegraphics[width=0.45\textwidth]{fig3b} % Sol1I.png \end{center} \caption{Case A: Riemann solutions in the phase I} \label{fig42} \end{figure} The solution begins with a liquid 3-wave from $V_L$ to $V_-$, followed by a 5-contact $V_+$, then it continues separately in each phase. In Phase I, the solution does not change across the 3-wave. It jumps by a 5-contact from $U_L$ to $U_+$, followed by a 1-wave from $U_+$ to $U_1$, then arrives at $U_R$ by a 2-wave, see Figure \ref{fig42} (right). In phase II, the solution arrives at $V_R$ from $V_+$ by a 4-wave. The whole Riemann solution has the configuration as in Figure \ref{fig43}, where the 1-, 2-, and 4-waves may interchange the order. \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig4} % Sol1xt.png \end{center} \caption{Case A: A whole Riemann solution where the contact is in the supersonic region $G_1$} \label{fig43} \end{figure} \begin{theorem}\label{theo41} Let $V_*=(\theta_*,v_*)$ be the intersection point of $\mathcal{W}_3(V_L)$ and $\mathcal{W}_4(V_R)$. Providing that the state $(\rho_L,u_L,v_*)$ is in the supersonic region $$ u_L-v_* > \sqrt{p'(\rho_L)}, $$ or, $\lambda_1(U_L)>\lambda_5(v_*)$, there exists an interval $I\ni \alpha_L$ such that whenever $\alpha_R\in I$, the states on both sides of the contact $U_\pm, V_\pm, U_-=U_L$, are well-defined and can be calculated. Let $U_1$ be the intersection point of the curves $\mathcal{W}_1(U_+)$ and $\mathcal{W}_2(U_R)$, and let $U_R$ be such that \begin{equation} \sigma_1(U_+,U_1) \ge v_+, \label{3.2} \end{equation} where $\sigma_1$ is the 1-shock speed. Then, the Riemann problem has a solution made up shocks, rarefaction waves, and a contact separated by states $U_+, U_1$ in phase I and by $V_\pm$ in phase II. Precisely, the solution can be described by: \begin{itemize} \item Solution in Phase I: $$ W_5(U_L,U_+)\oplus W_1(U_+,U_1)\oplus W_2(U_1,U_R). $$ \item Solution in Phase II: $$ W_3(V_L,V_-)\oplus W_5(V_-,V_+)\oplus W_4(V_+,V_R). $$ \end{itemize} \end{theorem} \begin{proof} Whenever the states are well-defined, the construction is clear. We will prove that these states exist. First, a strategy for computing these states can be done as follows. The solution in Phase I gives $v_\pm=v_0$ and 2 equations in \eqref{3.1}. The jump relations across the 5-contact give us 3 equations \begin{equation} \begin{gathered} \alpha_R\rho_+(u_+-v_0) =\alpha_L\rho_L(u_L-v_0):=m\,, \\ (u_+-v_0)^2+2h(\rho_+)=(u_L-v_0)^2+2h_L\,, \\ \beta_R q_+ - \beta_Lq_- +m(u_+-u_L) + (\alpha_Rp_+-\alpha_Lp_L)=0\,, \end{gathered}\label{3.3} \end{equation} where $p_+=p(\rho_+), q_\pm=q(\theta_\pm)=l\theta_\pm^\delta$ from the equation of state. The five equations \eqref{3.1}-\eqref{3.3} enable us to calculate the five quantities: $\theta_\pm, v_0,\rho_+,u_+$, and so give us the states $V_\pm, U_+$. Then, since $$ \{U_1\}=\mathcal{W}_1(U_+)\cap \mathcal{W}_2(U_R), $$ the state $U_1$ is determined by the equations $$ u_1=\omega_1(U_+;\rho_1) = \omega_2(U_R;\rho_1), $$ where $\omega_1$ and $\omega_2$ are defined by \eqref{2.11}, \eqref{2.12}, and \eqref{2.20}. We now prove that these states must exist. Indeed, eliminating $\theta_\pm$ from these five equations yields three equations for $\alpha_R,\rho_+,u_+,v_0$, as follows. Since the function $v=\omega_3(V_L;\theta)$ defined by \eqref{2.13} and \eqref{2.20} is decreasing as a function of $\theta$, it has an inverse function, denoted by $\theta=\omega_3^{-1}(V_L;v):=\theta_3(v)$, which is also decreasing in $v$. Similarly, since the function $v=\omega_4(V_R;\theta)$ defined by \eqref{2.14} and \eqref{2.20} is increasing as a function of $\theta$, it has an inverse function, denoted by $\theta=\omega_4^{-1}(V_R;v):=\theta_4(v)$, which is also increasing in $v$. From \eqref{3.1}, it holds that \begin{equation} \theta_{-}=\omega_3^{-1}(V_L;v_0),\quad \theta_{+}=\omega_4^{-1}(V_R;v_0). \label{3.4} \end{equation} Substituting \eqref{3.4} into \eqref{3.3} implies that $\alpha=\alpha_R$, $\rho=\rho_+$, $u=u_+,v=v_0$ satisfy \begin{equation} \begin{gathered} \alpha\rho(u-v) -\alpha_L\rho_L(u_L-v) =0\,, \\ (u-v)^2+2h(\rho)- (u_L-v)^2-2h_L=0\,,\\ \begin{aligned} &(1-\alpha) q_+(V_R;v) - \beta_Lq_-(V_L;v) +\alpha_L\rho_L(u_L-v)(u-u_L)\\ & + (\alpha p(\rho)-\alpha_Lp_L)=0, \end{aligned} \end{gathered}\label{3.5} \end{equation} and so $\rho_+,u_+,v_0 $ can be found in terms of $\alpha_R$, where \begin{equation} q_-(V_L;v)=q(\theta_3(V_L;v)), \quad q_+(V_R;v)=q(\theta_4(V_R;v)). \label{3.6} \end{equation} Since $q=q(\theta)$ is increasing in $\theta>0$, and $\theta=\theta_3(v)$ is decreasing and $\theta=\theta_4(v)$ is increasing, it holds that \begin{gather*} \frac{dq_-(V_L;v)}{dv}=\frac{dq}{d\theta}\frac{d\theta_3}{dv} <0,\\ \frac{dq_+(V_R;v)}{dv}=\frac{dq}{d\theta}\frac{d\theta_4}{dv} >0, \end{gather*} which mean that $q_-$ is decreasing in $v$, while $q_+$ is increasing in $v$. Let $\omega_3$ meet $\omega_4$ at $(\theta_*,v_*)$. That is, \begin{equation} v_*=\omega_3(V_L;\theta_{*}) =\omega_4(V_R;\theta_{*}). \label{3.7} \end{equation} This yields $q_+(V_R;v_*)=q_-(V_L;v_*)$. The three equations in \eqref{3.5} can be written as \begin{equation} F(X,Y)=0, \quad\text{where } X=\alpha,\; Y=(\rho,u,v).\label{3.8} \end{equation} It holds that $$ F(a,b)=0,\quad a=\alpha_L,\quad b=(\rho_L,u_L,v_*). $$ To see whether the implicit equation \eqref{3.8} can define a curve $Y=G(X)$ for $X$ near $a=\alpha_L$, we consider the matrix \begin{align*} &(\partial F_i(a,b)/\partial Y_j)_{i,j=\overline{1,3}}\\ &= \begin{pmatrix} \alpha_L(u_L-v_*)&\alpha_L\rho_L&0\\ 2h'(\rho_L)&2(u_L-v_*)&0\\ \alpha_L p'(\rho_L)&\alpha_L\rho_L(u_L-v_*) &(1-\alpha) \frac{dq_+}{dv}(v_*) - \beta_L \frac{dq_-}{dv}(v_*) \end{pmatrix}. \end{align*} As seen above, $q_+$ is increasing in $v$, and $q_-$ is decreasing in $v$, which yields $$ \frac{dq_+}{dv}(v_*) >0, \quad \frac{dq_-}{dv}(v_*)<0. $$ Thus, the determinant $$ |(\partial F_i(a,b)/\partial y_j)| = 2\alpha_L((1-\alpha)\frac{dq_+}{dv}(v_*) - \beta_L \frac{dq_-}{dv}(v_*) ) ((u_L-v_*)^2-h'(\rho_L)\rho_L) $$ has same sign as $((u_L-v_*)^2-h'(\rho_L)\rho_L)= (u_L-v_*)^2-p'(\rho_L) > 0$ when $U_L$ is in the supersonic region. So, the matrix $(\partial F_i(a,b)/\partial Y_j) $ is invertible. Applying the Implicit Function Theorem, we can see that there exists an interval $I$ containing $\alpha_L$ such that equation \eqref{3.8} determines a map $Y=Y(\alpha), \alpha\in I$. Thus, the states $U_\pm, V_\pm$ are well-defined for $\alpha_R\in I$. The inequality \eqref{3.2} implies that the contact wave $W_5(U_L,U_+)$ can be followed by the 1-wave $W_1(U_+,U_1)$. This completes the proof of Theorem \ref{theo41}. \end{proof} \subsection*{Case B: The contact belongs to the subsonic region $\lambda_1<\lambda_5<\lambda_2$} \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig5a} % Sol2xtI.png \quad \includegraphics[width=0.45\textwidth]{fig5b} % Sol2I.png \end{center} \caption{Case B: Riemann solutions in the phase I} \label{fig44} \end{figure} In Phase I, the solution begins with a 1-wave from $U_L$ to $U_-$, followed by a 5-contact from $U_-$ to $U_+$, then it arrives at $U_R$ by a 2-wave. In Phase II, the solution begins with a 3-wave from $V_L$ to $V_-$, followed by a 5-contact from $V_-$ to $V_+$, and then it arrives at $V_R$ by a 4-wave. The 1- and 3-waves may interchange the order, and the 2- and 4-waves may interchange the order, see Figure \ref{fig45}. \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig6} % Sol2xt.png \end{center} \caption{Case B: A whole Riemann solution where the contact is in the subsonic region} \label{fig45} \end{figure} \begin{theorem}\label{theo42} Let $U_*=(\rho_*,u_*)$ be the intersection point of $\mathcal{W}_1(U_L)$ and $\mathcal{W}_2(U_R)$ in the $(\rho,u)$-plane, and let $V_*=(\theta_*,v_*)$ be the intersection point of $\mathcal{W}_3(V_L)$ and $\mathcal{W}_4(V_R)$ in the $(\theta,v)$-plane. Providing that $(\rho_*,u_*,v_*)$ is in the subsonic region $$ (u_*-v_*)^20$, and $\theta=\theta_3(v)$ is decreasing and $\theta=\theta_4(v)$ is increasing, it holds that \begin{gather*} \frac{dq_-(V_L;v)}{dv}=\frac{dq}{d\theta}\frac{d\theta_3}{dv} <0,\\ \frac{dq_+(V_R;v)}{dv}=\frac{dq}{d\theta}\frac{d\theta_4}{dv} >0, \end{gather*} which means that $q_-$ is decreasing in $v$, while $q_+$ is increasing in $v$. Now, let $\mathcal{W}_1(U_L)$ meet $\mathcal{W}_2(U_R)$ at $(\rho_*,u_*)$. That is, \begin{equation} u_*=\omega_1(U_L;\rho_{*}) =\omega_2(U_R;\rho_{*}). \label{3.14} \end{equation} As above, let $\omega_3$ meet $\omega_4$ at $(\theta_*,v_*)$. That is, \begin{equation} v_*=\omega_3(V_L;\theta_{*}) =\omega_4(V_R;\theta_{*}). \label{3.15} \end{equation} This yields $$ q_+(V_R;v_*)=q_-(V_L;v_*). $$ For simplicity, we omit $U_L,U_R, V_L,V_R$ in \eqref{3.12} in the sequel whenever this does not cause any confusion. The three equations \eqref{3.12} have the form \begin{equation} F(X,Y)=0, \quad\text{where}\quad X=\alpha, Y=(\rho_-,\rho_+,v). \label{3.16} \end{equation} It holds that $$ F(a,b)=0,\quad a= \alpha_L, \quad b=(\rho_*,\rho_*,v_*). $$ To see whether the implicit equation \eqref{3.8} can define a curve $Y=G(X)$ for $X$ near $a=\alpha_L$, we consider the matrix $(\partial F_i(a,b)/\partial Y_j)_{i,j=\overline{1,3}}$, which is equal to $$ \begin{pmatrix} \alpha_L(u_*-v_*+\rho_*\omega_2'(\rho_*)) & -\alpha_L(u_*-v_*+\rho_*\omega_1'(\rho_*))& 0\\ 2(u_*-v_*)\omega'_2(\rho_*)+2h'(\rho_*)&-2(u_*-v_*)\omega'_1(\rho_*)-2h'(\rho_*)&0\\ \alpha_L\rho_*(u_*-v_*)\omega_2'(\rho_*) +\alpha_Lp'(\rho_*) & -\alpha_L\rho_*(u_*-v_*)\omega_1'(\rho_*) +\alpha_Lp'(\rho_*)& q_* \end{pmatrix}, $$ where $$ q_*=(1-\alpha_L)\frac{dq_+}{dv}(v_*) -\beta_L\frac{dq_-}{dv}(v_*). $$ As argued above, the function $q_+$ is increasing in $v$, and the function $q_-$ is decreasing in $v$. So, $$ \frac{dq_+}{dv}(v_*) >0, \quad \frac{dq_-}{dv}(v_*)<0. $$ This yields $q_*>0$. Thus, the determinant \begin{align*} &|(\partial F_i(a,b)/\partial y_j)| \\ &= q_* \begin{vmatrix} \alpha_L(u_*-v_*+\rho_*\omega_2'(\rho_*))& -\alpha_L(u_*-v_*+\rho_*\omega_1'(\rho_*))\\ 2(u_*-v_*)\omega'_2(\rho_*)+2h'(\rho_*)&-2(u_*-v_*)\omega'_1(\rho_*)-2h'(\rho_*)\\ \end{vmatrix}\\ &=q_*[\alpha_L(u_*-v_*+\rho_*\omega_2'(\rho_*))(-2(u_*-v_*)\omega'_1(\rho_*) -2h'(\rho_*))\\ &\quad +\alpha_L(u_*-v_*+\rho_*\omega_1'(\rho_*))(2(u_*-v_*)\omega'_2(\rho_*) +2h'(\rho_*))]\\ &=2q_*\alpha_L(u_*-v_*)^2(\omega_2'(\rho_*)-\omega_1'(\rho_*)) +\rho_*h'(\rho_*)(\omega_1'(\rho_*)-\omega_2'(\rho_*))\\ &=(\omega_2'(\rho_*)-\omega_1'(\rho_*))((u_*-v_*)^2-\rho_*h'(\rho_*))\\ &=(\omega_2'(\rho_*)-\omega_1'(\rho_*))((u_*-v_*)^2-p'(\rho_*)) <0, \end{align*} since $\omega_2'(\rho_*)>0$, $\omega_1'(\rho_*)<0$, and $U_*$ is in the subsonic region $(u_*-v_*)^2