# The Brinkman-Fourier System with Ideal Gas Equilibrium

Jan-Eric Sulzbach111Department of Applied Mathematics, Illinois Institute of
Technology, Chicago, IL 60616, United States Chun Liu∗222corresponding author:
<EMAIL_ADDRESS>
###### Abstract

In this work, we will introduce a general framework to derive the
thermodynamics of a fluid mechanical system, which guarantees the consistence
between the energetic variational approaches with the laws of thermodynamics.
In particular, we will focus on the coupling between the thermal and
mechanical forces. We follow the framework for a classical gas with ideal gas
equilibrium and present the existences of weak solutions to this thermodynamic
system coupled with the Brinkman-type equation to govern the velocity field.

##  1 Introduction

The Navier-Stokes-Fourier system, known for its power in modeling the
thermodynamics of fluids, is an extension of the classical Navier-Stokes
equations. Although a lot can be modeled in the classical setting, it has its
limitations encompassing of reality. More recent versions have included a
thermodynamic component in order to better describe more complex system.
Applications of this can be found in engineering, meteorology and even
astrophysics [Bir94], [HLB96], [CV02], [FN05]. This has opened up avenues for
researchers to update and expand upon previous applications of the Navier-
Stokes-Fourier system.

More recent advances and methods for the full Navier-Stokes-Fourier system can
be found in [FN09] and the references therein. Results in the field typically
focus on proving existence and uniqueness of solutions. Extensive theory on
weak solutions to the Navier-Stokes-Fourier system has been developed in
[FN05a], [Fei07], [FN12], and is continuously gaining interest among
researchers today. The existence of weak solutions can even be extended to
other domains, such as when Lipschitz boundaries are present [Pou07]. Many
mathematicians investigate long-term behavior of solutions, as in [FP07], in
order to gain information about equilibrium states, bounds on the energy, etc.

For the equations of viscous and heat-conductive gases, it can be shown that
unique and global strong solutions exist [Dan01]. Other results, focusing on
different aspects of the thermodynamics of fluids are the study of the gas
dynamics of thermal non-equilibrium models [Zen99], [Zen10]. Simplified models
with temperature dependent coefficients [Ann15] and other thermal effects such
as cross diffusion [BH15] and heat convection [NPT13] have been studied to see
the effect of the temperature dependence in the system of equations. Other
types of non-isothermal models are discussed in [ERS15], [AL19] and [LWL18].
Most recent results were obtained by [Tar19] and [LLT20] deriving better a
priori bounds and showing positivity for the absolute temperature and by
[LS21], where the well-posed of a non-isothermal reaction-diffusion sytem in a
critical Besov space is shown. In [GBY] the authors present a different free
energy formulation of the Navier-Stokes-Fourier system. Other works in this
field with similar approaches can be found in [Bul+20] and [Dre+20]. For the
Poisson-Nerst-Planck-Fourier system a global existence theorem could be proven
[Hsi+20].

The new and general approach of this paper is that we start from a given free
energy function and dissipation functional. Next, employing the basic laws of
thermodynamic in combination with the energetic variational approach and
appropriate constitutive relations, we derive a system of partial differential
equations describing the thermodynamics of the fluid. Due to the generality of
this approach it can be easily adapted to different regimes.

In this paper we will focus on the ideal gas equilibrium, i.e. the equation of
state for the internal energy and the pressure depends linearly on the density
and absolute temperature. For the momentum equation we choose a Brinkman-type
diffusion, i.e an interpolation between the Stokes and Darcy’s law. The
details can be found in Section 2.

The aim of this paper is to prove existence of local-in-time weak solutions of
the following system of partial differential equations

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0,$ (1.1)
$\displaystyle\nabla p=\mu\Delta u-\nu\rho u,$ (1.2)
$\displaystyle\partial_{t}s+\operatorname{div}(su)=\Delta+\operatorname{div}\big{(}\frac{\kappa\nabla\theta}{\theta}\big{)},$
(1.3)

where $\rho$ is the density, $u$ is the velocity and $\theta$ the absolute
temperature. The entropy production rate is denoted by $\Delta$ and depends on
$\rho$, $\theta$ and $u$. Moreover, the pressure is defined as
$p=k_{2}\rho\theta$ and corresponds to the ideal gas case and the dissipation
part in the momentum equation corresponds to the Brinkman-type equation. The
entropy is defined as $-s=k_{2}\rho(\log\rho-\frac{3}{2}(\log\theta+1))$. In
addition, we aim at a better understanding of the stability and dependence on
the initial data and second law of thermodynamics as expressed in [Daf79].

In Section 3 we obtain the weak formulation of this system consisting of
equations (1.1)-(1.3) and show that the absolute temperature is positive.
Moreover, we derive a priori bounds to the equations.

In Section 4 we state the existence theorem 4.1 and the higher regularity
theorem 4.11. In the proofs we follow the ideas and techniques presented in
[Fei04] and [FN09]. Key parts are the application of the div-curl lemma and
weak $L^{1}$-convergence results to control the absolute temperature in the
entropy equation.

##  2 Derivation of the System

In this section we derive the general model describing the thermal effects of
a fluid with the example of the ideal gas case. For the fluid model we
consider a Brinkman-type equation.
The unknown variables in the system are:

  1. 1.

a non-negative measurable function $\rho=\rho(t,x)$ the mass density;

  2. 2.

a vector field $u=u(t,x)$ denoting the velocity field of the fluid;

  3. 3.

a positive measurable function $\theta=\theta(t,x)$ the absolute temperature.

The new approach in the derivation of thermodynamic models focuses on the free
energy of the system $\psi(\rho,\theta)$ as a starting point and then applies
the laws of thermodynamics and the energetic variational approach to obtain
the complete model.

The notion of free energy is a useful concept in classical thermodynamics
[Bai99], [Sal01] because the change in the free energy is the maximum amount
of work that a thermodynamic system can perform in a process at constant
temperature, and its sign indicates whether a process is thermodynamically
favorable or forbidden. In the following derivation we use the definition of
the Helmholtz free energy. The second thermodynamic concept we use is the
entropy. Entropy can be understood as the measure of disorder of the system.
Note that for fixed density the entropy is a convex function in the
temperature and that for a fixed temperature the entropy is a convex function
in the density.

In the derivation we try to keep the statements as general as possible, since
we can apply this framework to different settings, i.e. the porous media case
or the Allen-Cahn and Cahn-Hilliard model.

From the thermodynamics of the ideal gas, [McQ76] and [BRR00], we know that
the relation between internal energy and the product of temperature and
density is linear and similar for the relation between the pressure and the
product of temperature and density. Working backwards from this observation
have the following definitions.

For the ideal gas we have the following definition of the free energy

$\displaystyle\psi(\rho,\theta)=k_{2}\theta\rho\log\rho-
k_{1}\rho\theta\log\theta.$ (2.1)

The entropy of a system is defined as follows

$\displaystyle s(\rho,\theta):=-\psi_{\theta}=-\rho\big{(}k_{2}\log\rho-
k_{1}(\log\theta+1)\big{)},$ (2.2)

where $\psi_{\theta}$ denotes the derivative of $\psi$ with respect to
$\theta$.

######  Remark 2.1.

We note that the $\theta\log\theta$ in the free energy is the weakest convex
function with faster than linear growth. Thus we have a non-vanishing
contribution of the temperature in the entropy.

Moreover, the convexity of the free energy with respect to the temperature
allows us to solve the equation of the entropy (2.2) for $\theta$. For the
given choice of the free energy this can be done explicitly

$\displaystyle\theta(\rho,s)=\frac{1}{e}\rho^{k_{2}/k_{1}}e^{s/k_{1}\rho}$

whereas in the general case this can be only done implicitly.
The internal energy is the Legendre transform in $\theta$ of the free energy,
i.e

$\displaystyle
e(\rho,\theta):=\psi-\psi_{\theta}\theta=\psi+s\theta=k_{1}\rho\theta.$ (2.3)

######  Remark 2.2.

Instead of having $(\rho,\theta)$ as state variables, we can also consider
$(\rho,s)$ as new state variables yielding

$\displaystyle
e_{1}(\rho,s)=e(\rho,\theta(\rho,s))=\frac{k_{1}}{e}\rho^{1+k_{2}/k_{1}}e^{s/k_{1}\rho}.$

This interpretation of the free energy is crucial because in the laws of
thermodynamics we have $\rho$ and $s$ as state variables.

We assume that we have a closed system, i.e. $\rho$ satisfies the following
continuity equation

$\displaystyle\rho_{t}+\operatorname{div}(\rho u)=0$ (2.4)

Thus integrating over the domain $\Omega$ and assuming that the microscopic
velocity $u$ satisfies $u\cdot n=0$ on the boundary we obtain the conservation
of mass

$\displaystyle\frac{d}{dt}\int_{\Omega}\rho(t,x)dx=0.$ (2.5)

In addition, we assume that $\theta$ moves with the particle and that the
kinematics of the temperature $\theta$ are governed by a transport equation
for the temperature. In other words, the temperature $\theta$ is transported
along the trajectory of the flow map with the velocity $u$.

######  Remark 2.3.

The assumption of the transport of the temperature $\theta$ with the particles
is suitable for situations like the ideal gas, where the temperature is not
determined by the environment or the case of a solution, where the temperature
is determined by the temperature of the solvent.

######  Remark 2.4.

We note that we have a weak duality of the time evolution of the temperature
and the total derivative of the entropy in the following way.
If $\theta$ evolves as $\frac{d}{dt}\theta=\theta_{t}+u\cdot\nabla\theta$ then
by testing this equation with the entropy $s$ in the weak form yields

$\displaystyle\int_{\Omega}\theta_{t}s+u\cdot\nabla\theta
s\,dx=-\int_{\Omega}s_{t}\theta+\operatorname{div}(su)\theta\,dx.$

Thus the kinematics of the entropy are

$\displaystyle\frac{d}{dt}s=s_{t}+\operatorname{div}(su)$ (2.6)

which are crucial in applying the laws of thermodynamics c.f. equations (2.15)
and (2.16).

Next, we choose the total energy and dissipation as follows

$\displaystyle
E^{tot}=\int_{\Omega}\psi(\rho,\theta),~{}~{}\mathcal{D}^{tot}=\int_{\Omega}\nu\rho
u^{2}+\mu|\nabla u|^{2}$

and employ the energetic variational approach to derive the forces [GKL17].

######  Remark 2.5.

Note that dissipation depends on both $u$ and $\nabla u$ and hence $\mu$ and
$\nu$ can be seen as interpolation parameters between the two pure cases.

Using the least action principle we have

$\displaystyle
A(x(t))=\int_{0}^{T}\mathcal{L}dt=-\int_{0}^{T}E^{tot}dt=-\int_{0}^{T}\int_{\Omega}\psi(\rho(x,t),\theta(x,t))dxdt,$

where $A$ denotes the action and $\mathcal{L}$ is the Lagrangian of the
action. Since we have no contribution of the kinetic part in the total energy
the variation of the action yields

$\displaystyle\delta_{x}A(x)=-\int_{0}^{T}\int_{\Omega}f_{cons}\delta_{x}dxdt.$

The next step is to compute the variation of the action, where we first
rewrite it in Lagrangian coordinates

$\displaystyle
A(x(X,t))=-\int_{0}^{T}\int_{\Omega_{0}^{X}}\psi\bigg{(}\frac{\rho_{0}(X)}{\det
F},\theta_{0}(X)\bigg{)}\det Fdxdt,$

where $F$ denotes the deformation gradient $F=\frac{\partial x}{\partial X}$.
We set $y(X,t)=\tilde{y}(x(X,t),t)$ and the variation $x^{\epsilon}=x+\epsilon
y$, $F^{\epsilon}=\frac{\partial x^{\epsilon}}{\partial X}$. Then

$\displaystyle\frac{d}{d\epsilon}\bigg{|}_{\epsilon=0}A(x^{\epsilon}(X,t))$
$\displaystyle=-\frac{d}{d\epsilon}\bigg{|}_{\epsilon=0}\int_{0}^{T}\int_{\Omega_{0}^{X}}\psi\bigg{(}\frac{\rho_{0}(X)}{\det
F^{\epsilon}}\theta_{0}(X)\bigg{)}\det F^{\epsilon}dXdt$
$\displaystyle=-\int_{0}^{T}\int_{\Omega_{0}^{X}}\frac{\partial}{\partial\phi}\psi\bigg{(}\frac{\rho_{0}(X)}{\det
F},\theta_{0}(X)\bigg{)}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\bigg{[}\frac{-\rho_{0}}{(\det F)^{2}}\det
F\operatorname{tr}\big{(}F^{-1}\frac{\partial y}{\partial
X}\big{)}\bigg{]}\det FdXdt$
$\displaystyle-\int_{0}^{T}\int_{\Omega_{0}^{X}}\psi\bigg{(}\frac{\rho_{0}(X)}{\det
F},\theta_{0}(X)\bigg{)}\bigg{[}\det
F\operatorname{tr}\big{(}F^{-1}\frac{\partial y}{\partial
X}\big{)}F^{-1}\bigg{]}dXdt.$ We transform the integral back to Eulerian
coordinates and obtain
$\displaystyle=-\int_{0}^{T}\int_{\Omega_{t}^{x}}-\frac{\partial}{\partial\rho}\psi\big{(}\rho,\theta\big{)}\rho\nabla_{x}\cdot\tilde{y}+\psi\big{(}\rho,\nabla_{x}\phi,\theta\big{)}\nabla_{x}\cdot\tilde{y}dxdt.$
Integration by parts yields
$\displaystyle=-\int_{0}^{T}\int_{\Omega_{t}^{x}}\bigg{[}\nabla_{x}\bigg{(}\rho\frac{\partial\psi}{\partial\rho}\bigg{)}-\nabla_{x}\psi\bigg{]}\cdot\tilde{y}dxdt$
$\displaystyle~{}~{}~{}-\int_{0}^{T}\int_{\partial\Omega_{t}^{x}}\bigg{(}-\rho\frac{\partial\psi}{\partial\rho}+\psi\bigg{)}\tilde{y}\cdot
ndS_{x}dt$ where $\tilde{y}\cdot n=0$ and thus the boundary terms equate to 0.
Putting everything together we have $\displaystyle\delta_{x}A(x)$
$\displaystyle=-\int_{0}^{T}\int_{\Omega_{t}^{x}}\nabla\big{(}\psi_{\rho}\rho-\psi\big{)}\cdot\tilde{y}\,dxdt.$

for an arbitrary smooth vector $\tilde{y}(x,t)$ satisfying $\tilde{y}\cdot
n=0$ on $\partial\Omega$.
Thus we have

$\displaystyle f_{cons}=\nabla\big{(}\rho\psi_{\rho}-\psi\big{)}=:\nabla p,$

where we define the gradient of the pressure to equal the conservative force.
By this definition the pressure law satisfies the following relation

$\displaystyle p(\rho,\theta)=\psi_{\rho}\rho-\psi=k_{2}\rho\theta.$ (2.7)

######  Lemma 2.6.

The pressure satisfies

$\displaystyle\nabla p=\rho\nabla\psi_{\rho}+s\nabla\theta.$

###### Proof.

From the definition of the pressure we have
$p(\rho,\theta)=\psi_{\rho}\rho-\psi$ and thus we compute

$\displaystyle\nabla p(\rho,\theta)$
$\displaystyle=\nabla(\psi_{\rho}\rho-\psi)=\rho\nabla\psi_{\rho}+\psi_{\rho}\nabla\rho-\nabla\psi$
$\displaystyle=\rho\nabla\psi_{\rho}+\psi_{\rho}\nabla\rho-\psi_{\rho}\nabla\rho-\psi_{\theta}\nabla\theta=\rho\nabla\psi_{\rho}+s\nabla\theta.$

∎

######  Remark 2.7.

From classical thermodynamics and statistical mechanics we know that the ideal
gas law is

$\displaystyle pV=NkT$

where $V$ is the volume, $N$ the number of particles and $k$ the Boltzmann
constant [Sal01]. Using tools from statistical mechanics, we can derive that

$\displaystyle S(N,V,E)\equiv S(N,V^{2/3}E)$

and as a consequence for a reversible adiabatic process we obtain

$\displaystyle V^{2/3}E=const.$

This relates the internal energy and the pressure as follows

$\displaystyle e(\rho,\theta)=\frac{3}{2}p(\rho,\theta).$

Thus we obtain a relation between the two constants $k_{1}$ and $k_{2}$

$\displaystyle k_{1}=\frac{3}{2}k_{2}.$

Following the maximum dissipation law, we now compute the variation of the
total dissipation

$\displaystyle\frac{1}{2}\delta_{u}\mathcal{D}=\int_{\Omega}\big{(}\nu\rho
u-\mu\Delta u\big{)}\tilde{u}dx=\int_{\Omega}f_{diss}\tilde{u}.$

Using the classical Newton’s force balance

$\displaystyle f_{cons}+f_{diss}=0$

yields a Brinkman-type equation [Bri47], [DB87], which interpolates between
the Darcy’s law and the Stokes equation

