Unconditional finite amplitude stability of a fluid in a mechanically isolated vessel with spatially non-uniform wall temperature

TL;DR

This paper proves that in a thermally isolated vessel with non-uniform wall temperature, the fluid's velocity and temperature perturbations diminish over time, ensuring stability of the steady state despite complex heat flux dynamics.

Contribution

It provides a thermodynamically grounded proof of the unconditional finite amplitude stability of the fluid's steady state in such a system.

Findings

Velocity field vanishes over time.

Temperature perturbations decay to steady state.

Stability holds regardless of initial conditions.

Abstract

A fluid occupying a mechanically isolated vessel with walls kept at spatially non-uniform temperature is in the long run expected to reach the spatially inhomogeneous steady state. Irrespective of the initial conditions the velocity field is expected to vanish, and the temperature field is expected to be fully determined by the steady heat equation. This simple observation is however difficult to prove using the corresponding governing equations. The main difficulties are the presence of the dissipative heating term in the evolution equation for temperature and the lack of control on the heat fluxes through the boundary. Using thermodynamically based arguments, it is shown that these difficulties in the proof can be overcome, and it is proved that the velocity and temperature perturbations to the steady state actually vanish as the time goes to infinity.

Equations

$${\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\hat{\theta }\left[\frac{1}{n}{\left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)}^n-\frac{1}{m}{\left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)}^m+\frac{n-m}{mn}\right]\mathrm{dv}\stackrel{t\to +\infty }{\to }0,$$

$$\rho c_{\mathrm{V},\mathrm{ref}}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }{\tilde{\theta }}^2\mathrm{dv}=-{\int }_{\Omega }{\kappa }_{\mathrm{ref}}\nabla \tilde{\theta }\bullet \nabla \tilde{\theta }\mathrm{dv}+{\int }_{\Omega }2\mu \tilde{\mathbb{D}}:\tilde{\mathbb{D}}\tilde{\theta }\mathrm{dv}+{\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\left(\tilde{\mathbf{v}}\bullet \nabla \hat{\theta }\right)\tilde{\theta }\mathrm{dv}.$$

$${\int }_{\tau =0}^{+\infty }\left({\int }_{\Omega }2\mu \tilde{\mathbb{D}}:\tilde{\mathbb{D}}\mathrm{dv}\right)\mathrm{d}\tau <+\infty ,$$

$$y:[0,+\infty )\mapsto {\mathbb{R}}^{+}$$

$${\int }_{\tau =0}^{+\infty }y(\tau )\mathrm{d}\tau \le C_1,$$

$$y(t)-y(s)\le {\int }_{\tau =s}^{t}f(y(\tau ))\mathrm{d}\tau +{\int }_{\tau =s}^{t}h(\tau )\mathrm{d}\tau$$

$$\underset{t\to +\infty }{\lim }y(t)=0.$$

$${\int }_{\tau =0}^{+\infty }y(\tau )\mathrm{d}\tau \le C_1,$$

$$\frac{\mathrm{d}y}{\mathrm{d}t}\le f(y)+h(t),$$

$$\underset{t\to +\infty }{\lim }y(t)=0.$$

$$g(t)=\{\begin{array}{ll}1,& t\in [N-\frac{1}{2N^2},N+\frac{1}{2N^2}],N\in \mathbb{N},\\ 0,& \text{otherwise},\end{array}$$

$${\int }_{t=0}^{+\infty }g(t)\mathrm{d}t=\sum _{N=1}^{+\infty }\frac{1}{N^2}<+\infty ,$$

$$\underset{t\to +\infty }{\lim }g(t)\ne 0,$$

$${\mathcal{V}}_{\mathrm{meq}}\left(\tilde{\mathbf{W}}\Vert \hat{\mathbf{W}}\right){=}_{\mathrm{def}}{\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\hat{\theta }\left[\frac{\tilde{\theta }}{\hat{\theta }}-\ln \left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)\right]\mathrm{dv}+{\int }_{\Omega }\frac{1}{2}\rho {\left|\tilde{\mathbf{v}}\right|}^2\mathrm{dv},$$

$$\begin{gathered} \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{V}}_{\mathrm{meq}}\left(\tilde{\mathbf{W}}\Vert \hat{\mathbf{W}}\right)=-{\int }_{\Omega }{\kappa }_{\mathrm{ref}}\hat{\theta }\nabla \ln \left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)\bullet \nabla \ln \left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)\mathrm{dv}-{\int }_{\Omega }\frac{2\mu \tilde{\mathbb{D}}:\tilde{\mathbb{D}}}{1+\frac{\tilde{\theta }}{\hat{\theta }}}\mathrm{dv} \\ +{\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\left(\nabla \hat{\theta }\bullet \tilde{\mathbf{v}}\right)\ln \left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)\mathrm{dv}, \end{gathered}$$

$${\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,m}\left(\tilde{\mathbf{W}}\Vert \hat{\mathbf{W}}\right){=}_{\mathrm{def}}{\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\hat{\theta }\left[\frac{\tilde{\theta }}{\hat{\theta }}-\frac{1}{m}\left({\left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)}^m-1\right)\right]\mathrm{dv}+{\int }_{\Omega }\frac{1}{2}\rho {\left|\tilde{\mathbf{v}}\right|}^2\mathrm{dv},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,m}-{\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,n}\right)\le -K^{m,n}\left({\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,m}-{\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,n}\right)+{\mathcal{H}}^{m,n},$$

$${\mathcal{Y}}^{m,n}{=}_{\mathrm{def}}{\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,m}-{\mathcal{V}}_{\mathrm{meq}}^{\vartheta ,n}={\int }_{\Omega }\rho c_{\mathrm{V},\mathrm{ref}}\hat{\theta }\left[\frac{1}{n}{\left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)}^n-\frac{1}{m}{\left(1+\frac{\tilde{\theta }}{\hat{\theta }}\right)}^m+\frac{n-m}{mn}\right]\mathrm{dv},$$

$$\frac{\mathrm{d}{\mathcal{Y}}^{m,n}}{\mathrm{d}t}\le K^{m,n}{\mathcal{Y}}^{m,n}+{\mathcal{H}}^{m,n}.$$

$$\mathrm{div}\mathbf{v}$$

$$\rho \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\rho c_{\mathrm{V},\mathrm{ref}}\frac{\mathrm{d}\theta }{\mathrm{d}t}$$

$$\psi {=}_{\mathrm{def}}-c_{\mathrm{V},\mathrm{ref}}\theta \left(\ln \left(\frac{\theta }{{\theta }_{\mathrm{ref}}}\right)-1\right).$$

$$\eta =c_{\mathrm{V},\mathrm{ref}}\ln \left(\frac{\theta }{{\theta }_{\mathrm{ref}}}\right),e=c_{\mathrm{V},\mathrm{ref}}\theta .$$

$${\mathbf{j}}_q{=}_{\mathrm{def}}-{\kappa }_{\mathrm{ref}}\nabla \theta ,$$

$$\rho \frac{\mathrm{d}\eta }{\mathrm{d}t}=\frac{{\zeta }_{\mathrm{th}}}{\theta }+\frac{{\zeta }_{\mathrm{mech}}}{\theta }-\mathrm{div}{\mathbf{j}}_{\eta },$$

$${\zeta }_{\mathrm{mech}}{=}_{\mathrm{def}}2\mu \mathbb{D}:\mathbb{D},{\zeta }_{\mathrm{th}}{=}_{\mathrm{def}}{\kappa }_{\mathrm{ref}}\frac{{\left|\nabla \theta \right|}^2}{\theta },{\mathbf{j}}_{\eta }{=}_{\mathrm{def}}\frac{{\mathbf{j}}_q}{\theta }.$$

$${\mathbf{v}|}_{\partial \Omega }=0,{\theta |}_{\partial \Omega }={\theta }_{\mathrm{bdr}},$$

$${\int }_{\partial \Omega }\mathbb{T}\mathbf{v}\bullet \mathbf{n}\mathrm{ds}=0.$$

$${\int }_{\partial \Omega }{\mathbf{j}}_q\bullet \mathbf{n}\mathrm{ds}\ne 0.$$

$$0=-\nabla \hat{p},0=\mathrm{div}\left({\kappa }_{\mathrm{ref}}\nabla \hat{\theta }\right),$$

$${\hat{\mathbf{v}}|}_{\partial \Omega }=0,{\hat{\theta }|}_{\partial \Omega }={\theta }_{\mathrm{bdr}},$$

$$\mathrm{div}\tilde{\mathbf{v}}$$

Full text

1 Introduction

The everyday experience is that if a fluid is put into a vessel, and if it is not allowed to substantially interact with the outside environment, then it eventually comes to the rest state. Moreover, the rest state is attained irrespective of the initial state of the fluid. The question is whether this behaviour can be deduced using the corresponding governing equations such as the Navier–Stokes–Fourier equations for an incompressible fluid.

In the case of a vessel that is completely isolated from the outside environment, see Figure 1a, the answer to this question is straightforward and positive. One can exploit basic ideas from continuum thermodynamics, see Coleman (1970), Gurtin (1973, 1975) and early considerations by Duhem (1911), and prove that both the velocity and temperature field eventually reach the spatially homogeneous equilibrium rest state. The same holds also for simple open systems such as a vessel immersed in a thermal bath, that is for the mechanically isolated vessel with walls kept at constant spatially uniform temperature, see Figure 1b. In a slightly more general case the answer to the question is much more complicated.

If we consider a mechanically isolated vessel with walls kept at a given spatially non-uniform temperature, see Figure 1c, that is if we change the temperature boundary condition, we are not anymore dealing with a thermodynamically isolated system or a system immersed in a thermal bath, and the concepts applicable in the stability analysis of these systems are of no use. Indeed, the change of temperature boundary condition, which on the formal level means the change from the zero Neumann boundary condition/constant Dirichlet boundary condition to the spatially inhomogeneous Dirichlet boundary condition, qualitatively changes both the steady state and the mode of interaction with the outside environment.

