The compressible Navier–Stokes equations with slip boundary conditions of friction type

Abstract

We study a mathematical model of a viscous compressible fluid obeying the slip boundary condition of friction type. We present a notion of weak solutions to this model, in which the momentum equation and the associated energy inequality are combined into a single relation. Moreover, the slip boundary condition of friction type is incorporated into this relation by the use of a boundary integral. Our main result proves the existence of such weak solutions. The proof of this result combines the classical existence theory for the compressible Navier–Stokes equations with an approximation of the aforementioned boundary integral via a convex regularization of the absolute value function.

1 Introduction

The subject of this article is a new model of a compressible viscous fluid satisfying the so-called slip boundary condition of friction type (SBCF). In this paper, we prove the existence of weak solutions to this model.

Classical (or strong) solutions to the compressible Navier–Stokes equations can be expected only for small data (or, more generally, for data close to an equilibrium). The first such result for the Cauchy problem for the Navier–Stokes Fourier system (when heat conductivity is included) goes back to the eighties [18]; for further developments, see e.g., [20] or [2]. However, classical solutions are not known to exist globally in time if the data are arbitrary. The concept of weak solutions was for the first time successfully used by Lions (see [17]) in the case of isentropic flow. In this book, several kinds of boundary conditions were considered: The no-slip boundary condition, which describes the vanishing of the fluid velocity on the boundary of the domain, periodic boundary conditions as well as the case of a fluid covering the whole space. A detailed proof of the existence of weak solutions in the case of the no-slip boundary condition can further be found in [21]. A weak solution in the case of heat-conducting fluids satisfying the no-slip boundary condition was for the first time constructed by Feireisl (see [4]) by combining the internal energy balance and the global energy balance. Another approach, presented by Feireisl and Novotný, is based on the entropy inequality (see [5]). In the latter book, the case of the complete-slip boundary condition, i.e., the case of fluids for which the normal component of the velocity vanishes on the boundary, is additionally taken into consideration. Moreover, the existence of weak solutions in the case of incompressible fluids is treated for example in [16] for the no-slip boundary condition, periodic boundary conditions as well as in the whole space ${\mathbb {R}}^N$.

The no-slip boundary condition has been the most widely used given its success in reproducing the standard velocity profiles for incompressible/compressible viscous fluids for many years. The no-slip hypothesis seems to be in good agreement with experiments but it can lead to certain rather surprising conclusions, e.g., the most striking one being the absence of collisions of rigid objects immersed in a linearly viscous fluid [11, 12].

The Navier–Stokes equations have also been studied in combination with more uncommon boundary conditions. The so-called Navier boundary condition, which allows for slip, offers more freedom and is likely to provide a physically acceptable solution at least to some of the paradoxical phenomena resulting from the no-slip boundary condition, see e.g., Moffat [19]. Recent developments in macrofluidic and nanofluidic technologies have renewed interest in the slip behavior that may become significant in the small spatial scales even for a relatively small Reynolds number (cf. Priezjev and Troian [22]). Mathematically, the behavior of the tangential component of the velocity is a delicate issue.

We further mention the Coulomb friction law boundary condition, which is used for the description of fluids that can slip on the boundary provided that the tangential component of the stress tensor is sufficiently large. In [1], the existence of weak solutions to the incompressible Navier–Stokes equations satisfying this boundary condition is proved in the case of two and three spatial dimensions. Another boundary condition modelling this phenomenon is the slip boundary condition of friction type introduced by Fujita [7] and Fujita [8] for the stationary Stokes and Navier Stokes equations. The same boundary condition was studied for the incompressible Navier–Stokes equations in [14], wherein the existence of solutions is proved globally in time in the 2D case and locally in time in the 3D case. A numerical analysis of the slip boundary condition of friction type can be found in [13]. Moreover, some applications to real-world problems with numerical simulations are given in [9, 10, 15].

In the present article, we combine, for the first time, the compressible Navier–Stokes equations with the slip boundary condition of friction type. We prove the existence of weak solutions in this setting. Since the slip boundary condition of friction type is particularly interesting for the modelling of fluids in moving domains or fluid–structure interaction, our result can be considered as a first stepping stone toward the study of these more sophisticated problems. From the mathematical point of view, the main novelty in our existence proof lies in the addition of one further approximation level to the classical approximation method used for the construction of weak solutions to the compressible Navier–Stokes equations with the no-slip boundary condition, c.f. for example [21]. This additional approximation level follows closely the approximation methods used in [1, 14] in the case of the incompressible Navier–Stokes equations with the Coulomb friction law boundary condition and the slip boundary condition of friction type, respectively. It consists of the addition of a boundary integral to the momentum equation which contains the gradient of a smooth and convex approximation of the absolute value of the velocity field, c.f. (28). Due to the convexity of the approximation, this boundary integral can later be replaced by the desired boundary integral which expresses the slip boundary condition of friction type in our weak formulation (12). Another novelty results from the fact that as in the incompressible case in [14], the weak formulation to our problem merges the momentum equation and the energy inequality into one single relation, c.f. (12). For technical reasons, however, we also need to study the momentum equation separately in order to deduce the same improved density estimates and the effective viscous flux identity as in the existence proof in the case of the no-slip boundary condition in [21], which are required for passing to the limit in the pressure term. As a consequence, we are forced to also pass to the limit in the momentum equation separately on every approximation level, which leads to a relation which we refer to as the alternative momentum equation, c.f. Remark 4.1.

The paper is organized as follows. In Sect. 2, we present the full model. A corresponding weak formulation of this model is presented in Sect. 3. In the same section, we further show that this weak formulation constitutes a suitable definition of weak solutions and present our main result. The full proof of the main result extends across Sects. 4.1–4.5.

2 Model

The model which we study in this paper is as follows. We consider a viscous compressible fluid occupying an open and bounded domain $\Omega \subset {\mathbb {R}}^3$ with locally Lipschitz boundary $\Gamma = \partial \Omega$ and with outward unit normal vector $\text {n}$ on $\Gamma$. The density $\rho :(0,T) \times \Omega \rightarrow {\mathbb {R}}$ and the velocity field $u :(0,T) \times \Omega \rightarrow {\mathbb {R}}^3$ of the fluid are determined via the system

$$\begin{aligned}&\partial _t \rho + \nabla \cdot (\rho u) = 0 \quad \text {in } (0,T)\times \Omega , \end{aligned}$$

(1)

$$\begin{aligned}&\partial _{t} (\rho u) + \nabla \cdot (\rho u \otimes u) = \nabla \cdot \sigma + \rho f \quad \text {in } (0,T)\times \Omega , \end{aligned}$$

(2)

$$\begin{aligned}&\rho (0) = \rho ^0,\quad (\rho u)(0) = q \quad \text {in } \Omega , \end{aligned}$$

(3)

$$\begin{aligned}&u\cdot \text {n} = 0 \quad \text {on } (0,T)\times \Gamma , \end{aligned}$$

(4)

$$\begin{aligned}&\left| \left( \sigma \text {n} \right) _\tau \right| \le g,\quad \left( \sigma \text {n} \right) _\tau \cdot u_\tau + g \left| u_\tau \right| = 0 \quad \text {on } (0,T)\times \Gamma , \end{aligned}$$

(5)

where the Cauchy stress tensor

$$\begin{aligned} \sigma = \sigma (u,p) := 2\nu {\mathbb {D}}(u) + \lambda (\nabla \cdot u)\text {Id} - p \text {Id},\quad {\mathbb {D}}(u) := \frac{1}{2} \nabla u + \frac{1}{2} (\nabla u)^T \end{aligned}$$

with the viscosity coefficients $\nu ,\lambda \in {\mathbb {R}}$, satisfying

$$\begin{aligned} \nu > 0,\quad \lambda + \nu \ge 0, \end{aligned}$$

can be split into its normal component $\sigma _\text {n} = \sigma \text {n} \cdot \text {n}$ and its tangential component $\sigma _\tau = \sigma \text {n} - \sigma _\text {n}\text {n}$. Moreover, the positive constant g in (5) is the threshold of slippage. Further, the pressure p is defined by the isentropic constitutive relation

$$\begin{aligned} p = a\rho ^\gamma ,\quad \gamma> \frac{3}{2},\ a > 0. \end{aligned}$$

In the considered model, Eqs. (1) and (2) represent the continuity equation and the momentum equation, respectively. The initial conditions are presented in (3). Finally, Eqs. (4) and (5) represent the slip boundary condition of friction type.

3 Weak formulation and main result

Here, we present the definition of a weak solution to the system (1)–(5) and state our main result. To this end, we denote by $H_\text {n}^1(\Omega ; {\mathbb {R}}^3)$ the Sobolev space of all functions in $H^1(\Omega ;{\mathbb {R}}^3)$ whose normal component vanishes on the boundary,

$$\begin{aligned} H_\text {n}^1\left( \Omega ;{\mathbb {R}}^3 \right) := \left\{ u \in H^1\left( \Omega ;{\mathbb {R}}^3 \right) :\ u\cdot \text {n} = 0 \ \text {on } \Gamma \right\} . \end{aligned}$$

Definition 1.1

Let $T>0$ and let $\Omega \subset {\mathbb {R}}^3$ be a bounded domain. Let $\nu , \lambda , a, \gamma \in {\mathbb {R}}$ be given constants which satisfy

$$\begin{aligned} \nu , a> 0,\quad \gamma > \frac{3}{2},\quad \nu + \lambda \ge 0. \end{aligned}$$

(6)

Further assume that $f \in L^\infty ((0,T)\times \Omega )$, $g \in L^2((0,T)\times \Gamma )$ and assume the initial data to satisfy the conditions

$$\begin{aligned} 0 \le \rho _0 \in L^\gamma (\Omega ),\quad q \in L^1(\Omega ),\quad \frac{\left| q\right| ^2}{\rho _0} \in L^1(\Omega ),\quad q = 0 \ \ \text {a.e. in } \left\{ x \in \Omega :\ \rho _0(x) = 0 \right\} \end{aligned}$$

(7)

Then, a pair of functions $(\rho , u)$, such that

$$\begin{aligned} 0 \le \rho \in L^\infty \left( 0,T;L^\gamma \left( \Omega ;{\mathbb {R}}\right) \right) \bigcap C\left( [0,T];L^1\left( \Omega ;{\mathbb {R}} \right) \right) \quad \text{ and } \quad u \in L^2\left( 0,T;H_{{\text {n}}}^1\left( \Omega , {\mathbb {R}}^3 \right) \right) , \end{aligned}$$

(8)

is said to be a weak solution to system (1)–(5) if it satisfies:

