# Riemann problem for constant flow with single-point heating source

Changsheng Yu Chengliang Feng Zhiqiang Zeng Tiegang Liu

###### Abstract

This work focuses on the Riemann problem of Euler equations with global
constant initial conditions and a single-point heating source, which comes
from the physical problem of heating one-dimensional inviscid compressible
constant flow. In order to deal with the source of Dirac delta-function, we
propose an analytical frame of double classic Riemann problems(CRPs) coupling,
which treats the fluids on both sides of the heating point as two separate
Riemann problems and then couples them. Under the double CRPs frame, the
solution is self-similar, and only three types of solution are found. The
theoretical analysis is also supported by the numerical simulation.
Furthermore, the uniqueness of the Riemann solution is established with some
restrictions on the Mach number of the initial condition.

Keywords: hyperbolic balance law, non-homogeneous Euler equations,
$\delta$-singularity, Riemann problem

AMS subject classifications: 35L81,80A20

##  1 Introduction

The Riemann problem of the one-dimensional inviscid compressible flow with
global constant initial conditions and a singe-point heating source is studied
in this paper. The heating point is located at $x=0$. The governing equations
is given by

$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=S,$ (1.1)

where

$U=\begin{pmatrix}\rho\\\ \rho u\\\ E\end{pmatrix},F=\begin{pmatrix}\rho u\\\
\rho u^{2}+p\\\ (E+p)u\end{pmatrix},S=\begin{pmatrix}0\\\ 0\\\
Q\delta(x)\end{pmatrix}.$

Here, $\rho$, $p$ and $E$ denote the thermodynamical variables: density,
pressure and total energy, respectively. $u$ is velocity. $Q>0$ is the heat
flux per unit time added to the flow. $\delta(x)$ is the Dirac delta-function.
The source means that $Q$ heat is added to the flow from the heating point per
unit time in the physical sense. We assume that the fluid is polytropic ideal
and the equaiotn of state is given by

$p=(\gamma-1)\rho e,\quad 1\leq\gamma\leq 3.$

where $\gamma$ is the ratio of specific heats and $e$ is the internal energy.
The initial condition is

$U(x,0)\equiv U_{1}=(\rho_{1},\rho_{1}u_{1},E_{1}),$ (1.2)

where $\rho_{1}$, $u_{1}$ and $E_{1}$ are constant. There is no loss of
generality in assuming $u_{1}>0$. The subscript ”1” means the initial state in
this paper. The physical problem described by the Riemann problem (1.1) and
(1.2) is that $Q$ heat is added to the one-dimensional constant flow per unit
time, as shown in Figure1.

Figure 1: one-dimensional constant flow with heat addition from a single point

The subscripts ”-” and ”+” denote the limiting states upstream and downstream
of the heating point, respectively. Write

$U_{-}(t)=U(0-,t),\quad U_{+}(t)=U(0+,t).$

The source term implies a jump of the energy flux across the heating point.

$\displaystyle\rho_{-}(t)u_{-}(t)$ $\displaystyle=\rho_{+}(t)u_{+}(t)$ (1.3a)
$\displaystyle\rho_{-}(t){u_{-}(t)}^{2}+p_{-}(t)$
$\displaystyle=\rho_{+}(t){u_{+}(t)}^{2}+p_{+}(t)$ (1.3b)
$\displaystyle(E_{-}(t)+p_{-}(t))u_{-}(t)+Q$
$\displaystyle=(E_{+}(t)+p_{+}(t))u_{+}(t).$ (1.3c)

If the velocity at the heating point is zero, thermal convection does not take
effect, then the heat addition has no effect on the flow. Thus in the
remainder of this paper we assume that

$u_{-}(t)\neq 0,\quad u_{+}(t)\neq 0.$

We only deal with the situation of constant $U_{-}(t)$ and $U_{+}(t)$.
Omitting the time parameter $t$, we can get the following equation from (1.3).

$\left(\frac{1}{2}{u_{-}}^{2}+h_{-}\right)(1+k)=\frac{1}{2}{u_{+}}^{2}+h_{+},$

where $h$ is enthalpy.
$k=\frac{Q}{\rho_{-}(t)u_{-}(t)\left(\frac{1}{2}{u_{-}}^{2}+h_{-}\right)}$ is
called the heating parameter and is assumed to be a constant parameter.

The heating equations are defined by

$\left\\{\begin{aligned} &\rho_{-}u_{-}=\rho_{+}u_{+}\\\
&\rho_{-}{u_{-}}^{2}+p_{-}=\rho_{+}{u_{+}}^{2}+p_{+}\\\
&\left(\frac{1}{2}{u_{-}}^{2}+h_{-}\right)(1+k)=\frac{1}{2}{u_{+}}^{2}+h_{+}\end{aligned}.\right.$
(1.4)

The solution of the heating equations (1.4) corresponds to a steady solution
of the equations (1.1).

The mathematical model of Figure1 can be applied to the study of one-
dimensional condensation problem (see [8, 9]). The condensation of the vapor
leads to the release of latent heat and therefore has heating effects on the
carrier fluid. Previous researches on condensation problems have revealed some
information about the solution of Riemann problem (1.1) and (1.2). The
solution of the heating equations (1.4) has already been available in [8, 7,
9, 20]. In [20] Schnerr made a good summary of the properties of the solution
to (1.4) at the case of $M_{-}>1$, which is typical for condensation problems.
Schnerr showed that there are two solutions of the heating equations (1.4),
one called the shock solution, which reduces to identity when $Q=0$, and the
other called the weak solution, which reduces to the adiabatic normal shock
solution when $Q=0$. For the heat addition of subsonic flow, Schnerr predicted
the appearance of the unsteady solution, but did not do a more detailed
analysis. For the unsteady solution of this Riemann problem, Dongen et al.[9]
studied the unsteady effects of the heat addition, and proposed three possible
solution structures. Those three structures are determined by the Mach numbers
of the fluid around the heating point. However, a complete and rigorous
theoretical proof is lacking at present. On the other hand, numerical
simulation as an effective tool has been employed to explore the wave patterns
of the heating problem. Chengwan et al.[6] used the ASCE method[17] to
numerically verify the above three structures by the simulation of the onset
of condensation in a slender Laval nozzle. In the wet nitrogen condensation
problem of the Laval nozzle, Chengwan showed the transition between the three
structures by adjusting the humidity. To our knowledge, there is no effective
theoretical method to study the exact solution of Riemann problem (1.1) and
(1.2).

Conservative hyperbolic equations with source terms can be transformed into
non-conservative hyperbolic equations without source terms. Adopting the
following form for the $\delta$-function

$\delta(x)=\frac{\partial H(x)}{\partial x},\quad H(x)=\begin{cases}0,\quad
x<0\\\ 1,\quad x>0\end{cases},$

where $H(x)$ is called Heaviside function, one can rewrite the Riemann problem
(1.1) and (1.2) into following non-conservative form.

$\frac{\partial\widetilde{U}}{\partial
t}+\widetilde{A}\frac{\partial\widetilde{U}}{\partial x}=0,$ (1.5)

where

$\widetilde{U}=\begin{pmatrix}\rho\\\ \rho u\\\ E\\\
h\end{pmatrix},\widetilde{A}=\begin{pmatrix}0&1&0&0\\\
\frac{1}{2}(\gamma-3)u^{2}&(3-\gamma)u&\gamma-1&0\\\
\frac{1}{2}(\gamma-2)u^{3}-\frac{c^{2}u}{\gamma-1}&\frac{3-2\gamma}{2}u^{2}+\frac{c^{2}}{\gamma-1}&\gamma
u&-Q\\\ 0&0&0&0\end{pmatrix},$

where $c$ is the speed of sound. The initial condition for the augmented
equations is

$\widetilde{U}(x,0)=\begin{pmatrix}\rho_{0}&\rho_{0}u_{0}&E_{0}&H(x)\end{pmatrix}^{T}.$

By defining proper entropy solution for non-conservative equations (1.5), one
can analyzed the solutions of the original Riemann problem and design
appropriate numerical methods to make numerical simulation (see [1, 10, 11]).
Applications of that method include the fluid in a nozzle with discontinuous
cross-sectional area(see [13, 14, 21]) and the shallow water equations with
discontinuous topography(see [2, 5, 15, 18, 22]). The augmented equations
(1.5) give rise to an additional linearly degenerated characteristic field and
an additional stationary discontinuity. Although the augmented equations do
not contain source terms, the solution process of the Riemann problem is very
complex and often problem-related due to the lack of conservation or strict
hyperbolicity.

The Riemann problem of homogeneous Euler equations with piecewise constant
initial conditions is called classical Riemann problem(CRP). Altough the
addition of the source term does not lead to the loss of self-similarity, the
singularity of the source term has a significant effect on the structure of
the solution. The Riemann solution of the Euler equations with smooth source
terms, which are often called generalized Riemann problem (GRP, see [3, 4, 24,
25]), has the same wave patterns as the Riemann solution of its corresponding
homogeneous Euler equations. However, the singularity source term affects the
wave patterns of the Riemann solution. In fact, as can be seen from [8], the
solution of Riemann problem (1.1) and (1.2) may be a four waves structure or a
five waves structure, which is far different from the wave structure of the
CRP.

In this work, an analytical frame, which is called the double CRPs frame, is
proposed to construct the exact solution of the Riemann problem (1.1) and
(1.2). It regards the fluids on both sides of the heating point as two
separate CRPs, and gives the upper limit of the number of waves first. The
solutions of the two CRPs are then coupled on the premise of maintaining the
physical properties of the heating point. Depending on the heating properties
and gasdynamics properties, the extra waves are deleted and the type of wave
is finally determined. One advantage of the double CRPs frame is that it is
independent of whether the fluid at the heating point is supersonic or
subsonic. Under this frame we demonstrate three possible structures of the
solution, which are verified by numerical tests in Section5.

The text is arranged as follows. In Section2, we will introduce the solution
of the heating equations and derive several useful properties. In Section3, an
analytical frame of double CRPs coupling will be introduced to constructively
solve the Riemann problem (1.1) and (1.2). Under this frame We will prove that
there are at most three types of the solution. In Section4, the structure of
the solution will be associated with the Mach number of the incoming flow to
illustrate the uniqueness of the solution. In Section5, we will give an
iterative method for the solution of each structure. Then we will give five
tests for the Riemann problem (1.1) and (1.2) with different initial
conditions and heating coefficients. Finally, a brief summary will be given in
Section6.

##  2 Solutions of the heating equations

According to [20], there are two branches of the solution to the heating
equations (1.4). We can choose the physical solution from these two branches
by the following property.

###### Property 2.1

If $M_{-}<1$, then $M_{+}\leq 1$. If $M_{-}>1$, then $M_{+}\geq 1$.

The physical solution corresponds to the weak solution in [20]. In this paper
the expression of Dongen et al.[9] is adopted. For the heat addition of
subsonic flow, the solution is

$\displaystyle
I\equiv\sqrt{\left(\gamma+\frac{1}{{M_{-}}^{2}}\right)^{2}-2(\gamma+1)\left(\frac{1}{{M_{-}}^{2}}+\frac{\gamma-1}{2}\right)(1+\kappa)},$
(2.1a)
$\displaystyle\frac{u_{+}}{u_{-}}=\frac{\rho_{-}}{\rho_{+}}=\frac{1}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}-I\right),$
(2.1b)
$\displaystyle\frac{p_{+}}{p_{-}}=\frac{{M_{-}}^{2}}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}+\gamma
I\right),$ (2.1c) $\displaystyle
M_{+}=\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}-I}{\gamma+\frac{1}{{M_{-}}^{2}}+\gamma
I}}.$ (2.1d)

For the heat addition of supersonic flow, the solution is

$\displaystyle
I\equiv\sqrt{\left(\gamma+\frac{1}{{M_{-}}^{2}}\right)^{2}-2(\gamma+1)\left(\frac{1}{{M_{-}}^{2}}+\frac{\gamma-1}{2}\right)(1+\kappa)},$
(2.2a)
$\displaystyle\frac{u_{+}}{u_{-}}=\frac{\rho_{-}}{\rho_{+}}=\frac{1}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}+I\right),$
(2.2b)
$\displaystyle\frac{p_{+}}{p_{-}}=\frac{{M_{-}}^{2}}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}-\gamma
I\right),$ (2.2c) $\displaystyle
M_{+}=\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}+I}{\gamma+\frac{1}{{M_{-}}^{2}}-\gamma
I}}.$ (2.2d)

$M$ is Mach number. The advantage of this expression is that the ratios of
variables before and after heat addition are only related to the upstream Mach
number of the heating point. In order to make the solution reasonable, there
is an upper bound on $k$ as follws.

$k\leq
k_{max}\mathop{=}\limits^{def}\frac{(1-{M_{-}}^{2})^{2}}{2(\gamma+1){M_{-}}^{2}\left(1+\frac{\gamma-1}{2}{M_{-}}^{2}\right)}.$
(2.3)

Figure 2: Relation between the maximum heating parameter and the upstream Mach
number of the heating point.

$k_{max}$ is called maximum heating parameter. The relation between $k_{max}$
and ${M_{-}}$ is shown in Figure2. When ${M_{-}}\to 0+$, $k_{max}$ approaches
asymptotically $+\infty$. When ${M_{-}}\to\infty$, $k_{max}$ approaches
asymptotically a finite number $\frac{1}{\gamma^{2}-1}$. The roots of the
equation

$k=\frac{(1-{M}^{2})^{2}}{2(\gamma+1){M}^{2}\left(1+\frac{\gamma-1}{2}{M}^{2}\right)}$
(2.4)

are

$M^{1,2}=\begin{cases}\sqrt{\frac{k(\gamma+1)+1\pm(\gamma+1)\sqrt{k(k+1)}}{1-k(\gamma^{2}-1)}},\quad\text{if}\
k(\gamma^{2}-1)\neq 1\\\ \sqrt{\frac{1}{2[k(\gamma+1)+1]}},\quad\text{if}\
k(\gamma^{2}-1)=1\end{cases}$

We denote

$M_{*}=\begin{cases}\sqrt{\frac{k(\gamma+1)+1-(\gamma+1)\sqrt{k(k+1)}}{1-k(\gamma^{2}-1)}},&\quad\text{if}\quad
k(\gamma^{2}-1)\neq 1\\\
\sqrt{\frac{1}{2[k(\gamma+1)+1]}},&\quad\text{if}\quad
k(\gamma^{2}-1)=1\end{cases}$
$M_{**}=\begin{cases}\sqrt{\frac{k(\gamma+1)+1+(\gamma+1)\sqrt{k(k+1)}}{1-k(\gamma^{2}-1)}},&\quad\text{if}\quad
k(\gamma^{2}-1)<1\\\ +\infty,&\quad\text{if}\quad k(\gamma^{2}-1)\geq
1\end{cases}$