$\displaystyle\nabla p=\mu\Delta u-\nu\rho u.$ (2.8)

######  Remark 2.8.

An overview over the energetic variational approach can be found in [LLZ05].
More recent approaches in adding the thermodynamics into the framework of the
energetic variational approach can be found in [LWL18] and [AL19].

The next step is to verify the physicality of the given free energy. Thus,
combining the above results, we can show that the Gibb’s equation is satisfied

$\displaystyle\theta
Ds(\rho,\theta)=De(\rho,\theta)+p(\rho,\theta)D\bigg{(}\frac{1}{\rho}\bigg{)},$
(2.9)

where we use the notation $Df$ to denote a total differential. Moreover we can
interpret this in terms of classical thermodynamics as follows

$\displaystyle E=TS-pV~{}~{}\text{ and }DE=TDS-pDV.$

For further details in classical thermodynamics we refer to [McQ76], [BRR00]
and [Sal01].

Next, we provide two useful basic Lemmas.

######  Lemma 2.9.

The internal energy as a function of the density and entropy satisfies

$\displaystyle{e_{1}}_{s}(\rho,s)=\theta$

###### Proof.

Note that $\psi_{1}(\rho,s)=\psi(\rho,\theta)$. Then,
$e_{1}(\rho,s)=e(\rho,\theta)=\psi+s\theta$ and thus

${e_{1}}_{s}={\psi_{1}}s+\theta+s\theta_{s}=\psi_{\theta}\theta_{s}+\theta+s\theta_{s}=\theta.$

∎

######  Lemma 2.10.

The internal energy in terms of density and entropy is related to the free
energy in terms of density and temperature

$\displaystyle\psi_{\rho}(\rho,\theta)={e_{1}}_{\rho}(\rho,s).$

###### Proof.

We compute the right hand side and obtain

$\displaystyle{e_{1}}_{\rho}(\rho,s)=\psi_{\rho}+\psi_{\theta}\theta_{\rho}+s\theta_{\rho}=\psi_{\rho}.$

∎

Now, we state the laws of thermodynamic. The first law of thermodynamics
relates the total derivative of the internal energy with work and heat

$\displaystyle De_{1}=\text{ work }+\text{ heat},$ (2.10)

where both the work and heat are in divergence form and more specific the heat
term can be expressed as $\nabla\cdot q$ with $q$ being the heat flux.
We recall that the total change in a quality can be expressed as a divergence
term plus some additional term. Applying this to the entropy we obtain

$\displaystyle s_{t}+\operatorname{div}(su)=\nabla\cdot j+\Delta,$ (2.11)

where the kinematics for the entropy are related to the transport of the
temperature as expressed in equation (2.6). And where we denote $j$ as the
entropy flux and $\Delta$ as the entropy production rate.

The second law of thermodynamics states that the entropy production is non-
negative:

$\displaystyle\Delta\geq 0.$ (2.12)

In order to derive a consistent model we need to supplement our equations by
constitutive relations. The relation between the heat flux $q$ and the entropy
flux $j$ is given by Durhem relation

$\displaystyle j\theta=q,$ (2.13)

and the relation between the heat flux $q$ and the absolute temperature
$\theta$ is given by the Fourier’s law

$\displaystyle q=\kappa\nabla\theta,$ (2.14)

where $\kappa$ is the heat conductivity.

In order to find the expression for the entropy production rate $\Delta$ we
compute

$\displaystyle\frac{d}{dt}\int_{\Omega}e_{1}(\rho,s)dx$
$\displaystyle=\int_{\Omega}\big{[}e_{1,\rho}\rho_{t}+e_{1,s}s_{t}\big{]}dx$
(2.15) Using the kinematics for the density $\rho$ from equation (2.4) we
obtain $\displaystyle=\int_{\Omega}\big{[}e_{1,\rho}\big{(}-\nabla\cdot(\rho
u)\big{)}+e_{1,s}s_{t}\big{]}dx$ Applying Lemma 2.10 yields
$\displaystyle=\int_{\Omega}\big{[}-\nabla\cdot\big{(}e_{1,\rho}\rho
u\big{)}+\big{(}\nabla\psi_{\rho}\rho\big{)}\cdot u+e_{1,s}s_{t}\big{]}dx$ In
order to have the full expression for the gradient of the pressure we have to
incorporate the term $s\nabla\theta$ which can only occur if the kinematics
for the entropy are as in equation (2.6). And by equation (2.11) we have
$\displaystyle=\int_{\Omega}\big{[}-\nabla\cdot\big{(}e_{1,\rho}\rho
u+e_{1,s}su\big{)}+\big{(}\nabla\psi_{\rho}\rho+s\nabla e_{1,s}\big{)}\cdot u$
$\displaystyle~{}~{}~{}~{}~{}+e_{1,s}\big{(}\nabla\cdot
j+\Delta\big{)}\big{]}dx$ By Lemma 2.9 and the Durhem equation (2.13) we have
$\displaystyle=\int_{\Omega}\big{[}-\nabla\cdot\big{(}(e_{1,\rho}\rho+\theta
s)u\big{)}+\big{(}\nabla\psi_{\rho}\rho+s\nabla\theta\big{)}\cdot u$
$\displaystyle~{}~{}~{}~{}~{}+\nabla\cdot
q-\frac{q\cdot\nabla\theta}{\theta}+\theta\Delta\big{]}dx$ Now, we can apply
Lemma 2.6 to obtain
$\displaystyle=\int_{\Omega}\big{[}-\nabla\cdot\big{(}e_{1,\rho}\rho u+\theta
su\big{)}+\nabla\cdot q+\nabla(\psi_{\rho}\rho-\psi)\cdot u$
$\displaystyle~{}~{}~{}~{}~{}-\frac{q\cdot\nabla\theta}{\theta}+\theta\Delta\big{]}dx$
From the definition of the pressure and the absence of external forces and
heat sources we have that $\displaystyle=\int_{\Omega}\big{[}\nabla p\cdot
u-\frac{q\cdot\nabla\theta}{\theta}+\theta\Delta\big{]}dx$ where we used that
the divergence terms equal to zero under the boundary conditions $u\cdot n=0$
and $\nabla\theta\cdot n=0$. Thus we have
$\displaystyle=\int_{\Omega}\big{[}\big{(}\mu\Delta u-\nu\rho u\big{)}\cdot
u-\frac{q\cdot\nabla\theta}{\theta}+\theta\Delta dx\big{]}$ and integration by
parts yields $\displaystyle=\int_{\Omega}\big{[}-\mu|\nabla u|^{2}-\nu\rho
u^{2}-\frac{q\cdot\nabla\theta}{\theta}+\theta\Delta\big{]}dx$ (2.16)

where we used that by equation (2.8) $\nabla p=\mu\Delta u-\nu\rho u$. Since
there are no external forces or heat sources in our system the total internal
energy must be conserved and we obtain that

$\displaystyle\Delta=\frac{1}{\theta}\bigg{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\frac{\kappa|\nabla\theta|^{2}}{\theta}\bigg{)}.$ (2.17)

We note that the second law of thermodynamics $\Delta\geq 0$ is satisfied as
long as $\theta>0$. Thus

$\displaystyle\frac{d}{dt}\int_{\Omega}e_{1}(\rho,s)dx=0$ (2.18)

and the total entropy is increasing

$\displaystyle\frac{d}{dt}\int_{\Omega}s(\rho,\theta)dx=\int_{\Omega}s_{t}+\operatorname{div}(su)dx=\int_{\Omega}\operatorname{div}j+\Delta\geq
0.$ (2.19)

Combining the above results we are now able to state our system of equations.
The model equations are

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)$ $\displaystyle=0$
$\displaystyle\partial_{t}s+\operatorname{div}(su)$
$\displaystyle=\operatorname{div}\bigg{(}\frac{\kappa\nabla\theta}{\theta}\bigg{)}+\Delta$
$\displaystyle\nabla p$ $\displaystyle=\mu\Delta u-\nu\rho u,$

where

$\displaystyle\Delta=\frac{1}{\theta}\bigg{(}|\nabla
u|^{2}+\nu\rho|u|^{2}+\frac{\kappa|\nabla\theta|^{2}}{\theta}\bigg{)}$

together with the state equations

$\displaystyle-s(\rho,\theta)$
$\displaystyle=k\rho\big{(}\log\rho-\frac{3}{2}(\log\theta+1)\big{)}$
$\displaystyle p(\rho,\theta)$ $\displaystyle=k\rho\theta.$

######  Remark 2.11.

The generality of this new approach in the derivation of the dynamics of
thermodynamic systems extends beyond the simple ideal gas case with
applications in the porous media equation or the Cahn-Hilliard equation, where
the only adaption to the new model takes place in the free energy.

##  3 Preliminaries

Starting from the derivation in the previous section we develop a
thermodynamically consistent mathematical model based on the unknown variables
$(\rho,\theta,u)$ satisfying the following properties

  1. 1.

the problem admits a local-in-time solution for any initial data of finite
energy;

  2. 2.

the total energy of the system remains constant in the absence of external
forces or heat sources;

  3. 3.

the entropy of the system is increasing, i.e. the system evolves to a state
maximizing the entropy;

  4. 4.

weak solutions coincide with classical solutions provided they are smooth
enough.

###  3.1 Weak Formulation

Now, let us summarize the weak formulation of the problem derived in the
previous chapter. Moreover, we also specify the minimal regularity of the
solutions required and interpret the weak formulation in terms of partial
differential equations provided that all quantities are smooth enough.

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded Lipschitz domain, where
$n=2,\,3$. Then we have the following equations.
The continuity equation

Weak formulation:
$\displaystyle\int_{0}^{T}\int_{\Omega}\rho\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}\,dxdt$
$\displaystyle=-\int_{\Omega}\rho_{0}\phi(0,\cdot)\,dx,$ (3.2) where $\phi\in
C_{c}^{1}([0,T)\times\overline{\Omega})$. Minimal regularity of solutions
required: $\displaystyle\rho\geq 0,~{}~{}\rho\in
L^{1}((0,T)\times\Omega),~{}~{}\rho u$ $\displaystyle\in
L^{1}((0,T)\times\Omega).$ (3.3) Formal interpretation:
$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)$
$\displaystyle=0~{}~{}\text{in }~{}(0,T)\times\Omega,$ (3.4)
$\displaystyle\rho(0,\cdot)$ $\displaystyle=\rho_{0},~{}~{}u\cdot
n|_{\partial\Omega}=0.$ (3.5)

Brinkman-type equation

Weak formulation:
$\displaystyle\int_{\Omega}k\rho\theta(t)\operatorname{div}\varphi\,dx$
$\displaystyle=\int_{\Omega}\mu\nabla u:\nabla\varphi+\nu\rho
u\cdot\varphi\,dx$ (3.7) where $\varphi\in C_{c}^{1}(\overline{\Omega})$ and
$\varphi|_{\partial\Omega}=0$. Minimal regularity of solutions required:
$\displaystyle\rho\theta\in L^{1}(\Omega),~{}~{}$ $\displaystyle\nabla u\in
L^{1}(\Omega),~{}~{}\rho u\in L^{1}(\Omega).$ (3.8) Formal interpretation:
$\displaystyle k\nabla(\rho\theta)$ $\displaystyle=\mu\Delta u-\nu\rho
u~{}~{}\text{in }~{}\Omega,$ (3.9) $\displaystyle u|_{\partial\Omega}$
$\displaystyle=0.$ (3.10)

Balance of internal energy

Weak formulation:
$\displaystyle\int_{0}^{T}\int_{\Omega}\rho\theta(t)\partial_{t}\psi(t)\,dxdt$
$\displaystyle=\psi(0)E_{0}$ (3.12) where $\psi\in C_{c}^{1}[0,T)$. Minimal
regularity of solutions required: $\displaystyle\rho\theta\in
L^{1}((0,T)\times$ $\displaystyle\Omega).$ (3.13) Formal interpretation:
$\displaystyle\frac{d}{dt}\int_{\Omega}\rho\theta\,dx=0~{}~{}\text{in
}(0,T),~{}~{}\int_{\Omega}\rho_{0}\theta_{0}\,dx$ $\displaystyle=E_{0}.$
(3.14)

Entropy production

Weak formulation:
$\displaystyle\begin{split}\int_{0}^{T}\int_{\Omega}s\big{(}\partial_{t}+u\cdot\nabla\phi\big{)}\,dxdt&-\int_{0}^{T}\int_{\Omega}\frac{\kappa\nabla\theta}{\theta}\cdot\nabla\phi\,dxdt\\\
&+\int_{0}^{T}\int_{\Omega}\sigma\phi\,dxdt=-\int_{\Omega}s_{0}\phi(0,\cdot)\,dx,\end{split}$
(3.16) where $\phi\in C_{c}^{1}([0,T)\times\overline{\Omega})$ and
$\displaystyle\sigma\geq\frac{1}{\theta}\bigg{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\frac{\kappa|\nabla\theta|^{2}}{\theta}\bigg{)}.$ (3.17)
Minimal regularity of solutions required: $\displaystyle\theta>0~{}\text{ a.a
on }(0,T)\times\Omega,~{}~{}\theta\in
L^{q}((0,T)\times\Omega),~{}~{}\nabla\theta\in L^{q}((0,T)\times\Omega)$
$\displaystyle,~{}~{}q>1$ (3.19) $\displaystyle s\in
L^{1}((0,T)\times\Omega),~{}~{}\frac{\nabla\theta}{\theta}\in
L^{1}((0,T)\times\Omega)$ (3.20) $\displaystyle\frac{|\nabla
u|^{2}}{\theta}\in
L^{1}((0,T)\times\Omega),~{}~{}\frac{\rho|u|^{2}}{\theta}\in
L^{1}((0,T)\times\Omega),~{}~{}\frac{|\nabla\theta|^{2}}{\theta^{2}}\in
L^{1}((0,T)\times\Omega)$ . (3.21) Formal interpretation:
$\displaystyle\partial_{t}s+\operatorname{div}(su)=\sigma+\operatorname{div}\bigg{(}\frac{\kappa\nabla\theta}{\theta}\bigg{)}~{}~{}\text{in
}(0,T)\times\Omega$ , (3.22) $\displaystyle
s(0,\cdot)=s_{0},~{}~{}\nabla\theta\cdot
n|_{\partial\Omega}=0,~{}~{}~{}~{}~{}~{}~{}~{}$ (3.23) where
$\displaystyle-s=k\rho(\log\rho-\frac{3}{2}(\log\theta+1)).~{}~{}~{}~{}~{}~{}~{}$
(3.24)

###  3.2 A Priori Estimates

The first observation we make is that the conservation of mass holds, i.e.

$\displaystyle\int_{\Omega}\rho(t,\cdot)\,dx=\int_{\Omega}\rho_{0}\,dx=M_{0}~{}~{}\text{for
a.a. }t\in(0,T).$ (3.25)

The second estimate that follows from the balance of energy is the energy
estimate

$\displaystyle ess\sup_{t\in(0,T)}\int_{\Omega}\rho\theta\,dx\leq
C(\rho_{0},\theta_{0}).$ (3.26)

From the second law of thermodynamics we observe that the total entropy of the
system $S$ defined as $S(t):=\int_{\Omega}s(t,x)dx$ is non-decreasing, i.e.

$\displaystyle\int_{\Omega}s(t,\cdot)\,dx\geq\int_{\Omega}s_{0}\,dx~{}~{}\text{for
a.a }t\in(0,T).$ (3.27)

In the next step we show the positivity of the absolute temperature.

First, we state a Lemma that for for a finite initial entropy the temperature
is positive on a set with positive Lebesgue measure. To this end, we define
two regions in the $(\rho,\theta)$-plane:

  * •

non-degenerate region: low density and/or sufficiently large temperature

$\displaystyle\rho\leq\overline{Z}\theta,~{}~{}\text{for some
}\overline{Z}>0;$

  * •

degenerate region: high density and/or very low temperature

$\displaystyle\rho>\overline{Z}\theta,~{}~{}\text{for some }\overline{Z}>0.$

Next, we set

$\displaystyle s_{\infty}=\lim_{\theta\to
0}s(\rho,\theta)\geq-\infty~{}~{}\text{for any fixed }\rho.$

######  Lemma 3.1.

Let $\Omega\subset\mathbb{R}^{n}$ with $n=2,\,3$ be a bounded Lipschitz
domain. Assume that the non-negative functions $\rho,\,\theta\in
L^{1}((0,T)\times\Omega)$ satisfy

$\displaystyle\int_{\Omega}\rho
dx=M_{0},~{}~{}\int_{\Omega}sdx>M_{0}s_{\infty}+\delta~{}~{}\text{for some
}\delta>0.$

Then there are $\underline{\theta}>0$ and $V_{0}>0$ such that

$\displaystyle\big{|}\\{x\in\Omega|~{}\theta(\cdot,x)>\underline{\theta}\\}\big{|}\geq
V_{0}.$

The Lemma and the idea of the proof are a slight modification of a result by
Feireisl and Novotný in [FN09, Lemma 2.2].

###### Proof.