However, even in this setting we again expect the fluid to reach a steady state. The expected steady state is the state with the zero velocity field, while the steady temperature field is given by the steady heat equation and it is spatially inhomogeneous. This means that from the thermodynamical perspective we are dealing with a non-equilibrium (entropy producing) steady state.

Despite the increased complexity of the steady state, it is still straightforward to prove that velocity field decays to zero, see for example Serrin (1959). This reflects the fact that the kinetic energy is spontaneously converted to the thermal energy via the dissipative heating, hence if there is no mechanical energy supply from outside, then we should observe the decay of the kinetic energy.

Concerning the temperature field, the situation is different. First, if the fluid is allowed to exchange heat with its surrounding, then it is a priori not clear, whether one will during the evolution of the system see heat influx or efflux. Second, the thermal energy is, in our simple case of incompressible Navier–Stokes–Fourier fluid, not converted to another type of energy. It is merely rearranged in space. This is in striking difference with the behaviour of the kinetic energy that is being dissipated and that naturally “disappears” from the system.

Consequently, the key question whether the governing equations actually predict the expected long time behaviour, becomes in the case of spatially inhomogeneous wall temperature rather difficult to answer. In particular, any characterisation of the long time behaviour of temperature perturbations seems to be difficult to obtain in such a setting. In fact, to the best of our knowledge, no proof regarding the decay of temperature perturbations is available in the literature so far.

1.1 Main result

We investigate the key question for the fluid described by the standard incompressible Navier–Stokes–Fourier model, and we provide a positive answer to the key question. (We consider the motion in the absence of external forces and we assume that the fluid density is constant. No Oberbeck–Boussinesq approximation is considered.) Although we consider stability of spatially inhomogeneous non-equilibrium steady state in a thermodynamically open system, we are able during the proof to exploit thermodynamical ideas, though these ideas go beyond the concepts introduced by Coleman (1970) in the stability analysis of spatially homogeneous steady states in isolated systems or systems immersed in a thermal bath.

In particular, assuming that the solution to the governing equations is the classical one, we prove that

[TABLE]

where $m,n\in(0,1)$, $n>m>\frac{n}{2}$, see Theorem 5.1 and Corollary 1. Here $\widehat{\theta}$ denotes the temperature field in the spatially inhomogeneous non-equilibrium steady state (solution to the steady heat equation), and $\widetilde{\theta}$ denotes the perturbation with respect to $\widehat{\theta}$, and the symbols $\rho$ and $c_{\mathrm{V},\mathrm{ref}}$ denote the density and the specific heat capacity at constant volume. The functional in (1.1) is non-negative and it vanishes if and only if the perturbation $\widetilde{\theta}$ vanishes everywhere in the domain $\Omega$ of interest, which means that (1.1) is the desired stability result. We note that (1.1) holds irrespective of the initial state of the system, and with minimal assumptions concerning the behaviour of the dissipative heating term in the evolution equation for the temperature.

1.2 Main issues

Before we proceed with the proof of (1.1), we briefly comment on the main issues in the analysis of long time behaviour of the temperature perturbations. (See Section 2 for the notation and precise formulation of the long time behaviour problem.)

The first issue is the presence of the dissipative heating term $2\mu\mathbb{D}:\mathbb{D}$ in the evolution equation for the temperature. If we were following the so-called energy method, see Joseph (1976a, b) or Straughan (2004), we would have investigated the behaviour of the temperature perturbations $\widetilde{\theta}$ by the means of functional $\int_{\Omega}\widetilde{\theta}^{2}\,\mathrm{d}\mathrm{v}$. The evolution equation for this functional is easy to obtain from the evolution equation for the perturbation $\widetilde{\theta}$. It suffices to take the product of the evolution equation for the perturbation with the perturbation itself and integrate by parts, and one arrives at the equation

[TABLE]

The first term on the right-hand side is a favourable one. It is negative and it pushes the temperature perturbation to zero. On the other hand the second term that comes from the dissipative heating does not have a sign, and it is hard to control/estimate unless one is willing to search for further information111Further information on the product $\widetilde{\mathbb{D}}:\widetilde{\mathbb{D}}$ can be formally obtained by the multiplication of the evolution equation (2.12b) by the Laplacian of the velocity perturbation, see Kagei et al. (2000) for a similar manipulation, and also the discussion of slightly compressible convection in Straughan (2004) and Richardson (1993). This approach would however work only for constant viscosity, and it would allow one to prove stability only for a restricted set of initial perturbations (small perturbations). We however aim at unconditional result – the size of initial perturbation must not be limited. This is what we expect intuitively from our physical system. concerning the behaviour of $2\mu\widetilde{\mathbb{D}}:\widetilde{\mathbb{D}}$. In particular, even if we know that the velocity perturbation $\widetilde{\bm{v}}$ decays, which is true, the fact that the velocity perturbation $\widetilde{\bm{v}}$ is small does not a priori imply that the velocity gradient $\widetilde{\mathbb{D}}$ is also small. (A rapidly oscillating function can be small, while its gradient can be arbitrarily large.) Further, one can not expect to obtain an explicit pointwise upper bound on the temperature perturbation $\widetilde{\theta}$. The source term $2\mu\widetilde{\mathbb{D}}:\widetilde{\mathbb{D}}$ in the heat equation may concentrate at a spatial point, and push the local value of the temperature above any a priori given threshold.

In fact, the only piece of information one is allowed to use in the analysis of the evolution of temperature perturbations is

[TABLE]

which means that the heat released by the dissipation of the kinetic energy is finite. However, it is not known when and where the heat is released. (We note that if the integral $\int_{\tau=0}^{\infty}g(\tau)\,\mathrm{d}\tau$ of a positive function $g$ is finite, it is not necessarily true that $g\to 0+$ as $\tau\to+\infty$. See also the discussion following Lemma 2 below.) In other words, the source term in the temperature evolution equation (2.12c) could be triggered at a priori unknown time instants and spatial locations.

The second issue is that the heat flux through the boundary is in the given setting beyond our control. We do not know a priori whether at the given time instant and spatial location the heat flows into the domain of interest or out of the domain of interest.

For these reasons the dissipative heating is conveniently neglected in most mathematical studies on the stability of fluid motion, see Joseph (1976a, b) or Straughan (2004).

1.3 Outline of the solution to the stability problem

If the only piece of information we are willing to use is (1.3), we can hardly expect any result concerning the decay rate of the perturbations. Consequently, if we ignore the spatial dependence, we need to work with very weak qualitative results concerning the asymptotic behaviour of time dependent functions. In our analysis we use the following lemma.

Lemma 1 (Decay of integrable functions)

[TABLE]

The lemma is taken from Zheng (2001, Lemma 1.3), and it originated in the work of Krejčí and Sprekels (1998, Theorem 9), see also Zheng (2004) for further discussion. Various generalisations of the lemma can be found in Ramm and Hoang (2012), and we also point out that the question on decay of integrable functions is also a subject of Barbălat lemma known in the optimal control theory, see Farkas and Wegner (2016). Lemma 1 can be also rewritten in the differential form, but in this case one has to assume the existence of the derivative. This modification of Lemma 1 reads as follows, see Zheng (2001, Lemma 1.3).

Lemma 2 (Decay of integrable functions – differential version)

Let $y:[0,+\infty)\mapsto{\mathbb{R}}^{+}$ be a continuous non-negative function such that

[TABLE]

where $C_{1}$ is a constant, and let us assume that $\frac{\mathrm{d}{y}}{\mathrm{d}{t}}$ is a locally integrable function on ${\mathbb{R}}^{+}$. Suppose that $y$ satisfies for all $t\in(0,+\infty)$ the condition

[TABLE]

where functions $f$ and $h$ have the same properties as in Lemma 1. Then

[TABLE]

In our case the existence of the time derivative of the corresponding quantity $y$ will be granted since we assume that we work with the classical solution of the governing equations. Consequently we can freely switch—whenever convenient—between the integral formulation used in Lemma 1 and the differential formulation used in Lemma 2.

We recall the well-known fact that the existence of the integral (1.4a) does not on its own imply (1.4c), the classical counterexample is the function

[TABLE]

see Figure 2 for a sketch of such a function. (If necessary, using the convolution with the standard smoothing kernel, the function can be made infinitely smooth.) For this function we have

[TABLE]

but $\lim_{N\to+\infty}g(N)=1$. Consequently we see that

[TABLE]

and the limit in fact does not exist at all. The condition (1.4b) or (1.6) in Lemma 1 and Lemma 2 respectively effectively prohibits the existence of the sufficiently fast thinning peaks as that shown in Figure 2. (The maximum slope of the graph that is the derivative is limited.) This provides an intuitive argument for the validity of the Lemma 1 and Lemma 2. The rigorous proof of the lemmas is given for example in Zheng (2004).

Once we decide that Lemma 1/Lemma 2 or its generalisation is the right tool for stability analysis, the only remaining task is to find the quantity $y(t)$ that satisfies the assumptions of the lemma, and that vanishes if and only if the temperature perturbation vanishes.

First, we try to solve this problem using a Lyapunov type functional constructed via the method proposed by Bulíček et al. (2019). The functional constructed by the method reads

[TABLE]

[TABLE]

and its time derivative reads

[TABLE]

see Section 3. We see that the first two terms on the right-hand side are negative, while last term can be positive or negative depending on the actual evolution of the perturbation. Consequently, it is not clear whether the time derivative of the proposed functional is negative. Despite these complications, the functional $\mathcal{V}_{\mathrm{meq}}$ is a promising one for the study of the stability problem, and it is almost the right one. We also note, that if we were dealing with a spatially homogeneous steady state $\widehat{\theta}$, which is the steady state in a thermodynamically isolated vessel or in a vessel immersed in a thermal bath, then $\nabla\widehat{\theta}=0$, and the last term would vanish, and the stability problem would be solved immediately.