It can easily be checked that $0<M_{*}<1$ and $M_{**}>1$, thus we arrive at
the following conclusion.

###### Lemma 2.1

  * (i)

If $k(\gamma^{2}-1)<1$, then $M_{-}\leq M_{*}$ or $M_{-}\geq M_{**}$.

  * (ii)

If $k(\gamma^{2}-1)\geq 1$, then $M_{-}\leq M_{*}$.

  * (iii)

If $M_{-}=M_{*}$ or $M_{-}=M_{**}$, then $I=0$ and $M_{+}=1$.

###### Lemma 2.2

If $k(\gamma^{2}-1)<1$, we have

$M_{*}=\sqrt{\frac{(\gamma-1)M_{**}^{2}+2}{2\gamma M_{**}^{2}-\gamma+1}},\quad
M_{**}=\sqrt{\frac{(\gamma-1)M_{*}^{2}+2}{2\gamma M_{*}^{2}-\gamma+1}}.$ (2.5)

The proof is trivial.

The downstream fluid is called thermal choked if $M_{+}=1$. If the upstream
fluid is sonic, then $k_{max}=0$, which means that any heat addition is not
allowed for the sonic flow. We now turn to the relations of the upstream and
downstream fluids.

###### Theorem 2.1

For the heat additon of subsonic flow we have $u_{+}>u_{-}$, $p_{+}<p_{-}$,
$\rho_{+}<\rho_{-}$ and $M_{+}>M_{-}$. For the heat additon of supersonic flow
we have $u_{+}<u_{-}$, $p_{+}>p_{-}$, $\rho_{+}>\rho_{-}$ and $M_{+}<M_{-}$.

###### Proof 2.1

For the heat addition of subsonic fluid, $k>0$ implies

$I<\sqrt{\left(\gamma+\frac{1}{{M_{-}}^{2}}\right)^{2}-2(\gamma+1)\left(\frac{1}{{M_{-}}^{2}}+\frac{\gamma-1}{2}\right)}=\frac{1}{{M_{-}}^{2}}-1,$
$\frac{u_{+}}{u_{-}}>\frac{1}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}+1-\frac{1}{{M_{-}}^{2}}\right)=1,$
$\frac{\rho_{+}}{\rho_{-}}=\frac{u_{-}}{u_{+}}<1,$
$\frac{p_{+}}{p_{-}}<\frac{{M_{-}}^{2}}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}-\gamma+\frac{\gamma}{{M_{-}}^{2}}\right)=1,$
$M_{+}=\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}-I}{\gamma+\frac{1}{{M_{-}}^{2}}+\gamma
I}}>\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}-\frac{1}{{M_{-}}^{2}}+1}{\gamma+\frac{1}{{M_{-}}^{2}}+\gamma(\frac{1}{{M_{-}}^{2}}-1)}}={M_{-}}.$

For the heat addition of supsonic fluid, $k>0$ implies

$I<\sqrt{\left(\gamma+\frac{1}{{M_{-}}^{2}}\right)^{2}-2(\gamma+1)\left(\frac{1}{{M_{-}}^{2}}+\frac{\gamma-1}{2}\right)}=1-\frac{1}{{M_{-}}^{2}},$
$\frac{u_{+}}{u_{-}}<\frac{1}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}+1-\frac{1}{{M_{-}}^{2}}\right)=1,$
$\frac{\rho_{+}}{\rho_{-}}=\frac{u_{-}}{u_{+}}>1,$
$\frac{p_{+}}{p_{-}}>\frac{{M_{-}}^{2}}{\gamma+1}\left(\gamma+\frac{1}{{M_{-}}^{2}}-\gamma+\frac{\gamma}{{M_{-}}^{2}}\right)=1,$
$M_{+}=\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}+I}{\gamma+\frac{1}{{M_{-}}^{2}}-\gamma
I}}<\sqrt{\frac{\gamma+\frac{1}{{M_{-}}^{2}}-\frac{1}{{M_{-}}^{2}}+1}{\gamma+\frac{1}{{M_{-}}^{2}}+\gamma(\frac{1}{{M_{-}}^{2}}-1)}}={M_{-}}.$

Given any $k>0$ and $1<\gamma<3$, we define

$\displaystyle\phi(M)=\frac{1}{\gamma+1}\left(\gamma+\frac{1}{{M}^{2}}-I\right),$
$\displaystyle\psi(M)=\frac{{M}^{2}}{\gamma+1}\left(\gamma+\frac{1}{{M}^{2}}+\gamma
I\right),$

where

$I=\sqrt{\left(\gamma+\frac{1}{{M}^{2}}\right)^{2}-2(\gamma+1)\left(\frac{1}{{M}^{2}}+\frac{\gamma-1}{2}\right)(1+\kappa)}.$

According to 2.1 and 2.2, $\frac{u_{+}}{u_{-}}$ and $\frac{p_{+}}{p_{-}}$ are
both functions of $M_{-}$. For the subsonic heat addition, we have

$\frac{u_{+}}{u_{-}}=\frac{\rho_{-}}{\rho_{+}}=\phi(M_{-}),\quad\frac{p_{+}}{p_{-}}=\psi(M_{-}).$
(2.6)

A tedious compution gives the following theorem.

###### Theorem 2.2

$\phi^{\prime}(M)>0,\quad\psi^{\prime}(M)<0$.

##  3 Structure of solution

A basic assumption of our double CRPs frame is that $U_{-}(t)$ and $U_{+}(t)$
are both constant vectors. The solution $U(x,t)$ of the Riemann problem (1.1)
and (1.2) satisfies the homogeneous Euler equations in both the left half
$\\{(x,t)|x<0,t\geq 0\\}$ and the right $\\{(x,t)|x>0,t\geq 0\\}$.
Consequently, $U(x,t)$ is self-similar. The exact solution in this paper
refers to the self-similar solution under this frame.

The exact solution $U(x,t)$ in the left half $\\{(x,t)|x<0,t\geq 0\\}$
satisfies the classical Euler equations, hence is the left half of the CRP
solution with $U_{1}$ and $U(0-,t)$ as the left and right initial conditions.
From the theory of CRP solution, we know that $U(x,t)$ in the left half
consists three discontinuities at most, which are two genuinely nonlinear
waves, namely shock waves or rarefaction waves, and a contact discontinuity
corresponding to the characteristic fields $u-a$, $u+a$ and $u$, respectively.
Similarly, $U(x,t)$ in the right half consists two genuinely nonlinear waves
and a contact discontinuity at most.

###### Theorem 3.1

The self-similar solution of the Riemann problem (1.1) and (1.2) consists of
seven discontinuities at most. They are a heating discontinuity at $x=0$, two
genuinely nonlinear waves and a contact discontinuity left to $x=0$, two
genuinely nonlinear waves and a contact discontinuity right to $x=0$
respectively, as shown in Figure3.

Figure 3: All possible waves for the Riemann problem (1.1) and (1.2)

The elementary waves on the left and right sides are denoted as $WL_{1}$,
$WL_{2}$, $WL_{3}$ and $WR_{1}$, $WR_{2}$, $WR_{3}$, respectively. The eight
constant regions are labeled $(1)\sim(8)$, respectively. $U_{1}\sim U_{8}$ are
the states in regoins $(1)\sim(8)$ and it is clear that $U_{1}=U_{8}$. Note
that the t-axis is a discontinuity and $U_{4}\neq U_{5}$. We adopt similiar
way to express the solution structure as in [21]. For examples,
$S(U_{1},U_{2})$ and $R(U_{1},U_{2})$ mean two states $U_{1}$ and $U_{2}$ are
connected by a shock and a rarefaction wave respectively, $C(U_{2},U_{3})$
means $U_{2}$ and $U_{3}$ are connected by a contact discontinuity, and
$H(U_{4},U_{5})$ means $U_{4}$ and $U_{5}$ are connected by a heating
discontinuity. All the six elementary waves in Figure3 can not exist at the
same time. The next step in our double CRPs coupling method is to eliminate
the redundant waves according to the heat addition properties and the
gasdynamics properties. We will give the main results in Theorem3.2.

###### Theorem 3.2

Under the double CRPs frame, there are three different structures of the exact
solution to the Riemann problem (1.1) and (1.2) as follows.

  * (1)

Type 1: if $M_{4}<M_{5}<1$, then the structure could be

$S(U_{1},U_{4})\oplus H(U_{4},U_{5})\oplus C(U_{5},U_{7})\oplus
S(U_{7},U_{8})$

.

  * (2)

Type 2: if $M_{4}<M_{5}=1$, then the structure could be

$S(U_{1},U_{4})\oplus H(U_{4},U_{5})\oplus R(U_{5},U_{6})\oplus
C(U_{6},U_{7})\oplus S(U_{7},U_{8})$

.

  * (3)

Type 3: if $M_{4}>M_{5}\geq 1$, then the structure could be

$H(U_{1},U_{5})\oplus R(U_{5},U_{6})\oplus C(U_{6},U_{7})\oplus
S(U_{7},U_{8})$

.

Figure4 depicts the three wave patterns in Theorem3.2. They will be denoted by
Type 1, Type 2 and Type 3, respectively. According to the Mach number at the
heating point and whether the thermal choked state appears, the proof of the
Theorem3.2 falls naturally into three parts: Lemma3.1, Lemma3.2 and Lemma3.3,
which correspond to Type 1, Type 2 and Type 3, respectively.

Figure 4: Three different structures of the exact solution to the Riemann
problem (1.1) and (1.2).

###### Lemma 3.1

Under the double CRPs frame, if $M_{4}<M_{5}<1$, then $WL_{2}$, $WL_{3}$ and
$WR_{1}$ do not exist, and then $u_{5}>u_{4}>0$; if $WL_{1}$ and $WR_{3}$
exist, they are both shock waves.

###### Proof 3.1

We first outline the proof. The problem is reduced to the determination of the
upstream fluid $U_{4}$ and the downstream fluid $U_{5}$ (either of which can
be determined by the other) at the heating point, such that $U_{4}$ and
$U_{1}$ can be connected by the fundamental waves with non-positive speed, and
that $U_{5}$ and $U_{8}=U_{1}$ can be connected by the fundamental waves with
non-negative speed.

According to (1.3), it follows that $u_{4}$ and $u_{5}$ have the same sign.
$WL_{3}$ and $WR_{1}$ do not exist since $M_{4}<1$ and $M_{5}<1$. There are
three cases for the sign of the velocity at heating point, as follows.

  * (1)

If $u_{4}>0$, then $u_{5}>0$ and $WL_{2}$ does not exist;

  * (2)

If $u_{4}<0$, then $u_{5}<0$ and $WR_{2}$ does not exist;

  * (3)

If $u_{4}=u_{5}=0$, then both the speeds of $WL_{2}$ and $WR_{2}$ are zero and
the widths of region $3$ and region $6$ are both zero.

These three structures are shown in Figure5.

Figure 5: Three possible structures when subsonic flow is heated without
thermal choke.

If $u_{4}<0,u_{5}<0$, then $WL_{1}$ is a shock wave and $WR_{3}$ is a
rarefaction wave. From Theorem2.1 we have $p_{4}<p_{5}$. The relations of the
elenmentary wave (see for instance [23]) imply $p_{1}\leq p_{2}=p_{4}$ and
$p_{5}\leq p_{8}$. It results $p_{5}\leq p_{1}$, which is contradictory with
$p_{5}>p_{4}\geq p_{1}$, hence the structure of $u_{4}<0,u_{5}<0$ does not
exist.

Using a similar method we can prove that the structure of $u_{4}=u_{5}=0$ dose
not exist. The structure of the solution is the left figure in Figure 5.

The next thing to do is to prove that $WL_{1}$ is a shock. If the assertion
would not hold, then $WL_{1}$ is a rarefaction wave, which implies
$p_{7}=p_{5}<p_{4}\leq p_{1}=p_{8}$. It follows that $WR_{3}$ is a rarefaction
wave and $u_{7}=u_{5}>u_{4}\geq u_{1}$. Applying the conservations of mass and
momentum at the control volume $[x_{0},x_{4}]\times[0,T]$ in the $x-t$ space,
as shown in the left figure of Figure6, we have

$\displaystyle\int_{x_{0}}^{x_{4}}\rho(x,T)dx$
$\displaystyle=\int_{x_{0}}^{x_{4}}\rho(x,0)dx,$
$\displaystyle\int_{x_{0}}^{x_{4}}\rho(x,T)u(x,T)dx$
$\displaystyle=\int_{x_{0}}^{x_{4}}\rho(x,0)u(x,0)dx.$

$U$ is constant at each region, thus