Assume there exist sequences $\rho_{n},\theta_{n}$ satisfying the assumptions
and such that

$\displaystyle\rho_{n}\to\rho~{}~{}\text{in
}L^{1}(\Omega),~{}~{}\int_{\Omega}\rho
dx=M_{0},~{}~{}\big{|}\\{x\in\Omega|~{}\theta_{n}(x)>\frac{1}{n}\\}\big{|}<\frac{1}{n}.$

Then $\theta_{n}\to 0$ in $L^{1}(\Omega)$. Moreover, by the definition of the
entropy we have

$\displaystyle|s(\rho_{n},\theta_{n})|\leq\rho_{n}(1+|\log\rho_{n}|+|\log\theta_{n}|).$

In the non-degenerate region we obtain

$\displaystyle\int_{\\{\rho_{n}\leq\overline{Z}\theta\\}}s(\rho_{n},\theta_{n})dx$
$\displaystyle\leq
c\int_{\\{\rho_{n}\leq\overline{Z}\theta\\}}\rho_{n}(1+|\log\rho_{n}|+|\log\theta_{n}|)dx$
$\displaystyle\leq
c(\overline{Z})\int_{\Omega}\theta_{n}(1+|\log\rho_{n}|+|\log\theta_{n}|)dx\to
0$

Thus,

$\displaystyle\limsup_{n\to\infty}\int_{\\{\rho_{n}\leq\overline{Z}\theta\\}}s(\rho_{n},\theta_{n})dx\leq
0.$

Moreover, in the degenerate region we have

$\displaystyle\int_{\\{\rho_{n}>\overline{Z}\theta_{n}\\}}s(\rho_{n},\theta_{n})dx=\int_{\\{Z\theta_{n}\geq\rho_{n}>\overline{Z}\theta_{n}\\}}s(\rho_{n},\theta_{n})dx+\int_{\\{\rho_{n}>Z\theta_{n}\\}}s(\rho_{n},\theta_{n})dx,$

where

$\displaystyle\int_{\\{Z\theta_{n}\geq\rho_{n}>\overline{Z}\theta_{n}\\}}s(\rho_{n},\theta_{n})dx\leq
S(\overline{Z})Z\int_{\Omega}\theta_{n}dx\to 0.$

Thus we conclude that

$\displaystyle\liminf_{n\to\infty}\int_{\\{\rho_{n}>Z\theta_{n}\\}}s(\rho_{n},\theta_{n})dx>M_{0}s_{\infty}~{}~{}\text{for
any }Z>\overline{Z}.$

However, this leads to a contradiction as

$\displaystyle\int_{\\{\rho_{n}>Z\theta_{n}\\}}s(\rho_{n},\theta_{n})dx\leq
S(Z)\int_{\\{\rho_{n}>Z\theta_{n}\\}}\rho_{n}dx\to S(Z)M_{0},$

where we used the notation

$\displaystyle s(\rho,\theta)=\rho
S(Z),~{}~{}Z=\frac{\rho}{\theta},~{}~{}\lim_{Z\to\infty}S(Z)=s_{\infty}.$

∎

Thus we have shown that $\int_{V}\log\theta\,dx$ is finite. Next, a version of
Poincare’s inequality provides the following.

######  Corollary 3.2 ([FN09, Proposition 2.2]).

Let $\Omega\subset\mathbb{R}^{n}$ with $n=2,\,3$ be a bounded Lipschitz
domain. Let $V\subset\Omega$ be a measurable set such that $|v|\geq V_{0}>0$.
Then there exists a constant $c(V_{0})$ such that

$\displaystyle\|v\|_{H^{1}(\Omega)}\leq c(V_{0})\bigg{(}\|\nabla
v\|_{L^{2}(\Omega)}+\int_{V}|v|dx\bigg{)}.$

Applying this result to $v=\log\theta$ we obtain the following estimate for
the temperature.

$\displaystyle\int_{0}^{T}\int_{\Omega}|\log\theta|^{2}+|\nabla\log\theta|^{2}dxdt\leq
c(data).$ (3.28)

This implies the positivity of the absolute temperature with a possible
exception on a set of measure zero.

##  4 Existence Theory

The ultimate goal of the forthcoming chapter is to show existence of weak
solutions to the system of equations (3.2)-(3.24).

######  Theorem 4.1 (Local-in-time Existence).

Let $\Omega\subset\mathbb{R}^{n}$ with $n=2,\,3$ be a bounded domain of class
$C^{2,\nu}$, $\nu\in(0,1)$. Assume that the data satisfies the initial
conditions

$\displaystyle\rho_{0}\in L^{p}(\Omega)~{}~{}\text{for some }p\geq
3,~{}~{}\int_{\Omega}\rho_{0}dx=M_{0}>0,$ $\displaystyle
E_{0}=\int_{\Omega}\rho_{0}\theta_{0}dx<\infty,~{}~{}s(\rho_{0},\theta_{0})\in
L^{1}(\Omega).$

In addition, let the initial density and temperature be positive, i.e.
$\rho_{0}(x)>0$ and $\theta_{0}(x)>0$ for all $x\in\Omega$. Then there exists
a time $T>0$ such that the system of equations for the thermal effects of an
ideal gas in the Brinkman model admits a weak solution $(\rho,u,\theta)$ on
$(0,T)\times\Omega)$ in the sense specified before, i.e. $(\rho,u,\theta)$
satisfy relations (3.2)-(3.24).

The proof of the main result consists of several steps outlined as follows:

  * •

The continuity equation is regularized with an artificial viscosity term and
the entropy production equation is replaced by the balance of internal energy.

  * •

Approximate solutions are obtained by a fixed point method.

  * •

Performing the first limit we convert the balance of internal energy to an
approximate entropy production equation containing an additional small
parameter.

  * •

We pass to the limit in the regularized continuity equation and finally pass
to the limit in the entropy production equation.

######  Remark 4.2.

When the energy law also contains a kinetic part the Brinkman-type equation
with additional inertial term becomes the compressible Navier-Stokes-Fourier
system for which the ideal gas case is still open [FN09].

######  Remark 4.3.

In the following subsections we set the parameter $k$ to equal one. This
allows to simplify the computations without changing the main theorem.

###  4.1 Approximate scheme

The first step in the proof of Theorem 4.1 is to regularize the equations
(3.2)-(3.24).

The continuity equation is regularized by means of an artificial viscosity
term

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho
u)=\epsilon\Delta\rho~{}~{}\text{in }(0,T)\times\Omega,$ (4.1)

with homogeneous Neumann boundary condition

$\displaystyle\nabla\rho\cdot n=0~{}~{}\text{on }(0,T)\times\partial\Omega,$
(4.2)

and the initial condition

$\displaystyle\rho(0,\cdot)=\rho_{0,\delta},$ (4.3)

where

$\displaystyle\rho_{0,\delta}\in
C^{2,\nu}(\overline{\Omega}),~{}~{}\inf_{x\in\Omega}\rho_{0,\delta}>0,~{}~{}\nabla\rho_{0,\delta}\cdot
n|_{\partial\Omega}=0.$ (4.4)

The Brinkman-type equation

$\displaystyle\int_{\Omega}\rho\theta\operatorname{div}\phi
dx=\int_{\Omega}\mu\nabla u:\nabla\phi dx+\int_{\Omega}\nu\rho u\phi dx,$
(4.5)

for any test function $\phi\in X_{n}$, where

$\displaystyle X_{n}\subset C^{2,\nu}(\Omega)\subset L^{2}(\Omega)$ (4.6)

is a finite dimensional vector space endowed with the Hilbert structure of the
Lebesgue space $L^{2}$ and the functions satisfy

$\displaystyle\phi=0~{}~{}\text{on }\partial\Omega~{}~{}\text{ no-slip
boundary conditions.}$ (4.7)

Instead of the entropy production, we consider a modified internal energy
equation of the form

$\displaystyle\begin{split}\partial_{t}e_{\delta}+\operatorname{div}(e_{\delta}u)-\kappa\Delta\theta&=\mu|\nabla
u|^{2}+\nu\rho|u|^{2}-\rho\theta\operatorname{div}u+\delta\frac{1}{\theta^{2}}\\\
&+\epsilon\delta(\rho^{\Gamma}+2)|\nabla\rho|^{2}-\delta\theta^{5}~{}~{}\text{in
}(0,T)\times\Omega,\end{split}$ (4.8)

with the Neumann boundary condition

$\displaystyle\nabla\theta\cdot n=0~{}~{}\text{on }(0,T)\times\partial\Omega,$
(4.9)

and the initial condition

$\displaystyle\theta(0,\cdot)=\theta_{0,\delta},$ (4.10)

where

$\displaystyle\theta_{0,\delta}\in H^{1}(\Omega)\cap
L^{\infty}(\Omega),~{}~{}ess\inf_{x\in\Omega}\theta_{0,\delta}>0.$ (4.11)

Here

$\displaystyle e_{\delta}(\rho,\theta)=\rho\theta.$ (4.12)

Moreover, the approximate internal energy balance holds

$\displaystyle\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0,\delta}\theta_{0,\delta}dx+\int_{0}^{t}\int_{\Omega}\frac{\delta}{\theta^{2}}-\delta\theta^{5}\,dxd\tau.$
(4.13)

for all $t\in[0,T]$.
The quantities $\epsilon$ and $\delta$ are small positive parameters, yielding
better estimates of the approximate scheme.

###  4.2 Solvability of the approximate system

The second step is to show existence of classical solutions to the approximate
system.

######  Theorem 4.4 (Global existence for the approximate system).

Let $\epsilon$, $\delta$ be given positive parameters. Under the hypotheses of
Theorem 4.1 there exists a $\Gamma_{0}>0$ such that for all
$\Gamma>\Gamma_{0}$ the approximate system (4.1)-(4.13) admits a strong
solution $(\rho,u,\theta)$ belonging to the following regularity class

$\displaystyle\rho\in
C([0,T];C^{2,\nu}(\overline{\Omega})),~{}~{}\partial_{t}\rho\in
C([0,T];C^{0,\nu}(\overline{\Omega})),~{}~{}\inf_{[0,T]\times\overline{\Omega}}\rho>0,$
$\displaystyle u\in C^{1}(X_{n}),$ $\displaystyle\theta\in
C([0,T];H^{2}(\Omega))\cap
L^{\infty}((0,T)\times\Omega),~{}~{}\partial_{t}\theta\in
L^{2}((0,T)\times\Omega),~{}~{}\operatorname*{ess\,inf}_{(0,T)\times\Omega}\theta>0.$

The idea of the proof follows standard arguments:

  * •

Given the velocity field $u$, the approximate continuity equation is solved
directly by standard parabolic theory;

  * •

After solving the continuity equation we determine the temperature $\theta$ of
the quasilinear parabolic problem, i.e. the internal energy equation, where
$u,~{}\rho$ play the role of given data;

  * •

To close the loop, the solution $u$ is looked for as the fixed point of an
integro-differential operator.

######  Lemma 4.5 ([FN09, Lemma 3.1]).

Let $\Omega\subset\mathbb{R}^{n}$ where $n=2,\,3$ be a bounded domain of class
$C^{2,\nu},~{}\nu\in(0,1)$ and let $u\in X_{n}$ be a given vector field.
Suppose the initial data has the regularity specified in Section 4.1.
Then the continuity equation posses a unique classical solution
$\rho=\rho_{u}$, more specifically

$\displaystyle\rho_{u}\in V\equiv\big{\\{}\rho\in
C([0,T];C^{2,\nu}(\overline{\Omega})),~{}\partial_{t}\rho\in
C([0,T];C^{0,\nu}(\overline{\Omega}))\big{\\}}.$ (4.14)

Moreover, the mapping $u\in X_{n}\to\rho_{u}$ maps bounded sets in $X_{n}$
into bounded sets in $V$ and is continuous with values in
$C^{1}([0,T]\times\overline{\Omega})$. Finally,

$\displaystyle\underline{\rho_{0}}\exp\bigg{(}-\int_{0}^{\tau}\|\operatorname{div}u\|_{L^{\infty}(\Omega)}dt\bigg{)}\leq\rho_{u}(\tau,x)\leq\overline{\rho_{0}}\exp\bigg{(}+\int_{0}^{\tau}\|\operatorname{div}u\|_{L^{\infty}(\Omega)}dt\bigg{)},$
(4.15)

for all $\tau\in[0,T],~{}x\in\Omega$, where
$\underline{\rho_{0}}=\inf_{x\in\Omega}\rho_{0,\delta},~{}\overline{\rho_{0}}=\sup_{x\in\Omega}\rho_{0,\delta}$.

In this part we focus on the quasilinear parabolic problem for the unknown
temperature $\theta$. First, we state a comparison principle.

######  Lemma 4.6 ([FN09, Lemma 3.2]).

Given the quantities $u,\,\rho_{u}$ satisfying the regularity $u\in X_{n}$,
$\rho\in C([0,T];C^{2}(\overline{\Omega}))$, $\partial_{t}\rho\in
C([0,T]\times\overline{\Omega})$, where $\inf_{(0,T)\times\Omega}\rho>0$, and
assume that $\underline{\theta}$ and $\overline{\theta}$ are a sub- and super-
solution to problem belonging to the regularity class

$\displaystyle\underline{\theta},\,\overline{\theta}\in
L^{2}(0,T;H^{2}(\Omega)),~{}\partial_{t}\underline{\theta},\,\partial_{t}\overline{\theta}\in
L^{2}((0,T)\times\Omega),$ (4.16) $\displaystyle
0<ess\inf_{(0,T)\times\Omega}\underline{\theta}\leq
ess\sup_{(0,T)\times\Omega}\underline{\theta}<\infty,~{}~{}0<ess\inf_{(0,T)\times\Omega}\overline{\theta}\leq
ess\sup_{(0,T)\times\Omega}\overline{\theta}<\infty,$ (4.17)

and satisfying

$\displaystyle\underline{\theta}(0,\cdot)\leq\overline{\theta}(0,\cdot)~{}~{}\text{a.e.
in }\Omega.$ (4.18)

Then

$\displaystyle\underline{\theta}(t,x)\leq\overline{\theta}(t,x)~{}~{}\text{a.e.
in }(0,T)\times\Omega.$

######  Remark 4.7.

If we assume in addition that

$\displaystyle 0<\underline{\theta}_{0}=ess\inf_{\Omega}\theta_{0,\delta}\leq
ess\sup_{\Omega}\theta_{0,\delta}=\overline{\theta}_{0}<\infty,$ (4.19)

the problem (4.8)-(4.12) admits at most one solution in the regularity class
specified above.

######  Corollary 4.8.

Let $u,\,\rho_{u}$ be in the regularity class as before and let the initial
data $\theta_{0,\delta}$ satisfy (4.19). Suppose that $\theta$ is a strong
solution of the problem belonging to the regularity class.

Then there exists two constants $\underline{\theta}$, $\overline{\theta}$
satisfying
$0<\underline{\theta}<\underline{\theta}_{0}<\overline{\theta_{0}}<\overline{\theta}$
and

$\displaystyle\underline{\theta}\leq\theta(t,x)\leq\overline{\theta}~{}~{}\text{for
a.a. }(t,x)\in(0,T)\times\Omega.$ (4.20)

###### Proof.

We check that the constant function $\underline{\theta}$ is a subsolution of
the problem

$\displaystyle\underline{\theta}\partial_{t}\rho+\underline{\theta}\operatorname{div}(\rho
u)=\mu|\nabla
u|^{2}+\nu\rho|u|^{2}-\underline{\theta}\rho\operatorname{div}u+\delta\frac{1}{\theta^{2}}+\epsilon\delta(\rho^{\Gamma}+2)|\nabla\rho|^{2}-\delta\theta^{5}.$

This is true as long as

$\displaystyle\frac{\delta}{\theta^{2}}\leq\underline{\theta}\rho\operatorname{div}u+\underline{\theta}(\partial_{t}\rho+\operatorname{div}(\rho
u))-\mu|\nabla
u|^{2}-\nu\rho|u|^{2}-\epsilon\delta(\rho^{\Gamma}+2)|\nabla\rho|^{2}.$

We note that all terms on the right-hand-side are bounded in terms of
$\|u\|_{X_{n}}$ and $\|\rho\|_{C^{1}}$ provided $0<\underline{\theta}<1$. Then
by the comparison principle (Lemma 4.4) the first inequality follows.
The upper bound can be established in a similar way by help of the dominating
term $-\delta\theta^{5}$ in equation (4.8). ∎

In addition, we observe the importance of the term
$\frac{\delta}{\theta^{2}}$. This term guarantees that the absolute
temperature stays bounded away from zero.

######  Lemma 4.9 ([FN09, Lemma 3.3]).

Let the data $\rho_{u},\,u$ belong to the regularity class as specified above
and let the initial data $\theta_{0,\delta}\in H^{1}(\Omega)$.

Then any strong solution $\theta$ of the problem belonging to the regularity
class (4.16) satisfies the estimate