• (i) the continuity equation in the distributional sense,

  $$\begin{aligned} \partial _t \rho + \nabla \cdot \left( \rho u \right) = 0 \quad&\text {in } {\mathcal {D}}'\left( (0,T) \times {\mathbb {R}}^3 \right) , \end{aligned}$$

  (9)

  and the renormalized sense,

  $$\begin{aligned}&\ \partial _t \zeta (\rho ) + \nabla \cdot \left( \zeta \left( \rho \right) u \right) + \left[ \zeta '\left( \rho \right) \rho - \zeta \left( \rho \right) \right] \nabla \cdot u = 0 \quad \text {in } {\mathcal {D}}' \left( (0,T) \times {\mathbb {R}}^3 \right) , \end{aligned}$$

  (10)
  $$\begin{aligned}&\quad \text {for all} \quad \zeta \in C^1\left( [0,\infty ) \right) : \quad \left| \zeta '(r)\right| \le cr^{\sigma } \quad \forall r \ge 1 \quad \text {for certain } c>0,\ \sigma > -1, \end{aligned}$$

  (11)

• (ii) the momentum and energy inequality

  $$\begin{aligned}&\int \limits _\Omega \frac{1}{2}\rho (0) \left| u(0)\right| ^2 + \frac{a\rho ^\gamma (0)}{\gamma - 1} \ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau ) \left| u(\tau ) \right| ^2 + \frac{a\rho ^\gamma (\tau )}{\gamma - 1} \ dx \nonumber \\ {}&\quad +\, \int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - (\rho u \otimes u) : \nabla \phi + \left( 2\nu {\mathbb {D}}(u) + \lambda \left( \nabla \cdot u \right) \, {\mathcal {I}} \right) : \left[ \nabla \phi - \nabla u \right] \nonumber \\ {}&\quad -\,p \, {\mathcal {I}} : \nabla \phi - \rho f \cdot \left[ \phi - u \right] \ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | - g |u| d\Gamma dt \ge 0 \end{aligned}$$

  (12)

  for almost all $\tau \in [0,T]$ and all $\phi \in {\mathcal {D}}((0,\tau )\times {\overline{\Omega }})$ with $\phi \cdot \text {n}|_{\partial \Omega } = 0$ and

• (iii) the initial conditions

  $$\begin{aligned} \rho (0) = \rho _0,\quad \lim _{\tau \rightarrow 0+} \int \limits _\Omega \rho (\tau ,x)u(\tau ,x) \cdot \phi (x) \ dx = \int \limits _\Omega q(x) \cdot \phi (x)\ dx \end{aligned}$$

  (13)

  for all $\phi \in {\mathcal {D}}({\overline{\Omega }})$ with $\phi \cdot {\text {n}}|_{\partial \Omega } = 0$.

In order to make sure that Definition (3.1) is a suitable definition of weak solutions, we show that any classical solution to the system (1)–(5) is also a weak solution and, vice versa, any weak solution with a sufficient amount of regularity solves problem (1)–(5) in the classical sense. In order to obtain the variational inequality (12) from system (1)–(5), we first pick an arbitrary time $\tau \in [0,T]$ and multiply the momentum equation (2) by an arbitrary function $\phi \in {\mathcal {D}}((0,\tau ) \times {\overline{\Omega }})$ with $\phi \cdot {\text {n}}|_{\partial \Omega } = 0$. Integrating (by parts) over $(0,\tau )\times \Omega$, we obtain the identity

$$\begin{aligned} \int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - (\rho u \otimes u) : \nabla \phi \ + \sigma : \nabla \phi \ dxdt =\int \limits _0^\tau \int \limits _\Omega \rho f \cdot \phi \ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } (\sigma \text {n})_\tau \cdot \phi \ d\Gamma dt. \end{aligned}$$

(14)

Similarly, we test the momentum equation (2) by u and subtract from it the continuity equation (1) tested by $\frac{1}{2}|u|^2$. Hence, we infer the energy inequality

$$\begin{aligned}&\frac{1}{2} \int \limits _\Omega \rho (\tau ) |u(\tau )|^2 + \frac{a\rho ^\gamma (\tau )}{\gamma - 1} \ dx + \int \limits _0^\tau 2\nu \left| {\mathbb {D}}(u) \right| ^2 + \lambda \left| \nabla \cdot u \right| ^2\ dxdt \\&\quad = \frac{1}{2} \int \limits _\Omega \rho (0) \left| u(0) \right| ^2 + \frac{a\rho ^\gamma (0)}{\gamma - 1} \ dx + \int \limits _{0}^\tau \int \limits _\Omega \rho f \cdot u\ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } \left( \sigma \text {n} \right) _{\tau } \cdot u \ d\Gamma dt. \end{aligned}$$

The last identity is subtracted from equation (14), and then we use the boundary condition (5). Finally, we obtain

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2}\rho (0) \left| u(0)\right| ^2 + \frac{a\rho ^\gamma (0)}{\gamma - 1} \ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau ) \left| u(\tau ) \right| ^2 + \frac{a\rho ^\gamma (\tau )}{\gamma - 1} \ dx + \int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - (\rho u \otimes u) : \nabla \phi \ \nonumber \\ {}&\qquad + \,\left( 2\nu {\mathbb {D}}(u) + \lambda \left( \nabla \cdot u \right) \, {\mathcal {I}} \right) : \left[ \nabla \phi - \nabla u \right] -p\, {\mathcal {I}} : \nabla \phi - \rho f \cdot \left[ \phi - u \right] \ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | - g |u| d\Gamma dt \nonumber \\ {}&\quad = \int \limits _0^\tau \int \limits _{\partial \Omega } (\sigma \text{ n})_\tau \cdot \left[ \phi - u \right] \ d\Gamma dt + \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | - g |u| d\Gamma dt \end{aligned}$$

(15)

$$\begin{aligned}&\quad = \int \limits _0^\tau \int \limits _{\partial \Omega } (\sigma \text{ n})_\tau \cdot \phi \ d\Gamma dt + \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | d\Gamma dt \ge 0, \end{aligned}$$

(16)

which is exactly the variational inequality (12). Conversely, we need to check that any sufficiently regular weak solution in the sense of Definition 3 also satisfies the system (1)–(5) in the classical sense. Hereof, continuity equation (1), the initial condition (3) and the boundary condition (4) are clear. For the derivation of the momentum equation, we test the variational inequality (12) by $\psi u \pm \phi$ for some arbitrary functions $\psi \in {\mathcal {D}}(0,T)$, $\phi \in {\mathcal {D}}((0,T)\times \Omega )$. Under exploitation of the assumed smoothness of $\rho$ and u, this yields the relation

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2}\rho (0) \left| u(0)\right| ^2 + \frac{a\rho ^\gamma (0)}{\gamma - 1} \ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau ) \left| u(\tau ) \right| ^2 + \frac{a\rho ^\gamma (\tau )}{\gamma - 1} \ dx + \int \limits _0^T \int \limits _\Omega \partial _t \left( \rho u \right) \cdot \left( \psi u \pm \phi \right) \ \\ {}&\quad + \,\left( \nabla \cdot \left( \rho u \otimes u \right) \right) \cdot \left( \psi u \pm \phi \right) + \left( \nabla \cdot \left( 2\nu {\mathbb {D}}(u) + \lambda \left( \nabla \cdot u \right) \, {\mathcal {I}} \right) \right) \cdot \left[ \left( 1 - \psi \right) u \pm \phi \right] + \nabla p \cdot \left( \psi u \pm \phi \right) \ \\ {}&\quad + \,\rho f \cdot \left[ \left( 1 - \psi \right) u \pm \phi \right] \ dxdt + \int \limits _0^T \int \limits _{\partial \Omega } g |\psi u| - g |u| d\Gamma dt \ge 0 \end{aligned}$$

Letting $\psi \rightarrow 1$, we see that the $\phi$-independent terms in this inequality cancel each other, and we are left with the equality

$$\begin{aligned}&\int \limits _0^T \int \limits _\Omega \partial _t \left( \rho u \right) \cdot \phi + \left( \nabla \cdot (\rho u \otimes u)\right) \cdot \phi \ - \left( \nabla \cdot \left( 2\nu {\mathbb {D}}(u) + \lambda \left( \nabla \cdot u \right) \, \mathcal {I} \right) \right) \cdot \phi + \nabla p \cdot \phi - \rho f \cdot \phi \ dxdt = 0. \end{aligned}$$

(17)

Hence, by the arbitrary choice of $\phi \in {\mathcal {D}}((0,T)\times \Omega )$, the classical formulation (2) of the momentum equation is satisfied. This in particular implies that identity (15) again holds true. Subtracting (15) from the given variational inequality (12), we find the estimate

$$\begin{aligned} -\int \limits _0^\tau \int \limits _{\partial \Omega } (\sigma \text {n})_\tau \cdot \left[ \phi - u \right] \ d\Gamma dt \le \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | - g |u| d\Gamma dt \end{aligned}$$

(18)

for any $\phi \in {\mathcal {D}}((0,\tau )\times {\overline{\Omega }})$ with $\phi \cdot \text {n}|_{\partial \Omega } = 0$ and, by a density argument, for any $\phi \in L^2(0,\tau ;H_{\text {n}}^1(\Omega ))$. We choose $s=T$ and test this inequality by $u \pm \phi$. Hence, replacing $\phi$ by $u \pm \phi$, we conclude, from the reverse triangle inequality, the estimate

$$\begin{aligned} \left| \int \limits _0^\tau \int \limits _{\partial \Omega } (\sigma \text {n})_\tau \cdot \phi \ d\Gamma dt \right| \le \int \limits _0^\tau \int \limits _{\partial \Omega } g |\phi | d\Gamma dt. \end{aligned}$$

It follows that $|(\sigma \text {n})_\tau |\le g$ on $\partial \Omega$ and further, from the estimate (18) with the choice $\phi = 0$ that $(\sigma \text {n})_\tau \cdot u + g |u| = 0$ on $\partial \Omega$. Hence, also the boundary condition (5) is satisfied.

We are now in the position to present the main result of our article, which is as follows:

Theorem 1.1

Let $T > 0$ and let $\Omega \subset {\mathbb {R}}^3$ be a bounded domain of class $C^{2,\eta } \bigcup C^{0,1}$ for some $\eta > 0$. Let the data $\nu , \lambda , a, \gamma \in {\mathbb {R}}$, $f \in L^\infty ((0,T)\times \Omega )$, $g \in L^2((0,T)\times \partial \Omega )$, $\rho _0 \in L^\gamma (\Omega )$ and $q \in L^1(\Omega )$ satisfy the conditions (6)–(7). Then, there exists a weak solution $(\rho , u)$, in the sense of Definition 3.1, to the system (1)–(5).

We remark that the $C^{2,\eta }$-regularity of $\Omega$ in Theorem 3.1 is necessary for the construction of the density in the approximate system in Sect. 4.1, c.f. [5, Lemma 3.1, Theorems 10.22, 10.23], [21, Proposition 7.39]. Moreover, the $C^1$-regularity of $\Omega$ is needed to extend u to an $L^2(0,T;H^1({\mathbb {R}}^3))$-function when showing that the couple $(\rho , u)$ satisfies the renormalized continuity equation in Sect. 4.4, c.f. [3, Sect. 5.4, Theorem 1].

4 Approximate system