$\displaystyle\int_{x_{0}}^{x_{1}}\rho(x,T)dx+\rho_{4}|x_{1}|+\rho_{5}|x_{2}|+\rho_{7}|x_{3}-x_{2}|+\int_{x_{3}}^{x_{4}}\rho(x,T)dx$
(3.1)
$\displaystyle=\int_{x_{0}}^{x_{1}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx+\rho_{4}|x_{1}|\frac{u_{4}}{u_{1}}+\rho_{5}|x_{2}|\frac{u_{5}}{u_{1}}+\rho_{7}|x_{3}-x_{2}|\frac{u_{7}}{u_{1}}+\int_{x_{3}}^{x_{4}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx.$

The velocity inside the rarefaction wave is monotonous, it follows that
$u(x,T)\geq max\\{u_{1},u_{4}\\}\geq u_{1}$ for $x_{0}\leq x\leq x_{1}$ and
$u(x,T)\geq max\\{u_{5},u_{8}\\}=u_{5}>u_{1}$ for $x_{3}\leq x\leq x_{4}$,
which imply

$\int_{x_{0}}^{x_{1}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx\geq\int_{x_{0}}^{x_{1}}\rho(x,T)dx,$
$\int_{x_{3}}^{x_{4}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx\geq\int_{x_{3}}^{x_{4}}\rho(x,T)dx.$

$\frac{u_{4}}{u_{1}}\geq 1$, $\frac{u_{5}}{u_{1}}>1$ and
$\frac{u_{7}}{u_{1}}>1$ hold for the right side of Lemma3.1. Substituting the
above inequalities to Lemma3.1, we have $x_{2}=x_{3}=0$ and $u_{5}=0$, which
leads to a contradiction.

Finally, we have to show that $WR_{3}$ is a shock. If the assertion is false,
then $WR_{3}$ is a non-degenerate rarefaction wave, which implies
$u_{4}<u_{5}=u_{7}\leq u_{8}=u_{1}$. Applying the conservations of mass and
momentum at the control volume $[x_{0},x_{3}]\times[0,T]$ in the $x-t$ space,
as shown in the right figure of Figure6, we have

Figure 6: The integral diagram in the proof of Lemma3.1. Left: $WL_{1}$ is a
rarefaction wave and $WR_{3}$ is a rarefaction wave. Right: $WL_{1}$ is a
shock and $WR_{3}$ is a rarefaction wave.

$\displaystyle\int_{x_{0}}^{x_{3}}\rho(x,T)dx$
$\displaystyle=\int_{x_{0}}^{x_{3}}\rho(x,0)dx,$
$\displaystyle\int_{x_{0}}^{x_{3}}\rho(x,T)u(x,T)dx$
$\displaystyle=\int_{x_{0}}^{x_{3}}\rho(x,0)u(x,0)dx.$

Thus

$\displaystyle\rho_{4}|x_{0}|+\rho_{5}|x_{1}|+\rho_{7}|x_{2}-x_{1}|+\int_{x_{2}}^{x_{3}}\rho(x,T)dx$
$\displaystyle=\rho_{4}|x_{0}|\frac{u_{4}}{u_{1}}+\rho_{5}|x_{1}|\frac{u_{5}}{u_{1}}+\rho_{7}|x_{2}-x_{1}|\frac{u_{7}}{u_{1}}+\int_{x_{2}}^{x_{3}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx.$

For the velocity $u$ insides the rarefaction wave $WR_{3}$, we have
$u(x,T)\leq max(u_{5},u_{8})=u_{8}=u_{1}$, which implies

$\int_{x_{2}}^{x_{3}}\rho(x,T)\frac{u(x,T)}{u_{1}}dx\leq\int_{x_{2}}^{x_{3}}\rho(x,T)dx.$

$\frac{u_{4}}{u_{1}}<1$, $\frac{u_{5}}{u_{1}}\leq 1$ and
$\frac{u_{7}}{u_{1}}\leq 1$ hold for the right side of Lemma3.1. We see at
once that $\frac{u_{7}}{u_{1}}=1$ and $U_{7}=U_{8}$, which lesds to a
contradiction.

The following two lemmas can be proved using a similar method.

###### Lemma 3.2

Under the double CRPs frame, if $M_{4}<M_{5}=1$, then $WL_{2}$ and $WL_{3}$ do
not exist, and then $u_{5}>u_{4}>0$; if $WL_{1}$, $WR_{1}$ and $WR_{3}$ exist,
they are a shock, a rarefaction wave and a shock, respectively.

###### Lemma 3.3

Under the double CRPs frame, if $M_{4}>M_{5}\geq 1$, then $WL_{1},$$WL_{2}$
and $WL_{3}$ do not exist, and then $u_{5}>u_{4}>0$; if $WR_{1}$ and $WR_{3}$
exist, they are a rarefaction wave and a shock, respectively.

Compared with Type 1 and Type 3, Type 2 has one more wave, and the condition
$M_{5}=1$ is used to match the number of conditions at this time. In the
remaindeer of this paper, each wave and each constant region are denoted as in
Figure4.

##  4 Uniqueness of solution

In this section, we will make a preliminary exploration of the uniqueness of
the solution.

A more refined analysis of the exact solution needs to involve the
quantitative relations of the fundamental waves and the heating discontinuity.
We begin by defining five functions to connect the states on the left and
right sides of the shock wave $WL_{1}$, which appears in the solution of Type
1 and Type 2.

$\displaystyle
f_{1}(M_{1},M_{2})\mathop{=}\limits^{def}\frac{(\gamma+1)(M_{1}-M_{2})^{2}}{(\gamma-1)(M_{1}-M_{2})^{2}+2},$
$\displaystyle
f_{2}(M_{1},M_{2})\mathop{=}\limits^{def}\frac{2\gamma(M_{1}-M_{2})^{2}-\gamma+1}{\gamma+1},$
$\displaystyle
f_{3}(M_{1},M_{2})\mathop{=}\limits^{def}\frac{(\gamma-1)(M_{1}-M_{2})^{2}+2}{(\gamma+1)(M_{1}-M_{2})^{2}},$
$\displaystyle
f_{4}(M_{1},M_{2})\mathop{=}\limits^{def}f_{3}(M_{1},M_{2})+\frac{M_{2}}{M_{1}}(1-f_{3}(M_{1},M_{2})),$
$\displaystyle
f_{5}(M_{1},M_{2})\mathop{=}\limits^{def}\frac{((\gamma-1)M_{1}+2M_{2})(M_{1}-M_{2})+2}{[2\gamma(\gamma-1)(M_{1}-M_{2})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-M_{2})^{2}-2(\gamma-1)]^{\frac{1}{2}}}.$

Applying the Rankine-Hugoniot conditions

$F(U_{4})-F(U_{1})=s_{L}(U_{4}-U_{1}),$

where $s_{L}$ is the speed of $WL_{1}$, we deduce that

$\displaystyle{\rho_{4}}/{\rho_{1}}=f_{1}(M_{1},M_{SL}),\quad$
$\displaystyle{p_{4}}/{p_{1}}=f_{2}(M_{1},M_{SL}),$ (4.1)
$\displaystyle({u_{4}-s_{L}})/({u_{1}-s_{L}})=f_{3}(M_{1},M_{SL}),\quad$
$\displaystyle{u_{4}}/{u_{1}}=f_{4}(M_{1},M_{SL}),$
$\displaystyle{a_{4}}/{a_{1}}=f_{2}(M_{1},M_{SL})^{\frac{1}{2}}f_{1}(M_{1},M_{SL})^{-\frac{1}{2}},\quad$
$\displaystyle M_{4}=f_{5}(M_{1},M_{SL}),$

where $M_{SL}=s_{L}/a_{1}$ is called the shock Mach number of $WL_{1}$.

We now turn to the overall relations of the solution. According to (2.6), the
relations of the stationary discontinuity at the origin are

$\displaystyle\frac{p_{5}}{p_{4}}=\psi(M_{4})=\psi(f_{5}(M_{1},M_{SL})),$
(4.2a)
$\displaystyle\frac{u_{5}}{u_{4}}=\phi(M_{4})=\phi(f_{5}(M_{1},M_{SL})).$
(4.2b)

It is obtained from the shock relation of $WR_{3}$ that

$u_{5}-u_{8}=\sqrt{\frac{\beta
p_{8}}{\rho_{8}}}\frac{p_{5}/p_{8}-1}{\sqrt{1+\tau p_{5}/p_{8}}}$

where $\beta=2/(\gamma-1),\tau=(\gamma+1)/(\gamma-1)$.

Substituting (4.1) and (4.2) into the above equation, we get

$M_{1}(f_{4}(M_{1},M_{SL})\phi(M_{4})-1)-\sqrt{\frac{\beta}{\gamma}}\frac{f_{2}(M_{1},M_{SL})\psi(M_{4})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{4})}}=0.$ (4.3)

The equation (4.1) forms a system of equations for $M_{SL}$ and $M_{4}$ as
follows.

$\left\\{\begin{aligned}
&M_{1}(f_{4}(M_{1},M_{SL})\phi(M_{4})-1)-\sqrt{\frac{\beta}{\gamma}}\frac{f_{2}(M_{1},M_{SL})\psi(M_{4})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{4})}}=0\\\
&M4=\frac{((\gamma-1)M_{1}+2M_{SL})(M_{1}-M_{SL})+2}{\sqrt{2\gamma(\gamma-1)(M_{1}-M_{SL})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-M_{SL})^{2}-2(\gamma-1)}}\end{aligned}\right.$
(4.4)

Another preliminary theorem is to conclude that $M_{SL}$ can be expressed as a
function of $M_{1}$ and $M_{4}$ in the solution of Type 1 or Type 2. From the
Lax entropy condition of the shock (see [19]), we have

$M_{1}-M_{SL}=\frac{u_{1}-s_{L}}{a_{1}}\geq 1.$ (4.5)

###### Lemma 4.1

Given any $M_{1}$ and $M_{4}$, the $M_{SL}$ satisfying
$M_{4}=f_{5}(M_{1},M_{SL})$ and $M_{1}-M_{SL}\geq 1$ exists and is unique.

###### Proof 4.1

The formula $M_{4}=f_{5}(M_{1},M_{SL})$ has another form.

$\displaystyle
M_{4}\sqrt{(\gamma-1)(M_{1}-M_{SL})^{2}+2}\sqrt{2\gamma(M_{1}-M_{SL})^{2}-\gamma+1}$
$\displaystyle=$
$\displaystyle-2(M_{1}-M_{SL})^{2}+(\gamma+1)M_{1}(M_{1}-M_{SL})+2.$

Dividing both sides of the above equation by $M_{1}-M_{SL}$, we have

$\displaystyle
M_{4}\sqrt{(\gamma-1)(M_{1}-M_{SL})^{2}+2}\sqrt{2\gamma-\frac{\gamma+1}{(M_{1}-M_{SL})^{2}}}$
$\displaystyle=$
$\displaystyle-2(M_{1}-M_{SL})+(\gamma+1)M_{1}(M_{1}-M_{SL})+\frac{2}{(M_{1}-M_{SL})^{2}}.$

We define

$\displaystyle\sigma_{1}(x)\mathop{=}\limits^{def}M_{4}\sqrt{(\gamma-1)x^{2}+2}\sqrt{2\gamma-\frac{\gamma+1}{x^{2}}},$
$\displaystyle\sigma_{2}(x)\mathop{=}\limits^{def}2x-(\gamma+1)M_{1}-\frac{2}{x},$
$\displaystyle\sigma(x)\mathop{=}\limits^{def}\sigma_{1}(x)+\sigma_{2}(x).$

The domains of $\sigma_{1}(x)$, $\sigma_{2}(x)$ and $\sigma(x)$ are
$\\{x|x\geq 1\\}$. Both $\sigma_{1}(x)$ and $\sigma_{2}(x)$ are monotone
increasing functions, hence $\sigma(x)$ is a monotone increasing function. It
is clear that

$\sigma(1)=(\gamma+1)(M_{4}-M_{1})\leq 0,$
$\lim_{x\to\infty}\sigma(x)=\lim_{x\to\infty}\sigma_{1}(x)+\lim_{x\to\infty}\sigma_{2}(x)=\infty.$

Consequently, the root of $\sigma(x)$ exists and is unique. $M_{1}-M_{SL}$ is
the root of $\sigma(x)$, and this completes the proof.

For the solutions of Type 1 and Type 2, it clear that

$M_{4}=\frac{((\gamma-1)M_{1}+2M_{SL})(M_{1}-M_{SL})+2}{\sqrt{2\gamma(\gamma-1)(M_{1}-M_{SL})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-M_{SL})^{2}-2(\gamma-1)}}.$

Treating $M_{1}-M_{SL}$ as the variable, we have

$A(M_{1}-M_{SL})^{4}+B(M_{1}-M_{SL})^{3}+C(M_{1}-M_{SL})^{2}+D(M_{1}-M_{SL})+E=0,$
(4.6)

where

$\displaystyle A=4-2\gamma(\gamma-1)M_{4}^{2},$ $\displaystyle
B=-4(\gamma+1)M_{1},$ $\displaystyle
C=(\gamma+1)^{2}M_{1}^{2}-8-(6\gamma-\gamma^{2}-1)M_{4}^{2},$ $\displaystyle
D=4(\gamma+1)M_{1},$ $\displaystyle E=4+2(\gamma-1)M_{4}^{2}.$

$M_{1}-M_{SL}$ is the root of a fourth-order polynomial function (4.6).
According to Lemma4.1, We have proved the following theorem.

###### Theorem 4.1

For the solutions of Type 1 and Type 2, $M_{SL}$ can be expressed as
$M_{SL}=f_{6}(M_{1},M_{4})$, which satisfies $M_{4}=f_{5}(M_{1},M_{SL})$ and
$M_{1}-M_{SL}\geq 1$.

Now we turn to the uniqueness of the solution. The proof consists of two
parts, one involes the states at the intermediate regions (Lemma4.2), and the
other is a division of these three structures (Lemma4.3). Both of these are
related to the Mach number $M_{1}$ of the initial flow.

We describe the strength of each nonlinear wave in terms of the ratio of
pressures on both sides of this wave.

###### Lemma 4.2

For these three types of solutions, the strength of each nonlinear wave is
determined by the Mach number $M_{1}$, independent of other variables.

A proof of this lemma will be given in AppendixProof of Lemma4.2.

We define

$X(M_{1},M_{4})\mathop{=}\limits^{def}M_{1}(f_{4}(M_{1},M_{SL})\phi(M_{4})-1)-\sqrt{\frac{\beta}{\gamma}}\frac{f_{2}(M_{1},M_{SL})\psi(M_{4})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{4})}},$ (4.7)

and

$Y(M_{1})\mathop{=}\limits^{def}X(M_{1},M_{*}),$ (4.8)

where $M_{SL}=f_{6}(M_{1},M_{4})$.

If $k(\gamma^{2}-1)\geq 1$, then $M_{**}$ does not exist. Therefore the
upstream flow of the heating point must be a subsonic flow. At this time,
there are only two possible structures: Type 1 and Type 2. In this case, the
downstream fluid must be subsonic. When $k(\gamma^{2}-1)<1$, all three
structures are possible. We associate the three structures with $M_{1}$ by the
following lemma.

###### Lemma 4.3

Under the double CRPs frame, the structures of the exact solution can only be
Type 1 or Type 2 for $k(\gamma^{2}-1)\geq 1$, and all three structures are
possible for $k(\gamma^{2}-1)<1$. This type of solution satisfies the
following conditions.

  * (i)

For the solution of Type 1, $Y(M_{1})\geq 0$ holds;

  * (ii)

For the solution of Type 2, $Y(M_{1})\leq 0$ and $M_{1}\geq M_{*}$ hold, and
$M_{1}\leq M_{**}$ holds at the case of $k(\gamma^{2}-1)<1$;

  * (iii)

For the solution of Type 3, $M_{1}\geq M_{**}$ holds.

A proof of this lemma will be given in AppendixProof of Lemma4.3.

It can be seen from Lemma4.3 that the structure of $M_{1}$ equaling to the
root of $Y(M_{1})$ is the demarcation structure of Type 1 and Type 2, which is

$S(U_{1},U_{4})\oplus H(U_{4},U_{5})\oplus C(U_{5},U_{7})\oplus
S(U_{7},U_{8})\quad\text{with}\quad M_{4}=M_{*}.$

And the structure of $M_{1}$ equaling to $M_{**}$ is the demarcation structure
of Type 2 and Type 3, which is

$S(U_{1},U_{4})\oplus H(U_{4},U_{5})\oplus R(U_{5},U_{6})\oplus
C(U_{6},U_{7})\oplus S(U_{7},U_{8})\quad\text{with}\quad s_{L}=0,$ (4.9)

or