$\displaystyle\begin{split}ess\sup_{t\in(0,T)}\|\theta\|_{H^{1}(\Omega)}^{2}&+\int_{0}^{T}\big{(}\|\partial_{t}\theta\|^{2}_{L^{2}(\Omega)}+\|\Delta\theta\|^{2}_{L^{2}(\Omega)}\big{)}dt\\\
&\leq
C\big{(}\Omega,\|\rho\|_{C^{1}([0,t]\times\overline{\Omega})},\|u\|_{X_{n}},\inf_{(0,t)\times\Omega}\rho,\|\theta_{0,\delta}\|_{H^{1}(\Omega)}\big{)}.\end{split}$
(4.21)

After establishing a priori bounds on the temperature $\theta$ we are able to
show existence of strong solutions to the approximate internal energy
equation. The key to this is that those bounds lead to the compactness of the
solutions in $L^{2}(0,T);H^{1}(\Omega))$. Note that we can rewrite the
approximate internal energy equation as a quasilinear parabolic equation in
the temperature $\theta$. For smooth enough data we can apply the results by
Ladyzhenskaya [LSU68] to obtain a unique strong solution.

######  Lemma 4.10.

Under the previous assumptions the problem (4.8)-(4.12) has a unique strong
solution $\theta=\theta_{u}$ belonging to the regularity class

$\displaystyle\begin{split}Y=\big{\\{}&\partial_{t}\theta\in
L^{2}((0,T)\times\Omega),~{}~{}\theta\in L^{\infty}(0,T;H^{2}(\Omega)\cap
L^{\infty}(\Omega)),\\\ &\frac{1}{\theta}\in
L^{\infty}((0,T)\times\Omega)\big{\\}}.\end{split}$ (4.22)

Moreover, the mapping $u\to\theta_{u}$ maps bounded sets in $X_{n}$ into
bounded sets in $Y$ and is continuous with values in
$L^{2}(0,T;H^{1}(\Omega))$.

Now, we are able to show existence of solutions to the approximate system. We
recall that $u\to(\rho_{u},\theta_{u})$ maps bounded sets in $X_{n}$ into
bounded sets in $C([0,T],C^{2,\nu}(\overline{\Omega}))\times
L^{\infty}(0,T;H^{2}(\Omega))$ and is continuous with values in
$C^{1}([0,T]\times\Omega)\times L^{2}(0,T;H^{1}(\Omega))$. Moreover from the
Brinkman-type equation we obtain using the positivity of $\rho$ that

$\displaystyle\|\rho u^{2}\|_{L^{1}(\Omega)}+\|\nabla
u\|^{2}_{L^{2}(\Omega)}\leq C\|\rho\theta\|^{2}_{L^{2}(\Omega)}\leq
C\big{(}\|\rho\|_{L^{4}(\Omega)}^{2}+\|\theta\|_{L^{4}(\Omega)}^{2}\big{)}.$

Applying the Poincare inequality and the Sobolev imbedding yields

$\displaystyle\|u\|_{H^{1}(\Omega)}\leq
C\big{(}\|\rho\|_{H^{1}(\Omega)}^{2}+\|\theta\|_{H^{1}(\Omega)}^{2}\big{)}.$
(4.23)

Moreover, we have

$\displaystyle\|\Delta
u\|_{L^{2}(\Omega)}\leq\|\nabla(\rho\theta)\|_{L^{2}(\Omega)}\leq\|\rho\|_{L^{\infty}(\Omega)}\|\nabla\theta\|_{L^{2}(\Omega)}+\|\theta\|_{L^{2}(\Omega)}\|\nabla\rho\|_{L^{\infty}(\Omega)}.$
(4.24)

Thus, for each $t\in(0,T)$ $(\rho,\theta)\to u$ maps bounded sets in
$C^{2,\nu}(\overline{\Omega})\times H^{2}(\Omega)$ into bounded sets in
$X_{n}$ and is continuous with values in $H^{1}(\Omega)$.

Now, a direct application of the Leray-Schauder fixed point theorem yields the
existence of a solution $(\rho,u,\theta)$ of the approximate system on a
(possibly short) time interval $(0,T_{n})$. Iterating this procedure yields
the existence of solutions on $(0,T)$ as long as the bounds are independent of
the time $T_{n}$.

###  4.3 Uniform estimates

In this section we establish uniform bounds, similar to those in chapter 3.
The existence of such uniform bounds guarantees the global existence of the
approximate solution in the desired spaces. Moreover, these estimates play a
crucial role in the limit passage in the following sections.

First, from the approximate continuity equation it follows that the total mass
of the system still remains constant in time, i.e.

$\displaystyle\int_{\Omega}\rho(t)dx=\int_{\Omega}\rho_{0,\delta}dx=M_{0,\delta}~{}~{}\text{for
all }t\in[0,T].$ (4.25)

Next, taking $u$ as a test function in the Brinkman-type equation (4.5) we
obtain

$\displaystyle\int_{\Omega}\rho\theta\operatorname{div}u-|\nabla
u|^{2}-\rho|u|^{2}dx=0.$ (4.26)

From the approximate internal energy equation (4.8) we have

$\displaystyle\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0,\delta}\theta_{0,\delta}dx+\int_{0}^{T}\int_{\Omega}\frac{\delta}{\theta^{2}}-\epsilon\delta\theta^{5}dxdt.$
(4.27)

Instead of working with the internal energy balance we manipulate the equation
(4.8) to get an approximate entropy production. To this end, we recall from
the classical theory of thermodynamics that

$\displaystyle\frac{de}{dt}=W+Q~{}~{}\text{and }T\frac{ds}{dt}=T\Delta+Q.$

Therefore, we compute

$\displaystyle~{}~{}~{}-\partial_{t}\big{(}\rho(\log\rho-\log\theta-1)\big{)}-\operatorname{div}\big{(}\rho(\log\rho-\log\theta-1)u\big{)}$
$\displaystyle=-\rho_{t}(\log\rho-\log\theta-1)-\rho_{t}+\frac{\rho\theta_{t}}{\theta}-(\log\rho-\log\theta-1)\operatorname{div}(\rho
u)$ $\displaystyle~{}~{}~{}\,-\rho u\cdot\nabla(\log\rho-\log\theta-1)$
$\displaystyle=(\log\theta-\log\rho)(\rho_{t}+\operatorname{div}(\rho
u))+\frac{\rho\theta_{t}}{\theta}+\operatorname{div}(\rho
u)-u\cdot\nabla\rho+\frac{\rho}{\theta}u\cdot\nabla\theta$
$\displaystyle=(\log\theta-\log\rho-1)\epsilon\Delta\rho+\rho\operatorname{div}u+\frac{\rho}{\theta}(\theta_{t}+u\cdot\nabla\theta)+\epsilon\Delta\rho$
$\displaystyle=(\log\theta-\log\rho-1)\epsilon\Delta\rho+\rho\operatorname{div}u+\frac{1}{\theta}\big{(}(\rho\theta)_{t}+\operatorname{div}(\rho\theta
u)\big{)}$

which gives the relation between the internal energy and the entropy of the
system.

Then, dividing the internal energy balance by $\theta$ yields

$\displaystyle\partial_{t}s+\operatorname{div}(su)-\frac{\Delta\theta}{\theta}$
$\displaystyle=\frac{1}{\theta}(\mu|\nabla
u|^{2}+\nu\rho|u|^{2})+\frac{\delta}{\theta^{3}}-\delta\theta^{4}$
$\displaystyle+\epsilon\Delta\rho(\log\theta-\log\rho-1).$ Rewriting this
expression yields the approximate entropy equation
$\displaystyle\partial_{t}s+\operatorname{div}(su)-\operatorname{div}(\frac{\nabla\theta}{\theta})$
$\displaystyle=\frac{1}{\theta}(\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\frac{|\nabla\theta|^{2}}{\theta})+\frac{\delta}{\theta^{3}}-\delta\theta^{4}$
$\displaystyle+\epsilon\Delta\rho(\log\theta-\log\rho-1).$

For higher regularity we modify the entropy production rate slightly and
obtain

$\displaystyle\begin{split}\partial_{t}s+\operatorname{div}(su)-\operatorname{div}(\frac{\nabla\theta}{\theta})&=\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}\\\
&-\delta\theta^{4}+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\Delta\rho(\log\theta-\log\rho-1).\end{split}$
(4.28)

The two approximate equations give rise to the following estimates, where we
multiplied equation (4.28) by an arbitrary positive constant
$\overline{\theta}$ and integrated over the space time domain.

$\displaystyle\int_{\Omega}H_{\delta,\overline{\theta}}(\rho,\theta)(\tau)dx+\overline{\theta}\int_{0}^{\tau}\int_{\Omega}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}dxdt$
$\displaystyle+\int_{0}^{\tau}\int_{\Omega}\frac{\delta}{\theta^{3}}+\delta\theta^{5}+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}dxdt$
$\displaystyle=\int_{\Omega}H_{\delta,\overline{\theta}}(\rho,\theta)(0)dx+\int_{0}^{\tau}\int_{\Omega}\frac{\delta}{\theta^{2}}+\delta\theta^{4}-\epsilon\overline{\theta}\Delta\rho(\log\theta-\log\rho-1)dxdt,$

where
$H_{\delta,\overline{\theta}}(\rho,\theta)=e-s=\rho\theta+\rho(\log\rho-\log\theta-1)$.
Integration by pats in the last term yields

$\displaystyle\begin{split}&\int_{\Omega}H_{\delta,\overline{\theta}}(\tau)dx+\overline{\theta}\int_{0}^{\tau}\int_{\Omega}\frac{1}{\theta}(\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}dxdt\\\
&+\int_{0}^{\tau}\int_{\Omega}\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho}dxdt\\\
&=\int_{\Omega}H_{\delta,\overline{\theta}}(0)dx+\int_{0}^{\tau}\int_{\Omega}\frac{\delta}{\theta^{2}}+\delta\theta^{4}+\epsilon\overline{\theta}\frac{\nabla\rho\cdot\nabla\theta}{\theta}dxdt~{}~{}\text{for
all }\tau\in[0,T].\end{split}$ (4.29)

We observe that we can control the terms on the right-hand side by the terms
on the left-hand side and the initial data.

Indeed, the quantity $\frac{\delta}{\theta^{2}}$ is dominated for low
temperatures by its counterpart $\frac{\delta}{\theta^{3}}$. Moreover, we
estimate

$\displaystyle\int_{0}^{\tau}\int_{\Omega}\epsilon\overline{\theta}\frac{\nabla\rho\cdot\nabla\theta}{\theta}dxdt$
$\displaystyle\leq\int_{0}^{\tau}\int_{\Omega}\bigg{|}\epsilon\overline{\theta}\frac{\nabla\rho\cdot\nabla\theta}{\theta}\bigg{|}dxdt$
$\displaystyle\leq\int_{0}^{\tau}\int_{\Omega}\epsilon\overline{\theta}\frac{|\nabla\theta|^{2}}{2\theta^{2}}dxdt+\int_{0}^{\tau}\int_{\Omega}\epsilon\overline{\theta}|\nabla\rho|^{2}dxdt,$

where we can bound both terms with terms on the left-hand side of equation
(4.29) for sufficiently small $\epsilon$.

Now, we can summarize the estimates as follows

$\displaystyle
ess\sup_{t\in(0,T)}\int_{\Omega}H_{\delta,\overline{\theta}}(\rho,\theta)(t)dx\leq
c$ (4.30)
$\displaystyle\int_{0}^{T}\int_{\Omega}\frac{1}{\theta}\bigg{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\bigg{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}dxdt\leq
c$ (4.31)
$\displaystyle\int_{0}^{T}\int_{\Omega}\epsilon\frac{|\nabla\rho|^{2}}{\rho}+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}dxdt\leq
c,$ (4.32)

where c is a positive constant depending on the initial data but is
independent of $T,\,\epsilon,\,\delta$ and $n$.

###  4.4 Limit Passage

The third step is to pass to the limit in the approximate/ regularized system.
This is done in three steps. First, we let $n\to\infty$. This is the
approximation step in the Fourier series for the Brinkman-type equation. Next,
we pass to the limit in $\epsilon$ as $\epsilon\to 0$, i.e. the additional
regularity term for the density $\rho$. And finally, we let $\delta\to 0$.
This is the most crucial step because only here we will have a smallness
condition on the time interval in which the weak solution exists.

Each step utilizes similar ideas: from finding uniform bounds to the div-curl
lemma. For the details of the theorems used in these steps we refer to the
appendix.

####  4.4.1 Limit $n\to\infty$

Let the vector space $X$ be defined as

$\displaystyle X\equiv\bigcup_{n}X_{n}~{}~{}\text{is dense in
}H^{1}_{0}(\Omega).$

For $\epsilon>0$ and $\delta>0$ fixed let $(\rho_{n},u_{n},\theta_{n})_{n}$ be
a sequence of approximate solutions.

In addition to the uniform estimates (4.30)-(4.32) we obtain the following
estimates:
From the quasilinear parabolic equation for the absolute temperature $\theta$
we obtain that
$0\leq\underline{\theta}\leq\theta_{n}(t,x)\leq\overline{\theta}$ for all
$(t,x)\in(0,T)\times\Omega$. Hence, $\|\nabla u_{n}\|^{2}_{L^{2}(\Omega)}\leq
c$ and together with the boundary condition $u|_{\partial\Omega}=0$ we have
$u_{n}\in H^{1}_{0}(\Omega)$. Then, by the Sobolev embedding $u_{n}\in
L^{6}(\Omega)$. The boundedness of the sequence then implies that

$\displaystyle u_{n}\to u~{}~{}\text{weakly in }H^{1}(\Omega).$ (4.33)

From the kinetic equation it follows that

$\displaystyle\frac{d}{dt}\frac{1}{2}\|\rho_{n}\|_{L^{2}(\Omega)}^{2}+\epsilon\|\nabla\rho_{n}\|^{2}_{L^{2}(\Omega)}\leq
C\|u_{n}\|_{L^{6}(\Omega)}^{2}\|\rho_{n}\|^{2}_{L^{4}(\Omega)}.$

By the interpolation inequality we have

$\displaystyle\frac{d}{dt}\|\rho_{n}\|_{L^{2}(\Omega)}^{2}+\epsilon\|\nabla\rho_{n}\|^{2}_{L^{2}(\Omega)}\leq
C(M_{0,\delta},\epsilon,\|u\|_{L^{6}(\Omega)},\Omega).$

Thus the sequence $\rho_{n}$ is bounded in $L^{\infty}(0,T;L^{2}(\Omega))\cap
L^{2}(0,T;H^{1}(\Omega))$. Therefore we can assume that

$\displaystyle\rho_{n}\to\rho~{}~{}\text{weakly-(*) in
}L^{\infty}(0,T;L^{2}(\Omega)).$ (4.34)

By applying Poincare inequality to the additional term in equation (4.32) we
get

$\displaystyle\rho_{n},~{}\rho_{n}^{\Gamma/2}~{}~{}\text{bounded in
}L^{2}(0,T;H^{1}(\Omega)),$ (4.35)

and by interpolation

$\displaystyle\rho_{n}~{}~{}\text{bounded in
}L^{\Gamma}(0,T;L^{3\Gamma}(\Omega)).$ (4.36)

By the boundedness of the entropy production rate we obtain that

$\displaystyle\nabla\theta_{n}^{\Gamma/2}~{}~{}\text{bounded in
}L^{2}(0,T;L^{2}(\Omega))$ (4.37) and
$\displaystyle\theta_{n}~{}~{}\text{bounded in }L^{5}((0,T)\times\Omega).$
(4.38)

In addition, $\nabla\log\theta_{n}$ is bounded in $L^{2}(0,T;L^{2}(\Omega))$.
This implies that $\theta_{n},~{}\theta_{n}^{\Gamma/2}\in
L^{2}(0,T;H^{1}(\Omega))$ and we obtain

$\displaystyle\theta_{n}\to\theta~{}~{}\text{weakly in
}L^{2}(0,T;H^{1}(\Omega)).$ (4.39)

Moreover, we have

$\displaystyle\int_{0}^{T}\int_{\Omega}\frac{1}{\theta^{3}}dxdt\leq\liminf_{n\to\infty}\int_{0}^{T}\int_{\Omega}\frac{1}{\theta_{n}^{3}}dxdt.$
(4.40)

By the standard Sobolev embedding we derive the higher integrability estimates
of $\theta_{n}$, i.e.

$\displaystyle\nabla\theta_{n}~{}~{}\text{bounded in
}L^{\Gamma}(0,T;L^{3\Gamma}(\Omega)).$ (4.41)

As a byproduct, we get that

$\displaystyle\log\theta_{n}~{}~{}\text{bounded in
}L^{q}((0,T)\times\Omega)~{}\text{for any finite }q\geq 1.$ (4.42)

Now, we are able to pass to the limit in the equations with respect to the
parameter $n$.
The limit in the Poisson equation is obtained via the standard Galerkin
argument, where we note that

$\displaystyle\int_{\Omega}\rho_{n}\theta_{n}\operatorname{div}u_{n}dx\leq\frac{1}{2}\|\nabla
u\|^{2}_{L^{2}(\Omega)}+C\big{(}\|\rho_{n}\|^{2}_{L^{4}(\Omega)}+\|\theta_{n}\|^{2}_{L^{4}(\Omega)}\big{)}.$
(4.43)