In this section, we present an approximate version of the problem introduced in Sect. 2, followed by a brief explanation of the individual approximation levels. We fix four parameters $n \in {\mathbb {N}}$, $\delta , \epsilon , \alpha >0$, each of them associated to one of these approximation levels. We further fix some parameter $\beta > \max \lbrace \gamma , 4 \rbrace$. By $V_n \subset C^2({\overline{\Omega }};{\mathbb {R}}^3) \subset L^2(\Omega ;{\mathbb {R}}^3)$ we denote an n-dimensional vector space equipped with the $L^2(\Omega )$-inner product, such that

$$\begin{aligned} \bigcup _{n\in {\mathbb {N}}}V_n \ \ \text {is dense in} \ \ W_{\text {n}}^{1,p}(\Omega ) := \left\{ \phi \in W^{1,p}(\Omega ):\ \left. \phi \cdot \text {n} \right| _{\Gamma } = 0 \right\} \quad \forall 1 \le p < \infty . \end{aligned}$$

(19)

For technical reasons (c.f. the deduction of the convergence (65)), we assume that without loss of generality the sequence of spaces $(V_n)_n$ contains a subsequence $(V_{0,n})_n$ of spaces $V_{0,n}\subset C_0^2({\overline{\Omega }};{\mathbb {R}}^3)$ such that

$$\begin{aligned} \bigcup _{n\in {\mathbb {N}}}V_{0,n} \ \ \text {is dense in} \ \ W_0^{1,p}(\Omega ) \quad \forall 1 \le p < \infty . \end{aligned}$$

(20)

Moreover, following the approximation method used for the proof of the existence of weak solutions to the Navier–Stokes equations with the Coulomb friction law boundary condition in [1, Sect. 3], we denote by

$$\begin{aligned} j_\delta (v) := \left\{ \begin{matrix} |v| &{} \text {for } |v| > \delta , \\ \frac{|v|^2}{2\delta } + \frac{\delta }{2} &{} \text {for } |v| \le \delta , \end{matrix}\right. \end{aligned}$$

a convex approximation $j_\delta \in C^1({\mathbb {R}}^3)\bigcap C^{1,1}_{\text {loc}}({\mathbb {R}}^3)$ of the absolute value function. We remark that while in [1] the local Lipschitz continuity of the gradient of $j_\delta$ is not required, it is necessary in our setting in order to achieve continuity of the operator ${\mathbb {T}}$ in the fixed point argument for the construction of an approximate solution in Sect. 4.1. The approximation $j_\delta$ further has the properties

$$\begin{aligned}&j_\delta (0) = 0, \end{aligned}$$

(21)

$$\begin{aligned}&{\text {grad}} j_\delta (v) \cdot v \ge 0 \quad \forall \ v \in {\mathbb {R}}^3, \end{aligned}$$

(22)

$$\begin{aligned}&\left| {\text {grad}} j_\delta (v) \right| \le 1 \quad \forall \ {v} \in {\mathbb {R}}^3, \end{aligned}$$

(23)

$$\begin{aligned}&\left| j_\delta (v) - \left| v \right| \right| \le \delta \quad \forall \ v \in {\mathbb {R}}^3, \end{aligned}$$

(24)

where ${\text {grad}} j_\delta$ denotes the gradient of $j_\delta$. Our approximate problem on the highest approximation level consists of finding functions

$$\begin{aligned} \rho _\delta \in&W:= \left\{ \psi \in C\left( [0,T]; C^{2,\eta } \left( {\overline{\Omega }} \right) \right) \bigcap C^1\left( [0,T]; C^{0,\eta } \left( {\overline{\Omega }} \right) \right) :\ \left. \nabla \psi \cdot \text {n} \right| _{\Gamma } = 0 \right\} , \end{aligned}$$

(25)

$$\begin{aligned} u_\delta \in&C\left( [0,T];V_n \right) \end{aligned}$$

(26)

which satisfy the approximate continuity equation

$$\begin{aligned} \partial _t \rho _\delta + \nabla \cdot \left( \rho _\delta u_\delta \right) =&\epsilon \Delta \rho _\delta , \quad \left. \nabla \rho _\delta \cdot \text {n} \right| _{\Gamma } = 0. \end{aligned}$$

(27)

in $(0,T)\times \Omega$ and the approximate momentum equation

$$\begin{aligned} \int \limits _\Omega \partial _t \left( \rho _\delta u_{\delta } \right) \cdot \phi \ dx&= \int \limits _\Omega \left( \rho _\delta u_\delta \otimes u_{\delta } \right) : \nabla \phi - 2\nu {\mathbb {D}}\left( u_{\delta }\right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u_\delta ) (\nabla \cdot \phi ) + a\rho _\delta ^\gamma \nabla \cdot \phi + \alpha \rho _\delta ^\beta \nabla \cdot \phi \nonumber \\&\quad +\, \rho _\delta f\cdot \phi - \epsilon \left( \nabla u_{\delta } \nabla \rho _\delta \right) \cdot \phi \ dx - \int \limits _{\partial \Omega } g {\text {grad}} j_\delta (u_\delta ) \cdot \phi \ d\Gamma \end{aligned}$$

(28)

in [0, T] for all $\phi \in C([0,T];V_n)$ as well as the initial conditions

$$\begin{aligned} \rho (0,x) = \rho _0(x),\quad u(0,x) = u_0(x) \quad \forall x \in \Omega . \end{aligned}$$

(29)

Here, the initial data $u_0$ for the velocity field is defined by

$$\begin{aligned} u_0 := P_n \left( \frac{q}{\rho _0} \right) \in V_n, \end{aligned}$$

(30)

where $P_n$ denotes the orthogonal projection from $L^2(\Omega )$ onto $V_n$, and the initial data $\rho _0$, q is assumed to satisfy the additional regularity criteria

$$\begin{aligned} \rho _0 \in C^{2,\eta }\left( {\overline{\Omega }}\right) ,\quad 0 < \alpha \le \rho _0 \le \alpha ^{-\frac{1}{2\beta }},\quad \left. \nabla \rho _0 \cdot {\text {n}} \right| _{\Gamma } = 0,\quad q \in C^2\left( {\overline{\Omega }}\right) . \end{aligned}$$

(31)

Having introduced the full approximate problem (25)–(29), we now give a short explanation of the individual approximation levels in the order, in which we will later pass to the limit in them, beginning with the $\delta$-level. On this level, following the proof of the existence of weak solutions to the incompressible Navier–Stokes equations with the Coulomb friction law boundary conditions in [1], we add a boundary integral containing the quantity ${\text {grad}} j_\delta (u_\delta )$ to the momentum equation. The convexity of $j_\delta$ then allows us to transform the approximate momentum equation (28) into an inequality in which the desired boundary condition is incorporated in the same way as in the momentum and energy inequality (12) in our weak formulation.

The remaining approximation levels coincide precisely with the corresponding approximation levels in the classical theory of the existence of weak solutions to the compressible Navier–Stokes equations, which can be found for example in [21, Chapter 7]. On the n-level, we carry out a Galerkin approximation, which allows us to find a solution to the approximate system. More precisely, this procedure reduces the problem to a finite dimensional problem in the spatial component, which can be solved via the classical theory of ordinary differential equations and a fixed point argument. The reason why we pass to the limit with respect to $\delta \rightarrow 0$ before passing to the limit in the Galerkin approximation lies in the high spatial regularity available on the Galerkin level. This regularity allows us to achieve uniform convergence of the velocity field $u_\delta$ when letting $\delta$ tend to zero, which is required for passing to the limit in the quantity $j_\delta (u_\delta )$.

On the $\epsilon$-level the additional quantity $\epsilon \Delta \rho _\delta$ is added to the continuity equation. This procedure (c.f. [21, Section 7.6]), known as the parabolic regularization of the continuity equation, is required to make sure that the density in our approximate system and consequently also in our final system is non-negative. For the sake of preserving an energy inequality under this modification of the continuity equation, the term $\epsilon (\nabla u_\delta \nabla \rho _\delta )$ is moreover added to the momentum equation.

Lastly, we have the $\alpha$-level, on which we add the artificial pressure $\alpha \rho ^\beta$ is added to the momentum equation. The choice $\beta > \max \lbrace \gamma , 4 \rbrace$ provides us with a higher regularity of the density, which in turn allows us to pass to the limit in the quantity $\epsilon (\nabla u_\delta \nabla \rho _\delta )$ in the limit passage with respect to $n \rightarrow \infty$, see [21, Sect. 7.8.2].

4.1 Solution to the approximate problem

Our proof of the existence of a solution to the approximate problem (25)–(29) mainly follows the classical existence theory for the compressible Navier–Stokes equations (c.f. for example [21, Sect. 7.7]) with the difference lying only in the consideration of the additional boundary integral in the momentum equation (28). We start by fixing an arbitrary function $w \in C([0,T];V_n)$. Then by the classical theory for the parabolic Neumann problem (see [5, Lemma 3.1, Theorems 10.22, 10.23], [21, Proposition 7.39]) there exists a unique function $\rho = \rho (w) \in W$ which solves the problem

$$\begin{aligned}&\partial _t \rho + \nabla \cdot \left( \rho w \right) = \epsilon \Delta \rho \quad \text {in } [0,T] \times \Omega , \end{aligned}$$

(32)

$$\begin{aligned}&\rho (0, \cdot ) = \rho _0(\cdot )\quad \text {in } \Omega ,\ \quad 0< {\underline{\rho }} \le \rho _0(\cdot ) \le {\overline{\rho }} < \infty , \quad \text {in } \Omega \end{aligned}$$

(33)

and which in addition satisfies the estimate

$$\begin{aligned} 0< {\underline{\rho }} \exp \left( -\left\| w \right\| _{L^1(0,t;V_n)} \right) \le \rho (t,\cdot ) \le {\overline{\rho }} \exp \left( \left\| w \right\| _{L^1(0,t;V_n)} \right) < \infty \quad \text {in } {\overline{\Omega }} \end{aligned}$$

(34)

for all $t \in [0,T]$. Further, this solution satisfies the estimates

$$\begin{aligned}&\left\| \rho (w) \right\| _{C \left( [0,T];C^{2,\eta }\left( {\overline{\Omega }}\right) \right) } + \left\| \rho (w) \right\| _{C^1 \left( [0,T];C^{0,\eta }\left( {\overline{\Omega }}\right) \right) } \le c(w), \end{aligned}$$

(35)

$$\begin{aligned}&\left\| \rho \left( w^1\right) - \rho \left( w^2\right) \right\| _{C([0,T];L^2(\Omega ))} \le c\left( w^1,w^2 \right) \left\| w^1 - w^2 \right\| _{C([0,T];W^{1,\infty }(\Omega ))} \end{aligned}$$

(36)

for all $w,w^1,w^2 \in C([0,T];V_n)$, where the constants $c(w), c(w^1,w^2)>0$ are bounded as long as $w,w^1,w^2$ are bounded in the norm on $C([0,T];V_n)$. Moreover, due to the bound (34) of $\rho (w)$ away from 0, it is easy to see from the classical theory of ordinary differential equations that there exists a unique solution $u = u(w) \in C([0,T];V_n)$ to the associated linearized problem

$$\begin{aligned} \int \limits _\Omega \partial _t \left( \rho (w) u \right) \cdot \phi \ dx&= \int \limits _\Omega \left( \rho (w) w \otimes u \right) : \nabla \phi - 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u ) (\nabla \cdot \phi ) \nonumber \\&\quad +\, \left( a\rho (w)^\gamma + \alpha \rho (w)^\beta \right) \nabla \cdot \phi + \rho (w) f\cdot \phi - \epsilon \left( \nabla u \nabla \rho (w) \right) \cdot \phi \ dx \nonumber \\&\quad -\, \int \limits _{\partial \Omega } g {\text {grad}} j_\delta (w) \cdot \phi \ d\Gamma \quad \text {in } [0,T], \end{aligned}$$