Second, in Section 4 we revisit the construction of the prospective Lyapunov type functional $\mathcal{V}_{\mathrm{meq}}$. We exploit the idea that the stability of the rest state should be an absolute fact independent of the choice of the description of the system. In particular, it should be independent of the choice of the temperature scale. In (1.11a) we however use the absolute temperature which is just a convenient measure of hotness. But it is not necessary to use the absolute temperature scale, other temperature scales can be used as well.

For the alternative temperature scale, we use the scale $\frac{\vartheta}{\vartheta_{\mathrm{ref}}}=_{\mathrm{def}}\left(\frac{\theta}{\theta_{\mathrm{ref}}}\right)^{1-m}$, where $\vartheta$ is the alternative temperature scale and $\vartheta_{\mathrm{ref}}$ and $\theta_{\mathrm{ref}}$ are some fixed reference temperatures and $m\in(0,1)$. It is straightforward to check that the alternative temperature scale preserves the ordering according to the hotness, which means that the heat still flows in the direction of the temperature difference and so forth, but quantitative characterisation of the flow is different.

With respect to this alternative choice of temperature scale, the fluid of interest virtually behaves like a fluid with a temperature dependent thermal conductivity and temperature dependent specific heat capacity at constant volume. For this virtual fluid we can again use the method proposed by Bulíček et al. (2019) and construct a candidate for a Lyapunov type functional $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$. If we rewrite the functional $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ using the original temperature scale $\theta$ it reads

[TABLE]

and its time derivative is again given by an explicit formula. Using this manipulation we obtain a whole class of functionals that are again almost appropriate for the stability analysis.

Third, we conclude our analysis by exploiting Lemma 1. Using the family of functionals (1.12) we derive, see Section 5, the inequality

[TABLE]

where $K^{m,n}$ is a positive constant, $\mathcal{H}^{m,n}$ is a non-negative function integrable with respect to time, and $m,n\in(0,1)$, $n>m>\frac{n}{2}$. In (1.13) it is straightforward to show that the quantity

[TABLE]

is a non-negative quantity that vanishes if and only if the temperature perturbation $\widetilde{\theta}$ vanishes. Using (1.13) we can conclude that $\int_{\tau=0}^{+\infty}\mathcal{Y}^{m,n}(\tau)\,\mathrm{d}\tau\leq C_{1}^{m,n}$, where $C_{1}^{m,n}$ is a positive constant. Estimating the right-hand side of (1.13) by its absolute value we also see that $\mathcal{Y}^{m,n}$ satisfies the inequality

[TABLE]

Inequality (1.13) and Lemma 1 then allow us to conclude that $\mathcal{Y}^{m,n}\xrightarrow{t\to+\infty}0$, which is in virtue of the definition of $\mathcal{Y}^{m,n}$ tantamount to (1.1). This is the desired stability result, see Theorem 5.1 and also Corollary 1 which gives the characterisation of the stability in terms of a Lebesgue space norm.

2 Preliminaries

Let us now introduce the notation and let us give a precise formulation of the problem. The notation regarding the governing equations and their thermodynamic basis is the standard one, we mainly follow Bulíček et al. (2019) and Málek and Průša (2018).

2.1 Governing equations and basic facts from thermodynamics

The fluid of interest is assumed to be the standard incompressible Navier–Stokes–Fourier fluid. In particular, we do not use any form of Oberbeck–Boussinesq type approximation, we just consider a fluid with a constant density. The governing equations in the absence of external forces read

[TABLE]

where $\theta$ denotes the absolute temperature, $\bm{v}$ denotes the velocity and $p$ stands for the pressure. Symbols $\rho$, $\mu$, $c_{\mathrm{V},\mathrm{ref}}$ and $\kappa_{\mathrm{ref}}$ denote the density, the viscosity, the specific heat at constant volume and the thermal conductivity, and all these material parameters are assumed to be positive constants. Furthermore, the symbol $\mathbb{D}=_{\mathrm{def}}\frac{1}{2}\left(\nabla\bm{v}+\nabla\bm{v}^{\top}\right)$ denotes the symmetric part of the velocity gradient, and $\frac{\mathrm{d}{}}{\mathrm{d}{t}}=\frac{\partial{}}{\partial{t}}+\bm{v}\bullet\nabla$ denotes the material time derivative.

Since we need some basic thermodynamics in the construction of Lyapunov type functionals via the method proposed by Bulíček et al. (2019), we recall some basic facts concerning the thermodynamic basis of the governing equations. The specific Helmholtz free energy for the incompressible Navier–Stokes–Fourier fluid is given by the formula

[TABLE]

Using the standard thermodynamical identities, see for example Callen (1985), the specific Helmholtz free energy $\psi$ in the form (2.2) leads to the following formulae for the specific entropy $\eta$ and the specific internal energy $e$,

[TABLE]

The heat flux vector $\bm{j}_{q}$ is given by the Fourier law,

[TABLE]

and the Cauchy stress tensor is given as $\mathbb{T}=_{\mathrm{def}}-p\mathbb{I}+2\mu\mathbb{D}$. The general evolution equation for the entropy reads

[TABLE]

where $\frac{\zeta_{\mathrm{th}}}{\theta}$ and $\frac{\zeta_{\mathrm{mech}}}{\theta}$ denote the entropy production terms related to the thermal and mechanical quantities, and $\bm{j}_{\eta}$ denotes the entropy flux. In our setting we have

[TABLE]

2.2 Boundary conditions

The boundary conditions are the standard boundary conditions for thermal convection. We consider a fluid in a rigid vessel $\Omega$ such that

[TABLE]

where $\theta_{\mathrm{bdr}}$ is a given function of position. In particular, we are interested in the case when $\theta_{\mathrm{bdr}}$ is a nontrivial function of position. Since the temperature evolution equation (2.1c) is a parabolic equation with a convective term and a positive source, we know that the temperature field $\theta=\widehat{\theta}+\widetilde{\theta}$ is under given boundary conditions (2.7) bounded from below uniformly in space and time, see for example Friedman (1964), Ladyzhenskaya et al. (1968) or Lieberman (1996). In particular, if the boundary temperature is positive $\theta_{\mathrm{bdr}}$, then $\theta$ is also always positive.

The no-slip boundary condition $\left.\bm{v}\right|_{\partial\Omega}=\bm{0}$ guarantees that the vessel is mechanically isolated. This means that the mechanical energy exchange between the vessel and the outside environment, which is given by the surface integral $\int_{\partial\Omega}\mathbb{T}\bm{v}\bullet\bm{n}\mathrm{d}\mathrm{s}$, is equal to zero,

[TABLE]

On the other hand the heat exchange between the vessel and the outside environment is given by the surface integral $-\int_{\partial\Omega}\bm{j}_{q}\bullet\bm{n}\mathrm{d}\mathrm{s}=\int_{\partial\Omega}\kappa_{\mathrm{ref}}\nabla\theta\bullet\bm{n}\mathrm{d}\mathrm{s}$. The temperature boundary condition $\left.\theta\right|_{\partial\Omega}=\theta_{\mathrm{bdr}}$ does not in general imply that the temperature gradient vanishes on the boundary as well, hence we in general expect that

[TABLE]

In this sense we characterise our system of interest as a mechanically isolated and thermally open system, and our choice of boundary conditions indeed corresponds to the system depicted in Figure 1c.

2.3 Spatially inhomogeneous non-equilibrium steady state

The non-equilibrium steady state $\widehat{\bm{W}}=_{\mathrm{def}}[\widehat{\bm{v}},\widehat{\theta},\widehat{p}]$ solves the system

[TABLE]

with boundary conditions

[TABLE]

which means that the fluid is at rest, $\widehat{\bm{v}}=\bm{0}$, and that the temperature field is given by the steady heat equation (2.10) with Dirichlet boundary condition (2.11). If $\theta_{\mathrm{bdr}}$ is a nontrivial function of position, then $\widehat{\theta}$ is a spatially inhomogeneous bounded function.

2.4 Evolution equations for a perturbation to the steady state and the stability problem

The velocity/temperature/pressure fields $\bm{W}=[\bm{v},\theta,p]$ are split into the steady part $\widehat{\bm{W}}=[\widehat{\bm{v}},\widehat{\theta},\widehat{p}]$ and the perturbation $\widetilde{\bm{W}}=_{\mathrm{def}}[\widetilde{\bm{v}},\widetilde{\theta},\widetilde{p}]$, that is $[\bm{v},\theta,\widehat{p}]=[\widehat{\bm{v}},\widehat{\theta},\widehat{p}]+[\widetilde{\bm{v}},\widetilde{\theta},\widetilde{p}]$. The governing equations for the triple $\widetilde{\bm{W}}=[\widetilde{\bm{v}},\widetilde{\theta},\widetilde{p}]$ are

[TABLE]

These equations are straightforward to obtain from the fact that $\bm{W}$ as well as $\widehat{\bm{W}}$ solve the governing equations (2.1). The initial conditions for the full evolution equations read $\left.\bm{v}\right|_{t=0}=\bm{v}_{\mathrm{init}}$, $\left.\theta\right|_{t=0}=\theta_{\mathrm{init}}$, which implies that the initial conditions for the perturbation are

[TABLE]

The boundary conditions for the perturbation are derived from (2.7), which upon using the fact that the steady state fulfills the boundary condition (2.11), leads to

[TABLE]

The task is to show that the perturbation vanishes as the time goes to infinity, that is

[TABLE]

The exact meaning of the convergence can be specified at will in the sense that it can be specifically tailored for the problem. In our setting and for the temperature perturbation the meaning of (2.15) is (1.1).

3 A candidate for a Lyapunov type functional

Bulíček et al. (2019) have proposed a thermodynamically based method for a systematic construction of a candidate for a Lyapunov type functional. We exploit this method, and we show that although the candidate for a Lyapunov type functional is not a genuine Lyapunov functional—its time derivative can not be shown to be negative—it is still useful in the stability analysis. The key findings of this section are the following.