Thus we can pass to the limit in the Brinkman equation by the bounds
established before

$\displaystyle-\rho u+\Delta u=\nabla(\rho\theta)~{}~{}\text{for a.e.
}x\in\Omega.$ (4.44)

From the kinetic equation we have

$\displaystyle\big{(}\partial_{t}-\epsilon\Delta\big{)}[\rho_{n}]=-\nabla\rho_{n}\cdot
u-\rho_{n}\operatorname{div}u_{n},$

where the terms on the right-hand side are bounded in
$L^{2}(0,T;L^{3/2}(\Omega))$. Thus

$\displaystyle\rho_{n}\to\rho~{}~{}\text{ a.e. in }(0,T)\times\Omega$ (4.45)

and we can let $n\to\infty$ in the approximate continuity equation to obtain

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho
u)=\epsilon\Delta\rho~{}~{}\text{a.e. in }(0,T)\times\Omega,$ (4.46)

where the density $\rho$ is a positive function satisfying

$\displaystyle\nabla\rho(t,\cdot)\cdot n|_{\partial\Omega}=0~{}~{}\text{for
a.a. }t\in(0,T)$ (4.47)

in the sense of traces, together with the initial data

$\displaystyle\rho(0,\cdot)=\rho_{0,\delta}.$ (4.48)

Moreover,

$\displaystyle\nabla\rho_{n}\to\nabla\rho~{}~{}\text{in
}L^{2}(0,T;L^{2}(\Omega)),$ (4.49)

where we used that

$\displaystyle\int_{\Omega}\rho_{n}^{2}(\tau)dx+2\epsilon\int_{0}^{\tau}\int_{\Omega}|\nabla\rho_{n}|^{2}dxdt\to\int_{\Omega}\rho_{\delta,0}^{2}dx-\int_{0}^{\tau}\int_{\Omega}\rho^{2}\operatorname{div}udxdt$

and

$\displaystyle\int_{0}^{\tau}\int_{\Omega}\rho^{2}\operatorname{div}udxdt=\int_{\Omega}\rho^{2}(\tau)dx+2\epsilon\int_{0}^{\tau}\int_{\Omega}|\nabla\rho|^{2}dxdt.$

Now, we can consider the limit in the entropy balance equation. The main
difficulty here is to show the strong convergence of the temperature. To this
end, we apply the Div-Curl Lemma discovered by Tartar [Tar75] to the function
$U$ and $V$, specified below. The details off the div-curl Lemma can be found
in the Appendix.

We rewrite the equation (4.28) in the following form

$\displaystyle\begin{split}\partial_{t}s&+\operatorname{div}(su)-\operatorname{div}(\frac{\nabla\theta}{\theta})-\epsilon\operatorname{div}(\nabla\rho(\log\theta-\log\rho-1))\\\
&=\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\\\
&+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\big{(}\frac{|\nabla\rho|^{2}}{\rho}+\frac{\nabla\rho\cdot\nabla\theta}{\theta}\big{)},\end{split}$
(4.50)

where we used that

$\displaystyle\epsilon\Delta\rho(\log\theta-\log\rho-1)=\epsilon\operatorname{div}(\nabla\rho(\log\theta-\log\rho-1))-\epsilon\big{(}\frac{|\nabla\rho|^{2}}{\rho}+\frac{\nabla\rho\cdot\nabla\theta}{\theta}\big{)}.$

Setting

$\displaystyle
U=[s,su-\frac{\nabla\theta}{\theta}-\epsilon\nabla\rho(\log\theta-\log\rho-1)]~{}~{}\text{and
}~{}V=[\theta,0,0,0]$

we can check the assumptions for the Div-Curl Lemma.

The temperature $\theta$ is bounded in $L^{2}((0,T)\times\Omega)$ and
$curl(V)$ yields only spatial partial derivatives and thus is bounded in
$L^{2}((0,T)\times\Omega)$ which is compact embedded into
$W^{-1,2}((0,T)\times\Omega)$. By the uniform estimates obtained before we
note that the right-hand side of equation (4.50) is bounded in
$L^{1}((0,T)\times\Omega)$ and therefore precompact in
$W^{-1,s}((0,T)\times\Omega)$ provided $s\in[1,\frac{4}{3})$. Thus it remains
to show that $U$ is bounded in a better space than $L^{1}$. To see this we
note that

$\displaystyle|s(\rho_{n},\theta_{n})|\leq
c(\rho_{n}+\rho_{n}|\log\rho_{n}|+\rho_{n}|\log\theta_{n}|)$

and by the uniform estimates $s$ is bounded in
$L^{\Gamma/3}((0,T)\times\Omega)$. In addition, $su$ is bounded in
$L^{p}((0,T)\times\Omega)$, where $\frac{1}{p}=\frac{1}{2}+\frac{3}{\Gamma}$
provided $\Gamma>6$.

For the other terms we have that $\nabla\log\theta_{n}$ is bounded in
$L^{2}((0,T)\times\Omega)$ and

$\displaystyle\epsilon\nabla\rho(\log\theta-\log\rho-1)~{}~{}\text{is bounded
in }L^{\frac{2\Gamma}{\Gamma+6}}((0,T)\times\Omega).$

Then the Div-Curl Lemma states that

$\displaystyle\overline{s(\rho,\theta)\,\theta}=\overline{s(\rho,\theta)}\theta,$
(4.51)

where the symbol $\overline{F(u)}$ denotes the weak $L^{1}$-limit of the
sequence $F(u_{n})$ of composed functions.

The goal is now, to conclude that we have almost everywhere convergence of
$\theta$. By the definition of the entropy we have

$\displaystyle s(\rho,\theta)=-\rho(\log\rho-\log\theta-1),$

where we note that the entropy $s$ is increasing in the temperature. Thus,

$\displaystyle\overline{\rho\log(\theta)\theta}\geq\overline{\rho\log\theta}\theta.$
(4.52)

Moreover, by the strong convergence of $\rho_{n}$ see equation (4.45), we have

$\displaystyle\overline{\rho\log(\theta)\theta}=\rho\overline{(\log\theta)\theta}.$

Combining the above equations we infer that

$\displaystyle\overline{\log(\theta)\theta}=\overline{\log\theta}\theta.$

By the strict convexity of the function $x\log x$ we have that

$\displaystyle\theta_{n}\to\theta~{}~{}\text{a.e in }(0,T)\times\Omega.$
(4.53)

For details in the argument we refer to Appendix.

Now, we can take the limit in the approximate entropy equation. To this end,
we first turn the equation into an inequality by applying Youngs inequality to
the $\nabla\rho\cdot\nabla\theta$ term.

$\displaystyle\begin{split}&\partial_{t}s+\operatorname{div}(su)-\operatorname{div}(\frac{\nabla\theta}{\theta})-\epsilon\operatorname{div}(\nabla\rho(\log\theta-\log\rho-1))\\\
\geq&\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\\\
&+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho},\end{split}$
(4.54)

As a consequence from the previous results we can identify the limits of the
individual terms.

$\displaystyle s(\rho_{n},\theta_{n})$ $\displaystyle\to
s(\rho,\theta)~{}~{}\text{in }L^{2}((0,T)\times\Omega),$ $\displaystyle
s(\rho_{n},\theta_{n})u_{n}$ $\displaystyle\to
s(\rho,\theta)u~{}~{}\text{weakly in }L^{1}((0,T)\times\Omega).$

The almost all convergence of $\theta_{n}$ implies that

$\displaystyle\frac{1}{\theta_{n}}$
$\displaystyle\to\frac{1}{\theta}~{}~{}\text{in }L^{2}((0,T)\times\Omega),$
$\displaystyle\nabla\log\theta_{n}$
$\displaystyle\to\nabla\log\theta~{}~{}\text{weakly in
}L^{1}((0,T)\times\Omega),$
$\displaystyle\theta_{n}^{\Gamma-1}\nabla\theta_{n}$
$\displaystyle\to\theta^{\Gamma-1}\nabla\theta~{}~{}\text{weakly in
}L^{p}((0,T)\times\Omega)~{}p>1.$

To control the $\epsilon$-term we note that

$\displaystyle|\epsilon(\log\theta_{n}-\log\rho_{n}-1)\nabla\rho_{n}|\leq
c|\nabla\rho_{n}|(|\log\theta_{n}|+|\log\rho_{n}|+1)$

and all terms on the right are bounded in $L^{p}((0,T)\times\Omega)$ for some
$p>1$. Thus we have

$\displaystyle\epsilon(\log\theta_{n}-\log\rho_{n}-1)\nabla\rho_{n}\to\epsilon(\log\theta-\log\rho-1)\nabla\rho~{}~{}\text{weakly
in }L^{1}((0,T)\times\Omega).$

Identifying the limit in the remaining terms of the entropy production rate
yields

$\displaystyle\frac{\nabla u_{n}}{\sqrt{\theta_{n}}}$
$\displaystyle\to\frac{\nabla u}{\sqrt{\theta}}~{}~{}\text{weakly in
}L^{2}((0,T)\times\Omega),$
$\displaystyle\frac{\nabla\rho_{n}}{\sqrt{\rho_{n}}}$
$\displaystyle\to\frac{\nabla\rho}{\sqrt{\rho}}~{}~{}\text{weakly in
}L^{2}((0,T)\times\Omega),$
$\displaystyle\sqrt{\Gamma\rho_{n}^{\Gamma-2}+2}\nabla\rho_{n}$
$\displaystyle\to\sqrt{\Gamma\rho^{\Gamma-2}+2}\nabla\rho~{}~{}\text{weakly in
}L^{2}((0,T)\times\Omega).$

These convergence results are sufficient to perform the weak limit in the
approximate entropy equation as $n\to\infty$. We note that the inequality is
preserved under the weak limit due to the lower semi-continuity of convex
superposition operators. This allows us to conclude that

$\displaystyle\begin{split}&\int_{0}^{T}\int_{\Omega}s(\rho,\theta)\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt+\int_{0}^{T}\int_{\Omega}\big{(}\frac{\nabla\theta}{\theta}-\epsilon(\log\theta-\log\rho-1)\nabla\rho\big{)}\cdot\nabla\phi
dxdt\\\ &+\int_{0}^{T}\int_{\Omega}\bigg{(}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\bigg{)}\phi
dxdt\\\
&+\int_{0}^{T}\int_{\Omega}\bigg{(}\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho}\bigg{)}\phi
dxdt\leq-\int_{\Omega}s(\rho_{\delta,0},\theta_{\delta,0})\phi(0,\cdot)dx\end{split}$
(4.55)

for all $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$.

Thus we can conclude that after the first limit the quantities satisfy the
following system of equations.

The approximate continuity equation

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho
u)=\epsilon\Delta\rho~{}~{}\text{in }(0,T)\times\Omega,$ (4.56)

with homogeneous Neumann boundary condition

$\displaystyle\nabla\rho\cdot n=0~{}~{}\text{on }(0,T)\times\partial\Omega,$
(4.57)

and the initial condition

$\displaystyle\rho(0,\cdot)=\rho_{0,\delta}.$ (4.58)

The Brinkman-type equation

$\displaystyle\int_{\Omega}\rho\theta\operatorname{div}\phi
dx=\int_{\Omega}\mu\nabla u:\nabla\phi+\nu\rho u\cdot\phi dx,$ (4.59)

for any test function $\phi\in C_{c}^{\infty}(\overline{\Omega})$ with

$\displaystyle\phi=0~{}~{}\text{on }\partial\Omega~{}~{}\text{ no-slip
boundary conditions.}$ (4.60)

The approximate internal energy balance

$\displaystyle\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0,\delta}\theta_{0,\delta}dx+\int_{0}^{t}\int_{\Omega}\frac{\delta}{\theta^{2}}-\delta\theta^{5}\,dxd\tau.$
(4.61)

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

The approximate entropy inequality

$\displaystyle\begin{split}&\int_{0}^{T}\int_{\Omega}s(\rho,\theta)\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt+\int_{0}^{T}\int_{\Omega}\big{(}\frac{\nabla\theta}{\theta}-\epsilon(\log\theta-\log\rho-1)\nabla\rho\big{)}\cdot\nabla\phi
dxdt\\\ &+\int_{0}^{T}\int_{\Omega}\bigg{(}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\bigg{)}\phi
dxdt\\\
&+\int_{0}^{T}\int_{\Omega}\bigg{(}\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho}\bigg{)}\phi
dxdt\leq-\int_{\Omega}s(\rho_{\delta,0},\theta_{\delta,0})\phi(0,\cdot)dx\end{split}$
(4.62)

for all $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$ with $\phi\geq
0$.

Rewriting the last equation

$\displaystyle\begin{split}&\int_{\Omega}s(\rho_{\delta,0},\theta_{\delta,0})\phi(0,\cdot)dx-\int_{0}^{T}\int_{\Omega}s(\rho,\theta)\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt\\\
&+\int_{0}^{T}\int_{\Omega}\big{(}\frac{\nabla\theta}{\theta}-\epsilon(\log\theta-\log\rho-1)\nabla\rho\big{)}\cdot\nabla\phi
dxdt\\\
&\geq\int_{0}^{T}\int_{\Omega}\bigg{(}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\bigg{)}\phi
dxdt\\\
&+\int_{0}^{T}\int_{\Omega}\bigg{(}\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho}\bigg{)}\phi
dxdt\end{split}$

for all $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$ with $\phi\geq
0$ we note that the left-hand side of the equation can be understood as a non-
negative linear form on the space of smooth function with compact support in
$[0,T)\times\overline{\Omega}$. By the Riesz representation theorem, there
exists a regular, non-negative Borel measure $\Sigma_{\epsilon,\delta}$ on
$[0,T)\times\overline{\Omega}$ that can be extended to
$[0,T]\times\overline{\Omega}$ such that

$\displaystyle\begin{split}&\int_{0}^{T}\int_{\Omega}s(\rho,\theta)\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt+\int_{0}^{T}\int_{\Omega}\big{(}\frac{\nabla\theta}{\theta}-\epsilon(\log\theta-\log\rho-1)\nabla\rho\big{)}\cdot\nabla\phi
dxdt\\\
&+\langle\Sigma_{\epsilon,\delta},\phi\rangle_{[\mathcal{M},C]([0,T]\times\overline{\Omega})}=-\int_{\Omega}s(\rho_{\delta,0},\theta_{\delta,0})\phi(0,\cdot)dx\end{split}$
(4.63)

for all $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$ with $\phi\geq
0$. Moreover,

$\displaystyle\Sigma_{\epsilon,\delta}\geq$
$\displaystyle\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}$
$\displaystyle+\epsilon\delta(\Gamma\rho^{\Gamma-2}+2)|\nabla\rho|^{2}+\epsilon\frac{|\nabla\rho|^{2}}{\rho}.$

####  4.4.2 Limit $\epsilon\to 0$

The next step is to let $\epsilon\to 0$ in the approximate system. To this
end, let $(\rho_{\epsilon},u_{\epsilon},\theta_{\epsilon})$ be a solution of
equations (4.56)-(4.63).

Then, similar to the previous part we obtain the following estimates
independent of $\epsilon$:

$\displaystyle\sup_{\epsilon>0}\bigg{\\{}\operatorname*{ess\,sup}_{t\in(0,T)}\int_{\Omega}H_{\delta,\overline{\theta}}(t)\,dx\bigg{\\}}<\infty$
(4.64)
$\displaystyle\sup_{\epsilon>0}\bigg{\\{}\Sigma_{\epsilon,\delta}\big{[}[0,T]\times\overline{\Omega}\big{]}\bigg{\\}}<\infty.$
(4.65)

This implies that

$\displaystyle\sup_{\epsilon>0}\bigg{\\{}\int_{0}^{T}\int_{\Omega}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+(\frac{\kappa}{\theta}+\delta\theta^{\Gamma-1})|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}\,dxdt\bigg{\\}}<\infty,$
(4.66)
$\displaystyle\sup_{\epsilon>0}\bigg{\\{}\epsilon\delta\int_{0}^{T}\int_{\Omega}\big{(}\Gamma\rho^{\Gamma-2}+2\big{)}|\nabla\rho|^{2}\,dxdt\bigg{\\}}<\infty,$
(4.67)
$\displaystyle\sup_{\epsilon>0}\bigg{\\{}\epsilon\int_{0}^{T}\int_{\Omega}\frac{|\nabla\rho|^{2}}{\rho}\,dxdt\bigg{\\}}<\infty.$
(4.68)

As in the previous part we conclude that

$\displaystyle\rho_{\epsilon}^{\Gamma/2}~{}~{}\text{bounded in
}L^{2}(0,T;H^{1}(\Omega)),$ (4.69)
$\displaystyle\theta_{\epsilon}^{\Gamma/2}~{}~{}\text{bounded in
}L^{2}(0,T;H^{1}(\Omega)).$ (4.70)

Thus from the Brinkman equation we obtain