$H(U_{1},U_{5})\oplus R(U_{5},U_{6})\oplus C(U_{6},U_{7})\oplus
S(U_{7},U_{8})\quad\text{with}\quad M_{5}=1.$ (4.10)

###### Remark 4.1

The structure of (4.9) is a limit structure of Type 2. $WL_{1}$ is a normal
shock and $M_{4}=M_{*}$ hold for this structure, and it is clear that
$M_{1}=M_{**}$. For the structure of (4.10), $M_{5}=1$ implies $M_{1}=M_{**}$.
Therefore the structure of (4.9) equals to the structure of (4.10).

According to Lemma4.2 and Lemma4.3, we can establish the following theorem.

###### Theorem 4.2

Under the double CRPs frame, we have

  * (i)

if $k(\gamma^{2}-1)\geq 1$, the solution of Riemann problem (1.1) and (1.2) is
unique;

  * (ii)

if $k(\gamma^{2}-1)<1$, the solution of Riemann problem (1.1) and (1.2) is
unique under the assumption that the root of $Y(M_{1})$ is not greater than
$M_{**}$.

If the root of $Y(M_{1})$ is greater than $M_{**}$, the uniqueness of the
solution has not been proved. To compare the size of the root of $Y(M1)$ and
$M_{**}$, we define

$R(\gamma,k)=\\{M|Y(M)=0\\},\quad T(\gamma,k)=M_{**}-R(\gamma,k).$

Figure 7: The contour of function $T(\gamma,k)$

Figure7 shows the contour of the function $T$. The curved boundary on the
upper right is the curve of $k(\gamma^{2}-1)=1$. From Figure7 it can be found
$T(\gamma,k)>0$, therefore the root of $Y(M_{1})$ is smaller than $M_{**}$ and
the assumption in the second part of Theorem4.2 is true. Thus we give the
following conjecture.

###### Conjecture 4.1

Under the double CRPs frame, the solution of Riemann problem (1.1) and (1.2)
is unique for any given $\gamma$ and $\kappa$, and the structure of the self-
similar solution is determined by $M_{1}$, as shown below.

  * (i)

If $Y(M_{1})\geq 0$, the structure of the solution is Type 1.

  * (ii)

If $Y(M_{1})\geq 0$ and $M_{1}\leq M_{**}(k(\gamma^{2}-1)<1)$, the structure
of the solution is Type 2.

  * (iii)

If $k(\gamma^{2}-1)<1$ and $M_{1}\geq M_{**}$, the structure of the solution
is Type 3.

##  5 Algorithm of Solution and Verification

Algorithm 1 The Construction Algorithm of the Exact Solution

the initial condition: $\rho_{1}$, $u_{1}$, $p_{1}$, $M_{1}$, $a_{1}$;

the intermediate states: $W_{4}=(\rho_{4},u_{4},p_{4})$,
$W_{5}=(\rho_{5},u_{5},p_{5})$, $W_{6}=(\rho_{6},u_{6},p_{6})$,
$W_{7}=(\rho_{7},u_{7},p_{7})$;

compute $M_{*}$ and $M_{**}$;

compute $Y(M_{1})$;

if $M_{1}\geq M_{**}$ then

$W_{4}\leftarrow W_{1}$;

compute $W_{5}$ by (2.2);

else if $Y(M_{1})\geq 0$ then

set the initial values of the bisection method: $x_{1}$, $x_{2}$;

$x_{0}\leftarrow(x_{1}+x_{2})/2$;

set the value of termination condition: $error$;

while $|x_{1}-x_{2}|>error$ do

$f_{1}\leftarrow((\gamma-1)M_{1}+2x_{1})(M_{1}-x_{1})+2-M_{*}(2\gamma(\gamma-1)(M_{1}-x_{1})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-x_{1})^{2}-2(\gamma-1))^{\frac{1}{2}}$;

$f_{2}\leftarrow((\gamma-1)M_{1}+2x_{2})(M_{1}-x_{1})+2-M_{*}(2\gamma(\gamma-1)(M_{1}-x_{2})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-x_{2})^{2}-2(\gamma-1))^{\frac{1}{2}}$;

if $f_{1}\times f_{2}<0$ then

$x_{0}\leftarrow x_{1}$;

else

$x_{1}\leftarrow x_{0}$;

end if

$x_{0}\leftarrow(x_{1}+x_{2})/2$;

end while

$\rho_{4}\leftarrow\rho_{1}f_{1}(M_{1},x_{0})$, $u_{4}\leftarrow
u_{1}f_{4}(M_{1},x_{0})$, $p_{4}\leftarrow p_{1}f_{2}(M_{1},x_{0})$,
$M_{4}\leftarrow M_{*}$;

compute $W_{5}$ by (2.1);

$p_{5}\leftarrow p_{4}\psi(M_{*})$, $u_{5}\leftarrow u_{4}\phi(M_{*})$,
$\rho_{5}\leftarrow(\rho_{4}u_{4})/u_{5}$;

else

//$f_{L}$ and $f_{R}$ are given by (5.2).

set the initial values of the bisection method: $y_{1}$, $y_{2}$;

$y_{0}\leftarrow(y_{1}+y_{2})/2$;

set the value of termination condition: $error$;

while $|y_{1}-y_{2}|>error$ do

$m_{1}\leftarrow(u_{1}-f_{L}(y_{1}))\left[\frac{\rho_{1}}{\gamma}\frac{\gamma+1+(\gamma-1)p_{1}/y_{1}}{(\gamma-1)y_{1}+(\gamma+1)p_{1}}\right]^{\frac{1}{2}}$,
$m_{2}\leftarrow(u_{1}-f_{L}(y_{2}))\left[\frac{\rho_{1}}{\gamma}\frac{\gamma+1+(\gamma-1)p_{1}/y_{2}}{(\gamma-1)y_{2}+(\gamma+1)p_{1}}\right]^{\frac{1}{2}}$;

$g_{1}\leftarrow
u_{1}+f_{R}(y_{1}\psi(m_{1}))-\phi(m_{1})(u_{1}-f_{L}(y_{1}))$,
$g_{2}\leftarrow
u_{1}+f_{R}(y_{2}\psi(m_{2}))-\phi(m_{1})(u_{1}-f_{L}(y_{2}))$;

if $g_{1}\times g_{2}<0$ then

$y_{0}\leftarrow y_{1}$;

else

$y_{1}\leftarrow y_{0}$;

end if

$y_{0}\leftarrow(y_{1}+y_{2})/2$;

end while

$\rho_{4}\leftarrow\rho_{1}\frac{(\gamma-1)p_{1}+(\gamma+1)y_{0}}{(\gamma-1)y_{0}+(\gamma+1)p_{1}}$,
$u_{4}\leftarrow u_{1}-f_{L}(y_{0})$, $p_{4}\leftarrow y_{0}$;

$M_{1}\leftarrow(u_{1}-f_{L}(y_{0}))\left[\frac{\rho_{1}}{\gamma}\frac{\gamma+1+(\gamma-1)p_{1}/y_{0}}{(\gamma-1)y_{0}+(\gamma+1)p_{1}}\right]^{\frac{1}{2}}$;

compute $W_{5}$ by (2.1);

end if

compute $W_{6}$ and $W_{7}$ by a Riemann solver $CRP(W_{5},W_{8})$.

Lemma4.3 makes it legitimate to apply $M_{1}$ to determine the structure of
the exact solution. Once the structure is determined, the solution becomes
easy. A algorithm for the solution is given in Algorithm1.

In the algorithm of Type 2, we use the bisection method to solve

$M_{*}=\frac{((\gamma-1)M_{1}+2M_{SL})(M_{1}-M_{SL})+2}{\sqrt{2\gamma(\gamma-1)(M_{1}-M_{SL})^{4}+(6\gamma-\gamma^{2}-1)(M_{1}-M_{SL})^{2}-2(\gamma-1)}}.$
(5.1)

In the algorithm of Type 1, we do not iteratively solve 4.4, but use the
pressure equations

$u_{1}+f_{R}(p_{4}\psi(M_{4}(p_{4})))-\phi(M_{4}(p_{4}))(u_{1}-f_{L}(p_{4}))=0,$

where

$M_{4}(p_{4})=[u_{1}-f_{L}(p_{1})]\left[\frac{\rho_{1}}{\gamma}\frac{\gamma+1+(\gamma-1)p_{1}/p_{4}}{(\gamma-1)p_{4}+(\gamma+1)p_{1}}\right]^{\frac{1}{2}},$
$\displaystyle
f_{L}(p)=\begin{cases}(p-p_{1})\left[\frac{A_{L}}{p+B_{L}}\right]^{\frac{1}{2}},\quad
if\ p>p_{1}\\\
\frac{2a_{1}}{\gamma-1}\left[\left(\frac{p}{p_{1}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right],\quad
if\ p\leq p_{1}\end{cases},$ (5.2) $\displaystyle
f_{R}(p)=\begin{cases}(p-p_{1})\left[\frac{A_{R}}{p+B_{R}}\right]^{\frac{1}{2}},\quad
if\ p>p_{1}\\\
\frac{2a_{1}}{\gamma-1}\left[\left(\frac{p}{p_{1}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right],\quad
if\ p\leq p_{1}\end{cases}.$

$A_{L},A_{R},B_{L},B_{R}$ are given by

$A_{L}=\frac{2}{(\gamma+1)\rho_{1}},B_{L}=\frac{\gamma-1}{\gamma+1}p_{1},A_{R}=\frac{2}{(\gamma+1)\rho_{1}},B_{R}=\frac{\gamma-1}{\gamma+1}p_{1}.$

Experiments show that this iterative method is less sensitive to the initial
value of the iteration. Lemma4.1 guarantees that the solution of the iterative
equation in the algorithm is unique. A classical Riemann solver
$CRP(W_{5},W_{8})$ of the Euler equations is needed in the algorithm. The wave
pattern is identical for this CRP, which consists a rarefaction wave
corresponding to the $u-a$ characteristic field and a shock wave corresponding
to the $u+a$ charatristic field.

We verify the existence of these three structures through numerical tests, and
compare the states of the constructed self-similar solution with the numerical
solution at intermediate regions. The following five tests are all for the
ideal gase with $\gamma=1.4$, and their exact solutions cover three proposed
structures. The initial conditions and heating parameters of the five tests
are shown in Table 1.

Table 1: Initial conditions and parameter setting in experiments | initial conditions | heating parameter | solution structure
---|---|---|---
| $\rho$ | $u$ | $p$
Test1 | 1.0 | 0.8 | 1.0 | 0.2 | Type 1
Test2 | 1.0 | 1.2 | 1.0 | 0.2 | Type 1
Test3 | 1.0 | 1.8 | 1.0 | 0.2 | Type 2
Test4 | 1.0 | 2.8 | 1.0 | 0.2 | Type 3
Test5 | 1.0 | 2.8 | 1.0 | 2.0 | Type 2

The comparison between the numerical solutions and the exact solutions are
shown in
FigureLABEL:figure_test1,figure_test2,figure_test3,figure_test4,figure_test5.
In the first four tests, the heating parameter is $k=0.2$, and
$k(\gamma^{2}-1)<1$ holds. $M_{*}$ and $M_{**}$ are 0.6136 and 1.8130,
respectively. The root of $Y(M_{1})$ is 1.0620, which is less than $M_{**}$,
hence the solution is unique. In the last test, the value of $k$ is 2.0.
$k(\gamma^{2}-1)>1$ holds at this time. According to Theorem4.2, the exact
solutions of these five Riemann problem are all unique.

(a) density (b) velocity (c) pressure

Figure 8: The numerical solution obtained by the RKDG and comparion to the
constructed self-similar solution for Test1 at $t=4.5s$.

(a) density (b) velocity (c) pressure

Figure 9: The numerical solution obtained by the RKDG and comparion to the
constructed self-similar solution for Test2 at $t=4.5s$.

(a) density (b) velocity (c) pressure

Figure 10: The numerical solution obtained by the RKDG and comparion to the
constructed self-similar solution for Test3 at $t=4.5s$.

(a) density (b) velocity (c) pressure

Figure 11: The numerical solution obtained by the RKDG and comparion to the
constructed self-similar solution for Test4 at $t=2.5s$.

(a) density (b) velocity (c) pressure

Figure 12: The numerical solution obtained by the RKDG and comparion to the
constructed self-similar solution for Test5 at $t=2.5s$.

We apply the Runge-Kutta discontinuous Galerkin method (RKDG) for numerical
simulation (see [26]). In the spatial direction, the solution is discretized
by piece-wise second-order ($P2$) polynomials. A total variation diminishing
(TVD) limiter is empolyed to avoid numerical oscillations. In the time
direction, we use the third-order Runge-Kutta method. The source is processed
through splitting. Note that the source term should be precessed at each time
step of Runge-Kutta method.

It can be observed from figures that the constructed exact solution and
numerical solution fit well in each intermediate regions. While inside several
cells near origin, the solutions of pressure are some different between the
exact solution and the numerical solution. The reason is that directly
splitting the source term has no ”well-balanced” property.

##  6 Conclusions

This paper focused on the Riemann problem of the Euler equations with a Dirac
delta-source in the energy conservation eqution. The double CRPs frame was
proposed to construct the self-similar solutions. We proved that there are
three types of the solution and studied the uniqueness of the solutions under
this frame. The present frame is completely different from the existing method
of dealing with the Riemann problem with discontinuous source. Based on the
solution of CRP, the strategy of the present frame is the elimination of
unconscionable waves on the general structure (Figure3). This frame can be
applied to the Riemann problem of other hyperbolic systems with Dirac delta-
function sources or other sources. Compared with some existing methods, it is
more simple and can naturally cover all possible structures. We verified the
double CRPs frame by comparing the constructed self-similar solution with the
numerical solution obtained by RKDG. Besides, the constructed solutions can be
used to evaluate the existing numerical methods for source terms, such as [11,
12, 13].

The uniqueness of the self-similar solutions under the double CRPs frame for
arbitrary initial conditions is an open question, which is the goal of our
future work. In addition, future work should focous on the heating addition of
unsteady flow, in which the transition of the three structures proposed in
this paper may occur.

## Appendix

### Proof of Lemma4.2

###### Proof 6.1

For the solution of Type 1, (4.1) and (4.2) imply

$\displaystyle p_{4}/p_{1}=f_{2}(M_{1},M_{SL}),$ $\displaystyle
p_{8}/p_{7}=p_{1}/p_{4}\times
p_{4}/p_{5}=\left[f_{2}(M_{1},M_{SL}\psi(f_{5}(M_{1},M_{SL})))\right]^{-1},$

where $M_{SL}$ is obtained from equation (4.4) and is determined by $M_{1}$.

For Structure 2, we have

$p_{4}/p_{1}=f_{2}(M_{1},M_{SL}),$

where $M_{SL}$ is obtained from equation (5.1) and is determined by $M_{1}$.

From the elementary wave equation in [23], we have

$\displaystyle
u_{6}=u_{5}-\frac{2a_{5}}{\gamma-1}\left[\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right],$
$\displaystyle
u_{7}=u_{8}+\sqrt{2}(p_{7}-p_{8})\left[(\gamma+1)p_{7}\rho_{8}+(\gamma-1)p_{8}\rho_{8}\right]^{-\frac{1}{2}}.$

According to $u_{6}=u_{7}$, $p_{6}=p_{7}$ and $U_{1}=U_{8}$, it follows that

$u_{5}-\frac{2a_{5}}{\gamma-1}\left[\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right]=u_{1}+\sqrt{2}(p_{6}-p_{1})\left[(\gamma+1)p_{6}\rho_{1}+(\gamma-1)p_{1}\rho_{1}\right]^{-\frac{1}{2}}.$

Divide both sides of the above equation by $a_{1}$, then

$\frac{u_{5}}{u_{1}}M_{1}-\frac{2}{\gamma-1}\frac{M_{1}}{M_{5}}\frac{u_{5}}{u_{1}}\left[\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right]=M_{1}+\frac{\sqrt{2}}{a_{1}}(p_{6}-p_{1})\left[(\gamma+1)p_{6}\rho_{1}+(\gamma-1)p_{1}\rho_{1}\right]^{-\frac{1}{2}}.$

Apply the equation of state to the above equation, we have

$\frac{u_{5}}{u_{1}}M_{1}-\frac{2}{\gamma-1}\frac{M_{1}}{M_{5}}\frac{u_{5}}{u_{1}}\left[\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}-1\right]=M_{1}+\sqrt{\frac{2}{\gamma}}(\frac{p_{6}}{p_{1}}-1)\left[(\gamma+1)\frac{p_{6}}{p_{1}}+(\gamma-1)\right]^{-\frac{1}{2}}.$

Therefore

$\left(A\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}+B\right)\sqrt{C\frac{p_{6}}{p_{5}}+\gamma-1}-\sqrt{\frac{2}{\gamma}}\left(\frac{C}{\gamma+1}\frac{p_{6}}{p_{5}}-1\right)=0,$

where $A=-\frac{2}{\gamma-1}M_{1}\phi(M_{*})f_{4}(M_{1},M_{SL})$,
$B=\frac{\gamma+1}{\gamma-1}M_{1}\phi(M_{*})f_{4}(M_{1},M_{SL})$,
$C=(\gamma+1)\psi(M_{*})f_{2}(M_{1},M_{SL})$ and $M_{SL}$ is obtained from
(5.1). Thus the strength of $WR_{1}$ in the solution of Type 2 is determined
by $M_{1}$.

For $WR_{3}$, we have

$\frac{p_{8}}{p_{7}}=\left[f_{2}(M_{1},M_{SL})\psi(M_{*})\frac{p_{6}}{p_{5}}\right]^{-1}.$

Therefore the strength of each nonlinear wave in the solution of Type 2 is
determined by $M_{1}$.

The analyse of Type 3 is similar, as follows.

$\displaystyle\left(A^{\prime}\left(\frac{p_{6}}{p_{5}}\right)^{\frac{\gamma-1}{2\gamma}}+B^{\prime}\right)\sqrt{C^{\prime}\frac{p_{6}}{p_{5}}+\gamma-1}-\sqrt{\frac{2}{\gamma}}\left(\frac{C^{\prime}}{\gamma+1}\frac{p_{6}}{p_{5}}-1\right)=0,$
$\displaystyle\frac{p_{8}}{p_{7}}=\left[\frac{p_{5}}{p_{1}}\frac{p_{6}}{p_{5}}\right]^{-1},$

where $A^{\prime}=-\frac{2}{\gamma-1}\frac{M_{1}}{M_{5}}\frac{u_{5}}{u_{1}}$,
$B^{\prime}=\frac{u_{5}}{u_{1}}M_{1}+\frac{2}{\gamma-1}\frac{M_{1}}{M_{5}}\frac{u_{5}}{u_{1}}$,
$C^{\prime}=(\gamma+1)\frac{p_{5}}{p_{1}}$. $\frac{u_{5}}{u_{1}}$,
$\frac{p_{5}}{p_{1}}$ and $M_{5}$ in the above equation are obtained from
Lemma3.1 and they are all detrmined by $M_{1}$.

### Proof of Lemma4.3

###### Proof 6.2

From (4.3), we have

$X(M_{1},M_{4})=0.$

Some tedious manipulation yields

$\frac{\partial f_{5}}{\partial M_{SL}}>0,$

and

$\frac{\partial f_{6}}{\partial M_{SL}}=\left(\frac{\partial f_{5}}{\partial
M_{SL}}\right)^{-1}>0.$

Then

$\frac{\partial X}{\partial M_{4}}=M_{1}\frac{\partial}{\partial
M_{4}}(f_{4}(M_{1},M_{SL})\phi(M_{4}))-\sqrt{\frac{\beta}{\gamma}}\frac{\partial}{\partial
M_{4}}\left(\frac{f_{2}(M_{1},M_{SL})\psi(M_{4})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{4})}}\right).$ (6.1)