Lemma 3 (Lyapunov type functional and its time derivative)

Let us consider the velocity/temperature perturbation $\widetilde{\theta}$, $\widetilde{\bm{v}}$ governed by equations (2.12). Let us define the functional

[TABLE]

where we use the notation introduced in Section 2. The functional remains non-negative for all possible temperature perturbations $\widetilde{\theta}\in(-\widehat{\theta},+\infty)$ and velocity perturbations $\widetilde{\bm{v}}$, and it vanishes if and only if the perturbations $\widetilde{\theta}$ and $\widetilde{\bm{v}}$ vanish everywhere in the domain of interest. Furthermore, the time derivative of the functional is given by the formula

[TABLE]

Lemma 4 (On boundedness of some integrals)

Let us consider the velocity/temperature perturbation $\widetilde{\theta}$, $\widetilde{\bm{v}}$ governed by equations (2.12). Then the dissipative heating term possesses a finite time-space integral,

[TABLE]

Note that Lemma 4 shows that all the terms on the right-hand side of (3.2) have bounded integrals if integrated with respect to time from zero to infinity. (Regarding the last term on the right-hand side see identity (3.21) and (3.3d).)

3.1 Construction of a candidate for a Lyapunov type functional

A candidate for a Lyapunov type functional is, according to Bulíček et al. (2019), defined as

[TABLE]

where

[TABLE]

and the functionals $S_{\widehat{\theta}}\left(\bm{W}\right)$ and $E_{\mathrm{tot}}\left(\bm{W}\right)$ are defined as

[TABLE]

The symbols $\left.\mathrm{D}S_{\widehat{\theta}}\left(\bm{W}\right)\right|_{\bm{W}=\widehat{\bm{W}}}\left[\widetilde{\bm{W}}\right]$ and $\left.\mathrm{D}E_{\mathrm{tot}}\left(\bm{W}\right)\right|_{\bm{W}=\widehat{\bm{W}}}\left[\widetilde{\bm{W}}\right]$ denote the Gâteaux derivative of the given functionals at point $\widehat{\bm{W}}$ in the direction $\widetilde{\bm{W}}$. In our case, the specific entropy $\eta$ and specific internal energy $e$ are given by formulae (2.3), and the Gâteaux derivatives of the functionals $S_{\widehat{\theta}}\left(\bm{W}\right)$ and $E_{\mathrm{tot}}\left(\bm{W}\right)$ read

[TABLE]

(Recall that we are investigating the stability of the non-equilibrium steady state $\widehat{\bm{W}}$ that does not evolve in time and where $\widehat{\bm{v}}=\bm{0}$.) Consequently, it is straightforward to see that the particular formulae for the functionals ${\mathcal{S}}_{\widehat{\theta}}(\left.\widetilde{\bm{W}}\right\|\widehat{\bm{W}})$ and ${\mathcal{E}}(\left.\widetilde{\bm{W}}\right\|\widehat{\bm{W}})$ read

[TABLE]

Using (3.8) in (3.4) we get the formula (3.1) for the candidate for a Lyapunov type functional. It is straightforward to check that the proposed functional vanishes if and only if the perturbation vanishes, and that the functional is non-negative.

3.2 Time derivative of our candidate for a Lyapunov type functional

It remains to investigate the time derivative of our candidate for a Lyapunov type functional $\mathcal{V}_{\mathrm{meq}}$. If we were dealing with a genuine Lyapunov functional we would need to show that it is negative. Using the explicit formulae for the entropy and internal energy (2.3) it follows that

[TABLE]

where we have introduced the relative entropy

[TABLE]

Note that in our case the explicit formula (2.3) for the entropy implies that

[TABLE]

Consequently, we can work with a very simple formula

[TABLE]

In order to get an explicit formula for the time derivative of the relative entropy $\eta_{\mathrm{diff}}$ in (3.12) we have to use the evolution equation for the entropy. The evolution equation for the entropy is (2.5), which implies that222Recall that the material time derivative in (3.13b) is taken with respect to the perturbed velocity field, that is $\frac{\mathrm{d}{\eta}}{\mathrm{d}{t}}(\widehat{\bm{W}}+\widetilde{\bm{W}})=\frac{\partial{\eta}}{\partial{t}}(\widehat{\bm{W}}+\widetilde{\bm{W}})+\left(\widehat{\bm{v}}+\widetilde{\bm{v}}\right)\bullet\nabla\eta(\widehat{\bm{W}}+\widetilde{\bm{W}})$, while the material time derivative in (3.13a) is taken with respect to the reference steady velocity field, that is $\frac{\mathrm{d}{\eta}}{\mathrm{d}{t}}(\widehat{\bm{W}})=\frac{\partial{\eta}}{\partial{t}}(\widehat{\bm{W}})+\widehat{\bm{v}}\bullet\nabla\eta(\widehat{\bm{W}})$.

[TABLE]

The relation between the temperature and the entropy is, for the model we are interested in, given by the formula (2.3). Furthermore, the formula for the thermal part of the entropy production $\zeta_{\mathrm{th}}$ (2.6), while the heat flux $\bm{j}_{q}$ is given by the Fourier law, which implies that

[TABLE]

Observations (3.14) and evolution equations (3.13) yield the evolution equation for the relative entropy $\eta_{\mathrm{diff}}$ in the form

[TABLE]

If we take into account that in our case we have $\widehat{\bm{v}}=\bm{0}$, and that the mechanical part of the entropy production $\zeta_{\mathrm{mech}}$ is given by (2.6), then (3.15) reduces to

[TABLE]

On the other hand, the evolution equations for the net total energy $E_{\mathrm{tot}}$ read

[TABLE]

Note that since we deal with a steady non-equilibrium state $\widehat{\bm{W}}$, then the left-hand side of (3.17b) is in fact zero. Subtracting (3.17a) and (3.17b) yields

[TABLE]

where we have used the specific constitutive relations for the incompressible Navier–Stokes–Fourier fluid, and the fact that $\widehat{\bm{v}}=\bm{0}$ and that $\widetilde{\bm{v}}$ vanishes on the boundary. Now we are in the position to substitute into (3.12). Using (3.18) and (3.16) in (3.12) yields

[TABLE]

Now we can utilise the lesson learned in the analysis of the simpler heat conduction problem in a rigid body, see Bulíček et al. (2019) for details, and we can rewrite (3.19) as

[TABLE]

where we have also used that fact that $\int_{\Omega}\rho\widehat{\theta}\widetilde{\bm{v}}\bullet\nabla\widehat{\eta}\,\mathrm{d}\mathrm{v}=0$. (This again follows from Stokes theorem and the fact that $\widehat{\theta}\nabla\widehat{\eta}$ can be rewritten as $\nabla F(\widehat{\theta})$.) This manipulation finishes the proof of Lemma 3. An alternative derivation of the formula for the time derivative of the proposed functional is given in Appendix A, where we in fact derive the pointwise evolution equation for the thermal part of the integrand in the functional $\mathcal{V}_{\mathrm{meq}}$.

Note that the last term in (3.20) can be also rewritten in a different form. Using the Stokes theorem and the fact that the perturbation velocity $\widetilde{\bm{v}}$ as well as $\eta_{\mathrm{diff}}$ vanish on the boundary, and that the perturbation velocity field has zero divergence, we see that

[TABLE]

which among other things shows that the term vanishes for spatially homogeneous field $\widehat{\theta}$.

All the terms in (3.20) are at least “quadratic” in the perturbation as expected. The right-hand side of (3.20) does not contain boundary terms, although we are working with a thermodynamically open system. This follows from the design of the proposed functional. We also see that the first two terms on the right-hand side of (3.20) have a sign, while the last term is indefinite. The presence of the last term prohibits one from showing that the time derivative of the functional $\mathcal{V}_{\mathrm{meq}}$ is, for a non-constant $\widehat{\theta}$, a non-positive quantity. Consequently, $\mathcal{V}_{\mathrm{meq}}$ can not directly serve as a Lyapunov type functional.

3.3 Estimate on the time derivative of the candidate for a Lyapunov type functional

The evolution equation for the perturbation velocity is (2.12b), which upon testing by the perturbation velocity $\widetilde{\bm{v}}$ yields

[TABLE]

(Recall that $\widetilde{\bm{v}}$ satisfies zero Dirichlet boundary conditions.) Integrating (3.22) with respect to time, we get

[TABLE]

Since $\frac{\rho}{2}\left\|\widetilde{\bm{v}}\right\|_{L^{2}\left(\Omega\right)}^{2}$ is a positive quantity, we see that (3.23) implies that the dissipated energy is finite, $\int_{\tau=0}^{t}\left(\int_{\Omega}2\mu\widetilde{\mathbb{D}}:\widetilde{\mathbb{D}}\,\mathrm{d}\mathrm{v}\right)\,\mathrm{d}\tau<+\infty.$ This is an important, yet trivial finding, and we will formulate it, for the sake of further reference, as Lemma 4, equation (3.3a).

Moreover, using the standard manipulation based on Poincaré and Korn (in)equality333Auxiliary tools from the theory of function spaces are for convenience summarised in Appendix C., we obtain the following estimate on the time derivative of the norm $\left\|\widetilde{\bm{v}}\right\|_{L^{2}\left(\Omega\right)}^{2}$

[TABLE]

which yields

[TABLE]

which is the standard result, see for example Serrin (1959). Note that if we were dealing for example with a temperature dependent viscosity, it would be sufficient to have the viscosity that is bounded from below, and we would still get an exponential decay of the net kinetic energy.

Using (3.25) and the $\varepsilon$–Young inequality $ab\leq\varepsilon a^{2}+\frac{b^{2}}{4\varepsilon}$ then gives us

[TABLE]

Furthermore, using again Poincaré inequality we get

[TABLE]

and utilising (3.27) and (3.26) in the formula (3.20) for the time derivative of the proposed Lyapunov type functional we arrive at the inequality