$\displaystyle u_{\epsilon}~{}~{}\text{bounded in }H^{1}(\Omega),$ (4.71)

and

$\displaystyle u_{\epsilon}\to u~{}~{}\text{weakly in }H^{1}(\Omega).$ (4.72)

Passing to the limit in the Brinkman equation yields

$\displaystyle\int_{\Omega}\mu\nabla u:\nabla\varphi dx+\int_{\Omega}\nu\rho
u\cdot\varphi dx-\int_{\Omega}\rho\theta\operatorname{div}\varphi dx=0,$
(4.73)

for any $\varphi\in C_{c}^{\infty}(\overline{\Omega})$ with
$\varphi|_{\partial\Omega}=0$.

Multiplying the approximate continuity equation by $\rho_{\epsilon}$ and
integrating by parts yields

$\displaystyle\frac{1}{2}\int_{\Omega}\rho_{\epsilon}^{2}(t)dx+\int_{0}^{t}\int_{\Omega}|\nabla\rho_{\epsilon}|^{2}dxdt=\frac{1}{2}\int_{\Omega}\rho_{0,\delta}^{2}dx-\frac{1}{2}\int_{0}^{t}\int_{\Omega}\rho_{\epsilon}^{2}\operatorname{div}u_{\epsilon}dxdt.$

Thus, we observe that

$\displaystyle\sqrt{\epsilon}\nabla\rho_{\epsilon}~{}~{}\text{is bounded in
}L^{2}(0,T;L^{2}(\Omega)),$ (4.74) $\displaystyle\rho_{\epsilon}~{}~{}\text{is
bounded in }L^{\infty}(0,T;L^{2}(\Omega)).$ (4.75)

and in particular

$\displaystyle\epsilon\nabla\rho_{\epsilon}\to 0~{}~{}\text{in
}L^{2}(0,T;L^{2}(\Omega)).$ (4.76)

This yields that

$\displaystyle\rho_{\epsilon}u_{\epsilon}\to\rho u~{}~{}\text{weakly in
}L^{2}(0,T;L^{2}(\Omega)).$ (4.77)

Thus, we can pass to the limit in the approximate continuity equation as
$\epsilon\to 0$ and the limit $\rho$ satisfies the integral identity

$\displaystyle\int_{0}^{T}\int_{\Omega}\big{(}\rho\partial_{t}\phi+\rho
u\cdot\nabla\phi\big{)}dxdt+\int_{\Omega}\rho_{0,\delta}dx=0$ (4.78)

for any test function $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$
i.e. $\rho,\,u$ satisfy in the following equation in the sense of
distributions

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0.$ (4.79)

It remains to pass to the limit in the approximate entropy equation. To this
end, we once more need to show strong convergence of the absolute temperature.
The idea is again to show uniform estimates for $\theta_{\epsilon}$ and then
apply the Div-Curl Lemma.

Taking a closer look at the bounds of the entropy production rate, we derive
that

$\displaystyle\theta_{\epsilon}^{\Gamma/2}~{}~{}\text{is bounded in
}L^{2}(0,T;H^{1}(\Omega)),$ $\displaystyle\theta_{\epsilon}^{-1}~{}~{}\text{is
bounded in }L^{3}((0,T)\times\Omega),$
$\displaystyle\theta_{\epsilon}~{}~{}\text{is bounded in
}L^{5}((0,T)\times\Omega),$ $\displaystyle\log\theta_{\epsilon}~{}~{}\text{is
bounded in }L^{2}(0,T:H^{1}(\Omega))\cap L^{q}((0,T)\times\Omega).$

Now, for the application of the Div-Curl Lemma we use the same idea as in the
previous limit case.

Setting

$\displaystyle
U_{\epsilon}=\bigg{[}s(\rho_{\epsilon},\theta_{\epsilon}),s(\rho_{\epsilon},\theta_{\epsilon})u_{\epsilon}+\frac{\nabla\theta_{\epsilon}}{\theta_{\epsilon}}+\epsilon\big{(}\log\theta_{\epsilon}-\log\rho_{\epsilon}-1\big{)}\nabla\rho_{\epsilon}\bigg{]},~{}~{}V_{\epsilon}=[\theta_{\epsilon},0,0,0]$
(4.80)

we observe that

$\displaystyle\operatorname{div}U_{\epsilon}=\Sigma_{\epsilon,\delta},~{}~{}\text{and
}\operatorname{curl}V_{\epsilon}$

are relatively precompact in $W^{-1,s}(\Omega)$ for $s\in[1,\frac{3}{2})$. The
boundedness of $U_{\epsilon}$ and $V_{\epsilon}$ in $L^{p}((0,T)\times\Omega)$
for some $p>1$ can be shown as follows.

The sequence $\theta_{\epsilon}$ is bounded in $L^{2}((0,T)\times\Omega)$ and
for the sequence $U_{\epsilon}$ we use the uniform estimates (4.64)-(4.68) and
the special structure of $s(\rho_{\epsilon},\theta_{\epsilon})$ to conclude
that it is bounded in $L^{p}((0,T)\times\Omega)$ for some $p>1$. Moreover,
$\epsilon(\log\theta_{\epsilon}-\log\rho_{\epsilon}-1)\to 0$ weakly in
$L^{p}((0,T)\times\Omega)$.

Hence, we obtain that

$\displaystyle\overline{s(\rho,\theta)\theta}=\overline{s(\rho,\theta)}\theta.$
(4.81)

By the monotonicity of the entropy and the weak convergence we conclude that
up to a subsequence

$\displaystyle\theta_{\epsilon}\to\theta~{}~{}\text{a.a. in
}(0,T)\times\Omega.$ (4.82)

In addition, we have that the limit temperature is positive a.a. on the set
$(0,T)\times\Omega$, more precisely

$\displaystyle\theta^{-3}\in L^{1}((0,T)\times\Omega).$ (4.83)

Now, we can let $\epsilon\to 0$ in the approximate entropy equation.
Using the previous relations we obtain that

$\displaystyle\frac{\kappa\nabla\theta_{\epsilon}}{\theta_{\epsilon}}\to\frac{\kappa\nabla\theta}{\theta}~{}~{}\text{weakly
in }L^{p}(0,T;\Omega),$

for some $p>1$.
Applying the Div-Curl Lemma once more with $V_{\epsilon}=[u_{\epsilon},0,0,0]$
and $U_{\epsilon}$ as before, we observe that

$\displaystyle
s(\rho_{\epsilon},\theta_{\epsilon})u_{\epsilon}\to\overline{s(\rho,\theta)}u~{}~{}\text{weakly
in }L^{p}((0,T)\times\Omega)$

for some $p>1$.
The terms appearing in $\sigma_{\epsilon,\delta}$ are weakly lower semi-
continuous as established in the previous part. Moreover, the
$\epsilon$-dependent terms are non-negative. Hence we can conclude that

$\displaystyle\Sigma_{\epsilon,\delta}\to\sigma_{\delta}~{}~{}\text{weakly in
}\mathcal{M}\big{(}[0,T]\times\overline{\Omega}\big{)},$

where $\sigma_{\delta}$ is a positive measure on
$[0,T]\times\overline{\Omega}$ satisfying

$\displaystyle\sigma_{\delta}\geq\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\big{(}\frac{\kappa}{\theta}+\frac{\delta}{2}\theta^{\Gamma-2}\big{)}|\nabla\theta|^{2}\big{)}+\frac{\delta}{\theta^{3}}+\delta\theta^{5}.$

It remains to show that $\overline{s(\rho,\theta)}=s(\rho,\theta)$. To this
end, we have to show the strong convergence in the densities. We follow the
ideas presented in [FN09, Chapter 3.6].

The first step is to introduce the test function $\varphi(x)$ for the momentum
equation, where

$\displaystyle\varphi(x)=\xi(x)\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho_{\epsilon}]~{}~{}\text{and
}\xi\in C_{c}^{\infty}(\Omega).$

Taking $\varphi$ as an admissible test function in the Brinkman-type equation
(4.59) yields

$\displaystyle\int_{\Omega}\rho\epsilon\theta_{\epsilon}\operatorname{div}\varphi
dx=\int_{\Omega}\mu\nabla
u_{\epsilon}:\nabla\varphi+\nu\rho_{\epsilon}u_{\epsilon}\cdot\varphi dx$

Taking into account the specific form of the test function and integration by
parts yields

$\displaystyle\int_{\Omega}\xi\bigg{(}\rho_{\epsilon}\theta_{\epsilon}\rho_{\epsilon}-\mu\nabla
u_{\epsilon}:R[\mathbb{1}_{\Omega}\rho_{\epsilon}]-\nu\rho_{\epsilon}u_{\epsilon}\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho_{\epsilon}]\bigg{)}dx=\sum_{i=1}^{2}I_{i,\epsilon},$
(4.84)

where

$\displaystyle I_{1,\epsilon}$
$\displaystyle=-\int_{\Omega}\rho_{\epsilon}\theta_{\epsilon}\nabla\xi\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho_{\epsilon}]dx$
$\displaystyle I_{2,\epsilon}$ $\displaystyle=\mu\int_{\Omega}\nabla
u_{\epsilon}:\nabla\xi\otimes\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho_{\epsilon}]dx$

and the symbol $R$ denotes the double Riesz transform.
Repeating the same argument for the limit of the Brinkman equation (4.73) with
the test function $\varphi$ being

$\displaystyle\varphi(x)=\xi(x)\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho]~{}~{}\text{and
}\xi\in C_{c}^{\infty}(\Omega).$

yields the following

$\displaystyle\int_{\Omega}\xi\bigg{(}\rho\theta\rho-\mu\nabla
u:R[\mathbb{1}_{\Omega}\rho]-\nu\rho
u\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho]\bigg{)}dx=\sum_{i=1}^{2}I_{i},$
(4.85)

where

$\displaystyle I_{1}$
$\displaystyle=-\int_{\Omega}\rho\theta\nabla\xi\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho]dx$
$\displaystyle I_{2}$ $\displaystyle=\mu\int_{\Omega}\nabla
u:\nabla\xi\otimes\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho]dx.$

From the previous estimates we recall that

$\displaystyle\rho_{\epsilon}\to\rho~{}~{}\text{ in
}C_{weak}([0,T],L^{2}(\Omega)).$

Taking into account relations (4.64) - (4.77) we observe that the integral
$I_{i,\epsilon}$ converges to its counterpart $I_{i}$ for $i=1,2$ and we infer

$\displaystyle\lim_{\epsilon\to 0}$
$\displaystyle\int_{\Omega}\xi\bigg{(}\rho_{\epsilon}\theta_{\epsilon}\rho_{\epsilon}-\mu\nabla
u_{\epsilon}:R[\mathbb{1}_{\Omega}\rho_{\epsilon}]-\nu\rho_{\epsilon}u_{\epsilon}\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho_{\epsilon}]\bigg{)}dx$
(4.86) $\displaystyle=\int_{\Omega}\xi\bigg{(}\rho\theta\rho-\mu\nabla
u:R[\mathbb{1}_{\Omega}\rho]-\nu\rho
u\cdot\nabla\Delta^{-1}[\mathbb{1}_{\Omega}\rho]\bigg{)}dx$ (4.87)

Moreover, the last terms on the left-hand side of the equality converges to
the last term on the right-hand side. The next step is to rewrite the
following term

$\displaystyle\int_{\Omega}\xi\mu\nabla
u:R[\mathbb{1}_{\Omega}\rho]dx=\int_{\Omega}\mu R:[\xi\nabla u]\rho dx$

where we used the properties of the double Riesz transform [FN09, Chapter 11]
and we observe that we can write

$\displaystyle\mu R:[\xi\nabla u]=\mu\xi\operatorname{div}u+\mu\omega(u),$

where $\omega(u)=R:[\xi\nabla u]-\xi R:[\nabla u]$ is the commutator. Applying
a result by Coifman and Meyer [CM75] and the previous bounds we obtain that

$\displaystyle\omega(u_{\epsilon})\rho_{\epsilon}\to\overline{\omega(u)}\rho~{}~{}\text{weakly
in }L^{1}((0,T)\times\Omega).$ (4.88)

This yields

$\displaystyle\overline{\omega(u)}=\omega(u).$

The proof of the convergence in (4.88) is shown by applying the Div-Curl Lemma
to

$U_{\epsilon}=[\rho_{\epsilon},\rho_{\epsilon}u_{\epsilon}]~{}~{}\text{and
}V_{\epsilon}=[\omega(u_{\epsilon}),0,0,0]$

Hence we obtain the following weak compactness identity for the effective
pressure

$\displaystyle\overline{\rho\theta\rho}-\mu\overline{\rho\operatorname{div}u}=\overline{\rho\theta}\rho-\mu\rho\operatorname{div}u.$
(4.89)

The final step is to multiply the continuity equation on
$G^{\prime}(\rho_{\epsilon})$, with $G$ being a smooth and convex function.
Then as $\epsilon\to 0$ we get

$\displaystyle\int_{\Omega}\overline{G(\rho)}(t)dx+\int_{0}^{t}\int_{\Omega}\overline{\big{(}G^{\prime}(\rho)\rho-G(\rho)\big{)}\operatorname{div}u}dxdt\leq\int_{\Omega}G(\rho_{0})dx$

for all $t\in(0,T)$ and we deduce that

$\displaystyle\int_{\Omega}\overline{\rho\log\rho}(t)dx+\int_{0}^{t}\int_{\Omega}\overline{\rho\operatorname{div}u}dxdt=\int_{\Omega}\rho_{0}\log\rho_{0}dx.$

Via the theory of renormalized solutions by DiPerna and Lions [DL89] we obtain

$\displaystyle\int_{\Omega}\rho\log\rho(t)dx+\int_{0}^{t}\int_{\Omega}\rho\operatorname{div}udxdt\leq\int_{\Omega}\rho_{0}\log\rho_{0}dx.$

Hence we obtain

$\displaystyle\overline{\rho\theta\rho}\geq\overline{\rho\theta}\rho$ (4.90)

and as a consequence of of equation (4.89)

$\displaystyle\overline{\rho\operatorname{div}u}\geq\rho\operatorname{div}u.$
(4.91)

Combining both estimates implies

$\displaystyle\overline{\rho\log\rho}=\rho\log\rho$ (4.92)

which yields the desired strong convergence of the density as the function
$x\log x$ is convex, i.e.

$\displaystyle\rho_{\epsilon}\to\rho~{}~{}\text{ a.a. in }(0,T)\times\Omega.$
(4.93)

This allows us to identify $\overline{s(\rho,\theta)}=s(\rho,\theta)$.

Having eliminated the $\epsilon$-dependent terms, we summarize the results.
For any $\delta>0$ we have constructed a trio $(\rho,u,\theta)$ solving the
following equations.

The continuity equation

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0~{}~{}\text{in
}(0,T)\times\Omega,$ (4.94)

with homogeneous Neumann boundary condition

$\displaystyle\nabla\rho\cdot n=0~{}~{}\text{on }(0,T)\times\partial\Omega,$
(4.95)

and the initial condition

$\displaystyle\rho(0,\cdot)=\rho_{0,\delta}.$ (4.96)

The Brinkman-type equation

$\displaystyle\int_{\Omega}\rho\theta\operatorname{div}\phi
dx=\int_{\Omega}\mu\nabla u:\nabla\phi dx+\int_{\Omega}\nu\rho u\cdot\phi dx,$
(4.97)

for any test function $\phi\in C_{c}^{\infty}(\overline{\Omega})$ with

$\displaystyle\phi=0~{}~{}\text{on }\partial\Omega~{}~{}\text{ no-slip
boundary conditions.}$ (4.98)

The approximate internal energy balance

$\displaystyle\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0,\delta}\theta_{0,\delta}dx+\int_{0}^{t}\int_{\Omega}\frac{\delta}{\theta^{2}}-\delta\theta^{5}\,dxd\tau.$
(4.99)

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

The approximate entropy inequality

$\displaystyle\begin{split}&\int_{0}^{T}\int_{\Omega}s(\rho,\theta)\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt+\int_{0}^{T}\int_{\Omega}\frac{\nabla\theta}{\theta}\cdot\nabla\phi
dxdt\\\
&+\langle\sigma_{\delta},\phi\rangle_{\mathcal{M}([0,T]\times\overline{\Omega}]}=-\int_{\Omega}s(\rho_{\delta,0},\theta_{\delta,0})\phi(0,\cdot)dx\end{split}$
(4.100)

for all $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$ where

$\displaystyle\sigma_{\delta}\geq\frac{1}{\theta}\bigg{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\big{(}\frac{1}{\theta}+\frac{\delta}{2}\theta^{\Gamma-1}\big{)}|\nabla\theta|^{2}+\frac{\delta}{\theta^{2}}+\delta\theta^{5}\bigg{)}.$
(4.101)

####  4.4.3 Limit $\delta\to 0$

The last step in this proof is to let $\delta\to 0$. To this end, let
$(\rho_{\delta},u_{\delta},\theta_{\delta})$ be a solution to the approximate
system (4.94)-(4.101). We recall that the total mass of the system is
conserved, i.e.