(37)

$$\begin{aligned} u(0,\cdot )&= u_0(\cdot ) \quad \text {in } \Omega . \end{aligned}$$

(38)

This allows us to consider the desired solution $u_\delta$ to the momentum equation (28) as a fixed point of the operator

$$\begin{aligned} {\mathbb {T}}: C\left( [0,T];V_n\right) \rightarrow C\left( [0,T];V_n\right) ,\quad {\mathbb {T}}(w) := u(w), \end{aligned}$$

mapping $w \in C([0,T];V_n)$ to the corresponding solution to the linearized problem (37). The existence of such fixed point follows from the version [3, Sect. 9.2.2, Theorem 4] of the Schauder fixed point theorem. We show that ${\mathbb {T}}$ is continuous, compact and fixed points of $s{\mathbb {T}}$ are bounded in $C([0,T];V_n)$ uniformly with respect to $s \in [0,1]$. To this end, we introduce the operator

$$\begin{aligned} {\mathcal {M}}_{\rho (w)(t)}: V_n \rightarrow V_n^*,\quad \left\langle {\mathcal {M}}_{\rho (w)(t)}v,\phi \right\rangle _{V_n^* \times V_n} := \int \limits _\Omega \rho (w)(t) v \cdot \phi \ dx \quad \forall \phi , v \in V_n. \end{aligned}$$

The bound (34) of $\rho (w)$ away from zero implies the existence of an inverse ${\mathcal {M}}_{\rho (w)(t)}^{-1}$ of ${\mathcal {M}}_{\rho (w)(t)}$, with the properties

$$\begin{aligned}&\left\| {\mathcal {M}}^{-1}_{\rho (w)(t)} \right\| _{{\mathcal {L}}(V_n^*,V_n)} \le \frac{1}{\inf _{(0,T)\times \Omega }\rho (w)}, \end{aligned}$$

(39)

$$\begin{aligned}&\left\| {\mathcal {M}}^{-1}_{\rho (w^1)(t)} - {\mathcal {M}}^{-1}_{\rho (w^2)(t)} \right\| _{{\mathcal {L}}(V_n^*,V_n)} \le \frac{c(n)}{\left( \inf _{(0,T)\times \Omega } \min \left\{ \rho \left( w^1\right) , \rho \left( w^2\right) \right\} \right) ^2} \left\| \rho \left( w^1\right) (t) - \rho \left( w^2\right) (t) \right\| _{L^1(\Omega )}, \end{aligned}$$

(40)

as well as

$$\begin{aligned}&\partial _t \left\langle {\mathcal {M}}_{\rho (w)(t)}v(t),\phi \right\rangle _{V_n^* \times V_n} \nonumber \\&\quad = \left\langle {\mathcal {M}}_{\rho (w)(t)}^{-1}{\mathcal {M}}_{\partial _t \rho (w)(t)}{\mathcal {M}}_{\rho (w)(t)}^{-1}v(t) + {\mathcal {M}}_{\rho (w)(t)}^{-1} \partial _t v(t),\phi \right\rangle _{V_n^* \times V_n} \quad \text {in } {\mathcal {D}}'(0,T), \end{aligned}$$

(41)

$$\begin{aligned}&\left\| {\mathcal {M}}_{\rho (w)(t)}^{-1} {\mathcal {M}}_{\partial _ t\rho (w)(t)}{\mathcal {M}}_{\rho (w)(t)}^{-1} \right\| _{{\mathcal {L}}(V_n^*,V_n)} \le \frac{c(n)}{\left( \inf _{(0,T)\times \Omega }\rho (w)\right) ^2} \left\| \partial _t \rho (w)(t) \right\| _{L^1(\Omega )} \end{aligned}$$

(42)

for any $t \in [0,T]$, any $w,w^1,w^2,v \in C([0,T];V_n)$ and any $\phi \in V_n$, c.f. [21, Sect. 7.7.1]. Denoting

$$\begin{aligned} \left\langle {\mathcal {N}}\left( w,\rho , u \right) , \phi \right\rangle _{V_n^*\times V_n}&:= \int \limits _\Omega \left( \rho (w) w \otimes u \right) : {\mathbb {D}} (\phi ) + a\rho ^\gamma (w) \nabla \cdot \phi + \alpha \rho ^\beta (w) \nabla \cdot \phi \\&\quad -\, 2\nu {\mathbb {D}}\left( u \right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u)(\nabla \cdot \phi ) + \rho (w) f\cdot \phi - \epsilon \left( \nabla u \nabla \rho (w) \right) \cdot \phi \ dx \\&\quad -\, \int \limits _{\partial \Omega } g {\text {grad}} j_\delta (w) \cdot \phi \ d\Gamma , \\&\qquad \left\langle \left( \rho _0 u_0 \right) ^*, \phi \right\rangle _{V_n^*\times V_n} := \int \limits _\Omega \rho _0 u_0 \cdot \phi \ dx, \end{aligned}$$

for any $\phi \in V_n$, the solution $u = {\mathbb {T}}(w)$ to the linearized problem (37), (38) can be expressed as

$$\begin{aligned} u(t) = {\mathcal {M}}_{\rho (w)(t)}^{-1} \left[ \left( \rho _0 u_0 \right) ^* + \int \limits _{0}^t {\mathcal {N}}\left( w,\rho (w), u \right) \ d\tau \right] . \end{aligned}$$

(43)

Combining this identity with the estimates (34), (35), (39), (40) and the local Lipschitz continuity of ${\text {grad}} j_\delta$, we deduce that the operator ${\mathbb {T}}$ is continuous. Further, the combination of the identity (43) with the identity (41) and the estimates (34), (36), (39) and (42) leads to the estimate

$$\begin{aligned} \left\| \partial _t u \right\| _{L^2(0,T;V_n)}^2 \le c(n,w) \end{aligned}$$

(44)

with a constant $c(n,w)>0$ which remains bounded as long as w is bounded in the norm of $C([0,T];V_n)$. From this estimate, we infer that the operator ${\mathbb {T}}$ is also compact. Finally, we consider an arbitrary number $s \in [0,1]$ and an arbitrary fixed point $u \in C([0,T];V_n)$ of the operator $s{\mathbb {T}}$. We test the corresponding linearized momentum equation (37) by u and subtract from it the corresponding continuity equation (32), tested by $\frac{1}{2}|u|^2$. This yields the energy equality

$$\begin{aligned}&\frac{d}{dt} \int \limits _\Omega \frac{1}{2} \rho (u)(t)|u(t)|^2 + s \frac{a\rho ^\gamma (u)(t)}{\gamma -1} + s \frac{\alpha \rho ^\beta (u)(t)}{\beta - 1}\ dx + \int \limits _\Omega 2\nu |{\mathbb {D}} (u)(t)|^2 + \lambda \left| \nabla \cdot u \right| ^2\ dx \nonumber \\&\qquad +\, s\epsilon \gamma \int \limits _\Omega \rho ^{\gamma -2}(u) \left| \nabla \rho (u) \right| ^2\ dx + s\epsilon \beta \int \limits _\Omega \rho ^{\beta -2}(u) \left| \nabla \rho (u) \right| ^2\ dx + \int \limits _{\partial \Omega } sg {\text {grad}} j_\delta (u) \cdot u\ d\Gamma \nonumber \\&\quad = s\int \limits _\Omega \rho (u)(t)f(t) \cdot u(t)\ dx \end{aligned}$$

(45)

for all $t \in [0,T]$. According to the property (22) of $j_\delta$, it holds that ${\text {grad}} j_\delta (u) \cdot u \ge 0$ and, by assumption, g is nonnegative. Hence, from the Gronwall lemma, we deduce that all fixed points u of $s{\mathbb {T}}$ are bounded in the norm of $C([0,T];V_n)$, independently of s. This and the continuity as well as the compactness of ${\mathbb {T}}$ provides the conditions for the fixed point theorem [3, Sect. 9.2.2, Theorem 4], which proves the existence of a fixed point $u \in C([0,T];V_n)$ of ${\mathbb {T}}$. Setting $(\rho _\delta ,u_\delta )=(\rho (u),u)$, the pair $(\rho _\delta ,u_\delta )$ constitutes the desired solution to our approximate problem (25)–(29). Integrating the energy inequality (45), which $\rho _\delta , u_\delta$ satisfy for $s=1$, we have shown the following proposition:

Proposition 1.1

Let the conditions of Theorem 3.1 be satisfied, let $n \in {\mathbb {N}}$, $\delta , \epsilon , \alpha > 0$ and let $\beta > \max \lbrace 4, \gamma \rbrace$. Moreover, let $u_0 \in V_n$ be defined by (30) and assume $\rho _0$, q, defined by (7), to satisfy the additional regularity conditions (31). Then, there exists a solution $(\rho _\delta , u_\delta ) \in W \times C([0,T];V_n)$ to the approximate problem (25)–(29) which in addition satisfies the energy equality

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2} \rho _\delta (\tau )|u_\delta (\tau )|^2 + \frac{a\rho _\delta ^\gamma (\tau )}{\gamma -1} + \frac{\alpha \rho _\delta ^\beta (\tau )}{\beta - 1}\ dx + \int \limits _0^\tau \int \limits _\Omega 2\nu |{\mathbb {D}} \left( u_\delta \right) (t)|^2 + \lambda \left| \nabla \cdot u_\delta \right| ^2 + \epsilon \gamma \rho _\delta ^{\gamma -2} \left| \nabla \rho _\delta \right| ^2 \nonumber \\&\qquad + \,\epsilon \beta \rho _\delta ^{\beta -2} \left| \nabla \rho _\delta \right| ^2\ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } g {\text {grad}} j_\delta \left( u_\delta \right) \cdot u_\delta \ d\Gamma dt \nonumber \\&\quad = \int \limits _0^\tau \int \limits _\Omega \rho _\delta f \cdot u_\delta \ dxdt + \int \limits _\Omega \frac{1}{2} \rho _0|u_0|^2 + \frac{a\rho _0 ^\gamma }{\gamma -1} + \frac{\alpha \rho _0 ^\beta }{\beta - 1}\ dx \end{aligned}$$

(46)

for all $\tau \in [0,T]$.

4.2 Limit passage with respect to $\delta \rightarrow 0$