[TABLE]

can be made positive by a suitable choice of the parameter $\varepsilon$. (Note that (3.28b) is always positive if the steady state temperature field $\widehat{\theta}$ is spatially homogeneous.) Therefore, the structure of the estimate on the time derivative is a favourable one. It contains non-positive terms and a positive term that decays in time to zero.

3.4 Boundedness of some integrals

Inequality (3.28a) can be gainfully exploited in the proof of the boundedness of some integrals, see Lemma 4. The last term on the right-hand side of (3.28a) is integrable with respect to time, which follows from the direct integration of the exponential. The integration of (3.28a) with respect to time yields

[TABLE]

and since $\mathcal{V}_{\mathrm{meq}}$ is a non-negative quantity, we see that the first two integrals on the right-hand side must be bounded. In particular, we get

[TABLE]

Now we can use (3.30) in (3.26), and we can conclude that

[TABLE]

which in virtue of (3.21) also means that $\int_{\tau=0}^{+\infty}\left|\int_{\Omega}\rho\widehat{\theta}\widetilde{\bm{v}}\bullet\nabla\eta_{\mathrm{diff}}\,\mathrm{d}\mathrm{v}\right|\,\mathrm{d}\tau<+\infty$. Finally, once we have obtained (3.31) we can go back to (3.20), and use the same argument as before, and show (3.3b) and (3.3c). Indeed, integrating (3.20) with respect to time we get

[TABLE]

and since $\mathcal{V}_{\mathrm{meq}}$ is a non-negative quantity, and the last term on the right-hand side of (3.32) is bounded, we see that the first two integrals on the right-hand side must be bounded. The manipulations described above give the proof of Lemma 4.

3.5 Remarks

Although we have not succeeded with the concept of Lyapunov functional we are close to the desired stability result. The time derivative of the proposed functional is negative up to an exponentially decaying positive term, see (3.28a). Moreover, the formula (3.28a) which reads

[TABLE]

is also almost the right one from the perspective of Lemma 1. Indeed, if $x\in(-1,x_{\mathrm{crit}})$ then

[TABLE]

where $x_{\mathrm{crit}}$ is a positive number. (See Figure 3a for a sketch of the functions.)

Consequently if the temperature ratio $\frac{\widetilde{\theta}}{\widehat{\theta}}$ is in the pointwise sense in the right interval, then we can rewrite (3.33) as

[TABLE]

which would allow us to directly exploit Lemma 1. (We would obtain the same integrand on the both sides of the inequality.) Unfortunately, we are not able to show that the temperature perturbation remains in the pointwise sense in the interval that guarantees the validity of (3.34). The temperature perturbation $\widetilde{\theta}$ is only known to be uniformly bounded from below, a uniform bound from above is not available.

4 Rethinking the development of the candidate for a Lyapunov type functional

The inequality (3.34) that would guarantee the validity of the manipulations described above in Section 3.5 seems to introduce a technical rather than a physical restriction. We will try to overcome this apparently unphysical restriction via a suitable choice of the temperature scale. The key findings of this section are the following.

Lemma 5 (Family of Lyapunov type functionals and their time derivatives)

Let us consider the velocity/temperature perturbation $\widetilde{\theta}$, $\widetilde{\bm{v}}$ governed by equations (2.12). Let us define the functionals

[TABLE]

where we use the notation introduced in Section 2, and where $m\in(0,1)$. For fixed $m\in(0,1)$ any functional in family (4.1) remains non-negative for all possible temperature perturbations $\widetilde{\theta}\in(-\widehat{\theta},+\infty)$, and it vanishes if and only if the perturbation $\widetilde{\theta}$ vanishes everywhere in the domain of interest. Furthermore, the time derivative of the functionals is given by the formula

[TABLE]

Lemma 6 (On boundedness of some integrals)

Let us consider the velocity/temperature perturbation $\widetilde{\theta}$, $\widetilde{\bm{v}}$ governed by equations (2.12), and let $m\in(0,1)$, then

[TABLE]

4.1 Description based on an alternative temperature scale

So far we have used the absolute temperature scale $\theta$, and we have found that the corresponding natural Lyapunov type functional is almost the right one. Now we can speculate about the following idea. The fact that the perturbed temperature field returns back to the non-equilibrium steady state can not depend on the choice of the temperature scale. The temperature scale is nothing that is given a priori, the absolute temperature scale is used because of its convenience. (See for example Pippard (1964) for a thorough discussion of the concept of temperature, and the theoretical benefits of the absolute temperature scale. See also Fosdick and Rajagopal (1983) for a discussion of fundamentals of thermometry.) In fact any strictly increasing function of the absolute temperature can be used as an alternative temperature scale. This alternative temperature scale will provide us essentially the same tool for qualitative comparison of hotness of two bodies in the sense it will preserve ordering with respect to the hotness. On the other hand the quantitative aspects of the ordering will differ, which might help us to overcome some technical difficulties in dealing with the governing equations. This idea motivates the following manipulation.

The non-equilibrium steady temperature field is given by the solution of steady heat equation $\operatorname{div}\left(\kappa_{\mathrm{ref}}\nabla\theta\right)=0$. If we rewrite the equation as

[TABLE]

then nothing changes. Now we can introduce a new temperature scale $\vartheta$ as

[TABLE]

where $m<1$ is an arbitrary number. (The restriction to $m<1$ means that the new temperature scale preserves the ordering according to hotness. Values $m>1$ would lead to the inversion of ordering, since the new temperature would measure coldness instead of hotness, and the heat flux would change its sign.) Using this newly defined temperature scale, we see that (4.4) reads

[TABLE]

Consequently, if we decide to use the alternative temperature scale (4.5), then the fluid of interest effectively behaves as a fluid with the temperature dependent thermal conductivity

[TABLE]

and the equation for the non-equilibrium steady state reads $\operatorname{div}\left(\kappa(\vartheta)\nabla\vartheta\right)=0.$ If we make the same substitution in the complete evolution equation for the temperature (2.1c), then we see that (2.1c) reads

[TABLE]

Consequently, the fluid of interest effectively behaves as a fluid with the temperature dependent specific heat at constant volume

[TABLE]

and the evolution equation for the newly defined temperature reads

[TABLE]

This means that if we use the alternative temperature scale then we are effectively dealing with a fluid with the temperature dependent heat conductivity and the temperature dependent specific heat capacity at constant volume. Concerning this system of governing equations, we know how to construct a candidate for the Lyapunov functional. All that needs to be done is to identify the specific Helmholtz free energy, and use the method proposed by Bulíček et al. (2019). The specific Helmholtz free energy is a solution to the equation

[TABLE]

which yields $\psi=-c_{\mathrm{V},\mathrm{ref}}\vartheta_{\mathrm{ref}}\frac{\theta_{\mathrm{ref}}}{\vartheta_{\mathrm{ref}}}\frac{1}{m}\left[(1-m)\left(\frac{\vartheta}{\vartheta_{\mathrm{ref}}}\right)^{\frac{1}{1-m}}-\frac{\vartheta}{\vartheta_{\mathrm{ref}}}\right]$ . (Note that this formula gives us an additional restriction $m>0$. One of the integration constants in (4.11) has been chosen in such a way that the entropy vanishes at $\vartheta=\vartheta_{\mathrm{ref}}$. The other integration constant is immaterial since we are finally interested only in the differences between the corresponding quantities.) The other quantities such as the entropy and the internal energy are then found using standard thermodynamical identities. This yields

[TABLE]

The formula for the candidate for a Lyapunov functional which is obtained by the same method as in Section 3 reads

[TABLE]

In virtue of (4.5) we also see that $1+\frac{\widetilde{\vartheta}}{\widehat{\vartheta}}=\left(1+\frac{\widetilde{\theta}}{\widehat{\theta}}\right)^{1-m}$, hence the functional in (4.13) can be rewritten in terms of the original absolute temperature $\theta$ as

[TABLE]

Furthermore, the formula for the time derivative $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ rewritten again in the terms of the original absolute temperature $\theta$ is the formula (4.2) reported in Lemma 5. The proof of formula (4.2) is given in Appendix A.

4.2 Estimate of the time derivative of the prospective Lyapunov type functionals and boundedness of some integrals

Using the whole family of functionals $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ we need to revisit the estimates established in Section 3. The objective is to get the same type of estimates as in Section 3.3 and Section 3.4. Namely we will be first dealing with the last term on the right-hand side of (4.2), and we will try to show that it has a bounded time integral. Let us denote

[TABLE]

First, we note that in virtue of the boundary condition (2.14), the function $d_{\mathrm{th},\mathrm{m}}$ vanishes on the boundary, $\left.d_{\mathrm{th},\,m}\right|_{\partial\Omega}=0$. Furthermore, we see that $\left(1+\frac{\widetilde{\theta}}{\widehat{\theta}}\right)^{m}-1=d_{\mathrm{th},\,m}^{2}+2d_{\mathrm{th},\,m}$. Using the just introduced notation, we can rewrite (4.2) as

[TABLE]

Let us see whether the last term without the sign can be “absorbed” in the first term. We see that the $\varepsilon$-Young inequality implies

[TABLE]

where we have used the decay of the velocity perturbation (3.25). Furthermore, we see that

[TABLE]

and we can focus on the term $\int_{\Omega}d_{\mathrm{th},\,m}^{2}\left|\widetilde{\bm{v}}\right|\,\mathrm{d}\mathrm{v}$ only.

Before we proceed with the estimates, we derive a simple inequality that will be useful in a moment. From the Sobolev embedding, see Lemma 12, and the Poincaré inequality, see Lemma 13, we know that $\left\|d_{\mathrm{th},\,m}\right\|_{L^{2s}\left(\Omega\right)}\leq C_{\mathrm{S}}\left\|d_{\mathrm{th},\,m}\right\|_{W^{1,2}\left(\Omega\right)}\leq C_{\mathrm{S}}\sqrt{1+C_{\mathrm{P}}}\left\|\nabla d_{\mathrm{th},\,m}\right\|_{L^{2}\left(\Omega\right)}$, where $2s\in[1,p^{\star}]$. (Recall that in the three dimensional case we have $p^{\star}=6$.) In particular, we see that for all $s\in\left[1,\frac{p^{\star}}{2}\right]$ we get