For the first part on the right of (6.1), we have

$\frac{\partial}{\partial
M_{4}}(f_{4}(M_{1},M_{SL})\phi(M_{4}))=\frac{\partial
f_{4}(M_{1},M_{4})}{\partial
M_{4}}\phi(M_{4})+\frac{d\phi(M_{4})}{dM_{4}}f_{4}(M_{1},M_{SL}),$
$\frac{\partial f_{4}(M_{1},M_{4})}{\partial M_{4}}=\frac{\partial
f_{4}(M_{1},M_{4})}{\partial M_{SL}}\frac{\partial M_{SL}}{\partial
M_{4}}=\frac{2(1+(M_{1}-M_{SL})^{-2})}{(\gamma+1)M_{1}}\frac{\partial
f_{6}}{\partial M_{4}}>0,$ $\phi^{\prime}(M_{4})>0.$

Thus

$\frac{\partial}{\partial M_{4}}(f_{4}(M_{1},M_{SL})\phi(M_{4}))>0.$

For the second part on the right of (6.1), we have

$\frac{\partial}{\partial
M_{4}}\left(\frac{f_{2}(M_{1},M_{SL})\psi(M_{4})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{4})}}\right)=\frac{\tau(f_{2}(M_{1},M_{SL})\psi(M_{4})+1)+2}{2(f_{2}(M_{1},M_{SL})\psi(M_{4}))^{3/2}}\frac{\partial}{\partial
M_{4}}(f_{2}(M_{1},M_{SL})\psi(M_{4})),$ $\frac{\partial}{\partial
M_{4}}(f_{2}(M_{1},M_{SL})\psi(M_{4}))=-\left(\frac{\partial f_{2}}{\partial
M_{1}}-\frac{\partial f_{2}}{\partial
M_{SL}}\right)\psi(M_{4})+\frac{d\psi(M_{4})}{dM_{SL}}f_{2}(M_{1},M_{SL}),$
$\frac{\partial f_{2}}{\partial M_{1}}-\frac{\partial f_{2}}{\partial
M_{SL}}>0,$ $\psi^{\prime}(M_{4})<0.$

Thus

$\frac{\partial}{\partial M_{4}}(f_{2}(M_{1},M_{SL})\psi(M_{4}))<0,$

and

$\frac{\partial}{\partial
M_{4}}\left(\frac{f_{2}(M_{1},M_{SL})g_{2}(M_{1},M_{SL})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})g_{2}(M_{1},M_{SL})}}\right)<0.$

Note that we have actually proved that

$\frac{\partial X}{\partial M_{4}}>0.$

According to Lemma2.1, it follows that $M_{4}<M_{*}$. Consequently

$Y(M_{1})=X(M_{1},M_{*})\geq X(M_{1},M_{4})=0.$

###### Proof 6.3

The solution right to the t-axis of Type 2 is the solution of
$CRP(U_{5},U_{8})$, whose wave pattern is $R(U_{5},U_{6})\oplus
C(U_{6},U_{7})\oplus S(U_{7},U_{8})$, it follows that (see [16])

$u_{5}-u_{8}\leq\sqrt{\frac{\beta
p_{8}}{\rho_{8}}}\frac{p_{5}/p_{8}-1}{\sqrt{1+\tau p_{5}/p_{8}}}.$

Represent the above equation by the Mach numbers, as follows.

$M_{1}(f_{4}(M_{1},M_{SL})\phi(M_{*})-1)-\sqrt{\frac{\beta}{\gamma}}\frac{f_{2}(M_{1},M_{SL})\psi(M_{*})-1}{\sqrt{1+\tau
f_{2}(M_{1},M_{SL})\psi(M_{*})}}\leq 0.$

Thus

$Y(M_{1})\leq 0.$

We denote

$M_{SL}^{\prime}=\frac{s_{L}}{a_{4}},$

where $s_{L}$ is the speed of $WL_{1}$ and $a_{4}$ is the speed of sound in
region 4. Similar to (4.1), it holds that

$M_{1}=\frac{((\gamma-1)M_{4}+2M_{SL}^{\prime})(M_{4}-M_{SL}^{\prime})+2}{\sqrt{2\gamma(\gamma-1)(M_{4}-M_{SL}^{\prime})^{4}+(6\gamma-\gamma^{2}-1)(M_{4}-M_{SL}^{\prime})^{2}-2(\gamma-1)}}.$

Define

$g(M_{4},M_{SL}^{\prime})=\frac{((\gamma-1)M_{4}+2M_{SL}^{\prime})(M_{4}-M_{SL}^{\prime})+2}{\sqrt{2\gamma(M_{4}-M_{SL}^{\prime})^{2}-\gamma+1}\sqrt{(\gamma-1)(M_{4}-M_{SL}^{\prime})^{2}+2}}.$

The domain of $g$ satisfies

$M_{4}\leq M_{*},\quad M_{SL}^{\prime}\leq 0.$

A routine computation gives rise to

$\frac{\partial g}{\partial M_{SL}^{\prime}}>0.$

Then

$g(M_{4},M_{SL}^{\prime})\leq
g(M_{4},0)=\sqrt{\frac{(\gamma-1)M_{4}^{2}+2}{2\gamma M_{4}^{2}-\gamma+1}},$

and

$M_{1}=g(M_{*},M_{SL}^{\prime})\leq
g(M_{*},0)=\sqrt{\frac{(\gamma-1)M_{*}^{2}+2}{2\gamma M_{*}^{2}-\gamma+1}}.$

According to Lemma2.2, it follows that

$M_{1}\leq M_{**}.$

The proof of (iii) is obvious.

## Acknowledgments

This work was supported by the NSFC-NSAF joint fund [No. U1730118]; and the
Science Challenge Project [No. JCKY2016212A502].

## References

  * [1] Rémi Abgrall and Smadar Karni.  A comment on the computation of non-conservative products.  Journal of Computational Physics, 229(8):2759–2763, 2010.
  * [2] Francisco Alcrudo and Fayssal Benkhaldoun.  Exact solutions to the riemann problem of the shallow water equations with a bottom step.  Computers & Fluids, 30(6):643 – 671, 2001.
  * [3] Matania Benartzi and Joseph Falcovitz.  A second-order godunov-type scheme for compressible fluid dynamics.  Journal of Computational Physics, 55(1):1–32, 1984.
  * [4] Matania Benartzi, Jiequan Li, and Gerald Warnecke.  A direct eulerian grp scheme for compressible fluid flows.  Journal of Computational Physics, 218(1):19–43, 2006.
  * [5] R Bernetti, V A Titarev, and Eleuterio F Toro.  Exact solution of the riemann problem for the shallow water equations with discontinuous bottom geometry.  Journal of Computational Physics, 227(6):3212–3243, 2008.
  * [6] Wan Cheng, Xisheng Luo, and Van Meh Rini Dongen.  On condensation-induced waves.  Journal of Fluid Mechanics, 651(1):145–164, 2010.
  * [7] Can F Delale, G H Schnerr, and Jurgen Zierep.  The mathematical theory of thermal choking in nozzle flows.  Zeitschrift für Angewandte Mathematik und Physik, 44(6):943–976, 1993.
  * [8] Can F. Delale, Günter H. Schnerr, and Marinus E. H. Van Dongen.  Condensation Discontinuities and Condensation Induced Shock Waves.  2007\.
  * [9] M. E. H. Van Dongen, X. Luo, G. Lamanna, and D. J. Van Kaathoven.  On condensation induced shock waves.  In Proc. 10th Chinese Symposium on Shock Waves, 2002.
  * [10] Laurent Gosse.  A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms.  Mathematical Models and Methods in Applied Sciences, 11(02):339–365, 2001.
  * [11] J M Greenberg, A Y Leroux, R Baraille, and A Noussair.  Analysis and approximation of conservation laws with source terms.  SIAM Journal on Numerical Analysis, 34(5):1980–2007, 1997.
  * [12] Shi Jin and Xin Wen.  Two interface-type numerical methods for computing hyperbolic systems with geometrical source terms having concentrations.  SIAM Journal on Scientific Computing, 26(6):2079–2101, 2005.
  * [13] Dietmar Kroner and Mai Duc Thanh.  Numerical solutions to compressible flows in a nozzle with variable cross-section.  SIAM Journal on Numerical Analysis, 43(2):796–824, 2005.
  * [14] Philippe G Lefloch and Mai Duc Thanh.  The riemann problem for fluid flows in a nozzle with discontinuous cross-section.  Communications in Mathematical Sciences, 1(4):763–797, 2003.
  * [15] Philippe G Lefloch and Mai Duc Thanh.  The riemann problem for the shallow water equations with discontinuous topography.  Communications in Mathematical Sciences, 5(4):865–885, 2007.
  * [16] T G Liu, Boo Cheong Khoo, and C W Wang.  The ghost fluid method for compressible gas-water simulation.  Journal of Computational Physics, 204(1):193–221, 2005.
  * [17] Xisheng Luo, B Bart Prast, Van Meh Rini Dongen, Hwm Harrie Hoeijmakers, and J Yang.  On phase transition in compressible flows: modelling and validation.  Journal of Fluid Mechanics, 548(1):403–430, 2006.
  * [18] Carlos Pares and Ernesto Pimentel.  The riemann problem for the shallow water equations with discontinuous topography: The wet–dry case.  Journal of Computational Physics, 378:344–365, 2019.
  * [19] D Peter.  Hyperbolic systems of conservation laws and the mathematical theory of shock waves /.  Society for Industrial and Applied Mathematics,, 1973.
  * [20] Gunter Schnerr.  Unsteadiness in condensing flow: Dynamics of internal flows with phase transition and application to turbomachinery.  Journal of Mechanical Engineering Science, 219, 2005.
  * [21] Mai Duc Thanh.  The riemann problem for a nonisentropic fluid in a nozzle with discontinuous cross-sectional area.  Siam Journal on Applied Mathematics, 69(6):1501–1519, 2009.
  * [22] Mai Duc Thanh.  Numerical treatment in resonant regime for shallow water equations with discontinuous topography.  Communications in Nonlinear Science and Numerical Simulation, 18(2):417–433, 2013.
  * [23] Eleuterio F. Toro.  Riemann solvers and numerical methods for fluid dynamics : a practical introduction.  Springer,.
  * [24] Eleuterio F Toro and Arturo Hidalgo.  Ader finite volume schemes for nonlinear reaction–diffusion equations.  Applied Numerical Mathematics, 59(1):73–100, 2009.
  * [25] Eleuterio F Toro and Gino I Montecinos.  Implicit, semi-analytical solution of the generalized riemann problem for stiff hyperbolic balance laws.  Journal of Computational Physics, 303:146–172, 2015.
  * [26] Yang Yang and Chi-Wang Shu.  Discontinuous galerkin method for hyperbolic equations involving $\delta$-singularities: negative-order norm error estimates and applications.  Numerische Mathematik, 124(4):753–781, 2013.