Our next goal is to pass to the limit in the regularization of the function j, i.e., the approximation parameter $\delta$ tend to zero. From the energy inequality (46), the equivalence of norms on the finite dimensional space $V_n$ and the estimates (34), (35) for the solution $\rho _\delta$ to the Neumann problem (32), (33) with $w = u_\delta$, we infer the uniform bounds

$$\begin{aligned} \left\| \rho _\delta \right\| _{C \left( [0,T];C^{2,\eta }\left( {\overline{\Omega }}\right) \right) } + \left\| \rho _\delta \right\| _{C^1 \left( [0,T];C^{0,\eta }\left( {\overline{\Omega }}\right) \right) } + \left\| \frac{1}{\rho _\delta } \right\| _{C([0,T]\times {\overline{\Omega }})} + \left\| u_{\delta } \right\| _{C([0,T];V_n)} \le c \end{aligned}$$

with a constant $c>0$ independent of $\delta$. In particular, the bound for $u_\delta$ implies that the bound (44) for $\partial _t u_\delta$ still holds true,

$$\begin{aligned} \left\| \partial _t u \right\| _{L^2(0,T;V_n)}^2 \le c \end{aligned}$$

(47)

with a constant $c>0$ independent of $\delta$. Consequently, making use of the Aubin–Lions Lemma, we may extract subsequences and find functions

$$\begin{aligned}&0 \le \rho \in \bigg \lbrace \psi \in C\left( [0,T];H^{1,2}(\Omega )\right) \bigcap C\left( [0,T];L^p(\Omega ) \right) \bigcap L^2 \left( 0,T; H^{2,2}(\Omega ) \right) : \nonumber \\&\partial _t \psi \in L^2 \left( (0,T)\times \Omega \right) ,\ \left. \nabla \psi \cdot {\text {n}} \right| _{\partial \Omega } = 0 \bigg \rbrace \quad \forall 1\le p < \infty , \nonumber \\&u \in \left\{ \phi \in C\left( [0,T];V_n \right) :\ \partial _t \phi \in L^2\left( 0,T;V_n \right) \right\} , \end{aligned}$$

(48)

such that

$$\begin{aligned}&\rho _{\delta } \rightarrow \rho \quad \text {in } C\left( [0,T];H^{1,2}(\Omega ) \right) \ \ \text {and} \ \ C\left( [0,T];L^p(\Omega ) \right) , \quad u_{\delta } \rightarrow u \ \ \quad \text {in } C\left( [0,T];V_n \right) , \end{aligned}$$

(49)

$$\begin{aligned}&\rho _{\delta } \rightharpoonup \rho \quad \text {in } L^2\left( 0,T;H^{2,2}(\Omega ) \right) ,\quad \partial _t \rho _{\delta } \rightharpoonup \partial _t \rho \quad \text {in } L^2\left( (0,T)\times \Omega \right) ,\quad \partial _t u_{\delta } \rightharpoonup \partial _t u \quad \text {in } L^2\left( 0,T;V_n \right) . \end{aligned}$$

(50)

Clearly, these convergences are sufficient to pass to the limit in the continuity equation (27) and infer that the limit functions $\rho$, u satisfy

$$\begin{aligned} \partial _t \rho + \nabla \cdot \left( \rho u \right) =&\epsilon \Delta \rho \quad \text {a.e. in } (0,T) \times \Omega . \end{aligned}$$

(51)

As in the weak formulation (12), we want to combine the momentum equation and the energy inequality into one single relation. To this end, we integrate the momentum equation (28) over $[0,\tau ]$ for some arbitrary $\tau \in [0,T]$ and subtract from it the energy equality (46). Further, we exploit the convexity and the $C^1$-regularity of $j_\delta$ to estimate

$$\begin{aligned} {\text {grad}} j_\delta \left( u_\delta \right) \cdot \left( \phi - u_\delta \right) \le j_\delta (\phi ) - j_\delta \left( u_\delta \right) , \end{aligned}$$

which allows us to bring the boundary integrals into the same form as in the weak formulation (12). Altogether we obtain the inequality

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2} \rho _0|u_0|^2 + a \frac{\rho _0^\gamma }{\gamma -1} + \frac{\alpha \rho _0 ^\beta }{\beta - 1}\ dx - \int \limits _\Omega \frac{1}{2} \rho _\delta (\tau )|u_\delta (\tau )|^2 + a \frac{\rho _\delta ^\gamma (\tau )}{\gamma -1} + \frac{\alpha \rho _\delta ^\beta (\tau )}{\beta - 1}\ dx \nonumber \\&\quad +\,\int \limits _0^\tau \int \limits _\Omega - \rho _\delta u_{\delta } \cdot \partial _t \phi - \left( \rho _\delta u_\delta \otimes u_{\delta } \right) : \nabla \phi + 2\nu {\mathbb {D}}\left( u_{\delta }\right) :{\mathbb {D}} (\phi - u_\delta ) + \lambda (\nabla \cdot u_\delta ) (\nabla \cdot (\phi - u_\delta )) \nonumber \\&\quad -\, a\rho _\delta ^\gamma \nabla \cdot \phi - \alpha \rho _\delta ^\beta \nabla \cdot \phi - \epsilon a \rho _\delta ^{\gamma - 2} \left| \nabla \rho _\delta \right| ^2 - \epsilon \beta \rho _\delta ^{\beta - 2} \left| \nabla \rho _\delta \right| ^2 + \epsilon \left( \nabla u_{\delta } \nabla \rho _\delta \right) \cdot \phi - \rho _\delta f\cdot (\phi - u_\delta )\ dxdt \nonumber \\&\quad +\, \int \limits _0^\tau \int \limits _{\partial \Omega } g j_\delta (\phi ) - gj_\delta (u_\delta ) \ d\Gamma dt \ge 0 \quad \forall \phi \in C_c^1 \left( (0,\tau );V_n \right) ,\ \tau \in [0,T]. \end{aligned}$$

(52)

Due to the uniform convergences (24) of $j_\delta$ and (49) of $u_\delta$, we can pass to the limit in the boundary integral,

$$\begin{aligned} \int \limits _0^\tau \int \limits _{\partial \Omega } g j_\delta (\phi ) - gj_\delta (u_\delta ) \ d\Gamma dt \rightarrow \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| \phi \right| - g \left| u \right| \ d\Gamma dt. \end{aligned}$$

The strong convergences (49) also allow us to pass to the limit in the remaining terms of the inequality (52). Hence, dropping the nonpositive quantity $-\epsilon a \rho _\delta ^{\gamma -2}|\nabla \rho _\delta |^2$ from the left-hand side of this inequality, we conclude that the limit functions $\rho$, u satisfy

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2} \rho _0|u_0|^2 + a \frac{\rho _0^\gamma }{\gamma -1} + \frac{\alpha \rho _0 ^\beta }{\beta - 1}\ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau )|u(\tau )|^2 + a \frac{\rho ^\gamma (\tau )}{\gamma -1} + \frac{\alpha \rho ^\beta (\tau )}{\beta - 1}\ dx \nonumber \\&\quad +\,\int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - \left( \rho u \otimes u \right) : \nabla \phi + 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi - u) + \lambda (\nabla \cdot u ) (\nabla \cdot (\phi - u)) \nonumber \\&\quad -\, a\rho ^\gamma \nabla \cdot \phi - \alpha \rho ^\beta \nabla \cdot \phi - \epsilon \beta \rho ^{\beta - 2} \left| \nabla \rho \right| ^2 + \epsilon \left( \nabla u \nabla \rho \right) \cdot \phi - \rho f\cdot (\phi - u)\ dxdt \nonumber \\&\quad +\, \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| \phi \right| - g \left| u \right| \ d\Gamma dt \ge 0 \quad \forall \phi \in C_c^1 \left( (0,\tau );V_n \right) ,\ \tau \in [0,T]. \end{aligned}$$

(53)

Remark 1.1

As a technical tool we will need in (65), (77) and (79) a limit version of the momentum equation (28) itself. In this limit equation, it will be sufficient to restrict ourselves to test functions vanishing on $\Gamma$. Using such test functions in the momentum equation (28), we see that the boundary integral vanishes, and we can pass to the limit to obtain the identity

$$\begin{aligned} -\int \limits _0^T\int \limits _\Omega \rho u \cdot \partial _t \phi \ dxdt&= \int \limits _0^T\int \limits _\Omega \left( \rho u \otimes u \right) : \nabla \phi - 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u ) (\nabla \cdot \phi ) + a\rho ^\gamma \nabla \cdot \phi \nonumber \\&\quad +\, \alpha \rho ^\beta \nabla \cdot \phi + \rho f\cdot \phi - \epsilon \left( \nabla u \nabla \rho \right) \cdot \phi \ dxdt \end{aligned}$$

(54)

for all $\phi \in C_c^1((0,T);V_n)$ such that $\phi |_{\Gamma } = 0$.

4.3 Limit passage with respect to $n \rightarrow \infty$

In this section, we pass to the limit in the Galerkin approximation, i.e., we let n tend to infinity. Choosing $\phi = 0$ in the momentum and energy inequality (53), we obtain a classical energy inequality. In combination with the classical regularity for the regularized continuity equation (51) (see for example [21, Lemmas 7.37, 7.38, Sect. 7.8.2]) this yields the uniform bounds

$$\begin{aligned}&\left\| \rho _n |u_n|^2 \right\| _{L^\infty (0,T;L^1(\Omega ))} + \left\| \rho _n \right\| _{L^\infty (0,T;L^\beta (\Omega ))} + \left\| u_n \right\| _{L^2(0,T;H^1(\Omega ))} \le c, \end{aligned}$$

(55)

$$\begin{aligned}&\epsilon ^\frac{1}{2} \left\| \rho _n^\frac{\beta }{2} \right\| _{L^2(0,T;H^1(\Omega ))} + \epsilon \left\| \nabla \rho _n \right\| _{L^{r}((0,T)\times \Omega )} + \epsilon \left\| \partial _t \rho _n \right\| _{L^{{\tilde{r}}}((0,T)\times \Omega )} + \epsilon ^2 \left\| \Delta \rho _n \right\| _{L^{{\tilde{r}}}((0,T)\times \Omega )} \le c \end{aligned}$$

(56)

for a constant $c>0$ independent of $n \in {\mathbb {N}}$, where we can choose

$$\begin{aligned} r := \frac{10\beta - 6}{3\beta + 3}> 2,\quad {\tilde{r}} := \frac{5\beta - 3}{4\beta } > 1, \end{aligned}$$

provided that $\beta > 4$. Interpolations between these bounds lead to the uniform bounds

$$\begin{aligned} \left\| \rho _n u_n \right\| _{L^\infty (0,T;L^\frac{2\beta }{\beta + 1}(\Omega ))} + \left\| \rho _n u_n \otimes u_n \right\| _{L^\frac{6\beta }{4\beta + 3}((0,T)\times \Omega )} + \epsilon ^{\frac{3}{5\beta }} \left\| \rho _n \right\| _{L^{\frac{5}{3}\beta }((0,T)\times \Omega )} \le c \end{aligned}$$