$\displaystyle\int_{\Omega}\rho_{\delta}(t,\cdot)dx=\int_{\Omega}\rho_{0,\delta}dx~{}~{}\text{for
any }t\in[0,T].$ (4.102)

We assume that

$\displaystyle\rho_{0,\delta}\to\rho_{0}~{}~{}\text{in }L^{1}(\Omega),$
(4.103)

and thus the bound is uniform for $\delta\to 0$.

Next, we apply the reverse Young’s inequality to the energy balance equation

$\displaystyle\int_{\Omega}\frac{1}{S(\rho/\theta)^{2}}\frac{1}{4}(\rho^{2}+\theta^{2})(t)dx\leq\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0,\delta}\theta_{0,\delta}dx+\int_{0}^{t}\int_{\Omega}\frac{\delta}{\theta^{2}}-\delta\theta^{5}dxd\tau,$
(4.104)

where $S(h)=S(1/h)>0$ for $h>0$ is the Specht radius [Tom02]. Since
$S(h)\to\infty$ as $h\to 0$ we need the ratio $\rho/\theta$ to be bounded from
below and above. This follows from the fact that we chose positive initial
data $\rho_{0}$ and $\theta_{0}$ and from the previous sections. As a
consequence the solution may exist only for a short time $T^{*}>0$.

Next, is the dissipation balance

$\displaystyle\int_{\Omega}H_{\overline{\theta}}(\rho,\theta)(t)dx+\overline{\theta}\sigma_{\delta}\big{[}[0,t]\times\overline{\Omega}\big{]}=\int_{\Omega}H_{\overline{\theta}}(\rho_{0,\delta},\theta_{0,\delta})dx+\int_{0}^{t}\int_{\Omega}\frac{\delta}{\theta^{2}}+\delta\theta^{4}dxd\tau$
(4.105)

satisfied for a.a. $t\in[0,T]$. Noting that the terms $\delta/\theta^{2}$ and
$\delta\theta^{4}$ are absorbed in the entropy production $\sigma_{\delta}$
and the uniform bounds

$\displaystyle\int_{\Omega}H_{\overline{\theta}}(\rho_{0,\delta},\theta_{0,\delta})dx\leq
c~{}~{}\text{uniformly for }\delta\to 0$ (4.106)

hold, we obtain the following uniform estimates depending only on the initial
data:

$\displaystyle\operatorname*{ess\,sup}_{t\in(0,T)}\|\rho_{\delta}(t)\|_{L^{2}(\Omega)}$
$\displaystyle\leq c,$ (4.107)
$\displaystyle\operatorname*{ess\,sup}_{t\in(0,T)}\|\theta_{\delta}(t)\|_{L^{2}(\Omega)}$
$\displaystyle\leq c,$ (4.108)
$\displaystyle\sigma_{\delta}\big{[}[0,t]\times\overline{\Omega}\big{]}$
$\displaystyle\leq c.$ (4.109)

have that

$\displaystyle\int_{0}^{T}\int_{\Omega}|\nabla\log\theta_{\delta}|^{2}dxdt\leq
c,$ (4.110) $\displaystyle\int_{0}^{T}\int_{\Omega}\frac{|\nabla
u|^{2}+\rho|u|^{2}}{\theta_{\delta}}dxdt\leq c,$ (4.111)
$\displaystyle\delta\int_{0}^{T}\int_{\Omega}\frac{1}{\theta_{\delta}^{3}}+\theta^{5}dxdt\leq
c,$ (4.112)
$\displaystyle\delta\int_{0}^{T}\int_{\Omega}\theta^{\Gamma-2}|\nabla\theta_{\delta}|^{2}dxdt\leq
c.$ (4.113)

Using that

$\displaystyle\nabla\log\theta_{\delta}~{}~{}\text{bounded in
}L^{2}((0,T)\times\Omega)$

we obtain that

$\displaystyle\nabla\theta_{\delta}~{}~{}\text{bounded in
}L^{2}((0,T)\times\Omega).$ (4.114)

Moreover applying Lemma 3.1 and Proposition 3.2 we have that

$\displaystyle\log\theta_{\delta}~{}~{}\text{bounded in
}L^{2}(0,T;H^{1}(\Omega)).$ (4.115)

Combining the above estimates, we see that

$\displaystyle\theta_{\delta}~{}~{}\text{is bounded in
}L^{2}(0,T;H^{1}(\Omega))$ (4.116)

and especially by the Sobolev embedding

$\displaystyle\theta_{\delta}~{}~{}\text{is bounded in
}L^{2}(0,T;L^{6}(\Omega)).$ (4.117)

Next, we estimate

$\displaystyle\int_{\Omega}|\nabla
u_{\delta}|^{3/2}=\int_{\Omega}\frac{|\nabla
u_{\delta}|^{3/2}}{\theta_{\delta}^{3/4}}\theta_{\delta}^{3/4}\leq
c\int_{\Omega}\frac{|\nabla
u_{\delta}|^{2}}{\theta}+c\int_{\Omega}\theta^{4}_{\delta}\leq c.$ (4.118)

Using the assumption that $\rho_{0}\in L^{3}(\Omega)$ and testing the
continuity equation with $\rho^{2}$ we obtain

$\displaystyle\int_{\Omega}\partial_{t}\rho\rho^{2}dx+\int_{\Omega}\operatorname{div}(\rho
u)\rho^{2}dx=0.$ Integrating by parts yields
$\displaystyle\frac{1}{3}\frac{d}{dt}\int_{\Omega}\rho^{3}dx\leq\frac{2}{3}\int_{\Omega}\rho^{3}|\operatorname{div}u|\,dx.$

Thus, we obtain

$\displaystyle\frac{d}{dt}\|\rho_{\delta}\|_{L^{3}(\Omega)}^{3}\leq
c\big{(}\|\rho_{\delta}\|_{L^{3}(\Omega)}^{3}\big{)}^{3}+c\|\nabla
u_{\delta}\|_{L^{3/2}(\Omega)}^{3/2},$ (4.119)

where we used an inverse type of the Jensen inequality [Tak+99]. Note, this
inequality requires lower and upper bounds on the density $\rho_{\delta}$
almost everywhere. Thus, we have an ordinary differential equation of the type

$\displaystyle x^{\prime}\leq c_{1}x^{3}+C_{2}.$

This is an ODE with Lipschitz right-hand side and by Picarc-Lindeloeff Theorem
the solution exists for a small time $T$. Hence,

$\displaystyle\rho_{\delta}~{}~{}\text{is bounded in
}L^{\infty}(0,T;L^{3}(\Omega).$ (4.120)

Now, we can conclude from the Brinkman equation that

$\displaystyle\|\rho u^{2}\|_{L^{1}(\Omega)}+\|\nabla
u_{\delta}\|^{2}_{L^{2}(\Omega)}\leq\int_{\Omega}\rho_{\delta}^{2}\theta^{2}_{\delta}\,dx\leq
c\|\rho_{\delta}\|_{L^{3}(\Omega)}^{3}+\|\theta_{\delta}\|^{6}_{L^{6}(\Omega)},$

where we used that $\rho$ is non-negative at least for a small time $T^{*}$.
And thus by the Poincare inequality

$\displaystyle u_{\delta}~{}~{}\text{is bounded in }H^{1}_{0}(\Omega).$
(4.121)

We remark that by repeating the previous two steps we obtain that

$\displaystyle\rho_{\delta}~{}~{}\text{is bounded in
}L^{\infty}(0,T;L^{p}(\Omega)$ (4.122)

for some $p>2$ at least for a small time $T>0$.

Now, we have all the necessary uniform estimates together in order to pass to
the limit in the equations.

For the continuity equation we get

$\displaystyle\rho_{\delta}\to\rho~{}~{}\text{weakly in
}L^{p}((0,T)\times\Omega)$ (4.123)

for some $p\geq 2$ and similarly

$\displaystyle\rho_{\delta}u_{\delta}\to\rho u~{}~{}\text{weakly in
}L^{p}((0,T)\times\Omega)$ (4.124)

for some $p>1$.

After passing to the limit in the continuity equation as $\delta\to 0$ the
limit satisfies the integral identity

$\displaystyle\int_{0}^{T}\int_{\Omega}\big{(}\rho\partial_{t}\phi+\rho
u\cdot\nabla\phi\big{)}dxdt+\int_{\Omega}\rho_{0}dx=0$ (4.125)

for any test function $\phi\in C_{c}^{\infty}([0,T)\times\overline{\Omega})$
i.e. $\rho,\,u$ satisfy in the following equation in the sense of
distributions

$\displaystyle\partial_{t}\rho+\operatorname{div}(\rho u)=0.$ (4.126)

For the Stokes equation we have

$\displaystyle u_{\delta}\to u~{}~{}\text{weakly in }H^{1}(\Omega)$ (4.127)

and

$\displaystyle\rho_{\delta}\theta_{\delta}\to\overline{\rho\theta}~{}~{}\text{weakly
in }L^{p}((0,T)\times\Omega)$ (4.128)

for some $p>1$. Thus the limit satisfies

$\displaystyle\int_{\Omega}\nu\rho u\cdot\varphi+\int_{\Omega}\mu\nabla
u:\nabla\varphi dx+\int_{\Omega}\rho\theta\operatorname{div}\varphi dx=0,$
(4.129)

for any $\varphi\in C_{c}^{\infty}(\overline{\Omega})$ with
$\varphi|_{\partial\Omega}=0$ or in the sense of distributions

$\displaystyle-\nu\rho u+\mu\Delta u=\nabla(\rho\theta).$ (4.130)

Before passing to the limit in the approximate entropy equation we need to
show the pointwise convergence of the temperature. Again as in the previous
sections the idea is to apply the Div-Curl Lemma. Setting

$\displaystyle
U_{\delta}=\bigg{[}s(\rho_{\delta},\theta_{\delta}),s(\rho_{\delta},\theta_{\delta})u_{\delta}+\frac{\nabla\theta_{\delta}}{\theta_{\delta}}],~{}~{}V_{\delta}=[\theta_{\delta},0,0,0]$
(4.131)

we observe that

$\displaystyle\operatorname{div}U_{\delta}=\sigma_{\delta},~{}~{}\text{and
}\operatorname{curl}V_{\delta}$

are relatively precompact in $W^{-1,s}(\Omega)$ for $s\in[1,\frac{3}{2})$.
Indeed,using equations (4.114)-(4.117) we see that

$\displaystyle\delta\theta^{\Gamma-1}_{\delta}\nabla\theta_{\delta}=\delta\frac{\Gamma}{2}\theta_{\delta}^{1/4}\theta_{\delta}^{\Gamma/2-1/4}\nabla\theta_{\delta}^{\Gamma/2}.$

Hence, we can conclude

$\displaystyle\delta\theta^{\Gamma-1}_{\delta}\nabla\theta_{\delta}\to
0~{}~{}\text{in }L^{p}((0,T)\times\Omega)$

as $\delta\to 0$ for a certain $p>1$. In addition, since $\theta_{\delta}>0$
for a.a. $(t,x)$ we obtain

$\displaystyle\delta\frac{1}{\theta_{\delta}^{3}}\to 0~{}~{}\text{in
}L^{1}((0,T)\times\Omega)$

as $\delta\to 0$.

The boundedness of $U_{\delta}$ in $L^{p}((0,T)\times\Omega)$ for some $p>1$
can be shown as follows.
For the sequence $U_{\delta}$ we use the uniform estimates (4.107)-(4.109) and
the special structure of $s(\rho_{\delta},\theta_{\delta})$ to conclude that
it is bounded in $L^{p}((0,T)\times\Omega)$ for some $p>1$. Hence, we obtain
that

$\displaystyle\overline{s(\rho,\theta)\theta}=\overline{s(\rho,\theta)}\theta.$
(4.132)

By the monotonicity of the entropy and the weak convergence we conclude that
up to a subsequence

$\displaystyle\theta_{\delta}\to\theta~{}~{}\text{a.a. in }(0,T)\times\Omega.$
(4.133)

In addition, we have that the limit temperature is positive a.a. on the set
$(0,T)\times\Omega$.

It remains to show that $\overline{s(\rho,\theta)}=s(\rho,\theta)$. We proceed
as in the previous section by showing the strong convergence of the density.
From the bounds obtained in equations (4.107)-(4.117) we see that the methods
from the $\epsilon$-limit can be applied in this setting too, cf. equations
(4.84)-(4.93).

Using the weak lower semi-continuity of convex functionals, we can let
$\delta\to 0$ in the approximate entropy balance to conclude that

$\displaystyle\begin{split}&\int_{0}^{T}\int_{\Omega}\overline{s(\rho,\theta)}\big{(}\partial_{t}\phi+u\cdot\nabla\phi\big{)}dxdt-\int_{0}^{T}\int_{\Omega}\frac{\kappa\nabla\theta}{\theta}\cdot\nabla\phi
dxdt\\\
&+\langle\sigma;\phi\rangle_{[\mathcal{M};C]([0,T]\times\overline{\Omega})}=-\int_{\Omega}s(\rho_{0},\theta_{0})dx\end{split}$
(4.134)

for any $\phi\in C_{c}^{\infty}([0,T]\times\overline{\Omega})$. Here
$\sigma\in\mathcal{M}^{+}([0,T]\times\overline{\Omega})$ is a weak-*-limit in
the space of measures $\mathcal{M}([0,T]\times\overline{\Omega})$ of the
sequence $\sigma_{\delta}$. Using the lower weak semi-continuity of convex
functionals and the fact that all $\delta$-dependent terms in the entropy
production rate are non-negative, we obtain that

$\displaystyle\sigma\geq\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}+\frac{|\nabla\theta|^{2}}{\theta}\big{)}.$ (4.135)

The last step is to take the limit in the internal energy balance. By equation
(4.128) we can pass to the limit and obtain

$\displaystyle\int_{\Omega}\rho\theta(t)dx=\int_{\Omega}\rho_{0}\theta_{0}dx~{}~{}\text{for
a.a. }t\in[0,T].$ (4.136)

This completes the proof of the theorem.

###  4.5 Higher regularity

In the proof of Theorem 4.1 we have noted that the weak solutions constructed
by the approximate scheme satisfy better regularity and integrability
properties.

######  Theorem 4.11 (Regularity of weak solutions).

Let $\Omega\subset\mathbb{R}^{n}$, where $n=2,\,3$, be a bounded Lipschitz
domain. Assume that the initial data $\rho_{0}$, $E_{0}$ and $s_{0}$ satisfy
the hypothesis of Theorem 4.1.
Then, in addition to the minimal regularity assumptions required in equations
(3.2), (3.6), (3.10), (3.14)-(3.16), there holds:

  1. i)

The weak solution satisfies

$\displaystyle\rho\in C_{weak}([0,T];L^{3}(\Omega))\cap
C([0,T];L^{1}(\Omega)),$ (4.137) $\displaystyle u\in H^{1}_{0}(\Omega),$
(4.138) $\displaystyle\theta\in L^{2}(0,T;H^{1}(\Omega))\cap
L^{\infty}(0,T:l^{2}(\Omega)),$ (4.139) $\displaystyle\log\theta\in
L^{2}(0,T;H^{1}(\Omega)),$ (4.140)

  2. ii)

The entropy satisfies

$\displaystyle ess\lim_{t\to
0^{+}}\int_{\Omega}s(\rho,\theta)(t,\cdot)\phi\,dx\geq\int_{\Omega}s(\rho_{0},\theta_{0})\phi\,dx~{}~{}\text{for
any }\phi\in C_{c}^{\infty}(\overline{\Omega}),~{}\phi\geq 0.$ (4.141)

If in addition, $\theta_{0}\in W^{1,\infty}(\Omega)$ then

$\displaystyle ess\lim_{t\to
0^{+}}\int_{\Omega}s(\rho,\theta)(t,\cdot)\phi\,dx=\int_{\Omega}s(\rho_{0},\theta_{0})\phi\,dx~{}~{}\text{for
all }\phi\in C_{c}^{\infty}(\overline{\Omega}).$ (4.142)

###### Proof.

The integrability properties follow directly from the proof of the existence
of weak solutions.

The strong continuity of the density is a general property of the transport
equation in the context of renormalized solutions.

The last part of the proof follows step 3 in the proof of Theorem 3.2 in
[FN09]. ∎

##  5 Conclusion and Remarks

In this section we conclude with several remarks.

The first one is that in the first part of this paper we showed how to apply
the energetic variational approach [LLZ05] and [HL+10] in the setting of fluid
mechanics and combine it with the temperature in a natural way. This leads to
the general framework of the free energy as starting point for the
thermodynamics of fluids. With the choice of free energy (in terms of
temperature and the state variable, the phase field variable here), and the
entropy production, as well as the kinematics/transport of these variables,
one should be able to uniquely determine the system.

Second, we observe the importance of the Laplacian in the velocity term. To
this end , we recall the following estimates obtained from the continuity and
momentum equation and the entropy production rate.