[TABLE]

Going back to the right-hand side of (4.18), and using the Hölder inequality with $p_{1}=3$ and $p_{2}=\frac{3}{2}$ we see that

[TABLE]

The first integral on the right-hand side of (4.20) can be estimated using (4.19) where we choose $s=3$. The second integral can be estimated using the trivial embedding of Lebesgue spaces, see Lemma 11, in particular for $p_{1}=\frac{3}{2}$ and $p_{2}=2$ we get $\left\|\widetilde{\bm{v}}\right\|_{L^{\frac{3}{2}}\left(\Omega\right)}\leq\left|\Omega\right|^{\frac{1}{6}}\left\|\widetilde{\bm{v}}\right\|_{L^{2}\left(\Omega\right)}$. Finally, we see that

[TABLE]

where we have again used the decay of the velocity perturbation (3.25).

Using the derived estimates (4.21) and (4.17) we revisit the expression for the time derivative of the proposed Lyapunov type functional (4.16), and we see that the right-hand side of (4.16) can be estimated as follows

[TABLE]

where we have split the first gradient term on the right-hand side of (4.16) to halves, and where we have used Poincaré inequality. Inequality (4.22) can be further rewritten as

[TABLE]

First we note that both the coefficients $C_{\nabla d_{\mathrm{th},\,m}}$ and $C_{d_{\mathrm{th},\,m}}$ are positive in the case of spatially homogeneous temperature distribution $\widehat{\theta}$, and that the last term on the right-hand side vanishes in this case as well. Consequently, $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ is a genuine Lyapunov type functional if we are interested in the steady state with spatially homogeneous temperature distribution $\widehat{\theta}$.

Moreover, if $\varepsilon$ in (4.17) is chosen sufficiently small we see that $C_{d_{\mathrm{th},\,m}}$ is positive even if the steady state temperature distribution $\widehat{\theta}$ is spatially inhomogeneous. Furthermore, the negative part of the coefficient $C_{\nabla d_{\mathrm{th},\,m}}$ is a decreasing function of time, while its positive part remains constant. The coefficient $C_{\nabla d_{\mathrm{th},\,m}}$ therefore becomes positive once we reach a (possibly large) time threshold, and then it remains positive.

Finally, the last term in (4.23a) is integrable with respect to time from zero to infinity. Having derived (4.23a) we are in fact in the same position as in Section 3, see especially formula (3.28). In particular, we can follow the manipulations in Section 3.4 and show that all the terms on the right-hand side of (4.2) are finite, which proves Lemma 6.

5 Application of the lemma on the decay of integrable functions

From Section 3 and Section 4 we know that for $m\in[0,1)$ we have a family of functionals $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ that have favourable properties. (The results obtained in Section 3 can be seen as a formal limit $m\to 0+$ of the results obtained in Section 4.) The functionals are non-negative for any perturbation, and they vanish if and only if the perturbation vanishes in the whole domain of interest. Furthermore, the key terms in the formulae for the time derivatives are negative. Unfortunately, none of these functionals alone are sufficient for establishing stability, namely none of them are suitable as a Lyapunov functional. However, we might still benefit from the fact that we have the whole family of functionals, and we can try to combine them in a suitable manner.

In order to proceed with this idea further, we need two simple lemmas concerning the behaviour of certain functions, see Lemma 8 and Lemma 9 in Appendix B. Note that Lemma 9 implies that for $m,n\in(0,1)$, $n>m>\frac{n}{2}$ and $x\in(-1,\infty)$ we have the inequality

[TABLE]

See also Figure 3b for a sketch of the graphs of the corresponding functions.

Let us consider the family of functionals $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ introduced in Lemma 5, that is

[TABLE]

We know, see Lemma 5, that the evolution equation for this functional reads

[TABLE]

where we have introduced the notation

[TABLE]

already known from Section 4.2. From Section 4, we also know that all the quantities on the right-hand side of (5.3) have finite time integrals, see Lemma 4 and Lemma 6. (In particular all the integrals in the remainder term $\mathcal{R}^{\vartheta,m}$ have finite integral with respect to time.) Now we can subtract the evolution equations for $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$ and $\mathcal{V}_{\mathrm{meq}}^{\vartheta,n}$, which yields

[TABLE]

Using the formulae for $v_{\mathrm{th},\,m}$ and $v_{\mathrm{th},\,n}$ we see that

[TABLE]

which upon using Lemma 8 implies that if $n>m$ then

[TABLE]

while the equality holds if and only if $\frac{\widetilde{\theta}}{\widehat{\theta}}$ vanishes. (Let us recall that $\frac{\widetilde{\theta}}{\widehat{\theta}}\in(-1,+\infty)$, see Section 2.2.) Consequently, we see that the difference

[TABLE]

is for $n,m\in(0,1)$, $n>m$ a non-negative functional that vanishes if and only it the perturbation $\widetilde{\theta}$ vanishes.

Furthermore the right-hand side of (5.5), which is the time derivative of $\mathcal{Y}^{m,n}$, can be rewritten as

[TABLE]

where we have added and subtracted the quantity $\int_{\Omega}4\frac{1-m}{m^{2}}\kappa_{\mathrm{ref}}\widehat{\theta}\left|\nabla d_{\mathrm{th},\,n}\right|^{2}\,\mathrm{d}\mathrm{v}$. (We again recall that all the terms on the right-hand side have a finite integral if we integrate them from zero to infinity with respect to time.) Concerning the first term on the right-hand side of (5.9), we can apply Poincaré inequality, and we get

[TABLE]

Now for $n,m\in(0,1)$, $n>m>\frac{n}{2}$ we have the pointwise inequality

[TABLE]

This is a straightforward consequence of Lemma 9, see (5.1), and the definitions of $d_{\mathrm{th},\,m}$ and $v_{\mathrm{th},\,m}$, see (5.4a) and (5.2b). Once we have (5.11) we can go back to (5.10), and we get the estimate

[TABLE]

Consequently, the inequality (5.9) reads

[TABLE]

This inequality implies that

[TABLE]

We can observe that in virtue of Lemma 6 we know that $\int_{\tau=0}^{+\infty}\mathcal{H}^{m,n}\left(\left.\widetilde{\bm{W}}\right\|\widehat{\bm{W}}\right)\,\mathrm{d}\tau<+\infty$.

Now we are finally in the position to use Lemma 2 on the decay of integrable functions. First, we integrate (5.14a) over the time interval $(s,t)$, which yields

[TABLE]

If we take $s=0$ and $t=+\infty$, we see that the term on the left-hand side is nonegative, while the first and the last term on the right-hand side are finite and positive. Consequently, the negative term $\int_{\tau=0}^{+\infty}\mathcal{Y}^{m,n}\left(\left.\widetilde{\bm{W}}\right\|\widehat{\bm{W}}\right)\,\mathrm{d}\tau$ must be bounded

[TABLE]

Once we have established the boundedness of the integral, we can go back to (5.14a) and estimate the right-hand side of (5.14a) by its norm. This yields

[TABLE]

and we see that the assumptions of Lemma 2 are satisfied. Consequently, we can conclude that for $n,m\in(0,1)$, $n>m>\frac{n}{2}$ we have $\mathcal{Y}^{m,n}\left(\left.\widetilde{\bm{W}}\right\|\widehat{\bm{W}}\right)\xrightarrow{t\to+\infty}0$, which in virtue of the definition of $\mathcal{Y}^{m,n}$, see (5.8), means that

[TABLE]

which gives us the decay of temperature perturbations. Now we recall that the decay of the velocity perturbation has been obtained already by a trivial manipulation described in Section 3.3, hence we know the temperature/velocity perturbation decays as desired. We can summarise our findings as a theorem.

Theorem 5.1 (Decay of temperature/velocity perturbations)

Let us consider the temperature/velocity perturbation $\widetilde{\theta}$, $\widetilde{\bm{v}}$ to the steady temperature/velocity field $\widehat{\theta}$, $\widehat{\bm{v}}=\bm{0}$. Let us assume that the perturbation is the classical solution to the governing equations (2.12) in domain $\Omega$ with boundary conditions (2.14) and arbitrary initial conditions (2.13), and that this solution exists for all times. Let $m,n\in(0,1)$, $n>m>\frac{n}{2}$, then the temperature perturbation $\widetilde{\theta}$ satisfies

[TABLE]

Furthermore, the velocity perturbation $\widetilde{\bm{v}}$ satisfies $\left\|\widetilde{\bm{v}}\right\|_{L^{2}\left(\Omega\right)}^{2}\leq\left.\left\|\widetilde{\bm{v}}\right\|_{L^{2}\left(\Omega\right)}^{2}\right|_{t=0}{\mathrm{e}}^{-\frac{2\mu}{C_{P}\rho}t}$, where $C_{P}$ is the Poincaré constant for the domain $\Omega$.

The integral (5.19) neither introduces a Lebesgue type norm of the temperature perturbation, nor does it provide an upper bound on a Lebesgue type norm of the temperature perturbation. This might be interpreted as a weakness of the result, since the functional does not measure the distance from the steady state in a proper sense of the word. Let us however see what happens if we rewrite (5.19) in terms of the relative entropy $\eta_{\mathrm{diff}}=c_{\mathrm{V},\mathrm{ref}}\ln\left(1+\frac{\widetilde{\theta}}{\widehat{\theta}}\right)$, see (3.10) and (3.11). We see that

[TABLE]