(57)

for another constant $c>0$ independent of n. Combining the bounds (55)–(57) with the Aubin–Lions Lemma, we may extract a subsequence and conclude the existence of functions $u \in L^2(0,T;H_{\text {n}}^1(\Omega ))$, $\overline{\rho u \otimes u} \in L^{\frac{6\beta }{4\beta +3}}((0,T)\times \Omega )$ and

$$\begin{aligned} 0 \le \rho \in&L^\infty \left( 0,T;L^\beta (\Omega )\right) \bigcap L^{2} \left( 0,T; H^1(\Omega ) \right) \bigcap L^{{\tilde{r}}}\left( 0,T; W^{2,{\tilde{r}}}(\Omega ) \right) \end{aligned}$$

with the properties

$$\begin{aligned} \partial _t \rho \in L^{{\tilde{r}}} \left( (0,T) \times \Omega \right) ,\quad \left. \nabla \rho \cdot {\text {n}} \right| _{\partial \Omega } = 0 \end{aligned}$$

(58)

such that

$$\begin{aligned}&\rho _n \rightarrow \rho \quad \text {in } L^\beta \left( Q \right) \ \ \text {and} \ \ L^2\left( 0,T;H^{1,2}(\Omega ) \right) ,\quad \rho _n \rightharpoonup \rho \quad \text {in } L^{{\tilde{r}}}\left( 0,T;W^{2,{\tilde{r}}}(\Omega ) \right) , \end{aligned}$$

(59)

$$\begin{aligned}&\ u_n \rightharpoonup u \quad \text {in } L^2 (0,T;H^{1,2} (\Omega )),\quad \partial _t \rho _n \rightharpoonup \partial _t \rho \quad \text {in } L^{{\tilde{r}}}\left( (0,T)\times \Omega \right) , \end{aligned}$$

(60)

$$\begin{aligned}&\rho _n u_n \buildrel *\over \rightharpoonup \rho u \quad \text {in } L^\infty \left( 0,T;L^\frac{2\beta }{\beta + 1}(\Omega )\right) ,\quad \rho _n u_n \otimes u_n \rightharpoonup \overline{\rho u \otimes u} \quad \text {in } L^\frac{6\beta }{4\beta + 3}\left( (0,T)\times \Omega \right) . \end{aligned}$$

(61)

With these convergences, we can directly pass to the limit in the continuity equation (51) and infer that

$$\begin{aligned} \partial _t \rho + \nabla \cdot \left( \rho u \right) =&\epsilon \Delta \rho \quad \text {a.e. in } (0,T) \times \Omega . \end{aligned}$$

(62)

This equation immediately implies that also the renormalized regularized continuity equation

$$\begin{aligned} \partial _t \zeta \left( \rho \right) + \nabla \cdot \left( \zeta \left( \rho \right) u \right) + \left[ \zeta ' \left( \rho \right) \rho - \zeta \left( \rho \right) \right] \nabla \cdot u - \epsilon \Delta \zeta \left( \rho \right) = - \epsilon \zeta '' \left( \rho \right) \left| \nabla \rho \right| ^2 \le 0 \quad \text {a.e. in } (0,T)\times \Omega \end{aligned}$$

(63)

holds true for all convex functions $\zeta \in C^2([0,+\infty ))$. In order to pass to the limit in the momentum and energy inequality (53) we need to identify the limit function $\overline{\rho u \otimes u}$ from the convergence (61). To this end, we test the alternative momentum equation (54) by $\psi \phi$ for some arbitrary functions $\psi \in {\mathcal {D}}(0,T)$ and $\phi \in V_N$, $N \le n$, such that $\phi |_\Gamma = 0$. Under exploitation of the uniform bounds (55)–(57), this leads us to the dual estimate

$$\begin{aligned}&\left\| \partial _t \int \limits _\Omega \rho _n u_n \cdot \phi \ dx \right\| _{L^{\min \left\{ \frac{6}{5}, \frac{2r}{2+r} \right\} }(0,T)} \le c \end{aligned}$$

for a constant $c>0$ depending on N but not on n. This allows us to infer from the Arzelà-Ascoli theorem that

$$\begin{aligned} \int \limits _\Omega \rho _n(\cdot , x) u_n (\cdot , x) \cdot \phi \ dx \rightarrow \int \limits _\Omega \rho (\cdot , x) u(\cdot , x) \cdot \phi \ dx \quad \text {in } C\left( [0,T] \right) \end{aligned}$$

(64)

for any fixed $\phi \in V_N$, $N \in {\mathbb {N}}$. Since the Galerkin spaces $V_N$ have been choosen such that the functions $\phi \in \bigcup _{N=1}^\infty V_N$ with $\phi |_\Gamma = 0$ are dense in $L^\frac{2\beta }{\beta +1}(\Omega )$ (c.f. (20)) and due to the continuity of the functions $\rho _n u_n$ with respect to the time variable (c.f. (49)) the convergence (64) suffices to infer that

$$\begin{aligned} \rho _n u_n \rightarrow \rho u \quad \text {in } C_{\text {weak}}\left( [0,T];L^\frac{5}{4}(\Omega ) \right) \ \ \text {and thus in}\ \ L^2\left( 0,T;H^{-1}(\Omega ) \right) , \end{aligned}$$

(65)

which is sufficient to identify, as desired,

$$\begin{aligned} \overline{\rho u \otimes u} = \rho u \otimes u \quad \text {a.e. in } (0,T) \times \Omega . \end{aligned}$$

(66)

For the limit passage in the boundary integrals, we note that by the weak convergence of $u_n$ in $L^2(0,T;H^1(\Omega ))$ and the trace theorem, $u_n$ also converges weakly in $L^2((0,T)\times \partial \Omega )$. Hence, the nonnegativity of $g \in L^2((0,T)\times \partial \Omega )$ and the weak lower semicontinuity of the $L^1((0,T)\times \partial \Omega )$-norm imply that

$$\begin{aligned} \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| u \right| \ d\Gamma dt \le \liminf _{n \rightarrow \infty } \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| u _n\right| \ d\Gamma dt \quad \forall \tau \in [0,T]. \end{aligned}$$

This, in combination with the convergences (59)–(61), the identification (66) of the weak limit of the convective term and the weak lower semicontinuity of norms, gives us all the ingredients required for passing to the limit in both the momentum and energy inequality (52) and the alternative momentum equation (54). Due to the density of the Galerkin functions in $W_{\text {n}}^{1,p}(\Omega )$ and $W_0^{1,p}(\Omega )$, $1 \le p < \infty$ [c.f. (19), (20)], we infer that

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2} \rho _0|u_0|^2 + a \frac{\rho _0^\gamma }{\gamma -1} + \frac{\alpha \rho _0 ^\beta }{\beta - 1}\ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau )|u(\tau )|^2 + a \frac{\rho ^\gamma (\tau )}{\gamma -1} + \frac{\alpha \rho ^\beta (\tau )}{\beta - 1}\ dx \nonumber \\&\quad +\,\int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - \left( \rho u \otimes u \right) : \nabla \phi + 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi - u) + \lambda (\nabla \cdot u ) (\nabla \cdot (\phi - u)) \nonumber \\&\quad - \,a\rho ^\gamma \nabla \cdot \phi - \alpha \rho ^\beta \nabla \cdot \phi - \epsilon \beta \rho ^{\beta - 2} \left| \nabla \rho \right| ^2 + \epsilon \left( \nabla u \nabla \rho \right) \cdot \phi - \rho f\cdot (\phi - u)\ dxdt \nonumber \\&\quad +\, \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| \phi \right| - g \left| u \right| \ d\Gamma dt \ge 0 \end{aligned}$$

(67)

holds true for almost all $\tau \in [0,T]$ and all $\phi \in {\mathcal {D}}((0,\tau )\times {\overline{\Omega }})$ with $\phi \cdot \text {n} |_{\partial \Omega } = 0$ and

$$\begin{aligned} -\int \limits _0^T\int \limits _\Omega \rho u \cdot \partial _t \phi \ dxdt&= \int \limits _0^T\int \limits _\Omega \left( \rho u \otimes u \right) : \nabla \phi - 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u ) (\nabla \cdot \phi ) + a\rho ^\gamma \nabla \cdot \phi \nonumber \\&\quad + \,\alpha \rho ^\beta \nabla \cdot \phi + \rho f\cdot \phi - \epsilon \left( \nabla u \nabla \rho \right) \cdot \phi \ dxdt, \end{aligned}$$

(68)

for all $\phi \in {\mathcal {D}}((0,T)\times \Omega )$.

4.4 Limit passage with respect to $\epsilon \rightarrow 0$

In this section, we consider the limit passage with respect to $\epsilon \rightarrow 0$ in order to get rid of the artificial regularization terms in the system. By setting $\phi = 0$ in the momentum and energy inequality (67) and a subsequent interpolation we infer, exactly as the corresponding bounds (55) and (57) in the previous limit passage, the uniform bounds

$$\begin{aligned}&\left\| \rho _\epsilon |u_\epsilon |^2 \right\| _{L^\infty (0,T;L^1(\Omega ))} + \left\| \rho _\epsilon \right\| _{L^\infty (0,T;L^\beta (\Omega ))} + \left\| u_\epsilon \right\| _{L^2(0,T;H^1(\Omega ))} \le c, \end{aligned}$$

(69)

$$\begin{aligned}&\left\| \rho _\epsilon u_\epsilon \right\| _{L^\infty (0,T;L^\frac{2\beta }{\beta + 1}(\Omega ))} + \left\| \rho _\epsilon u_\epsilon \otimes u_\epsilon \right\| _{L^\frac{6\beta }{4\beta + 3}((0,T)\times \Omega )} \le c \end{aligned}$$

(70)

for a constant $c>0$ independent of $\epsilon$. These bounds allow us to extract a subsequence and conclude the existence of functions

$$\begin{aligned} 0 \le \rho \in L^\infty \left( 0,T;L^\beta (\Omega ) \right) ,\quad u\in L^2\left( 0,T;H_{\text {n}}^1(\Omega )\right) \end{aligned}$$

(71)

such that

$$\begin{aligned} \rho _\epsilon \buildrel *\over \rightharpoonup \rho \ \ \ \text {in } L^\infty \left( 0,T;L^\beta (\Omega )\right) , \quad u_\epsilon \rightharpoonup u \ \ \ \text {in } L^2 \left( 0,T;H^{1,2}(\Omega )\right) . \end{aligned}$$

(72)

Under exploitation of continuity equation (62), the first one of these convergences further leads to

$$\begin{aligned} \rho _\epsilon \rightarrow \rho \quad \text {in } C_{\text {weak}} \left( [0,T];L^\beta (\Omega ) \right) \quad \text {and hence} \quad \rho _\epsilon u_\epsilon \rightharpoonup \rho u \quad \text {in } L^\infty \left( 0,T;L^\frac{2\beta }{\beta +1}(\Omega ) \right) . \end{aligned}$$

(73)

Moreover, we may test continuity equation (62) by $\rho _\epsilon$ to infer that