$\displaystyle\int_{\Omega}\nu\rho|u|^{2}+\int_{\Omega}\mu|\nabla u|^{2}$
$\displaystyle=\int_{\Omega}\rho\theta\operatorname{div}u$ (5.1)
$\displaystyle\frac{1}{\gamma+1}\frac{d}{dt}\int_{\Omega}\rho^{\gamma+1}+\frac{\gamma}{\gamma+1}\int_{\Omega}\rho^{\gamma+1}\operatorname{div}u$
$\displaystyle=0$ (5.2)
$\displaystyle\int_{\Omega}\frac{1}{\theta}\big{(}\mu|\nabla
u|^{2}+\nu\rho|u|^{2}\big{)}$ $\displaystyle\leq c$ (5.3)

We note that equation (5.3) gives a uniform bound on $\frac{\mu|\nabla
u|^{2}}{\theta}$ and this bound is then used in equation (5.1) and (5.2) to
obtain further estimates. If we were to let the parameter $\mu$ go to $0$ we
would loose these bounds, i.e. the control of the gradient of $u$ and thus we
cannot achieve the ultimate goal to consider the ideal gas under a Darcy-type
law.

Next, we note that adding a memory/ evolutionary term to the Brinkman-type
equation, i.e. $(\rho u)_{t}$ does not change analysis of the model. In
addition, as noted earlier in Remark 4.2, if we had a kinetic term in the
total energy we would obtain the incompressible Navier-Stokes-Fourier system
for which the existence of weak solutions for the ideal gas case is still open
[FN09] and [FNS11]. Thus our result gives in a sense a ”lower limit” existence
result of the full compressible Navier-Stokes-Fourier system, where the
difference and crucial aspect is the additional nonlinear term in the momentum
equation.

Finally, we want to remark that for a similar system of equations with only
Darcy-type dissipation we are able to show the well-posedness of the system in
a critical Besov space [LS21]. The difference in these two approaches is that
the first one uses energy methods for finding the weak solution, whereas the
second one uses scaling arguments and the algebra structure provided by the
critical Besov space.

## Acknowledgments

The authors would like to thank Prof. Anja Schlömerkemper for constructive
suggestions and discussions. This research was supported in part by the
National Science Foundation Grant NSF DMS–1714401 and the United States–Israel
Binational Science Foundation Grant BSF 2024246.

##  Appendix A Appendix

The following result stating the weak convergence of a product of functions is
due to [Tar75].

######  Theorem A.1 (Div-Curl Lemma).

Let $Q\subset\mathbb{R}^{n}$ be an open set. Assume

$\displaystyle U_{n}\to U~{}~{}\text{weakly in }L^{p}(Q),$ $\displaystyle
V_{n}\to V~{}~{}\text{weakly in }L^{q}(Q),$

where

$\displaystyle\frac{1}{p}+\frac{1}{q}=\frac{1}{r}<1.$

In addition, let

$\displaystyle\operatorname{div}U_{n}\equiv\nabla\cdot
U_{n}~{}~{}\text{precompact in }W^{-1,s}(Q),$
$\displaystyle\operatorname{curl}V_{n}\equiv\nabla
V_{n}-\nabla^{T}V_{n}~{}~{}\text{precompact in }W^{-1,s}(Q),$

for some $s>1$. Then

$\displaystyle U_{n}\cdot V_{n}\to U\cdot V~{}~{}\text{wekly in }L^{r}(Q).$

######  Theorem A.2 (Reverse Young’s inequality [Tom02]).

Let $a$, $b$ be positive real numbers and let $\nu\in[0,1]$. Then

$\displaystyle S\big{(}\frac{a}{b}\big{)}a^{1-\nu}b^{\nu}\geq(1-\nu)a+\nu b,$

where the constant $S(h)$ is called the Specht ratio and is defined by

$\displaystyle S(h)=S(1/h)=\frac{h^{\frac{1}{h-1}}}{e\log
h^{\frac{1}{h-1}}}~{}~{}h\neq 1,~{}\text{for }h>0.$

The next theorem stating the existence and uniqueness of ODEs is from [Rou13].

######  Theorem A.3 (Carathéodory).

Let $T$ be fixed and $f:~{}I\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a
Carathéodory mapping satisfying the growth condition
$|f(t,r)|\leq\gamma(t)+C|r|$ with some $\gamma\in L^{1}(I)$. Then:

  * i)

The initial-value problem

$\frac{du}{dt}=f(t,u(t))~{}\text{ for a.e. }~{}t\in I,~{}~{}u|_{t=0}=u_{0}$

has a solution $u\in W^{1,1}(I;\mathbb{R}^{n})$ on the interval $I=[0,T]$.

  * ii)

If $f(t,\cdot)$ is also Lipschitz continuous in the sense
$|f(t,r_{1})-f(t,r_{2})|\leq l(t)|r_{1}-r_{2}|$ with some $l\in L^{1}(I)$,
then the solution is unique.

The next theorem is an inverse of the Jensen’s inequality [Tak+99].

######  Theorem A.4 (Inverse Jensen’s inequality).

Let $f$ be a measurable function on a probability measure space
$(\Omega,\mathcal{F},\mu,)$ with $f(\Omega)\subset[m,M]$. Then for $p\geq 1$

$\displaystyle\int f^{p}d\mu\leq\alpha\bigg{(}\int fd\mu\bigg{)}^{p}+\beta$

holds for some $\alpha>0$ and $\beta=a(1-\frac{1}{p})x_{0}+b$ where

$\displaystyle
a=\frac{M^{p}-m^{p}}{M-m},~{}~{}b=\frac{Mm^{p}-mM^{p}}{M-m},~{}~{}m<x_{0}<M.$

The following two important theorems can be found in [FN09].

######  Theorem A.5.

Let $I\subset\mathbb{R}$ be an interval, $Q\subset\mathbb{R}^{n}$ a domain and
$(P,G)\in C(I)\times C(I)$ a couple of non-decreasing functions. Assume that
$\rho_{n}\in L^{1}(Q;I)$ is a sequence such that

$\displaystyle\begin{rcases}P(\rho_{n})\to\overline{P(\rho)}&\\\
G(\rho_{n})\to\overline{G(\rho)}&\\\
P(\rho_{n})G(\rho_{n})\to\overline{P(\rho)G(\rho)}&\end{rcases}~{}~{}\text{weakly
in }L^{1}(Q).$

  * i)

Then

$\displaystyle\overline{P(\rho)}\overline{G(\rho)}\leq\overline{P(\rho)G(\rho)}.$

  * ii)

If, in addition, $G\in C(\mathbb{R})$, $G(\mathbb{R})=\mathbb{R}$, $G$
strictly increasing and $P\in C(\mathbb{R})$, $P$ is non-decreasing then

$\displaystyle\overline{P(\rho)}=P\circ
G^{-1}\big{(}\overline{G(\rho)}\big{)}.$

  * iii)

In particular, if $G(z)=z$, then

$\displaystyle\overline{P(\rho)}=P(\rho).$

######  Theorem A.6.

Let $Q\subset\mathbb{R}^{n}$ be a measurable set and $\\{v_{n}\\}_{n}$ a
sequence of functions in $L^{1}(Q)$ such that

$\displaystyle v_{n}\to v~{}~{}\text{weakly in }L^{1}(Q).$

Let $\Phi:\mathbb{R}^{m}\to(-\infty,\infty]$ be a lower semi-continuous convex
function. Then

$\displaystyle\int_{Q}\Phi(v)dx\leq\liminf_{n\to\infty}\int_{Q}\Phi(v_{n})dx.$

Moreover, if

$\displaystyle\Phi(v_{n})\to\overline{\Phi(v)}~{}~{}\text{weakly in
}L^{1}(Q),$

then

$\displaystyle\Phi(v)\leq\overline{\Phi(v)}~{}~{}\text{a.a. on }Q.$

If, in addition, $\Phi$ is a strictly convex on an open convex set
$U\subset\mathbb{R}^{m})$, and

$\displaystyle\Phi(v)=\overline{\Phi(v)}~{}~{}\text{a.a. on }Q,$

then

$\displaystyle v_{n}(y)\to v(y)~{}~{}\text{for a.a. }y\in\\{y\in Q|v(y)\in
U\\}$

extracting a subsequence as the case may be.

## References

  * [Ann15] F. De Anna  “Global Weak Solutions for Boussinesq System with Temperature dependent Viscosity and bounded Temperature”  In _Adv. Differential Equations_ 21, 2015, pp. 1001–1048
  * [AL19] F. De Anna and C. Liu  “Non-isothermal General Ericksen–Leslie System: Derivation, Analysis and Thermodynamic Consistency”  In _Archive for Rational Mechanics and Analysis_ 231, 2019, pp. 637–717
  * [Bai99] R. Baierlein  “Thermal Physics”  Cambridge University Press, Cambridge, 1999
  * [BRR00] R.S. Berry, S.A. Rice and J. Ross  “Physical Chemistry”  Oxford University Press, Oxford, 2000
  * [Bir94] G. A. Bird  “Molecular Gas Dynamics And The Direct Simulation Of Gas Flows”  Clarendon Press, Oxford, 1994
  * [Bri47] H.C. Brinkman  “A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles”  In _Appl. Sci. Res., Sect. A_ 1, 1947, pp. 27
  * [BH15] M. Bulíček and J. Havrda  “On existence of weak solution to a model describing incompressible mixtures with thermal diffusion cross effects”  In _Z. Angew. Math. Mech._ 95, 2015, pp. 589–619
  * [Bul+20] M. Buliček, A. Jüngel, M. Pokornỳ and N. Zamponi  “Existence analysis of a stationary compressible fluid model for heat-conducting and chemically reacting mixtures”  In _arXiv preprint arXiv:2001.06082_ , 2020
  * [CV02] V.V Chepyzhov and M. I. Vishik  “Attractors for equations of mathematical physics”  American Mathematical Soc., Providence, 2002
  * [CM75] R. Coifman and Y. Meyer  “On commutators of singular integrals and bilinear singular integrals”  In _Transactions of the American Mathematical Society_ 212, 1975, pp. 315–331
  * [Daf79] C.M. Dafermos  “The second law of thermodynamics and stability”  In _Arch. Rational Mech. Anal._ 70, 1979
  * [Dan01] R. Danchin  “Global Existence in CriticalSpaces for Flows of Compressible Viscous and Heat-Conductive Gases”  In _Arch. Rational Mech. Anal._ 160, 2001, pp. 1–39
  * [DL89] R. DiPerna and P.-L. Lions  “Ordinary differential equations, transport theory and Sobolev spaces”  In _Inventiones mathematicae_ 98.3, 1989, pp. 511–547
  * [Dre+20] W. Dreyer, P.-E. Druet, P. Gajewski and C. Guhlke  “Analysis of improved Nernst–Planck–Poisson models of compressible isothermal electrolytes”  In _Zeitschrift für angewandte Mathematik und Physik_ 71.4, 2020, pp. 1–68
  * [DB87] L. Durlofsky and J.F. Brady  “Analysis of the Brinkman equation as a model for flow in porous media”  In _The Physics of fluids_ 30.11 American Institute of Physics, 1987, pp. 3329–3341
  * [ERS15] M. Eleuteri, E. Rocca and G. Schimperna  “On a non-isothermal diffuse interface model for two-phase flows of incompressible fluids”  In _Discrete & Continuous Dynamical Systems - A_ 35, 2015, pp. 2497–2522
  * [Fei04] E. Feireisl  “Dynamics of viscous compressible fluids”  Oxford Univ Press, Oxford, 2004
  * [Fei07] E. Feireisl  “Asymptotic analysis of the full Navier-Stokes-Fourier system: From compressible to incompressible fluid flows”  In _Russian Mathematical Surveys_ 62, 2007, pp. 511
  * [FN05] E. Feireisl and A. Novotný  “On a simple model of reacting compressible flows arising in astrophysics”  In _Proceedings of the Royal Society of Edinburgh: Section A Mathematics_ 135.6 Royal Society of Edinburgh Scotland Foundation, 2005, pp. 1169–1194
  * [FN05a] E. Feireisl and A. Novotný  “Weak Sequential Stability of the Set of Admissible Variational Solutions to the Navier-Stokes-Fourier System”  In _SIAM J. Math. Anal._ 37, 2005, pp. 619–650
  * [FN09] E. Feireisl and A. Novotný  “Singular Limits in Thermodynamics of Viscous Fluids”  Birkhäuser-Verlag, Basel, 2009
  * [FN12] E. Feireisl and A. Novotný  “Weak–Strong Uniqueness Property for the Full Navier–Stokes–Fourier System”  In _Arch. Rational Mech. Anal_ 204, 2012, pp. 683–706
  * [FNS11] E. Feireisl, A. Novotnỳ and Y. Sun  “Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids”  In _Indiana University Mathematics Journal_ , 2011, pp. 611–631
  * [FP07] E. Feireisl and H. Petzeltová  “On the Long-Time Behaviour of Solutions to the Navier–Stokes–Fourier System with a Time-Dependent Driving Force”  In _Journal of Dynamics and Differential Equations_ 19, 2007, pp. 685–707
  * [GBY] F. Gay-Balmaz and H. Yoshimura  “A free energy Lagrangian variational formulation of the Navier-Stokes-Fourier system”  In _International Journal of Geometric Methods in Modern Physics_ 16, pp. 1940006
  * [GKL17] M.-H. Giga, A. Kirshtein and C. Liu  “Variational modeling and complex fluids”  In _Handbook of mathematical analysis in mechanics of viscous fluids_ Springer International Publishing, Cham, 2017, pp. 1–41
  * [HLB96] P. Holmes, J. L. Lumley and G. Berkooz  “Turbulence, Coherent Structures, Dynamical Systems and Symmetry”, Cambridge Monographs on Mechanics  Cambridge University Press, Cambridge, 1996
  * [Hsi+20] C.-Y. Hsieh, T.-L. Lin, Chun C. Liu and P. Liu  “Global existence of the non-isothermal Poisson–Nernst–Planck–Fourier system”  In _Journal of Differential Equations_ 269.9, 2020, pp. 7287–7310
  * [HL+10] Y. Hyon and C. Liu  “Energetic variational approach in complex fluids: maximum dissipation principle”  In _Discrete & Continuous Dynamical Systems-A_ 26.4, 2010, pp. 1291
  * [LSU68] O.A. Ladyzhenskaya, V.A. Solonnikov and N.N. Uralceva  “Linear and Quasilinear Equations of Parabolic Type”  AMS, Providence, 1968
  * [LLT20] N.-A. Lai, C. Liu and A. Tarfulea  “Positivity of temperature for some non-isothermal fluid models”  In _arXiv preprint arXiv:2011.07192_ , 2020
  * [LLZ05] F.-H. Lin, C. Liu and P. Zhang  “On hydrodynamics of viscoelastic fluids”  In _Communications on Pure and Applied Mathematics_ 58.11, 2005, pp. 1437–1471
  * [LS21] C. Liu and J.-E. Sulzbach  “Well-Posedness for the Reaction-Diffusion Equation with Temperature in a critical Besov Space”  In _arXiv preprint arXiv:2101.10419_ , 2021
  * [LWL18] P. Liu, S. Wu and C. Liu  “Non-isothermal electrokinetics: energetic variational approach”  In _JournalCommunications in Mathematical Sciences_ 16, 2018, pp. 1451–1463
  * [McQ76] D.A. McQuarrie  “Statistical Mechanics”  Harper & Row, New York, 1976
  * [NPT13] T. Nishida, M. Padula and Y. Teramoto  “Heat Convection of Compressible Viscous Fluids: I”  In _J. Math. Fluid. Mech._ 15, 2013, pp. 525–536
  * [Pou07] L. Poul  “Existence of Weak Solutions to the Navier–Stokes–Fourier System on Lipschitz Domains”  In _Discrete and Continuous Dynamical Systems_ , 2007, pp. 834–843
  * [Rou13] T. Roubíček  “Nonlinear partial differential equations with applications”  Springer, Basel, 2013
  * [Sal01] S. Salinas  “Introduction to Statistical Physics”  Springer, New York, 2001
  * [Tak+99] S.-E. Takahasi, M. Tsukada, K. Tanahashi and T. Ogiwara  “An inverse type of Jensen’s Inequality”  In _Math. Japon._ 50, 1999, pp. 85–92
  * [Tar19] A. Tarfulea  “Improved a priori bounds for thermal fluid equations”  In _Trans. Amer. Math. Soc_ 371, 2019, pp. 2719–2737
  * [Tar75] L. Tartar  “Compensated compactnes and applications to partial differential equations”  Pitman, Boston, 1975
  * [Tom02] M. Tominaga  “Specht’s ratio in the Young inequality”  In _Sci. Math. Japon._ 55, 2002, pp. 583–588
  * [Zen99] Y. Zeng  “Gas Dynamics in Thermal Nonequilibrium and General Hyperbolic Systems with Relaxation”  In _Archive for Rational Mechanics and Analysis_ 155, 1999, pp. 225–279
  * [Zen10] Y. Zeng  “Gas Flows with Several Thermal Nonequilibrium Modes”  In _Arch. Rational Mech. Anal._ 196, 2010, pp. 191–225