Principle of structural analogy of solutions and

its application to nonlinear PDEs and delay PDEs

Andrei D. Polyanin

Ishlinsky Institute for Problems in Mechanics, Russian Academy of Sciences,

pr. Vernadskogo 101, bldg. 1, Moscow, 119526 Russia

e-mail<EMAIL_ADDRESS>

Using the principle of structural analogy of solutions, approaches have been
developed for constructing exact solutions of complex nonlinear PDEs,
including PDEs with delay, based on the use of special solutions to auxiliary
simpler related equations. It is shown that to obtain exact solutions of
nonlinear non-autonomous PDEs, the coefficients of which depend on time, it is
possible to use generalized and functional separable solutions of simpler
autonomous PDEs, the coefficients of which do not depend on time. Specific
examples of constructing exact solutions to nonlinear PDEs, the coefficients
of which depend arbitrarily on time, are considered. It has been discovered
that generalized and functional separable solutions of nonlinear PDEs with
constant delay can be used to construct exact solutions of more complex
nonlinear PDEs with variable delay of general form. A number of nonlinear
reaction-diffusion type PDEs with variable delay are described, which allow
exact solutions with generalized separation of variables.

Keywords: nonlinear PDEs, delay PDEs, reaction-diffusion equations, solution
methods, exact solutions, generalized separable solutions, functional
separable solutions

Mathematics Subject Classification: 35C05, 35K55, 35K57, 35R10

##  1 Introduction. Exact solutions

Exact solutions of mathematical physics equations and other partial
differential equations contribute to a better understanding of qualitative
features of many nonlinear phenomena and processes in various fields of
natural science. Since exact solutions are rigorous mathematical standards,
they can be used as test problems to assess the accuracy and verify the
adequacy of various numerical and approximate analytical methods for solving
nonlinear PDEs as well as more complex PDEs with delay.

Exact solutions of nonlinear PDEs are usually understood as solutions that are
expressed [1]:

(i) through elementary functions or in the form of quadratures (i.e., using
elementary functions and indefinite integrals);

(ii) through solutions of ordinary differential equations or systems of such
equations.

For more complex PDEs with delay to exact solutions from Items (i) and (ii) it
is necessary to add solutions that are expressed through solutions of ODEs
with delay or solutions of ODE systems with delay [2].

Exact solutions of nonlinear PDEs are most often constructed using the
classical method of symmetry reductions [3, 4, 5, 6], the direct method of
symmetry reductions [1, 7, 8, 9], the nonclassical symmetries methods [10, 11,
12, 13], methods of generalized separation of variables [1, 9, 14, 15, 16],
methods of functional separation of variables [1, 9, 17, 18], the method of
differential constraints [1, 9, 19, 20], and some other exact analytical
methods (see, for example, [21, 22, 23, 24, 25, 26, 27]). On methods for
constructing exact solutions of nonlinear delay PDEs and functional PDEs, see,
for example, [2, 28, 29, 30, 31, 32, 33, 34].

This article will describe new approaches for constructing exact solutions of
nonlinear non-autonomous equations of mathematical physics and other PDEs
(including PDEs with delay), based on the use of auxiliary solutions to
simpler equations. These approaches are based on the principle of structural
analogy of solutions, which is formulated as follows: exact solutions of
simpler equations can serve as the basis for constructing solutions to more
complex related equations. Note that this principle at the heuristic level was
successfully used in works [2, 26, 30, 31, 32, 35] to construct exact
solutions of nonlinear PDEs with delay, as well as functional PDEs, using
solutions of simpler related PDEs without delay. In this article, the results
obtained using the principle of analogies and illustrated with specific
examples, it was possible to generalize and strictly prove some constructive
statements.

##  2 General formulations of the problems under consideration

Let us give general formulations of two important problems discussed in this
article.

Problem 1. Let the exact solution $u=u(x,t)$ of some nonlinear PDE be known,
which depends on the free parameters $a_{1},\dots,a_{p}$. The problematic
question arises: in what cases can we say anything about the exact solutions
of a more complex PDE, which is obtained from the original one by replacing
the free parameters with arbitrary functions $a_{1}(t),\dots,a_{p}(t)$?

The next section will show how exact solutions with generalized or functional
separation of variables of nonlinear autonomous PDEs, the coefficients of
which do not depend on time, can be used to construct exact solutions of more
general nonlinear non-autonomous PDEs, the coefficients of which depend
arbitrarily on time $t$.

Problem 2. Let the exact solution $u=u(x,t)$ of some nonlinear PDE with a
constant delay $\tau$ be known. The problematic question arises: in what cases
can we say anything about the exact solutions of a more complex PDE, which is
obtained from the original one by replacing the constant delay $\tau$ with an
arbitrary variable delay $\tau(t)$?

Next, in Section 4 we will show how exact solutions with generalized or
functional separation of variables of nonlinear PDEs with constant delay can
be used to construct exact solutions of more general nonlinear PDEs with
variable delay of the general form.

##  3 Using solutions of autonomous PDEs to construct exact solutions of non-
autonomous PDEs

Let us describe a method for constructing exact solutions of non-autonomous
PDEs, whose coefficients depend on time $t$, which is based on the use of
solutions of simpler autonomous PDEs, the coefficients of which do not depend
on time. Following the principle ‘‘from simple to complex’’, we will first
demonstrate the sequence and logic of reasoning in such cases using several
specific illustrative examples, after analyzing which, we will then formulate
general conclusions in the form of rigorous statements.

Example 1. Let us consider the nonlinear reaction-diffusion type equation

$u_{t}=a_{1}uu_{xx}+a_{2},$ (1)

where $a_{1}$ and $a_{2}$ are free parameters. Equation (1) admits an exact
solution with generalized separation of variables, quadratic in the spatial
variable

$u=\psi_{1}(t)x^{2}+\psi_{2}(t),$ (2)

where the functions $\psi_{1}=\psi_{1}(t)$ and $\psi_{2}=\psi_{2}(t)$ satisfy
a nonlinear ODE system,

$\displaystyle\psi_{1}^{\prime}=2a_{1}\psi_{1}^{2},$ (3)
$\displaystyle\psi_{2}^{\prime}=2a_{1}\psi_{1}\psi_{2}+a_{2}.$

Here and below, the primes denote derivatives with respect to $t$. System (3)
is easily integrated, but we will not need this further.

Now formally replacing the parameters $a_{1}$ and $a_{2}$ in equation (1) with
arbitrary functions $a_{1}(t)$ and $a_{2}(t)$, we obtain

$u_{t}=a_{1}(t)uu_{xx}+a_{2}(t).$ (4)

It is easy to check that equation (4) also admits an exact solution of the
form (2), where the functions $\psi_{1}=\psi_{1}(t)$ and
$\psi_{2}=\psi_{2}(t)$ satisfy the nonlinear ODE system (3), in which the
parameters $a_{1}$ and $a_{2}$ must be respectively replaced by arbitrary
functions $a_{1}(t)$ and $a_{2}(t)$.

The considered illustrative example allows us to make a simple generalization,
which can be formulated in the form of the following statement.

Statement 1. Let the autonomous nonlinear partial differential equation

$u_{t}=F\bigl{(}x,u,u_{x},\dots,u^{(n)}_{x};a_{1},\dots,a_{p}\bigr{)},$ (5)

depending on the free parameters $a_{1},\dots,a_{p}$, has a generalized
separable solution in the form of a polynomial in the spatial variable $x$:

$\displaystyle u=\sum^{m}_{k=0}\psi_{k}(t)x^{k}.$ (6)

Then the more complex non-autonomous partial differential equation

$u_{t}=F\bigl{(}x,u,u_{x},\dots,u^{(n)}_{x};a_{1}(t),\dots,a_{p}(t)\bigr{)},$
(7)

which is obtained from equation (5) by formally replacing the parameters
$a_{1},\dots,a_{p}$ with arbitrary functions $a_{1}(t),\dots,a_{p}(t)$, also
admits an exact solution of the form (6) with other functions $\psi_{k}(t)$.

Statement 1 is a special case of the more general Statement 2 formulated
below.

Example 2. Let us consider another nonlinear equation of the reaction-
diffusion type

$u_{t}=a_{1}u_{xx}+u_{x}^{2}+bu^{2}+a_{2}\qquad(b>0),$ (8)

where $a_{1}$, $a_{2}$, and $b$ are free parameters. Equation (8) admits an
exact solution with generalized separation of variables containing a
trigonometric function $x$ of the form

$u=\psi_{1}(t)\cos(\sqrt{b}\,x)+\psi_{2}(t),$ (9)

where the functions $\psi_{1}=\psi_{1}(t)$ and $\psi_{2}=\psi_{2}(t)$ are
described by the following nonlinear ODE system:

$\displaystyle\psi_{1}^{\prime}=2b\psi_{1}\psi_{2}-a_{1}b\psi_{1},$ (10)
$\displaystyle\psi_{2}^{\prime}=b(\psi_{1}^{2}+\psi_{2}^{2})+a_{2}.$

Formally replacing the parameters $a_{1}$ and $a_{2}$ in equation (8) with
arbitrary functions $a_{1}(t)$ and $a_{2}(t)$, we obtain

$u_{t}=a_{1}(t)u_{xx}+u_{x}^{2}+bu^{2}+a_{2}(t)\qquad(b>0).$ (11)

It is easy to check that equation (11) also admits an exact solution of the
form (9), where the functions $\psi_{1}=\psi_{1}(t)$ and
$\psi_{2}=\psi_{2}(t)$ satisfy the nonlinear ODE system (10), in which the
parameters $a_{1}$ and $a_{2}$ must be replaced by arbitrary functions
$a_{1}(t)$ and $a_{2}(t)$.

An attempt to similarly replace the constant $b$ in equation (8) with an
arbitrary function $b(t)$ will be unsuccessful, since in this case the
function (9) with $b=b(t)$ will now not be an exact solution to such a
modified equation. This happened because the coordinate function
$\varphi(x)=\cos(\sqrt{b}\,x)$ in the solution (9) explicitly depends on the
parameter $b$.

The considered illustrative example allows us to make a fairly obvious
generalization, which we formulate in the form of the following statement.

Statement 2. Let the autonomous nonlinear partial differential equation (5),
depending on the free parameters $a_{1},\dots,a_{p}$, has a generalized
separable solution of the form

$\displaystyle u=\sum^{m}_{k=1}\psi_{k}(t)\varphi_{k}(x),$ (12)

in which all linearly independent coordinate functions $\varphi_{k}(x)$ do not
depend on the parameters $a_{1},\dots,a_{p}$, and the functions
$\psi_{k}=\psi_{k}(t)$ are described by an autonomous ODE system,

$\displaystyle\psi_{k}^{\prime}=f_{k}(\psi_{1},\dots,\psi_{m};a_{1},\dots,a_{p}),\quad\
k=1,\dots,m.$ (13)

Then a more complex non-autonomous partial differential equation (7), which is
obtained from equation (5) by formally replacing the parameters
$a_{1},\dots,a_{k}$ with arbitrary functions $a_{1}(t),\dots,a_{p}(t)$, admits
an exact solution of the form (12), where the coordinate functions
$\varphi_{k}(x)$ do not change, and the functions $\psi_{k}=\psi_{k}(t)$ are
described by a non-autonomous ODE system,

$\displaystyle\psi_{k}^{\prime}=f_{k}\bigl{(}\psi_{1},\dots,\psi_{m};a_{1}(t),\dots,a_{p}(t)\bigr{)},\quad\
k=1,\dots,m.$ (14)

This statement can be proven using the results of [14] (see also [9, 1]).
Indeed, let us substitute (12) into the autonomous equation (5), and then
eliminate the time derivatives using ODEs (13). Further replacing the
functions $\psi_{k}(t)$ with arbitrary constants $C_{k}$, we obtain

$F\,\Bigl{[}\sum^{m}_{k=1}C_{k}\varphi_{k}(x)\Bigr{]}=\sum^{m}_{k=1}f_{k}(C_{1},\ldots,C_{m};a_{1},\dots,a_{p})\varphi_{k}(x),$
(15)

where $F[u]$ denotes the right-hand side of equation (5). Relations (15) mean
(see [14, 1]) that a finite-dimensional linear subspace,

$\mathcal{L}_{\,m}=\bigl{\\{}\varphi_{1}(x),\ldots,\varphi_{m}(x)\bigr{\\}},$
(16)

the elements of which are all possible linear combinations of coordinate
functions $\varphi_{1}(x)$, …, $\varphi_{m}(x)$ included in the solution (12),
is invariant with respect to the nonlinear differential operator in the
spatial variable $x$, which is on the right side of the autonomous equation
(5). Nonlinear differential operator on the right side of the more complex
non-autonomous partial differential equation (7) obtained from the equation
(5) by formally replacing the free parameters $a_{1},\dots,a_{k}$ with
arbitrary functions $a_{1}(t),\dots,a_{p}(t)$, depends parametrically on time
$t$ (contains derivatives only with respect to the variable $x$). Therefore,
when substituting the expression (12) into the right side of the equation (7),
the functions $a_{1}(t),\dots,a_{p}(t)$ behave like constants [1], which
ultimately and leads to the non-autonomous ODE system (14) for functions
$\psi_{k}=\psi_{k}(t)$.