$$\begin{aligned} \epsilon ^\frac{1}{2} \left\| \nabla \rho _\epsilon \right\| _{L^2((0,T)\times \Omega )} \le c \quad \text {and thus} \quad \epsilon \nabla \rho _\epsilon \rightarrow 0 \quad \text {in } L^2\left( (0,T)\times \Omega \right) . \end{aligned}$$

(74)

Convergences (72)–(74) allow us to pass to the limit in the continuity equation (27) and infer that $\rho$ and u satisfy the continuity equation (9) in ${\mathcal {D}}'((0,T)\times \Omega )$. Due to the Lipschitz regularity of $\partial \Omega$, u can be extended continuously to a function $u \in L^2(0,T;H^1({\mathbb {R}}^3))$, see [3, Section 5.4, Theorem 1]. Extending also $\rho$ by 0 outside of $\Omega$, we infer that the continuity equation in fact holds true in ${\mathcal {D}}'((0,T)\times {\mathbb {R}}^3)$. From the regularization method by DiPerna and Lions (see [21, Theorem 6.9]), it follows that $\rho$ and u also satisfy the renormalized continuity equation (10), (11). Similar to the convergence of the convective term in the Galerkin limit [c.f. (61), (66)], we may deduce from the alternative momentum equation (68) that

$$\begin{aligned} \rho _\epsilon u_\epsilon \otimes u_\epsilon \rightharpoonup \rho u \otimes u \quad \text {in } L^\frac{6\beta }{4\beta + 3}\left( (0,T)\times \Omega \right) . \end{aligned}$$

(75)

In order to pass to the limit in the pressure terms, we need to find a uniform bound for the density in $L^q((0,T)\times \Omega )$ for some $q>\beta$. Thanks to the alternative momentum equation (68) such bound can be derived as in the case of the no-slip boundary condition, c.f. [6, Lemma 3.1]. Nameley, since the Bogovskii operator ${\mathcal {B}}_\Omega$ on $\Omega$ (c.f. [21, Section 3.3.1.2]) maps zero-mean functions in $L^p(\Omega )$, $1< p < \infty$, into $W_0^{1,p}(\Omega )$, we can test the alternative momentum equation (68) by functions of the form

$$\begin{aligned} \phi _\epsilon (t,x) := \psi (t) {\mathcal {B}}_\Omega \left[ \rho _\epsilon (t) - \frac{1}{|\Omega |} \int \limits _\Omega \rho _\epsilon (t,y)\ dy \right] (x),\quad 0 \le \psi \in {\mathcal {D}}\left( 0,T \right) . \end{aligned}$$

(76)

As ${\mathcal {B}}_\Omega$ can be understood as an inverse to the divergence operator, this allows us to find a constant $c>0$ independent of $\epsilon$ such that

$$\begin{aligned} \left\| \rho _\epsilon \right\| _{L^{\gamma + 1}((0,T) \times \Omega )} + \left\| \rho _\epsilon \right\| _{L^{\beta + 1}((0,T) \times \Omega )} \le c. \end{aligned}$$

(77)

Thus, we find subsequences and functions $\overline{\rho ^\gamma } \in L^\frac{\gamma + 1}{\gamma } ((0,T)\times \Omega )$, $\overline{\rho ^\beta } \in L^\frac{\beta + 1}{\beta } ((0,T)\times \Omega )$ such that

$$\begin{aligned} \rho _\epsilon ^\gamma \rightharpoonup \overline{\rho ^\gamma } \quad \text {in } L^\frac{\gamma + 1}{\gamma }\left( (0,T)\times \Omega )\right) ,\quad \rho ^\beta _\epsilon \rightharpoonup \overline{\rho ^\beta } \quad \text {in } L^\frac{\beta + 1}{\beta }\left( (0,T)\times \Omega \right) . \end{aligned}$$

(78)

Our next goal is to identify the limit functions $\overline{\rho ^\gamma }$ and $\overline{\rho ^\beta }$, for which we need the effective viscous flux identity

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} \int \limits _{0}^T \int \limits _{\Omega } \Phi (\lambda + 2\nu ) \left[ \rho _\epsilon \nabla \cdot u_\epsilon - \rho \nabla \cdot u \right] \ dxdt = \lim _{\epsilon \rightarrow 0} \int \limits _{0}^T \int \limits _{\Omega } \Phi \left( \left[ a\rho _\epsilon ^\gamma + \alpha \rho _\epsilon ^\beta \right] \rho _\epsilon - \left[ a\overline{\rho ^\gamma } + \alpha \overline{\rho ^\beta } \right] \rho \right) \ dxdt \end{aligned}$$

(79)

for all $0\le \Phi \in {\mathcal {D}}((0,T)\times \Omega )$. This identity can be proved by applying the method from Feireisl et al. [6, Lemma 3.2] to the alternative momentum equation (68). We test (68) and a corresponding limit identity, obtained from the convergences (72), (73) and (78), by functions of the form

$$\begin{aligned} \phi _\epsilon (t,x) := \Phi (t,x) \left( \nabla \Delta ^{-1} \right) \left[ \rho _\epsilon (t,\cdot ) \right] (t,x),\quad \phi (t,x) := \Phi (t,x) \left( \nabla \Delta ^{-1} \right) \left[ \rho (t,\cdot ) \right] (t,x), \end{aligned}$$

(80)

with $0 \le \Phi \in {\mathcal {D}}((0,T)\times \Omega )$, respectively. Subtracting the two resulting relations from each other we obtain the effective viscous flux identity exactly as in [6, Lemma 3.2]. Following the procedure in [6, Sect. 3.5], we consider—after a dominated convergence argument—both the renormalized continuity equation (63) on the $\epsilon$-level and the renormalized continuity equation (10) in the limit with the choice of the (strictly) convex function $\zeta (\xi ):= \xi \ln (\xi )$ in. This results in two relations which we subtract from each other to obtain the inequality

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} \int \limits _\Omega \rho (\tau ) \ln \left( \rho (\tau ) \right) \ - \rho _\epsilon (\tau ) \ln \left( \rho _\epsilon (\tau ) \right) \ dx \ge \lim _{\epsilon \rightarrow 0} \int \limits _0^\tau \int \limits _\Omega \rho _\epsilon \nabla \cdot u_\epsilon - \rho \nabla \cdot u \ dxdt \quad \forall \tau \in [0,T]. \end{aligned}$$

(81)

From the effective viscous flux identity (79) and the monotonicity of the (artificial) pressure function, it follows that the right-hand side of this relation is nonnegative. Further, since the mapping $\xi \rightarrow \xi \ln (\xi )$ is convex, we know that $\rho \ln (\rho ) \le \overline{\rho \ln (\rho )}$, where $\overline{\rho \ln (\rho )}$ denotes a weak limit of $\overline{\rho _\epsilon \ln (\rho _\epsilon )}$ in $L^1((0,T)\times \Omega )$. Combining these two facts, we conclude that

$$\begin{aligned} \rho \ln (\rho ) = \overline{\rho \ln (\rho )} \quad \text {a.e. in } (0,T) \times \Omega . \end{aligned}$$

By the relations between weakly convergent sequences and (strictly) convex functions (c.f. [5, Theorem 10.20]), this equation implies pointwise convergence of $\rho _\epsilon$, which in turn implies that, as desired,

$$\begin{aligned} \overline{\rho ^\gamma } = \rho ^\gamma \quad \text {a.e. in } (0,T)\times \Omega ,\quad \overline{\rho ^\beta } = \rho ^\beta \quad \text {a.e. in } (0,T)\times \Omega . \end{aligned}$$

(82)

Combining convergences (72)–(75), (78), the identification (82) of the limits of the pressure terms and the weak lower semicontinuity of norms, we can now pass to the limit in both the momentum and energy inequality (67) and the alternative momentum equation (68) and infer that

$$\begin{aligned}&\int \limits _\Omega \frac{1}{2} \rho _0|u_0|^2 + a \frac{\rho _0^\gamma }{\gamma -1} + \frac{\alpha \rho _0 ^\beta }{\beta - 1}\ dx - \int \limits _\Omega \frac{1}{2} \rho (\tau )|u(\tau )|^2 + a \frac{\rho ^\gamma (\tau )}{\gamma -1} + \frac{\alpha \rho ^\beta (\tau )}{\beta - 1}\ dx \nonumber \\&\quad +\,\int \limits _0^\tau \int \limits _\Omega - \rho u \cdot \partial _t \phi - \left( \rho u \otimes u \right) : \nabla \phi + 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi - u) + \lambda (\nabla \cdot u ) (\nabla \cdot (\phi - u)) \nonumber \\&\quad -\, a\rho ^\gamma \nabla \cdot \phi - \alpha \rho ^\beta \nabla \cdot \phi - \rho f\cdot (\phi - u)\ dxdt + \int \limits _0^\tau \int \limits _{\partial \Omega } g \left| \phi \right| - g \left| u \right| \ d\Gamma dt \ge 0 \end{aligned}$$

(83)

holds true for almost all $\tau \in [0,T]$ and all $\phi \in {\mathcal {D}}((0,\tau )\times {\overline{\Omega }})$ with $\phi \cdot \text {n}|_{\partial \Omega } = 0$ and

$$\begin{aligned} -\int \limits _0^T\int \limits _\Omega \rho u \cdot \partial _t \phi \ dxdt&= \int \limits _0^T\int \limits _\Omega \left( \rho u \otimes u \right) : \nabla \phi - 2\nu {\mathbb {D}}\left( u\right) :{\mathbb {D}} (\phi ) - \lambda (\nabla \cdot u ) (\nabla \cdot \phi ) + a\rho ^\gamma \nabla \cdot \phi \nonumber \\&\quad +\, \alpha \rho ^\beta \nabla \cdot \phi + \rho f\cdot \phi dxdt, \end{aligned}$$

(84)

holds true for all $\phi \in {\mathcal {D}}((0,T)\times \Omega )$.

4.5 Limit passage with respect to $\alpha \rightarrow 0$

Finally, it remains to get rid of the artificial pressure term in the momentum equation, i.e., to let $\alpha$ tend to zero. In addition, we return from the regularized initial data $\rho _{0,\alpha }$, $q_{\alpha }$ in the approximate problem [c.f. (31)] to the more general initial data $\rho _0$, q from the main result Theorem 3.1. More precisely, as in [6, Sect. 4], we choose $\rho _{0,\alpha }$, $q_{\alpha }$ satisfying the relations (31) for any fixed $\alpha > 0$ such that

$$\begin{aligned}&\rho _{0,\alpha } \rightarrow \rho _0 \quad \text {in } L^\gamma (\Omega ),\quad \alpha \rho _{0,\alpha }^\beta \rightarrow 0 \quad \text {in } L^1(\Omega ), \end{aligned}$$

(85)

$$\begin{aligned}&q_{\alpha } \rightarrow q \quad \text {in } L^1 (\Omega ),\quad \frac{\left| q_{\alpha }\right| ^2}{\rho _{0,\alpha }} \rightarrow \frac{\left| q\right| ^2}{\rho _{0}} \quad \ \ \ \text {in } L^1(\Omega ) \end{aligned}$$

