Applied Mathematics Letters 135 (2023) 108442

Contents lists available at ScienceDirect

Applied Mathematics Letters

www.elsevier.com/locate/aml

The Riemann problem for the inhomogeneous pressureless Euler

equations
Hongjun Cheng ∗, Hanchun Yang
Department of Mathematics, Yunnan University, Kunming, 650091, Yunnan, PR China

article info abstract

Article history: The Riemann problem for the homogeneous pressureless Euler equations has
Received 12 July 2022 been solved well. This paper constructively solves the Riemann problem for the
Received in revised form 8 September inhomogeneous pressureless Euler equations with space-dependent gravity and
2022 constant friction. An effective method to deal with the Riemann problems for
Accepted 8 September 2022 the inhomogeneous pressureless Euler equations is presented.
Available online 14 September 2022
© 2022 Elsevier Ltd. All rights reserved.
Keywords:
Pressureless Euler equations
Riemann problem
Vacuum state
Delta-shock

1. Introduction

The pressureless Euler equations [1–3] can describe the motion of free particles sticking together under
collision and explain the formation of large-scale structures in the universe. The distinct feature of this model
is that two blowup mechanisms occur simultaneously and delta-shocks [4–10] may develop in solutions. The
Riemann problem for the homogeneous pressureless Euler equations has been well studied [1–3]. We are
interested in the Riemann problem for the inhomogeneous case.
As for the inhomogeneous pressureless Euler equations, Shen [11] and Zhang & Zhang [12] have respec-
tively solved the 1-D Riemann problem with a Coulomb-like friction and a constant friction by using variable
substitutions, where how to find a suitable variable substitution is a key point. One can also refer to [13–16]
for some studies about the inhomogeneous pressureless Euler equations and [17] about an inhomogeneous
pressureless two-phase flow model.
In this paper, we are interested in the Riemann problem for the following inhomogeneous pressureless
Euler equations {
ρt + (ρu)x = 0,
( ) (1.1)
(ρu)t + ρu2 x = kxρ − αρu

∗ Corresponding author.
E-mail addresses: hjcheng@ynu.edu.cn (H. Cheng), hyang@ynu.edu.cn (H. Yang).

https://doi.org/10.1016/j.aml.2022.108442
0893-9659/© 2022 Elsevier Ltd. All rights reserved.
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

subject to the initial data {

(ρ− , u− ), x < 0,
(ρ, u)(x, t = 0) = (1.2)
(ρ+ , u+ ), x > 0,
where ρ ≥ 0 denotes the density, u the velocity, k > 0 and α > 0 are two constants. Recalling that the
source term for a gravity reads ρΦx with Φ being the gravitational potential, which can be a given function
in atmospheric modeling or determined by the Poisson equation in the case of self-gravity, the source term
kxρ in (1.1) characterizes the gravity with the quadratic static potential Φ = kx2 /2 while the usual constant
gravity kρ corresponds to a linear static potential Φ = kx. Besides the gravity, the frictions are also common
external forces. The source term αρu in (1.1) just models the friction with constant coefficient α. In fact,
the model (1.1) is also very excellent from the mathematical point of view. Especially, the equation for
the characteristic curve is a second order homogeneous linear ordinary differential equation with constant
coefficients (i.e., (2.3)), whose solvability enables us to derive the explicit solutions to the Riemann problem
(1.1) and (1.2). In the current paper, we also expect to present an effective method to investigate the
Riemann problems for the inhomogeneous pressureless Euler equations.
Firstly, in Section 2, we solve the general solution, the vacuum state and the contact discontinuity, with
which and by characteristic analysis, we solve the Riemann problem for u− < u+ by the pattern of waves
consisting of a vacuum state bounded by two contact discontinuities. Secondly, in Section 3, we prove the
simultaneous occurrence of two blowup mechanisms, introduce the delta-shocks and clarify the generalized
Rankine–Hugoniot relation and entropy condition, and finally solve the Riemann problem for u− > u+ by
the pattern of waves containing a single delta-shock. Finally, Section 4 presents the conclusions for this paper.
We find that the source terms do not lead to a change of the structures of the Riemann solutions but influence
the Riemann solutions in the way of exponential functions. It is observed that the characteristic analysis is
very handy to deal with the Riemann problems for the inhomogeneous pressureless Euler equations.

2. The vacuum solution to (1.1)–(1.2)

The system (1.1) has a double eigenvalue λ = u and only one right eigenvector r = (1, 0)T satisfying
∇λ · r ≡ 0. So it is extremely nonstrictly hyperbolic and λ is linearly degenerate. The characteristic curve
for (1.1) satisfies
dx
= u, (2.1)
dt
along which it holds that
du dρ
= kx − αu, = −ρux . (2.2)
dt dt
Therefore the characteristic curve satisfies the second order homogeneous linear ordinary differential
equations with constant coefficients
d2 x dx
2
+α − kx = 0, (2.3)
dt dt
whose characteristic equation is τ 2 + ατ − k = 0. By solving (2.3), we know that the characteristic curve
can be expressed as
x = c1 eλ1 t + c2 eλ2 t , (2.4)

along which
u = c1 λ1 eλ1 t + c2 λ2 eλ2 t (2.5)
√ √
−α− α2 +4k −α+ α2 +4k
where λ1 = 2 , λ2 = 2 , and c1 , c2 are two arbitrary constants.
2
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

Consider the system (1.1) with constant initial data (ρ, u)(x, t = 0) = (ρ0 , u0 ). For any point (x0 , 0),
from (2.4) and (2.5), we have x0 = c1 + c2 and u0 = c1 λ1 + c2 λ2 , which give c1 = uλ0 −λ 2 x0
1 −λ2
0 −u0
, c2 = λλ1 x1 −λ 2
.
Therefore, the characteristic curve across the point (x0 , 0) is

u0 − λ2 x0 λ1 t λ1 x0 − u0 λ2 t λ1 eλ2 t − λ2 eλ1 t eλ1 t − eλ2 t

x= e + e = x0 + u0 , (2.6)
λ1 − λ2 λ1 − λ2 λ1 − λ2 λ1 − λ2
along which
u0 − λ 2 x0 λ1 x0 − u0 λ1 λ2 λ1 eλ1 t − λ2 eλ2 t
u= λ1 eλ1 t + λ2 eλ2 t = (eλ2 t − eλ1 t )x0 + u0 . (2.7)
λ1 − λ2 λ1 − λ2 λ1 − λ2 λ1 − λ2
Eliminating x0 yields to
λ1 − λ2 λ1 λ2
u= e(λ1 +λ2 )t u0 + (eλ2 t − eλ1 t )x. (2.8)
λ1 e 2t
λ − λ2 eλ 1 t λ1 e t − λ2 eλ1 t
λ 2

Substituting (2.8) into the second equation in (2.2), one solves that
λ1 − λ2
ρ= ρ0 . (2.9)
λ1 eλ2 t − λ2 eλ1 t
Therefore this constant initial value problem admits a unique smooth solution (2.8)–(2.9), which are called
as the general solution to (1.1).
It is obvious that the vacuum state

ρ ≡ 0, u = arbitrary smooth function (2.10)

is the solution of (1.1).

Suppose that x = x(t) is a bounded discontinuity of solution of (1.1) and (ρ̄l , ūl ) and (ρ̄r , ūr ) are the left-
and right-side limits on the discontinuity of the solution (ρ, u)(x, t). Then the Rankine–Hugoniot relation
{
−x′ (t)[ρ] + [ρu] = 0,
(2.11)
−x′ (t)[ρu] + [ρu2 ] = 0

holds, where [g] = ḡl − ḡr . Remark that although a source term appears in (1.1), the Rankine–Hugoniot
relation has the same form as that for the homogeneous pressureless Euler equations. By solving (2.11), we
get that
x′ (t) = ūl = ūr . (2.12)
This is a contact discontinuity, which is just the characteristic curves for both sides in (x, t)-plane. Two states
can be connected by a contact discontinuity if and only if ūl = ūr . The contact discontinuity in (x, t)-plane
is characterized by x′ (t) = ūl = ūr .
With the general solution, the vacuum state and the contact discontinuity, we construct the solution to
(1.1)–(1.2) for u− < u+ . Similarly to (2.6), we know that across the point (x− < 0, 0) and (x+ > 0, 0), the
characteristic curves are respectively
λ1 eλ2 t − λ2 eλ1 t eλ1 t − eλ2 t
x= x− + u− (2.13)
λ1 − λ2 λ1 − λ2
and
λ1 eλ2 t − λ2 eλ1 t eλ1 t − eλ2 t
x= x+ + u+ , (2.14)
λ1 − λ2 λ1 − λ2
in which
eλ1 t − eλ2 t λ1 eλ2 t − λ2 eλ1 t
> 0, >0
λ1 − λ2 λ1 − λ2
3
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

Fig. 2.1. Vacuum solution (left) and delta-shock solution (right).

√ √
−α− α2 +4k −α+ α2 +4k
because λ1 = 2 < 0, λ 2 = 2 > 0. When u− < u+ , it can be seen that the characteristic
curves emitting from the x-axis never intersect with each other and there is no characteristic curve in the
region
eλ1 t − eλ2 t eλ1 t − eλ2 t
u− < x < u+ . (2.15)
λ1 − λ2 λ1 − λ2
Thus when u− < u+ , we can construct the Riemann solution consisting of two contact discontinuities and
a vacuum state as follows (see the left in Fig. 2.1)
(λ1 − λ2 )e(λ1 +λ2 )t λ1 λ2 (eλ2 t − eλ1 t )
⎧( )
λ1 − λ2
⎪
⎪ ρ − , u − + x , x < x1 (t),
λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t λ1 eλ)2 t − λ2 eλ1 t
⎪
⎪
⎪
⎪
(λ1 − λ2 )e(λ1 +λ2 )t λ1 λ2 (eλ2 t − eλ1 t )
⎨ (
(ρ, u)(x, t)= 0, U (x, t) + x , x1 (t) < x < x2 (t),
⎪
⎪ λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t
(λ +λ )t λ t λ t
( )
λ1 − λ2 (λ1 − λ2 )e 1 2 λ1 λ2 (e 2 − e 1 )
⎪
⎪
⎪
⎪
⎩ ρ+ , u+ + x , x > x2 (t),
λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t
(2.16)
where U (x, t) satisfies U (x1 (t) + 0, t) = u− , U (x2 (t) − 0, t) = u+ ; and
eλ1 t − eλ2 t eλ1 t − eλ2 t
x1 (t) = u− , x2 (t) = u+ .
λ1 − λ2 λ1 − λ2

3. The delta-shock solution to (1.1)–(1.2)

3.1. Simultaneous occurrence of two blowup mechanisms

Consider smooth solution for (1.1) with sufficiently smooth initial condition (ρ, u)(x, 0) = (ρ0 , u0 )(x)
satisfying u′0 (x) < 0. Similarly to (2.6)–(2.7), it is known that across any point (x0 , 0), the characteristic
curve is
λ1 eλ2 t − λ2 eλ1 t eλ1 t − eλ2 t
x= x0 + u0 (x0 ), (3.1)
λ1 − λ2 λ1 − λ2
along which
λ1 λ2 λ1 eλ1 t − λ2 eλ2 t
u= (eλ2 t − eλ1 t )x0 + u0 (x0 ). (3.2)
λ1 − λ2 λ1 − λ2
By using the implicit function theorem, it follows that when (λ2 eλ1 t − λ1 eλ2 t ) + (eλ2 t − eλ1 t )u′0 (x0 ) > 0,
there exists the function u = u(x, t) satisfying
∂u λ1 λ2 (eλ1 t − eλ2 t ) + (λ2 eλ2 t − λ1 eλ1 t )u′0 (x0 )
= . (3.3)
∂x (λ2 eλ1 t − λ1 eλ2 t ) + (eλ2 t − eλ1 t )u′0 (x0 )
4
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

Then solving the second equation in (2.2) leads to

(λ2 − λ1 )ρ0 (x0 )

ρ= . (3.4)
(λ2 eλ1 t − λ1 eλ2 t ) + (eλ2 t − eλ1 t )u′0 (x0 )

Then one obtains that as (λ2 eλ1 t −λ1 eλ2 t )+(eλ2 t −eλ1 t )u′0 (x0 ) → 0, it holds that (ρ, ∂u/∂x) → (∞, ∞) along
the characteristic curve (3.1) from below. This shows that the solution ρ and the gradient ∂u/∂x must blow
up simultaneously at a finite time. The above discussion shows the mechanism of occurrence of delta-shocks.

3.2. Definition of delta-shocks

A two-dimensional weighted Dirac delta function w(ζ)δL supported on a smooth curve L parameterized
as x = x(ζ), t = t(ζ)(c ≤ ζ ≤ d) is defined by
⟨ ⟩ ∫ d
w(ζ)δL , ϕ = w(ζ)ϕ(x(ζ), t(ζ))dζ (3.5)
c

for any ϕ ∈ C0∞ (R2 ).

A delta-shock is a discontinuity on which at least one state variable contains the Dirac delta function.
Let Ω is a region cut by a smooth curve S : x = x(t) into the left part Ωl and the right part Ωr . Seek the
delta-shock of (1.1) in the form
⎧
⎨(ρ , u )(x, t), (x, t) ∈ Ωl ,
( l l
⎪
)
(ρ, u)(x, t) = w(t)δS , uδ (t) , (x, t) ∈ S, (3.6)
⎪
(ρr , ur )(x, t), (x, t) ∈ Ωr ,
⎩

where (ρl , ul )(x, t) and (ρr , ur )(x, t) are the solution of (1.1); w(t) and uδ (t) are smooth.
We assert that if the relation, called as the generalized Rankine–Hugoniot relation,

dx(t)
⎧
⎪
⎪ = uδ (t),
⎨ dt
⎪
⎪
dw(t)
= −[ρ]uδ (t) + [ρu], (3.7)
⎪ dt
dw(t)u δ (t)
⎪
⎪
− kx(t)w(t) + αw(t)uδ (t) = −[ρu]uδ (t) + [ρu2 ]
⎪
⎩
dt
is satisfied, then (3.6) satisfies (1.1) in the sense of distributions, that is, it holds that
⎧∫∫ ∫∫ ⟨ ⟩ ⟨ ⟩
⎪
⎪
⎪ ρϕ t dxdt + ρuϕ x dxdt + w(t)δ S , ϕt + w(t)uδ (t)δ S , ϕx = 0,
⎨∫∫Ω/S Ω/S
⎪
⎪ ∫∫ ∫ ∫
ρuϕt dxdt + ρu2 ϕx dxdt + (kxρ − αρu)ϕdxdt (3.8)
⎪
⎪ Ω/S Ω/S Ω/S ⟩
⎪
⎪ ⟨ ⟩ ⟨ ⟨( ) ⟩
+ w(t)uδ (t)δS , ϕt + w(t)u2δ (t)δS , ϕx + kx(t)w(t) − αw(t)uδ (t) δS , ϕ = 0
⎪
⎩

for any ϕ(x, t) ∈ C0∞ (Ω ). The proof is similar to that in [2,3]; we omit it.
To guarantee uniqueness, we propose the entropy condition for the delta-shock

λ(ūl ) > dx(t)/dt > λ(ūr ), (3.9)

which means that all characteristic curves on both sides of the delta-shock are incoming.
5
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

3.3. Delta-shock solution

When u− > u+ , it is found that the characteristic curves (2.13) and (2.14) will overlap in the region
eλ1 t − eλ2 t eλ1 t − eλ2 t
u− > x > u+ . (3.10)
λ1 − λ2 λ1 − λ2
Therefore we suggest the delta-shock solution to (1.1) and (1.2) as follows (see the right in Fig. 2.1)
(λ1 − λ2 )e(λ1 +λ2 )t λ1 λ2 (eλ2 t − eλ1 t )
⎧( )
λ1 − λ2
⎪
⎪
λ t λ t −
ρ , u − + x , x < x(t),
⎪( λ1 e 2 − λ2 e 1 λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t )
⎪
⎪
⎪
(λ1 − λ2 )e(λ1 +λ2 )t λ1 λ2 (eλ2 t − eλ1 t )
⎨
(ρ, u)(x, t) = w(t)δ|x=x(t) , λ t λ t
Uδ + x , x = x(t), (3.11)
⎪
⎪ λ1 e 2 − λ2 e 1 λ1 eλ2 t − λ2 eλ1 t )
(λ1 − λ2 )e(λ1 +λ2 )t λ1 λ2 (eλ2 t − eλ1 t )
(
λ1 − λ2
⎪
⎪
⎪
⎪
⎩ ρ + , u+ + x , x > x(t),
λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t λ1 eλ2 t − λ2 eλ1 t
where x(0) = 0, w(0) = 0 and Uδ is constant.
Let us determine x(t), w(t) and Uδ by solving the generalized Rankine–Hugoniot relation (3.7) under the
entropy condition (3.9), at this moment, which is equivalent to

u− > Uδ > u+ . (3.12)

It can be calculated that u′δ (t) − kx(t) + αuδ (t) = 0. Then (3.7) becomes
dx(t)
⎧
⎪
⎪ = uδ (t),
⎨ dt
⎪
⎪
dw(t)
= −[ρ]uδ (t) + [ρu], (3.13)
⎪ dt
⎩u (t) dw(t) = −[ρu]u (t) + [ρu2 ].
⎪
⎪
⎪
δ δ
dt
Firstly, we multiply the second equation in (3.13) by uδ (t) and then subtract it from the third one to get

[ρ]u2δ (t) − 2[ρu]uδ (t) + [ρu2 ] = 0,

which gives
√ √ √ √
ρ̄l ūl ± ρ̄r ūr (λ1 − λ2 )e(λ1 +λ2 )t ρ− u− ± ρ+ u+ λ1 λ2 (eλ2 t − eλ1 t )
uδ (t) = √ √ = λ t λ t
· √ √ + · x(t).
ρ̄l ± ρ̄r λ1 e 2 − λ2 e 1 ρ− ± ρ+ λ1 eλ2 t − λ2 eλ1 t
Due to the entropy condition (3.12), we take the sign “+”, for which
√ √
ρ− u− + ρ+ u+
Uδ = √ √ . (3.14)
ρ− + ρ+
Then, because of
( )2
λ1 − λ2 √
−[ρ]uδ (t) + [ρu] = e(λ1 +λ2 )t ρ− ρ+ (u− − u+ ),
λ1 eλ2 t − λ2 eλ1 t
the second equation in (3.13) with w(0) = 0 gives
eλ1 t − eλ2 t √
w(t) = ρ− ρ+ (u− − u+ ). (3.15)
λ1 eλ2 t − λ2 eλ1 t
Finally, from the first equation in (3.13) with x(0) = 0, we can solve that
√ √
eλ1 t − eλ2 t ρ− u− + ρ+ u+ eλ1 t − eλ2 t
x(t) = Uδ · = √ √ · . (3.16)
λ1 − λ2 ρ− + ρ+ λ1 − λ2
6
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

4. Conclusion

In this paper, we have successfully solved the Riemann problem for the pressureless Euler equations
with space-dependent gravity and constant friction by the characteristic analysis method and obtained two
kinds of structures. In the first, two non-constant states are connected by two contact discontinuities and
an intermediate vacuum state. In the second, two non-constant states are connected by a single delta-shock.
It is found that the source terms do not lead to a change of the structures of the Riemann solutions but
influence the Riemann solutions in the way of exponential functions.
Letting k → 0 and/or α → 0 in the vacuum solution (2.16) and the delta-shock solution (3.11), we can
obtain the solutions to the Riemann problems for the corresponding limit systems. Besides, letting t → +∞,
it can be found that both the vacuum solution (2.16) and the delta-shock solution (3.11) tend to (0, λ2 x).
By a similar process, one can also solve the Riemann problems considered in [11,12], respectively. Indeed,
this paper provides a general method to deal with the Riemann problems for the inhomogeneous pressureless
Euler equations.

Data availability

No data was used for the research described in the article.

Acknowledgments

This paper is supported by the NSF of China (11861073, 12061084), NSF of Yunnan Province (2019FY00-
3007) and the high level talent training program of Yunnan Province (2019281).

References

[1] F. Bouchut, On zero pressure gas dynamics, in: B. Perthame (Ed.), in: Advances in Kinetic Theory and Computing:
Selected Papers, Series on Advances in Mathematics for Applied Sciences: Volume 22, World Scientific, River Edge, NJ,
1994, pp. 171–190.
[2] J. Li, T. Zhang, S. Yang, The Two-Dimensional Riemann Problem in Gas Dynamics, first ed., Longman Harlow, 1998.
[3] W. Sheng, T. Zhang, The Riemann problem for the transportation equations in gas dynamics, Mem. Amer. Math. Soc.
137 (654) 1999, http://dx.doi.org/10.1090/memo/0654.
[4] D. Tan, T. Zhang, Y. Zheng, Delta shock waves as limits of vanishing viscosity for hyperbolic systems of conservation
laws, J. Differential Equations 112 (1994) 1–32, http://dx.doi.org/10.1006/jdeq.1994.1093.
[5] H. Yang, Riemann problems for a class of coupled hyperbolic systems of conservation laws, J. Differential Equations
159 (1999) 447–484, http://dx.doi.org/10.1006/jdeq.1999.3629.
[6] J. Li, Note on the compressible Euler equations with zero temperature, Appl. Math. Lett. 14 (2001) 519–523, http:
//dx.doi.org/10.1016/S0893-9659(00)00187-7.
[7] G. Chen, H. Liu, Formation of delta-shocks and vacuum states in the vanishing pressure limit of solutions to the
isentropic Euler equations, SIAM J. Math. Anal. 34 (2003) 925–938, http://dx.doi.org/10.1137/S0036141001399350.
[8] E. Yu Panov, V.M. Shelkovich, δ ′ -Shock waves as a new type of solutions to systems of conservation laws, J. Differential
Equations 228 (2006) 49–86, http://dx.doi.org/10.1016/j.jde.2006.04.004.
[9] L. Guo, W. Sheng, T. Zhang, The two-dimensional Riemann problem for isentropic Chaplygin gas dynamic system,
Commun. Pure Appl. Anal. 9 (2009) 431–458, http://dx.doi.org/10.3934/cpaa.2010.9.431.
[10] H. Cheng, H. Yang, Riemann problem for the isentropic relativistic Chaplygin Euler equations, Z. Angew. Math. Phys.
63 (2012) 429–440, http://dx.doi.org/10.1007/s00033-012-0199-7.
[11] C. Shen, The Riemann problem for the pressureless Euler system with the Coulomb-like friction term, IMA J. Appl.
Math. 81 (2016) 76–99, http://dx.doi.org/10.1093/imamat/hxv028.
[12] Y. Zhang, R. Zhang, The Riemann problem for the equations of constant pressure fluid dynamics with nonlinear damping,
Int. J. Non Linear Mech. 133 (2021) 103712, http://dx.doi.org/10.1016/j.ijnonlinmec.2021.103712.
[13] Y. Ding, F. Huang, On a nonhomogeneous system of pressureless flow, Q. Appl. Math. 62 (2004) 509–528, http:
//dx.doi.org/10.1090/qam/2086043.
[14] I. Gallagher, L. Saint-raymond, On pressureless gases driven by a strong inhomogeneous magnetic field, SIAM J. Math.
Anal. 36 (2005) 1159–1176, http://dx.doi.org/10.1137/S0036141003435540.
7
H. Cheng and H. Yang Applied Mathematics Letters 135 (2023) 108442

[15] S. Ha, F. Huang, Y. Wang, A global unique solvability of entropic weak solution to the one-dimensional pressureless
Euler system with a flocking dissipation, J. Differential Equations 257 (2014) 1333–1371, http://dx.doi.org/10.1016/j.
jde.2014.05.007.
[16] S. Keita, Y. Bourgault, Eulerian droplet model: Delta-shock waves and solution of the Riemann problem, J. Math. Anal.
Appl. 472 (2019) 1001–1027, http://dx.doi.org/10.1016/j.jmaa.2018.11.061.
[17] C. Shen, M. Sun, Exact Riemann solutions for the drift-flux equations of two-phase flow under gravity, J. Differential
Equations 314 (2022) 1–55, http://dx.doi.org/10.1016/j.jde.2022.01.009.