Remark 1. On the left side of the nonlinear PDEs (5) and (7), instead of the
first derivative $\,u_{t}$ there can be the second derivative $\,u_{tt}$ or
any linear combination of derivatives with respect to $t$ of the form
$\,L[u]=\sum^{s}_{i=1}c_{i}(t)u^{(i)}_{t}$. In this case, in Statement 2 on
the left side of the ODEs (13) and (14), the first derivatives
$\,\psi_{k}^{\prime}$ must be replaced by the second derivatives
$\,\psi_{k}^{\prime\prime}$ or expressions $\,L[\psi_{k}]$.

Statement 2 allows for a simplified, more concise (but less informative)
formulation, which is given below.

Statement 2a. Let the autonomous partial differential equation (5), depending
on the free parameters $a_{1},\dots,a_{p}$, has solution (12), in which the
linearly independent coordinate functions $\varphi_{k}(x)$ are not depend on
the parameters $a_{1},\dots,a_{p}$. Then the more complex non-autonomous
partial differential equation (7), obtained from equation (5) by formally
replacing the parameters $a_{1},\dots,a_{k}$ with arbitrary functions
$a_{1}(t),\dots,a_{p}(t)$ also admits an exact solution of the form (12), in
which the coordinate functions $\varphi_{k}(x)$ do not change, and the
functions $\psi_{k}=\psi_{k}(t)$ are described by an appropriate system of
ODEs.

Statements 1 and 2 can also be generalized to the case of nonlinear PDEs,
which allow more complex exact solutions with functional separation of
variables. Let us illustrate this with a specific example.

Example 3. Let us consider a reaction-diffusion type equation with logarithmic
nonlinearity,

$u_{t}=a_{1}u_{xx}+a_{2}u\ln u,$ (17)

where $a_{1}$ and $a_{2}$ are free parameters. Equation (17) admits a
functional separable solution of the form [14, 1]:

$u=\exp[\psi_{1}(t)x^{2}+\psi_{2}(t)x+\psi_{3}(t)],$ (18)

where the functions $\psi_{n}=\psi_{n}(t)$ ($n=1,\,2,\,3$) satisfy the
nonlinear ODE system

$\displaystyle\psi_{1}^{\prime}$
$\displaystyle=4a_{1}\psi_{1}^{2}+a_{2}\psi_{1},$ (19)
$\displaystyle\psi_{2}^{\prime}$
$\displaystyle=4a_{1}\psi_{1}\psi_{2}+a_{2}\psi_{2},$
$\displaystyle\psi_{3}^{\prime}$
$\displaystyle=a_{2}\psi_{3}+2a_{1}\psi_{1}+a_{1}\psi_{2}^{2}.$

Now formally replacing the parameters $a_{1}$ and $a_{2}$ in equation (17)
with arbitrary functions $a_{1}(t)$ and $a_{2}(t)$, we obtain a more complex
equation,

$u_{t}=a_{1}(t)u_{xx}+a_{2}(t)u\ln u,$ (20)

It is easy to check that equation (20) also admits an exact solution of the
form (18), where the functions $\psi_{n}=\psi_{n}(t)$ ($n=1,\,2,\,3$) satisfy
the nonlinear ODE system (19), in which the parameters $a_{1}$ and $a_{2}$
must be respectively replaced by arbitrary functions $a_{1}(t)$ and
$a_{2}(t)$.

The considered illustrative example allows us to make a fairly simple
generalization given below.

Statement 3. Let the autonomous nonlinear partial differential equation (5),
depending on the free parameters $a_{1},\dots,a_{p}$, has a functional
separable solution of the form

$\displaystyle u=U(z),\quad\ z=\sum^{m}_{k=0}\psi_{k}(t)x^{k},$ (21)

where $U(z)$ is some function independent of the parameters
$a_{1},\dots,a_{p}$. Then the more complex non-autonomous partial differential
equation (7), which is obtained from equation (5) by formally replacing the
parameters $a_{1},\dots,a_{p}$ with arbitrary functions
$a_{1}(t),\dots,a_{p}(t)$, also admits an exact solution of the form (21) with
the same function $U(z)$, but with different functions $\psi_{k}(t)$.

This statement is proven by passing in the PDEs (5) and (7) to a new variable
$z$ using the substitution $u=U(z)$, after which the transformed equations
will have generalized separable solutions to which Statement 1 is already
applicable.

Statement 3 allows for further generalization, described below.

Statement 4. Let the autonomous nonlinear partial differential equation (5),
depending on the free parameters $a_{1},\dots,a_{p}$, has a functional
separable solution of the form

$\displaystyle u(x,t)=U(z),\quad\hbox{where}\quad
z=\sum^{m}_{k=1}\varphi_{k}(x)\psi_{k}(t),$ (22)

in which all coordinate functions $\varphi_{k}(x)$ and the external function
$U(z)$ do not depend on the parameters $a_{1},\dots,a_{p}$, and the functions
$\psi_{k}=\psi_{k}(t)$ are described by the autonomous ODE system (13). Then
the more complex non-autonomous partial differential equation (7), which is
obtained from the equation (5) by formally replacing the parameters
$a_{1},\dots,a_{p}$ with arbitrary functions $a_{1}(t),\dots,a_{p}(t)$, admits
an exact solution of the form (22), where the coordinate functions
$\varphi_{k}(x)$ and the function $U(z)$ do not change, and the functions
$\psi_{k}=\psi_{k}(t)$ are described by the non-autonomous ODE system (14).

Statement 4 is proven by passing to a new variable $z$ in the PDEs (5) and (7)
using the substitution $u=U(z)$, after which the transformed equations will
have generalized separable solutions to which Statement 2 applies.

Statement 4 allows for a simplified, shorter (but less informative)
formulation, which is given below.

Statement 4a. Let the autonomous nonlinear partial differential equation (5),
depending on the free parameters $a_{1},\dots,a_{p}$, has a functional
separable solution of the form (22), in which all coordinate functions
$\varphi_{k}(x)$ and the external function $U(z)$ do not depend on the
parameters $a_{1},\dots,a_{p}$. Then the more complex non-autonomous partial
differential equation (7), which is obtained from the equation (5) by formally
replacing the parameters $a_{1},\dots,a_{p}$ with arbitrary functions
$a_{1}(t),\dots,a_{p}(t)$, admits an exact solution of the form (22), where
the coordinate functions $\varphi_{k}(x)$ and the function $U(z)$ do not
change, and the functions $\psi_{k}=\psi_{k}(t)$ are described by an
appropriate ODE system.

##  4 Using solutions of PDEs with constant delay to construct exact solutions
of PDEs with variable delay

Let us describe a method for constructing exact solutions of PDEs with
variable delay of general form, based on the use of solutions of related
simpler PDEs with constant delay. We will demonstrate the sequence and logic
of useful preliminary reasoning in such cases using two specific examples,
after analyzing which we will formulate general conclusions in the form of
statements.

Example 4. Let us consider the reaction-diffusion equation with quadratic
nonlinearity and constant delay

$\displaystyle u_{t}=a(uu_{x})_{x}+b\bar{u},\quad\ \bar{u}=u(x,t-\tau),$ (23)

with $\tau=\text{const}>0$, which admits a generalized separable solution of
the form

$\displaystyle u=\psi_{1}(t)x^{2}+\psi_{2}(t),$ (24)

where the functions $\psi_{1}=\psi_{1}(t)$ and $\psi_{2}=\psi_{2}(t)$ are
described by a nonlinear ODE system with constant delay,

$\displaystyle\psi_{1}^{\prime}=6a\psi_{1}^{2}+b\bar{\psi}_{1},\quad\
\bar{\psi}_{1}=\psi_{1}(t-\tau),$ (25)
$\displaystyle\psi_{2}^{\prime}=2a\psi_{1}\psi_{2}+b\bar{\psi}_{2},\quad\
\bar{\psi}_{2}=\psi_{2}(t-\tau).$

Replacing the constant $\tau>0$ in (23) with an arbitrary positive time
function $\tau(t)$, we obtain a more complex nonlinear reaction-diffusion
equation with quadratic nonlinearity and variable delay,

$\displaystyle u_{t}=a(uu_{x})_{x}+b\bar{u},\quad\ \bar{u}=u(x,t-\tau(t)).$
(26)

Since the nonlinear terms containing derivatives with respect to $x$ on the
right-hand sides of equations (23) and (26) are the same, it is natural to
assume that the power structure of solutions in the spatial variable of both
equations will also be the same, and the time-dependent functional
coefficients at various powers of $x$ will change.

In other words, we are looking for exact solutions to the reaction-diffusion
equation with variable delay (26) in the same form (24) as solutions to the
original equation with constant delay (23). As a result, for the functions
$\psi_{1}=\psi_{1}(t)$, $\psi_{2}=\psi_{2}(t)$ we obtain the ODE system (25),
in which the constant $\tau$ must be replaced by the function $\tau(t)$.

The results described in this example allow for the generalization formulated
below.

Statement 5. Let the nonlinear PDE with constant delay

$u_{t}=F\bigl{(}x,u,u_{x},\dots,u^{(n)}_{x};\bar{u},\bar{u}_{x},\dots,\bar{u}^{(p)}_{x}\bigr{)},\quad\
\bar{u}=u(x,t-\tau),$ (27)

where $n>p$, have the generalized separable solution

$\displaystyle u=\sum^{m}_{k=1}\varphi_{k}(x)\psi_{k}(t),$ (28)

in which the linearly independent coordinate functions $\varphi_{k}(x)$ do not
depend on $\tau$, and the functions $\psi_{k}=\psi_{k}(t)$ are described by
the ODE system with constant delay

$\displaystyle\psi_{k}^{\prime}=f_{k}(\psi_{1},\dots,\psi_{m};\bar{\psi}_{1},\dots,\bar{\psi}_{m}),\quad\
\bar{\psi}_{k}=\psi_{k}(t-\tau),\quad\ k=1,\dots,m.$ (29)

Then a more complex PDE with variable delay,

$u_{t}=F\bigl{(}x,u,u_{x},\dots,u^{(n)}_{x};\bar{u},\bar{u}_{x},\dots,\bar{u}^{(p)}_{x}\bigr{)},\quad\
\bar{u}=u\bigl{(}x,t-\tau(t)\bigr{)},$ (30)

which is obtained from equation (27) by formally replacing the constant $\tau$
with an arbitrary time function $\tau(t)$, also admits an exact solution of
the form (28), where the coordinate functions $\varphi_{k}(x)$ do not change,
and the functions $\psi_{k}=\psi_{k}(t)$ are described by an ODE system with
variable delay,

$\displaystyle\psi_{k}^{\prime}=f_{k}(\psi_{1},\dots,\psi_{m};\bar{\psi}_{1},\dots,\bar{\psi}_{m}),\quad\
\bar{\psi}_{k}=\psi_{k}\bigl{(}t-\tau(t)\bigr{)},\quad\ k=1,\dots,m.$ (31)

The proof of Statement 5 is carried out in the same way as the proof of
Statement 2.

Remark 2. On the left side of the nonlinear PDEs with delay (27) and (30),
instead of the first derivative $\,u_{t}$ there can be the second derivative
$\,u_{tt}$ or any linear combination of derivatives with respect to $t$ of the
form $\ ,L[u]=\sum^{s}_{i=1}a_{i}(t)u^{(i)}_{t}$. In this case, in Statement 5
on the left side of the ODEs with delay (29) and (31), the first derivatives
$\,\psi_{k}^{\prime}$ must be replaced by the second derivatives
$\,\psi_{k}^{\prime\prime}$ or expressions $\,L[\psi_{k}]$.

Statement 5 can be generalized to the case of nonlinear PDEs with delay, which
allow more complex functional separable solutions. Let us illustrate this with
a specific example.

Example 5. Let us consider the reaction-diffusion type equation with two
logarithmic nonlinearities and constant delay

$u_{t}=u_{xx}+u(a\ln^{2}u+b\ln\bar{u}),\quad\ \bar{u}=u(x,t-\tau),$ (32)

which for $a>0$ admits a functional separable solution of the form [36]:

$u=\exp\bigl{[}\psi_{1}(t)\varphi(x)+\psi_{2}(t)\bigr{]},\quad\
\varphi(x)=C_{1}\cos\bigl{(}\sqrt{a}\,x\bigr{)}+C_{2}\sin\bigl{(}\sqrt{a}\,x\bigr{)},$
(33)

where $C_{1}$ and $C_{2}$ are arbitrary constants, and the functions
$\psi_{1}=\psi_{1}(t)$ and $\psi_{2}=\psi_{2}(t)$ satisfy the nonlinear ODE
system with constant delay

$\displaystyle\psi_{1}^{\prime}$
$\displaystyle=2a\psi_{1}\psi_{2}-a\psi_{1}+b\bar{\psi}_{1},\quad\
\bar{\psi}_{1}=\psi_{1}(t-\tau),$ (34) $\displaystyle\psi_{2}^{\prime}$
$\displaystyle=a(C_{1}^{2}+C_{2}^{2})\psi_{1}^{2}+a\psi_{2}^{2}+b\bar{\psi}_{2},\quad\
\bar{\psi}_{2}=\psi_{2}(t-\tau).$

Now formally replacing the constant $\tau$ in the equation (32) with an
arbitrary function $\tau(t)$, we obtain a more complex equation,

$u_{t}=u_{xx}+u(a\ln^{2}u+b\ln\bar{u}),\quad\
\bar{u}=u\bigl{(}x,t-\tau(t)\bigr{)}.$ (35)

It is easy to check that equation (35) also admits an exact solution of the
form (33), where the functions $\psi_{1}=\psi_{1}(t)$ and
$\psi_{2}=\psi_{2}(t)$ satisfy the nonlinear ODE system (34), in which the
constant $\tau$ must be replaced by an arbitrary function $\tau(t)$.

The considered illustrative example allows us to make a generalization, which
we will formulate in the form of a statement.

Statement 6. Let the nonlinear PDE with constant delay (27) have a functional
separable solution of the form

$\displaystyle u(x,t)=U(z),\quad\hbox{where}\quad
z=\sum^{m}_{k=1}\varphi_{k}(x)\psi_{k}(t),$ (36)

in which the linearly independent coordinate functions $\varphi_{k}(x)$ and
the external function $U(z)$ are independent on $\tau$, and the functions
$\psi_{k}=\psi_{k}(t)$ are described by the ODE system with constant delay
(29). Then the more complex PDE with variable delay (30), which is obtained
from equation (27) by formally replacing the constant $\tau$ with an arbitrary
function $\tau(t)$, admits an exact solution of the form (36), where the
coordinate functions $\varphi_{k}(x)$ and the function $U(z)$ are unchanged,
and the functions $\psi_{k}=\psi_{k}(t)$ are described by the ODE system with
variable delay (31).

Table 1 lists some nonlinear reaction-diffusion equations with variable delay,
allowing generalized and functional separable solutions. It is believed that
$\tau=\tau(t)$ is an arbitrary positive continuously differentiable function
(in particular, in the case of proportional delay in the equations one should
put $\tau=(1-p)t$, i.e. $t-\tau=pt$ , where $0<p<1$). Note that the functional
separable solutions of the last two equations in Table 1 are presented in
implicit form.