Since the temperature $\theta=\widehat{\theta}+\widetilde{\theta}$ is bounded from below uniformly in space and time, see Section 2.2, we know that the relative entropy $\eta_{\mathrm{diff}}$ is also bounded from below uniformly in space and time. Concerning the behaviour of the function $f(x)=_{\mathrm{def}}\frac{1}{n}{\mathrm{e}}^{nx}-\frac{1}{m}{\mathrm{e}}^{mx}+\frac{n-m}{mn}$, that appears on the right-hand side of (5.20) we can use Lemma 10, and show that for any $l\in{\mathbb{N}}$, $l\geq 3$ there exists a constant $L$ such that

[TABLE]

(See Figure 4 for a visualisation of this inequality.) The inequality (5.21) then implies that

[TABLE]

where the integral on the right-hand side converges to zero as $t\to+\infty$. Consequently, we have proved the following corollary.

Corollary 1 (Decay of relative entropy)

Let the assumptions of Theorem 5.1 be fulfilled, and let $l\in[3,+\infty)$, then

[TABLE]

Note that the convergence to zero in Corrolary 1 is not uniform with respect to $l$, since the constant $L$ in (5.21) depends on $l$. The corollary shows that the relative entropy $\eta_{\mathrm{diff}}$ converges to zero in any Lebesgue space norm $L^{p}\left(\Omega\right)$, $p\in[1,+\infty)$, hence the relative entropy unlike the temperature perturbation exhibits the decay in a well established normed space.

6 Conclusion

We have investigated the stability of a spatially inhomogeneous non-equilibrium steady state in a thermodynamically open system. Namely we have considered an incompressible heat conducting viscous fluid occupying a vessel with walls kept at spatially non-uniform temperature. The steady state in this system is the rest state (zero velocity) with the temperature field given by the solution of the steady heat equation. This steady state is expected to be stable with respect to arbitrary perturbation.

This trivial fact is difficult to prove using the corresponding evolution equations. The main difficulty is twofold. First, the system of interest is a thermodynamically open system, which means that we have a priori no control on the fluxes through the boundary. Second, we consider the full system of governing equations including the dissipative heating term. This term is neglected in most of the mathematical works, see Joseph (1976a, b), Straughan (2004) and references therein. (Albeit these works aim at a different and more complex problem of thermal convection within the context of Oberbeck–Boussinesq approximation, the key difficulties are the same.) This holds also for works that deal with the temperature dependent thermal conductivity or temperature dependent viscosity, see for example Flavin and Rionero (1999) and Diaz and Straughan (2004). If the stability is considered including the dissipative heating, the available works see Richardson (1993) and Kagei et al. (2000) typically lead to conditional results, include ad hoc constructed functionals, require high regularity of the velocity field, and are limited to constant transport coefficients.

In our analysis we have assumed that the governing equations—the Navier–Stokes–Fourier equations for an incompressible fluid—have the classical solution that exists for all times. Using the thermodynamically motivated method for the construction of prospective Lyapunov type functionals, see Bulíček et al. (2019), we have then been able to address all these difficulties. The main theorem, see Theorem 5.1, then states that the temperature and velocity perturbations vanish as time goes to infinity, and that this behaviour holds for any choice of initial perturbation and for any shape of the vessel. This implies the recovery of the spatially inhomogeneous non-equilibrium steady state for any initial perturbation.

We note that the only piece of information regarding the behaviour of the dissipative heating term we have used in the stability analysis has been rather weak, which is a very welcome feature. In fact, we have used just the positivity of the term and its integrability in space and time. Consequently, it might be speculated that a similar thermodynamically based analysis is feasible not only for the incompressible Navier–Stokes–Fourier fluid but also for more complex models, provided that one knows a thermodynamic basis for such models444For a thermodynamic basis for some popular models see for example Rajagopal and Srinivasa (2000), Heida and Málek (2010), Málek et al. (2015), Hron et al. (2017), Málek et al. (2018), Málek et al. (2018), and Dostalík et al. (2019), where the thermodynamic basis is developed in the form directly applicable in the outlined stability analysis. Note that the basic thermodynamic building blocks are also provided for the models constructed within other thermodynamic frameworks such as the GENERIC framework, see Pavelka et al. (2018).. In particular, if the given model properly characterises the degradation of the mechanical energy to the thermal energy, then we would again get the integrability and positivity of the dissipative heating term, and we could essentially follow the same argument as in the present contribution. The only difference would be a different formula for the dissipative heating, and a more involved treatment of the term arising from the convective derivative.

Once it is clear how to handle the dissipative heating term in the stability analysis of a (physically) simple system, one might speculate about more ambitious applications. In particular, one can again consider a heat conducting viscous fluid occupying a vessel with walls kept at spatially non-uniform temperature. This time one can however consider a fluid with a temperature dependent density, which leads to the Rayleigh–Bénard convection problem. In such a setting the presence of the buoyancy term makes the stability analysis more challenging, since the evolution equation for the thermal and mechanical variables are tightly coupled. The expected result in the Rayleigh–Bénard setting would be that the classical results, see Joseph (1976b), are not substantially altered by the presence of the dissipative heating term. Such an analysis is however beyond the scope of this work, and it will be carried out in the future.

Appendix A Derivation of the formulae for the time derivatives of the functionals

The derivation presented in the main text can be greatly simplified if we restrict ourselves to formal manipulations. We recall that the main objective is to find the time derivative of the functional

[TABLE]

where we have used the identity $\widehat{\theta}\left.\frac{\partial{\eta}}{\partial{\theta}}\right|_{\theta=\widehat{\theta}}=\left.\frac{\partial{e}}{\partial{\theta}}\right|_{\theta=\widehat{\theta}}$, which holds universally for any material. In the particular case of the fluid described by the Helmholtz free energy in the form (2.2), we get the explicit formula for the functional $\mathcal{V}_{\mathrm{meq}}$, see Lemma 3, equation (3.1). It is convenient to rewrite (3.1) in terms of the relative entropy, $\eta_{\mathrm{diff}}=c_{\mathrm{V},\mathrm{ref}}\ln\left(1+\frac{\widetilde{\theta}}{\widehat{\theta}}\right)$, see (3.10) and (3.11). If we do so, the counterpart of (3.1) reads

[TABLE]

We know that the evolution equation for the entropy of the perturbed state reads

[TABLE]

which in our case of constant heat conductivity $\kappa$ and vanishing steady state velocity $\widehat{\bm{v}}=\bm{0}$ simplifies to

[TABLE]

The last equation can be again rewritten in terms of the relative entropy,

[TABLE]

Note that at this point we heavily exploit the fact that in our case the relative entropy is given by a formula that nicely matches with the structure of the heat flux term. Using the evolution equation for the relative entropy, one can derive the evolution equation for any quantity of the form $c_{\mathrm{V},\mathrm{ref}}\widehat{\theta}f({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}})$ , where $f$ is a given function, see Lemma 7, equation (A.6). Using the lemma for particular choices of $f$ then allows us to seamlessly derive formulae for the time derivatives of the functionals introduced in Section 3 and Section 4.

Lemma 7 (Pointwise evolution equation for functions of the exponential of the relative entropy)

Let $\frac{\mathrm{d}{{}_{\widetilde{\bm{v}}}}}{\mathrm{d}{t}}=\frac{\partial{}}{\partial{t}}+\widetilde{\bm{v}}\bullet\nabla$ denote the material time derivative with respect to the perturbed velocity field $\widetilde{\bm{v}}$, and let $f$ denote a given function. The evolution equation for the quantity $c_{\mathrm{V},\mathrm{ref}}\widehat{\theta}f({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}})$ reads

[TABLE]

Proof

The proof is based on the direct computation. Using the evolution equation for the relative entropy, see (A.5), we see that

[TABLE]

where $f^{\prime}$ denotes the derivative of function $f$ with respect to its argument. The first term on the right-hand side can be rewritten as a divergence term plus a source term,

[TABLE]

Now we use the identity ${\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}}\nabla f^{\prime}\left({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}}\right)=\nabla\left({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}}f^{\prime}\left({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}}\right)-f\left({\mathrm{e}}^{\frac{\eta_{\mathrm{diff}}}{c_{\mathrm{V},\mathrm{ref}}}}\right)\right)$, which allows us to rewrite the corresponding term on right-hand side of (A.8) as

[TABLE]

where we have used the fact that $\operatorname{div}\left(\kappa_{\mathrm{ref}}\nabla\widehat{\theta}\right)=0$. Using (A.9) on the right-hand side of (A.8) reveals that

[TABLE]

This finishes the manipulation with the first term on the right-hand side of (A.7). Now we focus on the prospective convective term on the right-hand side of (A.7) . We see that

[TABLE]

Next we observe that $\nabla\widehat{\eta}=c_{\mathrm{V},\mathrm{ref}}\frac{\nabla\widehat{}\theta}{\widehat{\theta}}$, and using (A.10) and (A.11) in (A.7) we finally obtain the evolution equation (A.6).

Now we are in the position to exploit Lemma 7 in the derivation of formulae for the time derivative of the functionals introduced in Lemma 3 and Lemma 5.

Indeed, if we set

[TABLE]

then $f^{\prime}(y)=1-\frac{1}{y}$, $f^{\prime\prime}(y)=\frac{1}{y^{2}}$, and using Lemma 7 we obtain the pointwise evolution equation

[TABLE]

The last equation can be upon straightforward manipulations and upon using the fact that $\zeta_{\mathrm{mech}}\left(\widehat{\bm{W}}+\widetilde{\bm{W}}\right)=2\mu\widetilde{\mathbb{D}}:\widetilde{\mathbb{D}}$ rewritten in terms of the temperature

[TABLE]

This is the sought pointwise evolution equation for the thermal part of the integrand in the functional $\mathcal{V}_{\mathrm{meq}}$. Integrating (A.14) over the (material) domain $\Omega$ leads directly to the result reported in Lemma 3. (Recall that the boundary condition $\left.\widetilde{\theta}\right|_{\partial\Omega}=0$ guarantees that the boundary terms vanish.)