(86)

for $\alpha \rightarrow 0$. As in the previous limit passages, we infer, from the choice $\phi = 0$ in the momentum and energy inequality (83) and an ensuing interpolation, the uniform bounds

$$\begin{aligned}&\left\| \rho _\alpha |u_\alpha |^2 \right\| _{L^\infty (0,T;L^1(\Omega ))} + \left\| \rho _\alpha \right\| _{L^\infty (0,T;L^\gamma (\Omega ))} + \left\| u_\alpha \right\| _{L^2(0,T;H^1(\Omega ))} \le c, \end{aligned}$$

(87)

$$\begin{aligned}&\left\| \rho _\alpha u_\alpha \right\| _{L^\infty (0,T;L^\frac{2\gamma }{\gamma + 1}(\Omega ))} + \left\| \rho _\alpha u_\alpha \otimes u_\alpha \right\| _{L^\frac{6\gamma }{4\gamma + 3}((0,T)\times \Omega )} \le c \end{aligned}$$

(88)

for a constant $c>0$ independent of $\alpha$. This allows us to find a subsequence as well as functions

$$\begin{aligned} 0 \le \rho \in L^\infty \left( 0,T;L^\gamma (\Omega ) \right) ,\quad u\in L^2\left( 0,T;H_{\text {n}}^1(\Omega )\right) \end{aligned}$$

(89)

such that

$$\begin{aligned} \rho _\alpha \buildrel *\over \rightharpoonup \rho \quad \text {in } L^\infty \left( 0,T;L^\gamma (\Omega )\right) , \quad u_\alpha \rightharpoonup u \quad \text {in } L^2 \left( 0,T;H^{1,2}(\Omega )\right) \end{aligned}$$

(90)

as well as, under exploitation of the continuity equation (9) and the alternative momentum equation (84),

$$\begin{aligned} \rho _\alpha \rightarrow \rho \quad \text {in } C_{\text {weak}} \left( [0,T];L^\gamma (\Omega ) \right) ,\quad \rho _\alpha u_\alpha \rightarrow \rho u \quad \text {in } C_{\text {weak}} \left( [0,T];L^\frac{2\gamma }{\gamma + 1} (\Omega ) \right) \end{aligned}$$

(91)

and consequently

$$\begin{aligned} \rho _\alpha u_\alpha \otimes u_\alpha \rightharpoonup \rho u \otimes u \quad \text {in } L^\frac{6\gamma }{4\gamma + 3}((0,T)\times \Omega ). \end{aligned}$$

(92)

Due to the convergences (91) and the continuity equation (9) on the $\alpha$-level, the limit functions $\rho$ and u satisfy the same continuity equation (9) in ${\mathcal {D}}'((0,T)\times \Omega )$. For the limit passage in the pressure terms the derivation of the improved uniform bounds (77) of the density on the $\epsilon$-level needs to be modified. More specifically, the derivation of these bounds relies on the fact that on the $\epsilon$-level the density is bounded uniformly in $L^\infty (0,T;L^\beta (\Omega ))$, which is not the case in our current situation. As a compensation, the density in the test functions (76) needs to be replaced by a suitable smooth approximation of $\rho _\alpha ^\theta$ for some sufficiently small value $\theta > 0$. This procedure, which is described in detail in [6, Section 4.1], leads, by a use of the resulting test functions in the alternative momentum equation (84), to the desired improved pressure estimates

$$\begin{aligned} \left\| \rho _\alpha \right\| _{L^{\gamma + \theta }((0,T) \times \Omega )} + \alpha ^\frac{1}{\beta + \theta }\left\| \rho _\alpha \right\| _{L^{\beta + \theta }((0,T) \times \Omega )} \le c \end{aligned}$$

with a constant $c>0$ independent of $\alpha$. In particular, we may extract a subsequence and find a function $\overline{\rho ^\gamma } \in L^\frac{\gamma + 1}{\gamma }((0,T)\times \Omega )$ such that

$$\begin{aligned} \rho _\alpha ^\gamma \rightharpoonup \overline{\rho ^\gamma } \quad \text {in } L^\frac{\gamma + \theta }{\gamma }\left( (0,T)\times \Omega \right) ,\quad \alpha \rho ^\beta _\alpha \rightarrow 0 \quad \text {in } L^\frac{\beta + \theta }{\beta }\left( (0,T)\times \Omega \right) . \end{aligned}$$

(93)

For the identification of the limit function $\overline{\rho ^\gamma }$, we further follow the procedure in [6, Sect. 4.3] and deduce the following modified version of the effective viscous flux identity (79) on the $\epsilon$-level,

$$\begin{aligned}&\lim _{\alpha \rightarrow 0} \int \limits _{0}^T \int \limits _{\Omega } \Phi (\lambda + 2\nu ) \left[ T_k\left( \rho _\alpha \right) \nabla \cdot u_\alpha - \overline{T_k\left( \rho \right) } \nabla \cdot u \right] \ dxdt \nonumber \\&\quad = \lim _{\alpha \rightarrow 0} \int \limits _{0}^T \int \limits _{\Omega } \Phi \left( a\rho _\alpha ^\gamma T_k \left( \rho _\alpha \right) - a\overline{\rho ^\gamma }\ \overline{T_k\left( \rho \right) } \right) \ dxdt \end{aligned}$$

(94)

for all $\Phi \in {\mathcal {D}}((0,T)\times \Omega )$, where $T_k \le 2k$, $k \in {\mathbb {N}}$, constitutes a suitable smooth and concave cut-off version of the identity function on $[0,\infty )$ and $\overline{T_k(\rho )}$ denotes a weak limit of $T_k(\rho _\alpha )$ in $L^1((0,T)\times \Omega )$. In the derivation of this identity, we again have to make up for the lower integrability of the density as compared to the density on the $\epsilon$-level. This is achieved by replacing the test functions (80) used on the $\epsilon$-level by test functions of the form

$$\begin{aligned} \phi _\alpha (t,x) := \Phi (x) \left( \nabla \Delta ^{-1} \right) \left[ T_k\left( \rho _\alpha \right) (t,\cdot ) \right] (x),\quad \phi (t,x) := \Phi (x) \left( \nabla \Delta ^{-1} \right) \left[ \overline{ T_k\left( \rho \right) }(t,\cdot ) \right] (x). \end{aligned}$$

where $\Phi \in {\mathcal {D}}((0,T)\times \Omega )$. In our case, we use these test functions in the alternative momentum equation (84) on the $\alpha$-level and a corresponding limit identity, respectively. Comparing the resulting identities, we arrive at the desired effective viscous flux identity (94), exactly as in the proof of [6, Lemma 4.2]. From this identity, the concavity of $T_k$ and the convexity of $\xi \mapsto \xi ^\gamma$ we deduce, exactly as in the proof of [6, Lemma 4.3], boundedness of the oscillation defect measure,

$$\begin{aligned} {\textbf {osc}}_{\gamma + 1}\left[ \rho _\alpha \rightarrow \rho \right] \left( (0,T)\times {\mathbb {R}}^3 \right) := \sup _{k \ge 1} \left[ \limsup _{\alpha \rightarrow 0} \int \limits _0^T \int \limits _{{\mathbb {R}}^3} \left| T_k \left( \rho _\alpha \right) - T_k\left( \rho \right) \right| ^{\gamma + 1}\ dxdt \right] < \infty . \end{aligned}$$

(95)

Next, we choose $\zeta = T_k$ in the renormalized continuity equation (10) on the $\alpha$-level and pass to the limit with respect to $\alpha$. Since $\overline{T_k(\rho )}$ is bounded uniformly, the resulting limit identity can be renormalized under exploitation of the regularization technique by DiPerna and Lions, see [21, Lemma 6.9]. Subsequently, using the bound (95) of the oscillation defect measure, we may let k tend to infinity to infer that $\rho$ and u also satisfy the renormalized continuity equation (10). For the details of this procedure, we refer to the proof of [6, Lemma 4.4]. Under exploitation of the dominated convergence theorem, we may use the choice

$$\begin{aligned} \zeta (\xi ) := \zeta _k(\xi ) := \xi \int \limits _1^\xi \frac{T_k(s)}{s^2}\ ds \end{aligned}$$

in both the renormalized continuity equation (10) on the $\alpha$-level and in the limit. Comparing the resulting equations to each other and passing to the limit with respect to $\alpha \rightarrow 0$, we infer that

$$\begin{aligned}&\lim _{\alpha \rightarrow 0}\int \limits _{\Omega } \left( \zeta _k(\rho _\alpha ) - \zeta _k(\rho ) \right) (\tau )\ dx + \lim _{\alpha \rightarrow 0} \int \limits _0^\tau \int \limits _{\Omega } T_k\left( \rho _\alpha \right) \nabla \cdot u_\alpha - \overline{T_k(\rho )} \nabla \cdot u \ dxdt \nonumber \\&\quad =\int \limits _0^\tau \int \limits _{\Omega } T_k(\rho ) \nabla \cdot u - \overline{T_k(\rho )} \nabla \cdot u \ dxdt \end{aligned}$$

(96)

for all $\tau \in [0,T]$. Here, the second term on the left-hand side is nonnegative, which follows from the effective viscous flux identity (94), the fact that both the mappings $\xi \mapsto \xi ^\gamma$ and $\xi \mapsto T_k(\xi )$ are nondecreasing and the classical relations between weakly convergent sequences and monotone functions (c.f. [5, Theorem 10.19]). Moreover, the right-hand side of the equation (96) vanishes for $k \rightarrow \infty$ as can be seen from the bound (95) of the oscillation defect measure. Consequently, letting k tend to infinity also on the left-hand side of this identity, we infer that

$$\begin{aligned} \lim _{\alpha \rightarrow 0} \int \limits _\Omega \rho (\tau ) \ln \left( \rho (\tau ) \right) \ - \rho _\alpha (\tau ) \ln \left( \rho _\alpha (\tau ) \right) \ dx \ge 0. \end{aligned}$$

Exactly as in the limit passage with respect to $\epsilon \rightarrow 0$ [c.f. (82)], this estimate yields pointwise convergence of $\rho _\alpha$ and consequently the identity $\overline{\rho ^\gamma } = \rho ^\gamma$ almost everywhere in $(0,T)\times \Omega$. Therefore, using the convergences (90)–(92) and (93) as well as the weak lower semicontinuity of norms we may pass to the limit in the momentum and energy inequality (83) and infer that $\rho$ and u satisfy the momentum and energy inequality (12). Finally, we note that $\rho$, as a solution to the renormalized continuity equation (10), is an element of the space $C([0,T];L^1(\Omega ))$, c.f. [4, Proposition 4.3]. Due to this continuity in the time variable and the convergence (85) of the initial data it satisfies the initial condition $\rho = \rho _0$ in the classical sense. The initial condition stated for $\rho u$ in (13) follows from the convergences (86) and (91) This concludes the proof of Theorem 3.1.