Table 1. Reaction-diffusion type equations with variable delay of a general
type, allowing solutions with generalized and functional separation of
variables. Notation: $\bar{u}=u(x,t-\tau(t))$, $\bar{\psi}=\psi(t-\tau(t))$;
$f(z)$, $g(z)$, and $\tau(t)$ are arbitrary functions; $C_{1}$, $C_{2}$, and
$C_{3}$ are arbitrary constants.

|  |  |
---|---|---|---
|  Original equation |  Form of exact solution |  Determining equations
|  |  |
|  |  |
|  |  |
|  $u_{t}=au_{xx}+uf(\bar{u}/u)$ |  $u=[C_{1}\cos(\beta x)+C_{2}\sin(\beta x)]\psi(t)$; $u=[C_{1}\exp(-\beta x)+C_{2}\exp(\beta x)]\psi(t)$; $u=(C_{1}x+C_{2})\psi(t)$ |  $\psi^{\prime}_{t}=-a\beta^{2}\psi+\psi f(\bar{\psi}/\psi)$; $\psi^{\prime}_{t}=a\beta^{2}\psi+\psi f(\bar{\psi}/\psi)$; $\psi^{\prime}_{t}=\psi f(\bar{\psi}/\psi)$
|  |  |
|  |  |
|  $u_{t}=au_{xx}+bu\ln u$ $+\,uf(\bar{u}/u)$ |  $u=\varphi(x)\psi(t)$ |  $a\varphi^{\prime\prime}_{xx}=C_{1}\varphi-b\varphi\ln\varphi$, $\psi^{\prime}_{t}=C_{1}\psi+b\psi\ln\psi$ $+\,\psi f(\bar{\psi}/\psi)$
|  |  |
|  |  |
|  $u_{t}=au_{xx}+f(u-\bar{u})$ |  $u=C_{2}x^{2}+C_{1}x+\psi(t)$ |  $\psi^{\prime}_{t}=2C_{2}a+f(\psi-\bar{\psi})$
|  |  |
|  |  |
|  $u_{t}=au_{xx}+bu$ $+\,f(u-\bar{u})$ |  $u=C_{1}\cos(\lambda x)+C_{2}\sin(\lambda x)+\psi(t)$, where $\lambda=\sqrt{b/a}$ (if $ab>0$); $u=C_{1}\exp(-\lambda x)+C_{2}\exp(\lambda x)+\psi(t)$, where $\lambda=\sqrt{-b/a}$ (if $ab<0$) |  $\psi^{\prime}_{t}=b\psi+f(\psi-\bar{\psi})$; $\psi^{\prime}_{t}=b\psi+f(\psi-\bar{\psi})$
|  |  |
|  |  |
|  $u_{t}=a(u^{k}u_{x})_{x}+uf(\bar{u}/u)$ |  $u=\varphi(x)\psi(t)$ |  $a(\varphi^{k}\varphi^{\prime}_{x})^{\prime}_{x}=C_{1}\varphi$, $\psi^{\prime}_{t}=C_{1}\psi^{k+1}+\psi f(\bar{\psi}/\psi)$
|  |  |
|  |  |
|  $u_{t}=a(u^{k}u_{x})_{x}+bu^{k+1}$ $+\,uf(\bar{u}/u)$ |  $u=[C_{1}\cos(\beta x)+C_{2}\sin(\beta x)]^{\frac{1}{k+1}}\psi(t)$, where $\beta=\sqrt{b(k+1)/a\\!}$, $b(k+1)>0$; $u=(C_{1}e^{-\beta x}+C_{2}e^{\beta x})^{\frac{1}{k+1}}\psi(t)$, where $\beta=\sqrt{-b(k+1)/a\\!}$, $\,b(k+1)<\\!0$; $u=C_{1}\exp\bigl{(}-\frac{b}{2a}x^{2}+C_{2}x\bigr{)}\psi(t)$ for $k=-1$; $u=\varphi(x)\psi(t)$ (generalizes previous solutions) |  $\psi^{\prime}_{t}=\psi f(\bar{\psi}/\psi)$; $\psi^{\prime}_{t}=\psi f(\bar{\psi}/\psi)$; $\psi^{\prime}_{t}=\psi f(\bar{\psi}/\psi)$; $a(\varphi^{k}\varphi^{\prime}_{x})^{\prime}_{x}+b\varphi^{k+1}=C_{1}\varphi$, $\psi^{\prime}_{t}=C_{1}\psi^{k+1}+\psi f(\bar{\psi}/\psi)$
|  |  |
|  |  |
|  $u_{t}=a(e^{\lambda u}u_{x})_{x}+f(u-\bar{u})$ |  $\displaystyle u=\frac{1}{\lambda}\ln(C_{1}\lambda x^{2}+C_{2}x+C_{3})+\psi(t)$ |  $\psi^{\prime}_{t}=2aC_{1}e^{\lambda\psi}+f(\psi-\bar{\psi})$
|  |  |
|  |  |
|  $u_{t}=a(e^{\lambda u}u_{x})_{x}+be^{\lambda u}$ $+f(u-\bar{u})$ |  $u=\frac{1}{\lambda}\ln[C_{1}\\!\cos(\beta x)\\!+\\!C_{2}\sin(\beta x)]\\!+\\!\psi(t)$, where $\beta=\sqrt{b\lambda/a}$, $b\lambda>0$; $u=\frac{1}{\lambda}\ln(C_{1}e^{-\beta x}+C_{2}e^{\beta x})+\psi(t)$, where $\beta=\sqrt{-b\lambda/a}$, $b\lambda<0$; $u=\varphi(x)+\psi(t)$ (generalizes previous solutions) |  $\psi^{\prime}_{t}=f(\psi-\bar{\psi})$; $\psi^{\prime}_{t}=f(\psi-\bar{\psi})$; $a(e^{\lambda\varphi}\varphi^{\prime}_{x})^{\prime}_{x}+be^{\lambda\varphi}=C_{1}$, $\psi^{\prime}_{t}=C_{1}e^{\lambda\psi}+f(\psi-\bar{\psi})$
|  |  |
|  |  |
|  $u_{t}=[(a\ln u+b)u_{x}]_{x}$ $-\,cu\ln u+uf(\bar{u}/u)$ |  $\displaystyle u=\exp(\pm\lambda x)\psi(t)$, $\lambda=\sqrt{c/a}$ |  $\psi^{\prime}_{t}\\!=\lambda^{2}(a\\!+b)\psi+\psi f(\bar{\psi}/\psi)$
|  |  |
|  |  |
|  $u_{t}=a[f^{\prime}(u)u_{x}]_{x}+b$ $\quad\ +\frac{1}{f^{\prime}(u)}g(f(u)-f(\bar{u}))$ |  $\displaystyle f(u)=\psi(t)-\frac{b}{2a}x^{2}+C_{1}x+C_{2}$ |  $\psi^{\prime}_{t}=g(\psi-\bar{\psi})$
|  |  |
|  |  |
|  $u_{t}=a[uf^{\prime}(u)u_{x}]_{x}$ $\quad\ +\frac{1}{f^{\prime}(u)}[bf(u)+cf(\bar{u})]$ |  $\displaystyle f(u)=\varphi(t)x+\psi(t)$ |  $\varphi^{\prime}_{t}=b\varphi+c\bar{\varphi}$, $\psi^{\prime}_{t}=b\psi+c\bar{\psi}+a\varphi^{2}$
|  |  |

##  5 Brief conclusions

It has been established that generalized and functional separable solutions of
nonlinear autonomous PDEs, the coefficients of which do not depend on time
$t$, can be used to construct exact solutions of more general nonlinear non-
autonomous PDEs, the coefficients of which depend arbitrarily on time.
Examples of exact solutions of nonlinear non-autonomous PDEs are given.

It is shown that generalized and functional separable solutions of nonlinear
PDEs with constant delay make it possible to find exact solutions of more
complex nonlinear PDEs with variable delay of general form. Some nonlinear
reaction-diffusion equations with variable delay are described, allowing exact
solutions with a generalized separation of variables.

## References

  * 1. Polyanin, A.D., Zhurov, A.I.: Separation of Variables and Exact Solutions to Nonlinear PDEs. CRC Press, Boca Raton–London (2022)
  * 2. Polyanin, A.D., Sorokin, A.D., Zhurov, A.I.: Delay Ordinary and Partial Differential Equations. CRC Press, Boca Raton–London (2023)
  * 3. Ovsiannikov, L.V.: Group Analysis of Differential Equations. Academic Press, New York (1982)
  * 4. Bluman, G.W., Kumei, S.: Symmetries and Differential Equations. Springer, New York (1989)
  * 5. Ibragimov, N.H. (ed.): CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 1. Symmetries, Exact Solutions and Conservation Laws. CRC Press, Boca Raton (1994)
  * 6. Olver, P.J.: Application of Lie Groups to Differential Equations, 2nd ed. Springer, New York (2000)
  * 7. Clarkson, P.A., Kruskal, M.D.: New similarity reductions of the Boussinesq equation. J. Math. Phys. 30(10), 2201–2213 (1989)
  * 8. Hood, S.: On direct, implicit reductions of a nonlinear diffusion equation with an arbitrary function — generalizations of Clarkson’s and Kruskal’s method. IMA J. Appl. Math. 64(3), 223–244 (2000)
  * 9. Polyanin, A.D., Zaitsev, V.F.: Handbook of Nonlinear Partial Differential Equations, 2nd ed. CRC Press, Boca Raton (2012)
  * 10. Bluman, G.W., Cole, J.D.: The general similarity solution of the heat equation. J. Math. Mech. 18, 1025–1042 (1969)
  * 11. Levi, D., Winternitz, P.: Nonclassical symmetry reduction: Example of the Boussinesq equation. J. Phys. A. 22, 2915–2924 (1989)
  * 12. Arrigo, D., Broadbridge, P., Hill, J.M.: Nonclassical symmetry solutions and the methods of Bluman–Cole and Clarkson–Kruskal. J. Math. Phys. 34, 4692–4703 (1993)
  * 13. Cherniha, R., Davydovych, V.: Nonlinear Reaction-Diffusion Systems: Conditional Symmetry, Exact Solutions and Their Applications in Biology. Springer, Cham (2017)
  * 14. Galaktionov, V.A., Svirshchevskii, S.R.: Exact Solutions and Invariant Subspaces of Nonlinear Partial Differential Equations in Mechanics and Physics. Chapman & Hall/CRC Press, Boca Raton (2007)
  * 15. Kosov, A.A., Semenov, E.I.: Exact solutions of the generalized Richards equation with power-law nonlinearities. Diff. Equat. 56(9), 1119–1129 (2020)
  * 16. Kosov, A.A., Semenov, E.I.: Reduction method and new exact solutions of the multidimensional nonlinear heat equation. Diff. Equat. 58(2), 187–194 (2022)
  * 17. Doyle, P.W., Vassiliou, P.J.: Separation of variables for the 1-dimensional non-linear diffusion equation. Int. J. Non-Linear Mech. 33(2), 315–326 (1998)
  * 18. Estevez, P.G., Qu, C., Zhang, S.: Separation of variables of a generalized porous medium equation with nonlinear source. J. Math. Anal. Appl. 275, 44–59 (2002)
  * 19. Kaptsov, O.V., Verevkin, I.V.: Differential constraints and exact solutions of nonlinear diffusion equations. J. Phys. A: Math. Gen. 36, 1401–1414 (2003)
  * 20. Meleshko, S.V.: Methods for Constructing Exact Solutions of Partial Differential Equations. Springer, New York (2005)
  * 21. Kudryashov, N.A.: Simplest equation method to look for exact solutions of nonlinear differential equations. Chaos, Solitons and Fractals. 24(5), 1217–1231 (2005)
  * 22. Kudryashov, N.A.: Methods of Nonlinear Mathematical Physics (in Russian). Izd. Dom Intellekt Publ., Dolgoprudnyi (2010)
  * 23. Conte, R., Musette, M. The Painlevé Handbook, 2nd ed. Springer, Cham (2020)
  * 24. Calogero, F., Degasperis, A.: Spectral Transform and Solitons: Tolls to Solve and Investigate Nonlinear Evolution Equations. North Holland, Amsterdam (1982)
  * 25. Ablowitz, M.J., Clarkson, P.A.: Solitons, Nonlinear Evolution Equations and Inverse Scattering. Cambridge Univ. Press, Cambridge (1991)
  * 26. Aksenov, A.V., Polyanin, A.D.: Methods for constructing complex solutions of nonlinear PDEs using simpler solutions. Mathematics 9(4), 345 (2021)
  * 27. Bedrikovetsky, P., Borazjani, S.: Exact solutions for gravity-segregated flows in porous media. Mathematics 10, 2455 (2022).
  * 28. Meleshko, S.V., Moyo S.: On the complete group classification of the reaction–diffusion equation with a delay. J. Math. Anal. Appl. 338, 448–466 (2008)
  * 29. Long, F.-S., Meleshko, S.V. On the complete group classification of the one-dimensional nonlinear Klein–Gordon equation with a delay. Math. Methods Appl. Sci. 39(12), 3255–3270 (2016)
  * 30. Polyanin, A.D., Sorokin, V.G.: A method for constructing exact solutions of nonlinear delay PDEs. J. Math. Anal. Appl. 494, 124619 (2021)
  * 31. Polyanin, A.D., Sorokin, V.G.: Construction of exact solutions to nonlinear PDEs with delay using solutions of simpler PDEs without delay. Comm. Nonlinear Sci. Numer. Simul. 95, 105634 (2021)
  * 32. Polyanin, A.D., Sorokin, V.G.: Nonlinear pantograph-type diffusion PDEs: Exact solutions and the principle of analogy. Mathematics. 9(5), 511 (2021)
  * 33. Aibinu, M.O., Thakur, S.C., Moyo, S.: Exact solutions of nonlinear delay reaction-diffusion equations with variable coefficients. Partial Diff. Equat. Appl. Math. 4, 100170 (2021)
  * 34. Sorokin, V.G., Vyazmin, A.V.: Nonlinear reaction-diffusion equations with delay: Partial survey, exact solutions, test problems, and numerical integration. Mathematics. 10(11), 1886 (2022)
  * 35. Polyanin, A.D., Sorokin, V.G.: Exact solutions of reaction-diffusion PDEs with anisotropic time delay. Mathematics. 11(14), 3111 (2023)
  * 36. Polyanin, A.D., Zhurov, A.I. Exact solutions of linear and non-linear differential-difference heat and diffusion equations with finite relaxation time. Int. J. Non-Linear Mech. 54, 115–126 (2013)