On the other hand, if we choose $m\in{\mathbb{R}}$, and if we set

[TABLE]

then $f^{\prime}(y)=1-y^{m-1}$, $f^{\prime\prime}(y)=-(m-1)y^{m-2}$, and using Lemma 7 we obtain the pointwise evolution equation

[TABLE]

Using the specific formula for the mechanical dissipation, and rewriting (A.13) in terms of the temperature we get

[TABLE]

This is the sought pointwise evolution equation for the thermal part of the integrand in the functional $\mathcal{V}_{\mathrm{meq}}^{\vartheta,\,m}$. Integrating (A.14) over the (material) domain $\Omega$ leads directly to the result reported in Lemma 5.

Appendix B Auxiliary tools

Lemma 8

Let $m,n\in(0,1)$ and let $n>m$. Let $x\in(-1,+\infty)$ and let us define the function

[TABLE]

The function $f(x,m,n)$ is for $x\in(-1,+\infty)$ a non-negative function which vanishes if and only if $x=0$.

Proof

Taking the derivative of (B.1) with respect to $x$ yields $\frac{\mathrm{d}{}}{\mathrm{d}{x}}f(x,m,n)=\left(1+x\right)^{m-1}\left[\left(1+x\right)^{n-m}-1\right]$. If $n>m$, then the derivative of $f(x,m,n)$ is positive for $x>0$ and negative for $x<0$. Further the value of $f(x,m,n)$ at $x=0$ is $f(0,m,n)=0$. Consequently, function $f(x,m,n)$ is for all $x\in(-1,+\infty)$ a non-negative function which vanishes if and only if $x=0$.

Lemma 9

Let $m,n\in(0,1)$ and let $n>m>\frac{n}{2}$. Let $x\in(-1,+\infty)$ and let us define the function

[TABLE]

The function $g(x,m,n)$ has the limits $\lim_{x\to-1+}g(x,m,n)=-3+\frac{n}{m}<0,$ and $\lim_{x\to+\infty}g(x,m,n)=-\infty,$ and the function $g(x,m,n)$ is non-positive for all $x\in(-1,+\infty)$, and it vanishes if and only if $x=0$.

Proof

Apparently, the function vanishes for $x=0$, and the computation of the limits is straightforward. If we take the derivative of this function with respect to $x$ we get

[TABLE]

and we see that the derivative is negative for $x>0$, and it is positive for $x<0$. Consequently, we see that $g(x,m,n)\leq 0$ for all $x\in(-1,+\infty)$, and that it vanishes if and only if $x=0$.

Lemma 10

Let $m,n\in(0,1)$ and $n>m>\frac{n}{2}$, and let $f(x)=_{\mathrm{def}}\frac{1}{n}{\mathrm{e}}^{nx}-\frac{1}{m}{\mathrm{e}}^{mx}+\frac{n-m}{mn}$. Let $x_{\mathrm{crit}}$ is an arbitrary negative number, and let $l$ is a given natural number, $l\geq 3$, then there exists a constant $L$ (possibly large) such that $\frac{\left|x\right|^{l}}{L}\leq f(x)$ holds for all $x\in[x_{\mathrm{crit}},+\infty)$.

Proof

Using the series expansion for the exponential we can observe that $\frac{1}{n}{\mathrm{e}}^{nx}-\frac{1}{m}{\mathrm{e}}^{mx}+\frac{n-m}{mn}=\sum_{k=2}^{+\infty}\frac{n^{k-1}-m^{k-1}}{k!}x^{k}$, and since $n>m$, we see that the inequality

[TABLE]

holds for any non-negative $x$ and arbitrary $l\in{\mathbb{N}}$, $l\geq 2$, and the equality occurs if and only if $x=0$. Hence if we define function $h$ as

[TABLE]

we see that $h(x)<f(x)$ for all positive $x$ and $h(x)=f(x)$ for $x=0$.

Moreover, a straightforward calculation also shows that the function $f(x)-h(x)$ has for $l\geq 3$ strict local minimum at zero, hence $h(x)<f(x)$ even in some left neighborhood of zero, that is for some $x\in(-\varepsilon,0)$. On the other hand function $h(x)$ is for $x<0$ unbounded, while $f(x)$ is bounded, hence they must intersect, $h(x_{\mathrm{int}})=f(x_{\mathrm{int}})$, at some point $x_{\mathrm{int}}$, where $x_{\mathrm{int}}<0$. (Among possible intersecting points we chose the one with the smallest magnitude.) If $x_{\mathrm{int}}\leq x_{\mathrm{crit}}$, then we are done, and we have found the desired function. In this case $L$ is given by the formula $\frac{1}{L}=\frac{n^{l-1}-m^{l-1}}{l!}$.

If $x_{\mathrm{crit}}<x_{\mathrm{int}}$ it remains to flatten the graph of $h(x)$ such that the intersection point is moved sufficiently far to the left. This can be done by the means of the transformation $g(x)=_{\mathrm{def}}\frac{1}{K}h\left(x\right),$ where $K$ is a sufficiently large positive constant. If we choose $\frac{1}{K}=_{\mathrm{def}}\frac{h(x_{\mathrm{int}})}{h(x_{\mathrm{crit}})}$, then $\frac{1}{K}<1$. (Function $h$ is for $x<0$ a strictly decreasing function.) We can observe that $g(x)<f(x)$ for $x\in[x_{\mathrm{int}},+\infty)$. On the other hand, if $x\in[x_{\mathrm{crit}},x_{\mathrm{int}})$ then

[TABLE]

where we have exploited the fact that $f$ is for $x<0$ a strictly decreasing function. Consequently $g(x)\leq f(x)$ holds for all $x\in[x_{\mathrm{crit}},+\infty)$ as desired. Moreover, the constant $L$ is given by the formula $\frac{1}{L}=\frac{f(x_{\mathrm{int}})}{\left|x_{\mathrm{crit}}\right|^{l}}$.

Appendix C Auxiliary tools – inequalities in function spaces

For further reference we recall the standard embedding theorems for Lebesgue and Sobolev spaces, as well as other standard inequalities, for proofs see for example Evans (1998) or Adams and Fournier (2003).

Lemma 11 (Trivial embedding of Lebesgue spaces)

Let $\Omega$ be a domain with a finite volume, and let $p_{1},p_{2}\in[1,+\infty]$ such that $p_{2}\geq p_{1}$, then $\left\|f\right\|_{L^{p_{1}}\left(\Omega\right)}\leq\left|\Omega\right|^{\frac{p_{2}-p_{1}}{p_{1}p_{2}}}\left\|f\right\|_{L^{p_{2}}\left(\Omega\right)}$.

Lemma 12 (Continuous embedding of Sobolev spaces into Lebesgue spaces)

Let $\Omega\subset{\mathbb{R}}^{n}$ be a domain with $\mathcal{C}^{0,1}$ boundary and let $p\in[1,n)$. Then for all $q\in[1,p^{\star}]$ where $p^{\star}=_{\mathrm{def}}\frac{np}{n-p}$ there exists a constant $C_{\mathrm{S}}$ that depends on $\Omega$, $p$, $n$ and $q$ such that $\left\|f\right\|_{L^{q}\left(\Omega\right)}\leq C_{\mathrm{S}}\left\|f\right\|_{W^{1,p}\left(\Omega\right)}$. The constant $C_{\mathrm{S}}$ is referred to as the Sobolev embedding constant.

Lemma 13 (Poincaré inequality)

Let $f\in W_{0}^{1,p}\left(\Omega\right)$, where $p\in[1,+\infty)$ and $\Omega\in{\mathcal{C}}^{0,1}$. Then there exists a constant $C$ that depends only on $\Omega$ and $p$ such that $\left\|f\right\|_{L^{p}\left(\Omega\right)}\leq C(\Omega,p)\left\|\nabla f\right\|_{L^{p}\left(\Omega\right)}$. In particular, if $p=2$ then $\left\|f\right\|_{L^{2}\left(\Omega\right)}^{2}\leq C_{\mathrm{P}}\left\|\nabla f\right\|_{L^{2}\left(\Omega\right)}^{2}$, where $C_{\mathrm{P}}$ is referred to as the Poincaré constant.

Lemma 14 (Hölder inequality)

Let $\left\{p_{i}\right\}_{i=1}^{k}$ be a sequence of exponents such that $p_{i}\in[1,+\infty]$ and $\sum_{k=1}^{+\infty}\frac{1}{p_{k}}=1$. Let us consider functions $\left\{f_{i}\right\}_{i=1}^{k}$, such that $f_{i}\in L^{p_{i}}\left(\Omega\right)$, then $\left\|\prod_{i=1}^{k}f_{i}\right\|_{L^{1}\left(\Omega\right)}\leq\prod_{i=1}^{k}\left\|f_{i}\right\|_{L^{p_{i}}\left(\Omega\right)}$.

Lemma 15 ($\varepsilon$-Young inequality)

Let $a,b\in{\mathbb{R}}^{+}$, $\varepsilon\in{\mathbb{R}}^{+}$ and $p,q\in(1,+\infty)$ such that $\frac{1}{p}+\frac{1}{q}=1$. Then $ab\leq\varepsilon a^{p}+\frac{1}{\left(\varepsilon p\right)^{\frac{q}{p}}q}b^{q}$.

Lemma 16 (Korn equality)

Let $\Omega$ be a bounded domain with smooth boundary. Let $\bm{v}$ be a smooth vector field that vanishes on the boundary $\left.\bm{v}\right|_{\partial\Omega}=0$, then $2\int_{\Omega}{\mathbb{D}:\mathbb{D}}\,\mathrm{d}\mathrm{v}=\int_{\Omega}{\nabla\bm{v}:\nabla\bm{v}}\,\mathrm{d}\mathrm{v}+\int_{\Omega}{\left(\operatorname{div}\bm{v}\right)^{2}}\,\mathrm{d}\mathrm{v}$.