Generalized Boussinesq System with Energy Dissipation: Existence of Stationary Solutions

Abstract

In this paper, we investigate the solvability of a boundary value problem for a heat and mass transfer model with the spatially averaged Rayleigh function. The considered model describes the 3D steady-state non-isothermal flow of a generalized Newtonian fluid (with shear-dependent viscosity) in a bounded domain with Lipschitz boundary. The main novelty of our work is that we do not neglect the viscous dissipation effect in contrast to the classical Boussinesq approximation, and hence, deal with a system of strongly nonlinear partial differential equations. Using the properties of the averaging operation and d-monotone operators as well as the Leray–Schauder alternative for completely continuous mappings, we prove the existence of weak solutions without any smallness assumptions for model data. Moreover, it is shown that the set of all weak solutions is compact, and each solution from this set satisfies some energy equalities.

1. Introduction and Problem Statement

The classical Boussinesq approximation system [1,2] was proposed as a simplified model for heat and mass transfer in a linear viscous fluid and is widely used in studying non-isothermal flows (see, for example, [3,4,5,6,7,8,9,10] and the references cited therein). This system of partial differential equations (PDEs) includes the motion equations, the incompressibility condition, and the energy balance equation. Typically, the nonlinear term characterizing the conversion of kinetic energy into thermal energy due to the viscous friction effect is ignored in the energy balance equation. This term is commonly referred to as the Rayleigh dissipation function or simply the dissipation function.

Clearly, the replacement of the dissipation function by zero is reasonable only under some conditions, namely in situations when its values are small compared to other terms in the energy balance equation, allowing it to be neglected. The main motivation for this simplification lies in the fact that ignoring the effect of viscous dissipation greatly facilitates mathematical analysis and finding of solutions to heat and mass transfer models that are based on the Boussinesq system. However, from a physical point of view as well as for certain practical applications, it is interesting to consider the “full equations”, that is, the equations that include all the nonlinear terms [11,12,13,14,15,16]. Studying the Boussinesq system with energy dissipation is important as it offers insights into complex dynamics of fluid motion and energy transfer, contributing to advancements in environmental science, engineering, and climate modelling. On the other hand, the mathematical difficulties carried by these PDEs are interesting by themselves.

The primary challenges in constructing solutions are related to deriving a priori estimates for a temperature distribution and establishing the convergence of Galerkin’s approximate solutions. While the compactness method based on standard a priori estimates for a velocity field successfully handles the limit passage in all terms of the motion equation, these estimates are insufficient for the limit passage in the strongly nonlinear energy balance equation. To overcome these difficulties, in the present work, we will use the following approach: replacing the Rayleigh dissipation function with its spatially averaged (regularized) version [17].

We consider the following model for the 3D steady-state non-isothermal flow of a generalized Newtonian fluid (with variable viscosity) through a porous media:

(1)$\left\{\begin{array}{cc}\hfill (\mathit{u}·\nabla )\mathit{u}=\mathrm{div}\left[\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right)\right]-\nabla \pi -\alpha \mathit{u}+\mathit{f}(\mathit{x},\theta )& \mathrm{in}\Omega ,\hfill \\ \hfill \mathrm{div}\mathit{u}=0& \mathrm{in}\Omega ,\hfill \\ \hfill (\mathit{u}·\nabla )\theta =k\Delta \theta +{\Phi }_{\rho }\left(\mathit{u}\right)+g(\mathit{x},\theta )& \mathrm{in}\Omega ,\hfill \\ \hfill \mathit{u}=\mathbf{0}& \mathrm{on}\partial \Omega ,\hfill \\ \hfill k{\displaystyle \frac{\partial \theta }{\partial \nu }}+\beta \theta =0& \mathrm{on}\partial \Omega .\hfill \end{array}\right.$

Here,

• $\Omega$ is the flow domain, $\Omega \subset {\mathbb{R}}^3$;

• $\partial \Omega$ denotes the boundary of $\Omega$;

• $\mathit{u}$ is the flow velocity;

• $\pi$ is the pressure;

• $\mathcal{D}\left(\mathit{u}\right)$ denotes the rate-of-deformation tensor with the components

  ${\mathcal{D}}_{ij}\left(\mathit{u}\right):=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),i,j\in \{1,2,3\};$

• $\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)>0$ is the effective viscosity, dependent on the Euclidean norm of the rate-of-deformation tensor;

• $\theta$ is the temperature;

• $\mathit{f}(\mathit{x},\theta )$ is the body force acting on the fluid;

• k stands for the thermal conductivity coefficient, $k>0$;

• $\alpha$ is the Darcy (permeability of porous medium) coefficient, $\alpha >0$;

• ${\Phi }_{\rho }\left(\mathit{u}\right):={\displaystyle \frac{1}{c_p}}\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2$ represents the Rayleigh dissipation function with the spatially averaged rate-of-deformation tensor ${\mathcal{D}}_{\rho }\left(\mathit{u}\right)$ (see [18]) that is defined as follows:

  (2)${\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left(\mathit{x}\right):={\int }_{{\mathbb{R}}^3}\rho (\mathit{x}-\mathit{y})\mathcal{D}\left(\tilde{\mathit{u}}\right)\left(\mathit{y}\right)\mathit{dy},$

  where

  $\tilde{\mathit{u}}\left(\mathit{x}\right):=\left\{\begin{array}{c}\mathit{u}\left(\mathit{x}\right)\mathrm{if} \mathit{x}\in \Omega ,\hfill \\ \mathbf{0}\mathrm{if} \mathit{x}\in {\mathbb{R}}^3\setminus \Omega ,\hfill \end{array}\right.$

  and $\rho :{\mathbb{R}}^3\to \mathbb{R}$ is a smooth function with compact support such that

  ${\int }_{{\mathbb{R}}^3}\rho \left(\mathit{y}\right)\mathit{dy}=1$

  and $\rho \left(\mathit{x}\right)=\rho \left(\mathit{y}\right)$ for $\mathit{x}$ and $\mathit{y}$ with the same Euclidean norm;

• $c_p$ signifies the specific heat capacity of the fluid, $c_p>0$;

• $\beta$ stands for the heat exchange coefficient at the walls of the vessel $\Omega$, $\beta >0$;

• $g(\mathit{x},\theta )$ is the heat source intensity;

• $\mathit{\nu }$ denotes the unit outward normal vector to the surface $\partial \Omega$.

It is worth noting that, as the diameter of $\mathrm{supp}\left(\rho \right)$ tends to zero, the averaged tensor ${\mathcal{D}}_{\rho }\left(\mathit{u}\right)$ converges to the original tensor $\mathcal{D}\left(\mathit{u}\right)$ in the ${\mathit{L}}^q$-norm. Therefore, we can achieve high quality of approximation. The averaging operation defined by (2) also possesses other useful properties [19] and is widely used in mathematical hydrodynamics. For instance, this operation is employed to establish differentiability properties of generalized solutions to the Navier–Stokes equations [20], as well as to investigate the solvability of boundary and initial-boundary value problems for models of non-Newtonian fluids with objective derivatives in rheological equations [18,21,22].

Rheological models with shear-dependent viscosity are applied in studying multiphase flows of such media as polymer solutions, concentrated suspensions, clayey mixtures, corn flour in water, oil paints, molten steel, and others (see [23,24]).

Within the scope of the present paper, our main aim is to establish the solvability of problem (1) in the weak formulation. The proof of the existence theorem proceeds as follows. We first introduce an auxiliary system of nonlinear equations parametrized by $\lambda \in [0,1]$. A priori bounds for solutions of this system are derived, which do not depend on $\lambda$. Next, we interpret the new problem as a one-parameter operator equation and study some properties of the corresponding operators. Specifically, we demonstrate that, under reasonable assumptions regarding the model data, one of the operators is d-monotone, while the others are completely continuous. By applying the Leray–Schauder fixed-point principle, we conclude the existence of weak solutions to problem (1) in the Cartesian product of some subspaces of the Sobolev space $H^1(\Omega )$. Moreover, we establish that the set of all weak solutions is compact in these subspaces. Note that this work is a continuation of [25], where a similar heat and mass transfer model is considered without taking into account dissipative heating in the energy balance equation.

The mathematical analysis of flow models with shear-dependent viscosity started with the works by Ladyzhenskaya [26,27] and Litvinov [23]. These authors exclusively considered isothermal flows (for both the steady-state and unsteady cases) and, by the Ritz, Galerkin, and Faedo–Galerkin schemes, constructed solutions to the corresponding boundary and initial-boundary value problems. Their results have been extended to a flow model with nonlocal boundary conditions, where the fluid slips along an impermeable solid boundary of the flow domain [28], as well as to an optimal flow control problem [29] and a model for a rigid viscoplastic media of the Bingham kind with threshold slippage [30]. It also should be mentioned at this point that there exists extensive literature devoted to studying optimization and control problems for equations governing heat and mass transfer in a fluid with constant viscosity [31,32,33,34,35,36,37,38]. Finally, we mention the works [39,40,41], in which the unique solvability of the evolutionary Boussinesq system with constant viscosity and energy dissipation is proved under suitable smallness assumptions on data. In contrast to the case of the Navier–Stokes system, for such models (even in the weak formulation of associated initial-boundary value problems), there are to date no results on the existence of large-data solutions that are global in time.

The remainder of the present paper is organized as follows. The next section is completely devoted to necessary preliminaries. In Section 3, we provide a description of assumptions on model data. Moreover, this section contains the rigorous definition of weak solutions to problem (1). In Section 4, we formulate and prove the main results of this work—Theorem 1 on the existence of weak solutions to problem (1) and their properties.

2. Preliminaries

In this section, we collect notions and statements that will be needed for establishing our main results.

2.1. Notations and Function Spaces

Let

$\mathit{X}:\mathit{Y}:=\sum _{i,j=1}^3{X}_{ij}{Y}_{ij},\left|\mathit{X}\right|:=\sqrt{\mathit{X}:\mathit{X}},$

for $3\times 3$-matrices $\mathit{X}$ and $\mathit{Y}$.

As usual, $\mathbb{R}:=(-\infty ,+\infty )$, ${\mathbb{R}}_{+}:=[0,+\infty )$, and ${\mathbb{R}}^n$ stands for the n-dimensional Euclidean space.

Let $\mathbf{E}$ and $\mathbf{F}$ be normed linear spaces. By $\mathcal{L}(\mathbf{E},\mathbf{F})$, we denote the space of all continuous linear operators from $\mathbf{E}$ into $\mathbf{F}$.

By ⇀ (→, ⇉, resp.) we denote weak (strong, uniform, resp.) convergence. The symbol ↪ stands for a continuous embedding, while $\hookrightarrow \hookrightarrow$ denotes a compact embedding.

Let $\Omega$ be a bounded domain in ${\mathbb{R}}^3$ with the Lipschitz-continuous boundary $\partial \Omega$. We will use the Lebesgue spaces $L^p(\Omega )$ with $p\ge 1$ and the Sobolev space $H^1(\Omega ):=W^{1,2}(\Omega )$. The corresponding classes of vector-valued functions $\mathit{w}:\Omega \to {\mathbb{R}}^3$ are denoted by the same symbols but in bold font, that is, ${\mathit{L}}^p(\Omega ):=L^p{(\Omega )}^3$ and ${\mathit{H}}^1(\Omega ):=H^1{(\Omega )}^3.$ The definitions of these spaces and analysis of their properties can be found in [42,43,44].

When solving problem (1) in the space ${\mathit{H}}^1(\Omega )$, it is convenient to use the norm

${\parallel \omega \parallel }_{H^1(\Omega )}:={\left({\int }_{\Omega }{|\nabla \omega |}^2\mathit{dx}+{\int }_{\partial \Omega }{\left|\omega \right|}^{2}dS\right)}^{\frac{1}{2}},$

which is equivalent to the standard $H^1$-norm, computed as follows:

${\parallel \omega \parallel }_1:={\left({\int }_{\Omega }{\left|\omega \right|}^2+{|\nabla \omega |}^2\mathit{dx}\right)}^{\frac{1}{2}}.$

We now provide two function spaces developed in the theory of Navier–Stokes equations:

$\begin{array}{cc}\hfill {\mathit{C}}_{0,\sigma }^{\infty }(\Omega )& :=\left\{\mathit{v}\in {\mathit{C}}^{\infty }(\Omega ):\mathrm{supp}\left(\mathit{v}\right)\subset \Omega \mathrm{and}\mathrm{div}\mathit{v}\equiv 0\mathrm{in}\mathrm{the}\mathrm{domain}\Omega \right\},\hfill \\ \hfill \mathit{V}(\Omega )& :=\mathrm{the}\mathrm{closure}\mathrm{of}\mathrm{the}\mathrm{set}{\mathit{C}}_{0,\sigma }^{\infty }(\Omega )\mathrm{in}\mathrm{the}\mathrm{space}{\mathit{H}}^1(\Omega ).\hfill \end{array}$

We equip $\mathit{V}(\Omega )$ with the scalar product defined as follows:

${(\mathit{v},\mathit{w})}_{\mathit{V}(\Omega )}:={\int }_{\Omega }\mathcal{D}\left(\mathit{v}\right):\mathcal{D}\left(\mathit{w}\right)\mathit{dx},$

for arbitrary vector functions $\mathit{v},\mathit{w}\in \mathit{V}(\Omega )$. By applying the Korn inequality (see, for example, [23], Chapter I, § 2), one can verify that this scalar product is well defined. Moreover, the Euclidean norm ${\parallel \mathit{v}\parallel }_{\mathit{V}(\Omega )}:={(\mathit{v},\mathit{v})}_{\mathit{V}(\Omega )}^{1/2}$ is equivalent to the norm naturally induced from the space ${\mathit{H}}^1(\Omega )$.

As the main space in which we will find solutions to problem (1), we consider the Cartesian product $\mathit{V}(\Omega )\times H^1(\Omega )$ equipped with the max-norm:

${\parallel (\mathit{v},\omega )\parallel }_{\mathit{V}(\Omega )\times H^1(\Omega )}:=max\left\{{\parallel \mathit{v}\parallel }_{\mathit{V}(\Omega )},{\parallel \omega \parallel }_{H^1(\Omega )}\right\}.$

2.2. Some Properties of Spatially Averaged Rate-of-Deformation Tensor

Lemma 1.

Suppose Ω is a bounded Lipschitz-continuous domain of ${\mathbb{R}}^3$, the tensor ${\mathcal{D}}_{\rho }\left(\mathit{u}\right)$ is defined by (2), and

$\begin{array}{cc}\hfill M_i(\rho ,\Omega )& :=\underset{\mathit{x}\in {\mathbb{R}}^3}{max}\left|\frac{\partial \rho }{\partial x_i}\right|\sqrt{\mathrm{meas}(\Omega )},i\in \{1,2,3\},\hfill \\ \hfill M(\rho ,\Omega )& :=max\left\{M_1(\rho ,\Omega ),M_2(\rho ,\Omega ),M_3(\rho ,\Omega )\right\}.\hfill \end{array}$

Then the operator ${\mathcal{D}}_{\rho }:\mathit{V}(\Omega )\to {\mathit{C}}^{\infty }\left(\overline{\Omega }\right)$ has the following properties.

• (i). For any vector function $\mathit{u}\in \mathit{V}(\Omega )$, the estimate

  (3) $\parallel {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}\le M(\rho ,\Omega ){\parallel \mathit{u}\parallel }_{{\mathit{L}}^2(\Omega )}$

  is valid.

• (ii). If ${\left\{{\mathit{u}}_{\ell }\right\}}_{\ell =1}^{\infty }$ is a sequence such that

  (4) ${\mathit{u}}_{\ell }\rightharpoonup {\mathit{u}}_{0}inthespace\mathit{V}(\Omega )as\ell \to \infty ,$

  then

  (5) ${\mathcal{D}}_{\rho }\left({\mathit{u}}_{\ell }\right)\rightrightarrows {\mathcal{D}}_{\rho }\left({\mathit{u}}_0\right)on\overline{\Omega }as\ell \to \infty .$

Proof. 

Let $\mathit{u}\in \mathit{V}(\Omega )$. Applying integration by parts, we obtain

(6)$\begin{array}{cc}\hfill {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)}_{ij}\left(\mathit{x}\right)& =\frac{1}{2}{\int }_{\Omega }\rho (\mathit{x}-\mathit{y})\left(\frac{\partial u_i\left(\mathit{y}\right)}{\partial y_j}+\frac{\partial u_j\left(\mathit{y}\right)}{\partial y_i}\right)\mathit{dy}\hfill \\ \hfill & =-\frac{1}{2}{\int }_{\Omega }{u}_i\left(\mathit{y}\right)\frac{\partial }{\partial y_j}\rho (\mathit{x}-\mathit{y})\mathit{dy}\hfill \\ \hfill & -\frac{1}{2}{\int }_{\Omega }{u}_j\left(\mathit{y}\right)\frac{\partial }{\partial y_i}\rho (\mathit{x}-\mathit{y})\mathit{dy},\forall \mathit{x}\in \overline{\Omega },i,j\in \{1,2,3\},\hfill \end{array}$

and, by using the Cauchy–Schwarz inequality, we derive inequality (3).

From the Rellich–Kondrachov theorem (see, for example, [45], Chapter II, § 2.6.1) it follows that

$\mathit{V}(\Omega )\hookrightarrow \hookrightarrow {\mathit{L}}^p(\Omega ),\forall p\in [1,6),$

which together with (4) yield

(7)${\mathit{u}}_{\ell }\to {\mathit{u}}_0\mathrm{in}\mathrm{the}\mathrm{space}{\mathit{L}}^2(\Omega )\mathrm{as}\ell \to \infty .$

Moreover, in view of (3), we have

(8)$\parallel {\mathcal{D}}_{\rho }({\mathit{u}}_{\ell }-{\mathit{u}}_0){\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}\le M(\rho ,\Omega ){\parallel {\mathit{u}}_{\ell }-{\mathit{u}}_0\parallel }_{{\mathit{L}}^2(\Omega )}.$

Combining (7) and (8), we obtain

$\underset{\ell \to \infty }{lim}{\parallel {\mathcal{D}}_{\rho }({\mathit{u}}_{\ell }-{\mathit{u}}_0)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}=0.$

This implies (5). Thus, Lemma 1 is proved. □

Taking into account relations (3) and (6), we arrive at the following statement.

Lemma 2.Under the conditions of Lemma 1, the operator ${\mathcal{D}}_{\rho }:{\mathit{L}}^2(\Omega )\to \mathit{C}\left(\overline{\Omega }\right)$ is well defined and continuous.

2.3. Continuous Invertibility of Monotone-Type Operators in Banach Spaces

When studying flow models for fluids with shear-dependent viscosity, it is convenient to use the notion of so-called d-monotone operators. Let us recall the definition of the d-monotonicity property.

Let $\mathbf{E}$ be a reflexive separable Banach space with the dual space ${\mathbf{E}}^{*}$. By ${⟨·,·⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}$ denote a duality pairing between ${\mathbf{E}}^{*}$ and $\mathbf{E}$.

Definition 1.An operator $\mathbf{A}:\mathbf{E}\to {\mathbf{E}}^{*}$ is called d-monotone if there exists a strictly increasing function $\varkappa :{\mathbb{R}}_{+}\to \mathbb{R}$ such that

${⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\ge {\left(\varkappa \right(\parallel \mathbf{u}\parallel }_{\mathbf{E}}{)-\varkappa (\parallel \mathbf{v}\parallel }_{\mathbf{E}}{\left)\right)(\parallel \mathbf{u}\parallel }_{\mathbf{E}}-{\parallel \mathbf{v}\parallel }_{\mathbf{E}}),$

for any $\mathbf{v},\mathbf{u}\in \mathbf{E}$.

It is easy to see that this condition is less restrictive than the requirement of the strong monotonicity:

${⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\ge c{\parallel \mathbf{u}-\mathbf{v}\parallel }_{\mathbf{E}}^2,c=\mathrm{const},$

for any $\mathbf{u},\mathbf{v}\in \mathbf{E}$.

Proposition 1.Let $\mathbf{A}:\mathbf{E}\to {\mathbf{E}}^{*}$ be a continuous d-monotone operator with $\varkappa \left(s\right)\equiv c_{0}s$, where $c_0$ is a positive constant. Then:

• the operator $\mathbf{A}:\mathbf{E}\to {\mathbf{E}}^{*}$ is invertible;

• the operator ${\mathbf{A}}^{-1}:{\mathbf{E}}^{*}\to \mathbf{E}$ is continuous.

Proof. 

First, we show that the operator $\mathbf{A}$ is strictly monotone. Clearly, we have

(9)${⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\ge c_0{(\parallel \mathbf{u}\parallel }_{\mathbf{E}}-{\parallel \mathbf{v}\parallel }_{\mathbf{E}}{)}^2,$

for any $\mathbf{u},\mathbf{v}\in \mathbf{E}$. Taking into account this inequality and $c_0>0$, it is easy to see that

${⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\ge 0.$

Furthermore, if

${⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}=0,$

then we have

$\begin{array}{cc}\hfill 0=& {⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\mathbf{v}\right),\mathbf{u}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\hfill \\ \hfill =& 2{⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\frac{\mathbf{u}+\mathbf{v}}{2}\right),\frac{\mathbf{u}-\mathbf{v}}{2}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\hfill \\ \hfill & +2{⟨\mathbf{A}\left(\frac{\mathbf{u}+\mathbf{v}}{2}\right)-\mathbf{A}\left(\mathbf{v}\right),\frac{\mathbf{u}-\mathbf{v}}{2}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\hfill \\ \hfill =& 2{⟨\mathbf{A}\left(\mathbf{u}\right)-\mathbf{A}\left(\frac{\mathbf{u}+\mathbf{v}}{2}\right),\mathbf{u}-\frac{\mathbf{u}+\mathbf{v}}{2}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\hfill \\ \hfill & +2{⟨\mathbf{A}\left(\frac{\mathbf{u}+\mathbf{v}}{2}\right)-\mathbf{A}\left(\mathbf{v}\right),\frac{\mathbf{u}+\mathbf{v}}{2}-\mathbf{v}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\hfill \\ \hfill \ge & 2c_0{\left({\parallel \mathbf{u}\parallel }_{\mathbf{E}}-{∥\frac{\mathbf{u}+\mathbf{v}}{2}∥}_{\mathbf{E}}\right)}^2+2c_0{\left({∥\frac{\mathbf{u}+\mathbf{v}}{2}∥}_{\mathbf{E}}-{\parallel \mathbf{v}\parallel }_{\mathbf{E}}\right)}^2,\hfill \end{array}$

and hence

${∥\frac{\mathbf{u}+\mathbf{v}}{2}∥}_{\mathbf{E}}={\parallel \mathbf{u}\parallel }_{\mathbf{E}}={\parallel \mathbf{v}\parallel }_{\mathbf{E}}.$

In the Hilbert space $\mathbf{E}$, the last relations imply the equality $\mathbf{u}=\mathbf{v}$, which can be verified, for instance, by using the parallelogram law.

Next, we observe that the mapping $\mathbf{A}$ is coercive. Specifically, by setting $\mathbf{v}=\mathbf{0}$ in (9), we obtain

$\begin{array}{cc}\hfill {⟨\mathbf{A}\left(\mathbf{u}\right),\mathbf{u}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}\ge & {⟨\mathbf{A}\left(\mathbf{0}\right),\mathbf{u}⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}+c_0{\parallel \mathbf{u}\parallel }_{\mathbf{E}}^2\hfill \\ \hfill \ge & -{\parallel \mathbf{A}\left(\mathbf{0}\right)\parallel }_{{\mathbf{E}}^{*}}{\parallel \mathbf{u}\parallel }_{\mathbf{E}}+c_0{\parallel \mathbf{u}\parallel }_{\mathbf{E}}^2\hfill \\ \hfill =& {\xi (\parallel \mathbf{u}\parallel }_{\mathbf{E}}{)\parallel \mathbf{u}\parallel }_{\mathbf{E}},\hfill \end{array}$

with

$\xi \left(s\right):=-{\parallel \mathbf{A}\left(\mathbf{0}\right)\parallel }_{{\mathbf{E}}^{*}}+c_{0}s,\underset{s\to +\infty }{lim}\xi \left(s\right)=+\infty .$

Now, consider a sequence ${\left\{{\mathbf{u}}_n\right\}}_{n=1}^{\infty }\subset \mathbf{E}$ such that

${\mathbf{u}}_n\rightharpoonup {\mathbf{u}}_0\mathrm{in}\mathrm{the}\mathrm{space}\mathbf{E}\mathrm{as}n\to \infty$

and

(10)$\underset{n\to \infty }{lim}{⟨\mathbf{A}\left({\mathbf{u}}_n\right)-\mathbf{A}\left({\mathbf{u}}_0\right),{\mathbf{u}}_n-{\mathbf{u}}_0⟩}_{{\mathbf{E}}^{*}\times \mathbf{E}}=0.$

Due to inequality (9), relation (10) implies that

(11)$\underset{n\to \infty }{lim}\parallel {\mathbf{u}}_n{\parallel }_{\mathbf{E}}={\parallel {\mathbf{u}}_0\parallel }_{\mathbf{E}}.$

Since $\mathbf{E}$ is a Hilbert space, the weak convergence

${\mathbf{u}}_n\rightharpoonup {\mathbf{u}}_0\mathrm{in}\mathbf{E}\mathrm{as}n\to \infty$

together with (11) imply the strong convergence

${\mathbf{u}}_n\to {\mathbf{u}}_0\mathrm{in}\mathbf{E}\mathrm{as}n\to \infty .$

Thus, we deduce that the operator $\mathbf{A}$ satisfies the so-called $\left(S\right)$-property (see [46], Chapter III, § 1).

Taking into account the properties of the operator $\mathbf{A}$ established above, proving Proposition 1 is sufficient by applying Theorem 2.2 from [46], Chapter III, § 2. □

2.4. The Leray–Schauder Alternative for Completely Continuous Mappings

Proposition 2.

Let us suppose that $\mathbf{E}$ is a real Banach space, $\mathcal{B}[\mathbf{0},r]:=\left\{\mathbf{v}\in {\mathbf{E}:\parallel \mathbf{v}\parallel }_{\mathbf{E}}\le r\right\}$, and $\mathbf{K}:\mathcal{B}[\mathbf{0},r]\times [0,1]\to \mathbf{E}$ is a completely continuous mapping such that $\mathbf{K}(·,0):\mathcal{B}[\mathbf{0},r]\to \mathbf{E}$ is the zero operator. Then:

• (i). either the operator $\mathbf{K}(·,1):\mathcal{B}[\mathbf{0},r]\to \mathbf{E}$ has a fixed point in the ball $\mathcal{B}[\mathbf{0},r],$ or

• (ii). there exists a pair $({\mathbf{v}}_{*},{\lambda }_{*})\in \partial \mathcal{B}[\mathbf{0},r]\times [0,1]$ such that ${\mathbf{v}}_{*}=\mathbf{K}({\mathbf{v}}_{*},{\lambda }_{*})$.

For more details on fixed-point results and their applications in nonlinear analysis, we refer to the monographs [47,48].

2.5. Continuity of Superposition Operator in Lebesgue Spaces

Proposition 3

(Krasnoselskii’s theorem). Let $\mathcal{G}$ be a bounded domain in the real n-space ${\mathbb{R}}^n$ and let $Q:\mathcal{G}\times \mathbb{R}\to \mathbb{R}$ be a function such that

• there exist constants $p_1\ge 1,$$p_2\ge 1,$$b>0$, and a function $a\in L^{p_2}\left(\mathcal{G}\right)$ such that

  $\left|Q(\mathit{x},y)\right|\le a\left(\mathit{x}\right)+{b\left|y\right|}^{p_1/p_2},$

  for any $y\in \mathbb{R}$ and almost every $\mathit{x}\in \mathcal{G};$

• the function $Q(·,y):\mathcal{G}\to \mathbb{R}$ is measurable for any $y\in \mathbb{R};$

• the function $Q(\mathit{x},·):\mathbb{R}\to \mathbb{R}$ is continuous for almost every $\mathit{x}\in \mathcal{G}.$

Then the superposition operator ${\mathcal{N}}_Q$ defined by ${\mathcal{N}}_Q\left[w\right](·):=Q\left(·,w(·)\right)$ is continuous and bounded as a mapping from $L^{p_1}\left(\mathcal{G}\right)$ into $L^{p_2}\left(\mathcal{G}\right)$.

A detailed proof of this proposition can be found in the monograph [49], Chapter I, § 2.

3. Description of Assumptions on Model Data and Weak Formulation of Problem

We assume that the model under consideration is subject to the five conditions:

• (H.1). the function $\mu :{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ is continuous, and there exist constants ${\mu }_{\min }$ and ${\mu }_{\max }$ such that $0<{\mu }_{\min }\le \mu \left(\tau \right)\le {\mu }_{\max }$ for all $\tau \in {\mathbb{R}}_{+}$;

• (H.2). with a positive constant ${\mu }_{★}$, the inequality

  $\mu \left(t\right)t-\mu \left(s\right)s\ge {\mu }_{★}(t-s)$

  holds for any $t,s\in {\mathbb{R}}_{+}$ such that $t\ge s$;

• (H.3). the functions $f_1(·,y),f_2(·,y),f_3(·,y):\Omega \to \mathbb{R}$ and $g(·,y):\Omega \to \mathbb{R}$ are measurable for any $y\in \mathbb{R}$;

• (H.4). the functions $f_1(\mathit{x},·),f_2(\mathit{x},·),f_3(\mathit{x},·):\mathbb{R}\to \mathbb{R}$ and $g(\mathit{x},·):\mathbb{R}\to \mathbb{R}$ are continuous for almost every $\mathit{x}\in \Omega$;

• (H.5). there exist functions $f_{\max }\in L^2(\Omega )$ and $g_{\max }\in L^2(\Omega )$ such that

  $\left|\mathit{f}(\mathit{x},y)\right|\le f_{\max }\left(\mathit{x}\right),\left|g(\mathit{x},y)\right|\le g_{\max }\left(\mathit{x}\right)$

  for any $y\in \mathbb{R}$ and almost every $\mathit{x}\in \Omega$.

It can easily be checked that condition (H.2) holds for any $C^1$-smooth function $\mu :{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ satisfying

$\mu \left(\tau \right)\ge {\mu }_{\min }>0,{\mu }^{\prime }\left(\tau \right)\ge 0,\forall \tau \in {\mathbb{R}}_{+},$

where ${\mu }^{\prime }:=d\mu /d\tau$. Indeed, by the Mean Value Theorem (Lagrange’s theorem) and condition (H.1), we derive

$\begin{array}{cc}\hfill \mu \left(t\right)t-\mu \left(s\right)s& =\left[{\mu }^{\prime }\left(\zeta (t,s)\right)\zeta (t,s)+\mu \left(\zeta (t,s)\right)\right](t-s)\hfill \\ \hfill & \ge {\mu }_{\min }(t-s),\hfill \end{array}$

for arbitrary $t,s\in {\mathbb{R}}_{+}$, $t>s$, and some $\zeta (t,s)\in [s,t]$.

Note that the requirements (H.1) and (H.2) imposed on the viscosity function $\mu$ are less restrictive compared to the conditions used in [23] (see page 53) for the variational formulation of the problem concerning isothermal flows of non-Newtonian fluids. Notably, a Newtonian fluid with constant viscosity $\mu \equiv {\mu }_{\min }$ satisfies (H.1) and (H.2) with ${\mu }_{★}={\mu }_{\min }$.

Now, let us provide the weak formulation of problem (1).

Definition 2.

We shall say that a pair $(\mathit{u},\theta )$ is a weak solution to problem (1) if the following two conditions hold:

• the pair $(\mathit{u},\theta )$ belongs to the space $\mathit{V}(\Omega )\times H^1(\Omega );$

• the equalities

  $\begin{array}{cc}\hfill & -\sum _{i=1}^3{\int }_{\Omega }{u}_i\mathit{u}·\frac{\partial \mathit{v}}{\partial x_i}\mathit{dx}+{\int }_{\Omega }\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right):\mathcal{D}\left(\mathit{v}\right)\mathit{dx}\hfill \\ \hfill & +\alpha {\int }_{\Omega }\mathit{u}·\mathit{v}\mathit{dx}={\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{v}\mathit{dx},\hfill \end{array}$

  $\begin{array}{cc}\hfill & -\sum _{i=1}^3{\int }_{\Omega }{u}_i\theta \frac{\partial \eta }{\partial x_i}\mathit{dx}+k{\int }_{\Omega }\nabla \theta ·\nabla \eta \mathit{dx}+\beta {\int }_{\partial \Omega }\theta \eta dS\hfill \\ \hfill & ={\displaystyle \frac{1}{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\eta \mathit{dx}+{\int }_{\Omega }g(\mathit{x},\theta )\eta \mathit{dx}\hfill \end{array}$

  are valid for all test functions $\mathit{v}\in \mathit{V}(\Omega )$ and $\eta \in H^1(\Omega )$.

By $\mathfrak{M}$ we denote the set of all weak solutions to problem (1).

Remark 1.Having a weak solution $(\mathit{u},\theta )$ in hand, the corresponding pressure π can be obtained by the well-known de Rham theorem (see, for example, [50], Chapter I, § 1). The triplet $(\mathit{u},\theta ,\pi )$ is called a full weak solution.

4. Main Results and Their Proof

The main results of this article are formulated in the following theorem.

Theorem 1.Suppose Ω is a bounded Lipschitz domain in ${\mathbb{R}}^3$ and conditions(H.1)–(H.5)are satisfied. Then:

• (a). problem (1) has at least one weak solution, that is, $\mathfrak{M}\ne ⌀;$

• (b). any pair $(\mathit{u},\theta )\in \mathfrak{M}$ satisfies the energy equalities:

  (12)${\int }_{\Omega }{\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2\mathit{dx}+\alpha {\int }_{\Omega }{\left|\mathit{u}\right|}^2\mathit{dx}={\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{u}\mathit{dx},$

  (13)$\begin{array}{cc}\hfill & k{\int }_{\Omega }{|\nabla \theta |}^2\mathit{dx}+\beta {\int }_{\partial \Omega }{\left|\theta \right|}^{2}dS\hfill \\ \hfill & ={\displaystyle \frac{1}{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\theta \mathit{dx}+{\int }_{\Omega }g(\mathit{x},\theta )\theta \mathit{dx};\hfill \end{array}$

• (c). the set $\mathfrak{M}$ is compact in the space $\mathit{V}(\Omega )\times H^1(\Omega )$ as well as in the space ${\mathit{L}}^p(\Omega )\times L^p(\Omega )$, where the exponent p can be chosen arbitrarily from the closed interval $[1,6]$.

Proof. 

We perform the proof of statement (a) in five steps.

Step 1: Auxiliary problem. Let us consider a one-parameter problem:

Find a pair of functions $(\mathit{u},\theta )$ from the space $\mathit{V}(\Omega )\times H^1(\Omega )$ that satisfy the following two equations:

(14)$\begin{array}{cc}\hfill & -\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\mathit{u}·\frac{\partial \mathit{v}}{\partial x_i}\mathit{dx}+{\int }_{\Omega }\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right):\mathcal{D}\left(\mathit{v}\right)\mathit{dx}\hfill \\ \hfill & +\lambda \alpha {\int }_{\Omega }\mathit{u}·\mathit{v}\mathit{dx}=\lambda {\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{v}\mathit{dx},\forall \mathit{v}\in \mathit{V}(\Omega ),\hfill \end{array}$

(15)$\begin{array}{cc}\hfill & -\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\theta \frac{\partial \eta }{\partial x_i}\mathit{dx}+k{\int }_{\Omega }\nabla \theta ·\nabla \eta \mathit{dx}+\beta {\int }_{\partial \Omega }\theta \eta dS\hfill \\ \hfill & ={\displaystyle \frac{\lambda }{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\eta \mathit{dx}+\lambda {\int }_{\Omega }g(\mathit{x},\theta )\eta \mathit{dx},\forall \eta \in H^1(\Omega ),\hfill \end{array}$

where $\lambda$ is a numerical parameter belonging to the interval $[0,1]$.

Step 2: A priori bounds. Let us derive a priori bounds for solutions of the above problem. Suppose that a pair $(\mathit{u},\theta )\in \mathit{V}(\Omega )\times H^1(\Omega )$ satisfies (14) and (15) for some $\lambda \in [0,1]$.

By setting $\mathit{v}=\mathit{u}$ in equality (14), we obtain

(16)$\begin{array}{cc}\hfill & -\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\mathit{u}·\frac{\partial \mathit{u}}{\partial x_i}\mathit{dx}+{\int }_{\Omega }{\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2\mathit{dx}\hfill \\ \hfill & +\lambda \alpha {\int }_{\Omega }{\left|\mathit{u}\right|}^2\mathit{dx}=\lambda {\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{u}\mathit{dx}.\hfill \end{array}$

The first term on the left-hand side of (16) is zero. Indeed, applying integration by parts, we obtain

(17)$\begin{array}{cc}\hfill \sum _{i=1}^3{\int }_{\Omega }{u}_i\mathit{u}·\frac{\partial \mathit{u}}{\partial x_i}\mathit{dx}& =\frac{1}{2}\sum _{i=1}^3{\int }_{\Omega }{u}_i\frac{{\partial \left|\mathit{u}\right|}^2}{\partial x_i}\mathit{dx}\hfill \\ \hfill & =-\frac{1}{2}\sum _{i=1}^3{\int }_{\Omega }\frac{\partial u_i}{\partial x_i}{\left|\mathit{u}\right|}^2\mathit{dx}\hfill \\ \hfill & =-\frac{1}{2}\sum _{i=1}^3{\int }_{\Omega }\left(\mathrm{div}\mathit{u}\right){\left|\mathit{u}\right|}^2\mathit{dx}=0.\hfill \end{array}$

Of course, here we used the two identities: $\mathrm{div}\mathit{u}\equiv 0$ in $\Omega$ and $\mathit{u}\equiv 0$ on $\partial \Omega$.

In view of (17), equality (16) reduces to

(18)${\int }_{\Omega }{\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2\mathit{dx}+\lambda \alpha {\int }_{\Omega }{\left|\mathit{u}\right|}^2\mathit{dx}=\lambda {\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{u}\mathit{dx}.$

From this equality, using assumptions (H.1) and (H.5) as well as the Cauchy–Schwarz inequality, we deduce that

(19)$\begin{array}{cc}\hfill {\mu }_{\min }{\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )}^2& \le {\int }_{\Omega }{\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2\mathit{dx}\le \lambda {\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{u}\mathit{dx}\le {\int }_{\Omega }{f}_{\max }\left|\mathit{u}\right|\mathit{dx}\hfill \\ \hfill & \le {\left({\int }_{\Omega }{\left|f_{\max }\right|}^2\mathit{dx}\right)}^{\frac{1}{2}}{\left({\int }_{\Omega }{\left|\mathit{u}\right|}^2\mathit{dx}\right)}^{\frac{1}{2}}\hfill \\ \hfill & =\parallel f_{\max }{\parallel }_{L^2(\Omega )}{\parallel \mathit{u}\parallel }_{{\mathit{L}}^2(\Omega )}.\hfill \end{array}$

Note that

(20)${\parallel \mathit{u}\parallel }_{{\mathit{L}}^2(\Omega )}\le {\parallel \mathbb{I}\parallel }_{\mathcal{L}(\mathit{V}(\Omega ),{\mathit{L}}^2(\Omega ))}{\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )}.$

Here and in the succeeding discussion, the symbol $\mathbb{I}$ denotes the identity operator.

Taking into account (20), we can continue (19) as follows:

${\mu }_{\min }{\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )}^2\le {\parallel \mathbb{I}\parallel }_{\mathcal{L}(\mathit{V}(\Omega ),{\mathit{L}}^2(\Omega ))}\parallel f_{\max }{\parallel }_{L^2(\Omega )}{\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )},$

whence

(21)${\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )}\le \frac{1}{{\mu }_{\min }}{\parallel \mathbb{I}\parallel }_{\mathcal{L}(\mathit{V}(\Omega ),{\mathit{L}}^2(\Omega ))}{\parallel f_{\max }\parallel }_{L^2(\Omega )}=:\phi (\Omega ,{\mu }_{\min },f_{\max }).$

Next, on setting $\eta =\theta$ in (15), we find

(22)$\begin{array}{cc}\hfill & -\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\theta \frac{\partial \theta }{\partial x_i}\mathit{dx}+k{\int }_{\Omega }{|\nabla \theta |}^2\mathit{dx}+\beta {\int }_{\partial \Omega }{\left|\theta \right|}^{2}dS\hfill \\ \hfill & ={\displaystyle \frac{\lambda }{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\theta \mathit{dx}+\lambda {\int }_{\Omega }g(\mathit{x},\theta )\theta \mathit{dx}.\hfill \end{array}$

The first term on the left-hand side of the last equality is zero. Indeed, using the integration by parts formula and taking into account the identities $\mathrm{div}\mathit{u}\equiv 0$ in $\Omega$ and $\mathit{u}\equiv 0$ on $\partial \Omega$, we obtain

$\begin{array}{cc}\hfill \sum _{i=1}^3{\int }_{\Omega }{u}_i\theta \frac{\partial \theta }{\partial x_i}\mathit{dx}& =\frac{1}{2}\sum _{i=1}^3{\int }_{\Omega }{u}_i\frac{{\partial \left|\theta \right|}^2}{\partial x_i}\mathit{dx}\hfill \\ \hfill & =-\frac{1}{2}\sum _{i=1}^3{\int }_{\Omega }\frac{\partial u_i}{\partial x_i}{\left|\theta \right|}^2\mathit{dx}\hfill \\ \hfill & =-\frac{1}{2}{\int }_{\Omega }\left(\mathrm{div}\mathit{u}\right){\left|\theta \right|}^2\mathit{dx}=0.\hfill \end{array}$

Therefore, equality (22) simplifies to

(23)$k{\int }_{\Omega }{|\nabla \theta |}^2\mathit{dx}+\beta {\int }_{\partial \Omega }{\left|\theta \right|}^{2}dS={\displaystyle \frac{\lambda }{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\theta \mathit{dx}+\lambda {\int }_{\Omega }g(\mathit{x},\theta )\theta \mathit{dx}.$

From (23), using assumption (H.5) and the Cauchy–Schwarz inequality, we derive

(24)$\begin{array}{cc}\hfill {min\{k,\beta \}\parallel \theta \parallel }_{H^1(\Omega )}^2& \le k{\int }_{\Omega }{|\nabla \theta |}^2\mathit{dx}+\beta {\int }_{\partial \Omega }{\left|\theta \right|}^{2}dS\hfill \\ \hfill & ={\displaystyle \frac{\lambda }{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\theta \mathit{dx}+\lambda {\int }_{\Omega }g(\mathit{x},\theta )\theta \mathit{dx}\hfill \\ \hfill & \le {\displaystyle \frac{{\mu }_{\max }}{c_p}}\parallel {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}^2{\int }_{\Omega }\left|\theta \right|\mathit{dx}+{\int }_{\Omega }|g_{\max }\left|\right|\theta |\mathit{dx}\hfill \\ \hfill & \le \left({\displaystyle \frac{{\mu }_{\max }}{c_p}}\parallel {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}^2\sqrt{\mathrm{meas}(\Omega )}+{\parallel g_{\max }\parallel }_{L^2(\Omega )}\right){\parallel \theta \parallel }_{L^2(\Omega )}.\hfill \end{array}$

Notice that

$\begin{array}{cc}\hfill {\parallel \theta \parallel }_{L^2(\Omega )}& \le {\parallel \mathbb{I}\parallel }_{\mathcal{L}(H^1(\Omega ),L^2(\Omega ))}{\parallel \theta \parallel }_{H^1(\Omega )}.\hfill \end{array}$

Therefore, (24) can be continued as follows:

$\begin{array}{cc}\hfill {min\{k,\beta \}\parallel \theta \parallel }_{H^1(\Omega )}^2& \le \left({\displaystyle \frac{{\mu }_{\max }}{c_p}}\parallel {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}^2\sqrt{\mathrm{meas}(\Omega )}+{\parallel g_{\max }\parallel }_{L^2(\Omega )}\right)\hfill \\ \hfill & \times {\parallel \mathbb{I}\parallel }_{\mathcal{L}(H^1(\Omega ),L^2(\Omega ))}{\parallel \theta \parallel }_{H^1(\Omega )},\hfill \end{array}$

whence

(25)$\begin{array}{cc}\hfill {\parallel \theta \parallel }_{H^1(\Omega )}& \le \left({\displaystyle \frac{{\mu }_{\max }}{c_p}}\parallel {\mathcal{D}}_{\rho }{\left(\mathit{u}\right)\parallel }_{\mathit{C}\left(\overline{\Omega }\right)}^2\sqrt{\mathrm{meas}(\Omega )}+{\parallel g_{\max }\parallel }_{L^2(\Omega )}\right)\frac{{\parallel \mathbb{I}\parallel }_{\mathcal{L}(H^1(\Omega ),L^2(\Omega ))}}{min\{k,\beta \}}.\hfill \end{array}$

Taking into account inequalities (3), (20), and (21), from (25) we can easily derive the following bound for $\theta$:

(26)$\begin{array}{cc}\hfill {\parallel \theta \parallel }_{H^1(\Omega )}& \le ({\displaystyle \frac{{\mu }_{\max }}{c_p}}{M}^2(\rho ,\Omega ){\parallel \mathbb{I}\parallel }_{\mathcal{L}(\mathit{V}(\Omega ),{\mathit{L}}^2(\Omega ))}^2{\phi }^2(\Omega ,{\mu }_{\min },f_{\max })\sqrt{\mathrm{meas}(\Omega )}\hfill \\ \hfill & +\parallel g_{\max }{\parallel }_{L^2(\Omega )})\frac{{\parallel \mathbb{I}\parallel }_{\mathcal{L}(H^1(\Omega ),L^2(\Omega ))}}{min\{k,\beta \}}\hfill \\ \hfill & =:\psi (\Omega ,{\mu }_{\min },{\mu }_{\max },f_{\max },g_{\max },c_p,\rho ,k,\beta ).\hfill \end{array}$

Step 3: Operator formulation. Now, let us proceed with an operator interpretation of system (14)–(15). To achieve this, we introduce the eight operators:

$\begin{array}{cc}\hfill & {\mathbb{A}}_1:\mathit{V}(\Omega )\to {\left[\mathit{V}(\Omega )\right]}^{*},\hfill \\ \hfill & {⟨{\mathbb{A}}_1\left(\mathit{u}\right),\mathit{v}⟩}_{{\left[\mathit{V}(\Omega )\right]}^{*}\times \mathit{V}(\Omega )}:={\int }_{\Omega }\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right):\mathcal{D}\left(\mathit{v}\right)\mathit{dx};\hfill \end{array}$

$\begin{array}{cc}\hfill & {\mathbb{A}}_2:H^1(\Omega )\to {\left[H^1(\Omega )\right]}^{*},\hfill \\ \hfill & {⟨{\mathbb{A}}_2\left(\theta \right),\eta ⟩}_{{\left[H^1(\Omega )\right]}^{*}\times H^1(\Omega )}:=k{\int }_{\Omega }\nabla \theta ·\nabla \eta \mathit{dx}+\beta {\int }_{\partial \Omega }\theta \eta dS;\hfill \end{array}$

$\begin{array}{cc}\hfill & \mathbb{A}:\mathit{V}(\Omega )\times H^1(\Omega )\to {\left[\mathit{V}(\Omega )\right]}^{*}\times {\left[H^1(\Omega )\right]}^{*},\hfill \\ \hfill & \mathbb{A}(\mathit{u},\theta ):=\left({\mathbb{A}}_1\left(\mathit{u}\right),{\mathbb{A}}_2\left(\theta \right)\right);\hfill \end{array}$

$\begin{array}{cc}\hfill & {\tilde{\mathbb{K}}}_1:{\mathit{L}}^4(\Omega )\times L^4(\Omega )\times [0,1]\to {\left[\mathit{V}(\Omega )\right]}^{*},\hfill \\ \hfill & {⟨{\tilde{\mathbb{K}}}_1(\mathit{u},\theta ,\lambda ),\mathit{v}⟩}_{{\left[\mathit{V}(\Omega )\right]}^{*}\times \mathit{V}(\Omega )}:=\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\mathit{u}·\frac{\partial \mathit{v}}{\partial x_i}\mathit{dx}-\lambda \alpha {\int }_{\Omega }\mathit{u}·\mathit{v}\mathit{dx}+\lambda {\int }_{\Omega }\mathit{f}(\mathit{x},\theta )·\mathit{v}\mathit{dx};\hfill \end{array}$

$\begin{array}{cc}\hfill & {\tilde{\mathbb{K}}}_2:{\mathit{L}}^4(\Omega )\times L^4(\Omega )\times [0,1]\to {\left[H^1(\Omega )\right]}^{*},\hfill \\ \hfill & {⟨{\tilde{\mathbb{K}}}_2(\mathit{u},\theta ,\lambda ),\eta ⟩}_{{\left[H^1(\Omega )\right]}^{*}\times H^1(\Omega )}:=\lambda \sum _{i=1}^3{\int }_{\Omega }{u}_i\theta \frac{\partial \eta }{\partial x_i}\mathit{dx}+{\displaystyle \frac{\lambda }{c_p}}{\int }_{\Omega }\mu \left(\right|{\mathcal{D}}_{\rho }\left(\mathit{u}\right)\left|\right)|{\mathcal{D}}_{\rho }{\left(\mathit{u}\right)|}^2\eta \mathit{dx}\hfill \\ \hfill & +\lambda {\int }_{\Omega }g(\mathit{x},\theta )\eta \mathit{dx};\hfill \end{array}$

$\begin{array}{cc}\hfill & \tilde{\mathbb{K}}:{\mathit{L}}^4(\Omega )\times L^4(\Omega )\times [0,1]\to {\left[\mathit{V}(\Omega )\right]}^{*}\times {\left[H^1(\Omega )\right]}^{*},\hfill \\ \hfill & \tilde{\mathbb{K}}(\mathit{u},\theta ,\lambda ):=\left({\tilde{\mathbb{K}}}_1(\mathit{u},\theta ,\lambda ),{\tilde{\mathbb{K}}}_2(\mathit{u},\theta ,\lambda )\right);\hfill \end{array}$

$\begin{array}{cc}\hfill & \mathbb{J}:\mathit{V}(\Omega )\times H^1(\Omega )\times [0,1]\to {\mathit{L}}^4(\Omega )\times L^4(\Omega )\times [0,1],\hfill \\ \hfill & \mathbb{J}(\mathit{u},\theta ,\lambda ):=(\mathit{u},\theta ,\lambda );\hfill \end{array}$

$\begin{array}{cc}\hfill & \mathbb{K}:\mathit{V}(\Omega )\times H^1(\Omega )\times [0,1]\to {\left[\mathit{V}(\Omega )\right]}^{*}\times {\left[H^1(\Omega )\right]}^{*},\hfill \\ \hfill & \mathbb{K}(\mathit{u},\theta ,\lambda ):=(\tilde{\mathbb{K}}\circ \mathbb{J})(\mathit{u},\theta ,\lambda ).\hfill \end{array}$

It is easy to see that system (14)–(15) is equivalent to the operator equation

(27)$\mathbb{A}(\mathit{u},\theta )=\mathbb{K}(\mathit{u},\theta ,\lambda ).$

Step 4: Properties of operators. Let us establish some properties of the operators $\mathbb{A}$ and $\mathbb{K}$.

First, using assumption (H.2), we obtain

$\begin{array}{cc}\hfill (\mu & \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right)-\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\mathcal{D}\left(\mathit{w}\right)):\left(\mathcal{D}\left(\mathit{u}\right)-\mathcal{D}\left(\mathit{w}\right)\right)\hfill \\ \hfill & ={\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2-\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right):\mathcal{D}\left(\mathit{w}\right)\hfill \\ \hfill & -\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\mathcal{D}\left(\mathit{u}\right):\mathcal{D}\left(\mathit{w}\right)+{\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\left|\mathcal{D}\left(\mathit{w}\right)\right|}^2\hfill \\ \hfill & \ge {\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|}^2-\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|\left|\mathcal{D}\left(\mathit{w}\right)\right|\hfill \\ \hfill & -\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\left|\mathcal{D}\left(\mathit{u}\right)\right|\left|\mathcal{D}\left(\mathit{w}\right)\right|+{\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\left|\mathcal{D}\left(\mathit{w}\right)\right|}^2\hfill \\ \hfill & =\left\{\mu \left(\right|\mathcal{D}\left(\mathit{u}\right)\left|\right)\left|\mathcal{D}\right(\mathit{u}\left)\right|-\mu \left(\right|\mathcal{D}\left(\mathit{w}\right)\left|\right)\left|\mathcal{D}\right(\mathit{w}\left)\right|\right\}\left(\right|\mathcal{D}\left(\mathit{u}\right)|-\left|\mathcal{D}\left(\mathit{w}\right)\right|)\hfill \\ \hfill & \ge {\mu }_{★}\left(\right|\mathcal{D}\left(\mathit{u}\right)|-{\left|\mathcal{D}\left(\mathit{w}\right)\right|)}^2,\hfill \end{array}$

for any $\mathit{u},\mathit{w}\in \mathit{V}(\Omega )$. From this, we deduce the inequality

${⟨{\mathbb{A}}_1\left(\mathit{u}\right)-{\mathbb{A}}_1\left(\mathit{w}\right),\mathit{u}-\mathit{w}⟩}_{{\left[\mathit{V}(\Omega )\right]}^{*}\times \mathit{V}(\Omega )}\ge {\mu }_{★}{\left({\parallel \mathit{u}\parallel }_{\mathit{V}(\Omega )}-{\parallel \mathit{w}\parallel }_{\mathit{V}(\Omega )}\right)}^2.$

This means that the operator ${\mathbb{A}}_1$ is d-monotone with $\varkappa \left(s\right):={\mu }_{★}s$, for $s\in {\mathbb{R}}_{+}$. Applying Proposition 1, we obtain that the operator ${\mathbb{A}}_1$ is continuously invertible.

The same holds true for the operator ${\mathbb{A}}_2$. In order to prove the continuous invertibility of this operator, we observe that

${⟨{\mathbb{A}}_2\left(\theta \right),\theta ⟩}_{{\left[H^1(\Omega )\right]}^{*}\times H^1(\Omega )}\ge min\{k,\beta \}{\parallel \theta \parallel }_{H^1(\Omega )}^2$

and apply the classical Lax–Milgram lemma (see [45], Chapter I, § 1.3.1).

Thus, we see that the operator $\mathbb{A}$ is continuously invertible, and the following relation

${\mathbb{A}}^{-1}({\mathit{h}}_1,{\mathit{h}}_2)=\left({\mathbb{A}}_1^{-1}\left({\mathit{h}}_1\right),{\mathbb{A}}_2^{-1}\left({\mathit{h}}_2\right)\right)$

holds for any $({\mathit{h}}_1,{\mathit{h}}_2)\in {\left[\mathit{V}(\Omega )\right]}^{*}\times {\left[H^1(\Omega )\right]}^{*}$.

Next, note that, due to compact embedding results concerning Sobolev spaces (see [45], Chapter II, § 2.6.1), we have

$H^1(\Omega )\hookrightarrow \hookrightarrow L^p(\Omega ),\forall p\in [1,6).$

Therefore, the operator $\mathbb{J}$ is completely continuous. Moreover, by conditions (H.3)–(H.5), Proposition 3 and Lemma 2, it can be directly checked that the operator $\tilde{\mathbb{K}}$ is continuous. Thus, we see that the operator $\mathbb{K}=\tilde{\mathbb{K}}\circ \mathbb{J}$ (and consequently, the operator ${\mathbb{A}}^{-1}\circ \mathbb{K}$) is completely continuous.

Step 5: Solvability. Applying the operator ${\mathbb{A}}^{-1}$ to both sides of Equation (27), we obtain

$(\mathit{u},\theta )=({\mathbb{A}}^{-1}\circ \mathbb{K})(\mathit{u},\theta ,\lambda ).$

Taking into account the previously derived a priori bounds (21) and (26), which hold for any solutions of (27), we apply Proposition 2 to conclude that the equation

$(\mathit{u},\theta )=({\mathbb{A}}^{-1}\circ \mathbb{K})(\mathit{u},\theta ,1)$

has a solution $({\mathit{u}}_0,{\theta }_0)$ in a ball

$\mathit{B}[(\mathbf{0},0),r]:=\{(\mathit{v},\omega ):(\mathit{v},\omega )\in \mathit{V}(\Omega )\times H^1(\Omega )and\parallel (\mathit{v},\omega ){\parallel }_{\mathit{V}(\Omega )\times H^1(\Omega )}<r\}$

with sufficiently large radius r.

For example, we can take $r=r_{★}$, where

$r_{★}:=max\left\{\phi (\Omega ,{\mu }_{\min },f_{\max }),\psi (\Omega ,{\mu }_{\min },{\mu }_{\max },f_{\max },g_{\max },c_p,\rho ,k,\beta )\right\}+1.$

It is evident that

$\mathbb{A}({\mathit{u}}_0,{\theta }_0)=\mathbb{K}({\mathit{u}}_0,{\theta }_0,1)$

and the pair of functions $({\mathit{u}}_0,{\theta }_0)$ is a weak solution to problem (1). Thus, we establish that $\mathfrak{M}\ne ⌀$, and hence part (a) is proved.

In order to establish (b), we just substitute $\lambda =1$ into (18) and (23) and obtain the required energy equalities (12) and (13).

Finally, let us prove statement (c). It is readily apparent from estimates (21) and (26) that set $\mathfrak{M}$ is bounded in the Cartesian product $\mathit{V}(\Omega )\times H^1(\Omega )$. Moreover, this set coincides with the fixed-points set of the completely continuous operator $({\mathbb{A}}^{-1}\circ \mathbb{K})(·,·,1)$. Taking into account these facts, it is easy to show the compactness of $\mathfrak{M}$ in $\mathit{V}(\Omega )\times H^1(\Omega )$.

Further, using the Sobolev embedding theorem (see, for example, [45], Chapter II, § 2.3.4), we obtain

$\mathit{V}(\Omega )\hookrightarrow {\mathit{L}}^p(\Omega ),H^1(\Omega )\hookrightarrow L^p(\Omega ),\forall p\in [1,6].$

Therefore, the set $\mathfrak{M}$ belongs to the space ${\mathit{L}}^p(\Omega )\times L^p(\Omega )$ with $p\in [1,6]$ and is compact in this space, since a continuous image of a compact set is compact. Thus, we have completed the proof of Theorem 1. □

Remark 2.

In the particular case when $\mu =\mathrm{const}$ and $\alpha =0$, the results of Theorem 1 coincide with the results from the work [17] by the first author.

5. Conclusions

In this paper, we have proved the existence of stationary weak solutions to the generalized Boussinesq system with the spatially averaged Rayleigh function in the energy balance equation. Unlike conventional approaches that overlook the viscous dissipation effect, our model incorporates this factor, departing from the classical Boussinesq approximation. Notably, our existence theorem does not necessitate smallness model data assumptions. The proof of this theorem is based on an operator interpretation of the boundary value problem under consideration and applying the Leray–Schauder alternative. The key points in our proof strategy are deriving suitable a priori estimates for weak solutions and establishing the strong continuity of the used averaging operation and the continuous invertibility of d-monotone operators. As a result, we overcame the mathematical challenges associated with solving the corresponding strongly nonlinear PDEs. Furthermore, our analysis revealed that the set of weak solutions is compact in the Cartesian product of Sobolev spaces, and each solution from this set satisfies some energy equalities. The approach proposed in the present work advances the insight of non-Newtonian fluid dynamics and provides new ways for exploring intricate heat and mass transfer models with energy dissipation. The authors suggest the following directions for future research:

• existence and uniqueness of other types of solutions (strong, classical, etc.);

• continuous dependence of solutions on model data;

• global solvability of time-dependent problems;

• stability/instability issues;

• flow control and optimization problems.

Author Contributions

Conceptualization, E.S.B.; investigation, E.S.B. and O.Y.S.; writing—original draft preparation, E.S.B. and O.Y.S.; writing—review and editing, O.Y.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

References

1. Oberbeck, A. Über die Wärmeleitung der Flüssigkeiten bei der Berücksichtigung der Strömungen infolge von Temperaturdifferenzen. Ann. Phys. Chem.; 1879; 7, pp. 271-292. [DOI: https://dx.doi.org/10.1002/andp.18792430606]

2. Boussinesq, J. Théorie Analytique de la Chaleur; Gauthier-Villars: Paris, France, 1903.

3. Shinbrot, M.; Kotorynski, W.P. The initial value problem for a viscous heat-conducting fluid. J. Math. Anal. Appl.; 1974; 45, pp. 1-22. [DOI: https://dx.doi.org/10.1016/0022-247X(74)90115-2]

4. Alekseev, G.V.; Smishlyaev, A.B.; Tereshko, D.A. The solvability of a boundary value problem for time-independent equations of heat and mass transfer under mixed boundary conditions. Comput. Math. Math. Phys.; 2003; 43, pp. 63-77.

5. Hmidi, T.; Rousset, F. Global well-posedness for the Navier–Stokes–Boussinesq system with axisymmetric data. Ann. Henri Poincare; 2010; 27, pp. 1227-1246. [DOI: https://dx.doi.org/10.1016/j.anihpc.2010.06.001]

6. Liu, X. Global existence and uniqueness of solutions to the three-dimensional Boussinesq equations. Bound. Value Probl.; 2016; 2016, 85. [DOI: https://dx.doi.org/10.1186/s13661-016-0592-9]

7. Jiu, Q.; Yu, H. Global well-posedness for 3D generalized Navier–Stokes–Boussinesq equations. Acta Math. Appl. Sin. Engl. Ser.; 2016; 32, pp. 1-16. [DOI: https://dx.doi.org/10.1007/s10255-016-0539-z]

8. Khor, C.; Xu, X. Temperature patches for the subcritical Boussinesq–Navier–Stokes system with no diffusion. J. Funct. Anal.; 2022; 283, 109501. [DOI: https://dx.doi.org/10.1016/j.jfa.2022.109501]

9. Ershkov, S.; Burmasheva, N.; Leshchenko, D.D.; Prosviryakov, E.Y. Exact solutions of the Oberbeck–Bussinesk equations for the description of shear thermal diffusion of Newtonian fluid flows. Symmetry; 2023; 15, 1730. [DOI: https://dx.doi.org/10.3390/sym15091730]

10. Ershkov, S.V.; Prosviryakov, E.Y.; Burmasheva, N.V.; Christianto, V. Solving the hydrodynamical system of equations of inhomogeneous fluid flows with thermal diffusion: A review. Symmetry; 2023; 15, 1825. [DOI: https://dx.doi.org/10.3390/sym15101825]

11. Kagei, Y.; Růžička, M.; Thäter, G. Natural convection with dissipative heating. Commun. Math. Phys.; 2000; 214, pp. 287-313. [DOI: https://dx.doi.org/10.1007/s002200000275]

12. Palani, G.; Kim, K.Y. Viscous dissipation effects on heat transfer in flow over an inclined plate. J. Appl. Mech. Tech. Phys.; 2010; 51, pp. 241-248. [DOI: https://dx.doi.org/10.1007/s10808-010-0034-x]

13. Moslemi, M.; Javaherdeh, K. Viscous dissipation effect in the free convection of non-Newtonian fluid with heat generation or absorption effect on the vertical wavy surface. J. Appl. Math.; 2021; 2021, 7567981. [DOI: https://dx.doi.org/10.1155/2021/7567981]

14. Goruleva, L.S.; Prosviryakov, E.Y. A new class of exact solutions to the Navier–Stokes equations with allowance for internal heat release. Opt. Spectrosc.; 2022; 130, pp. 365-370. [DOI: https://dx.doi.org/10.1134/S0030400X22070037]

15. Privalova, V.V.; Prosviryakov, E.Y. A new class of exact solutions of the Oberbeck–Boussinesq equations describing an incompressible fluid. Theor. Found. Chem. Eng.; 2022; 56, pp. 331-338. [DOI: https://dx.doi.org/10.1134/S0040579522030113]

16. Baranovskii, E.S. Exact solutions for non-isothermal flows of second grade fluid between parallel plates. Nanomaterials; 2023; 13, 1409. [DOI: https://dx.doi.org/10.3390/nano13081409]

17. Baranovskii, E.S. The stationary Navier–Stokes–Boussinesq system with a regularized dissipation function. Math. Notes; 2024; 115, in press

18. Vorotnikov, D.A. An objective model of viscoelastic fluid: Solvability of motion equations and attractors. Proceedings of the Fluid DTU Summer School on Complex Motion in Fluids, Krogerup Hojskole; Copenhagen, Denmark, 19–25 August 2007; 23.

19. Sobolev, S.L. Some Applications of Functional Analysis in Mathematical Physics; 3rd ed. AMS: Providence, RI, USA, 1991.

20. Ladyzhenskaya, O.A. The Mathematical Theory of Viscous Incompressible Flow; Gordon and Breach: New York, NY, USA, 1969.

21. Baranovskii, E.S.; Artemov, M.A. Global existence results for Oldroyd fluids with wall slip. Acta Appl. Math.; 2017; 147, pp. 197-210. [DOI: https://dx.doi.org/10.1007/s10440-016-0076-z]

22. Baranovskii, E.S. Steady flows of an Oldroyd fluid with threshold slip. Commun. Pure Appl. Anal.; 2019; 18, pp. 735-750. [DOI: https://dx.doi.org/10.3934/cpaa.2019036]

23. Litvinov, W.G. Motion of Nonlinear-Viscous Fluid; Nauka: Moscow, Russia, 1982.

24. Chhabra, R.P.; Richardson, J.F. Non-Newtonian Flow and Applied Rheology; 2nd ed. Butterworth-Heinemann: Oxford, UK, 2008; [DOI: https://dx.doi.org/10.1016/B978-0-7506-8532-0.X0001-7]

25. Domnich, A.A.; Artemov, M.A.; Shishkina, O.Y. On the boundary value problem for a model of nonisothermal flows of a non-Newtonian fluid. J. Appl. Ind. Math.; 2020; 14, pp. 37-45. [DOI: https://dx.doi.org/10.1134/S1990478920010056]

26. Ladyzhenskaya, O.A. New equations for the description of the motions of viscous incompressible fluids, and global solvability for their boundary value problems. Proc. Steklov Inst. Math.; 1967; 102, pp. 95-118.

27. Ladyzhenskaya, O.A. On some nonlinear problems in the theory of continuous media. Amer. Math. Soc. Transl. Ser. 2; 1968; 70, pp. 73-88. [DOI: https://dx.doi.org/10.1090/trans2/070]

28. Kuzmin, M.Y. A mathematical model of the motion of a nonlinear viscous fluid with the condition of slip on the boundary. Russ. Math.; 2007; 51, pp. 51-60. [DOI: https://dx.doi.org/10.3103/S1066369X07050064]

29. Baranovskii, E.S.; Artemov, M.A. Existence of optimal control for a nonlinear-viscous fluid model. Int. J. Differ. Equ.; 2016; 2016, 9428128. [DOI: https://dx.doi.org/10.1155/2016/9428128]

30. Baranovskii, E.S. On flows of Bingham-type fluids with threshold slippage. Adv. Math. Phys.; 2017; 2017, 7548328. [DOI: https://dx.doi.org/10.1155/2017/7548328]

31. Fursikov, A.V.; Imanuvilov, O.Y. Local exact boundary controllability of the Boussinesque equations. SIAM J. Control Optim.; 1998; 36, pp. 391-421. [DOI: https://dx.doi.org/10.1137/S0363012996296796]

32. Alekseev, G.V. Solvability of stationary boundary control problems for heat convection equations. Sib. Math. J.; 1998; 39, pp. 844-858. [DOI: https://dx.doi.org/10.1007/BF02672906]

33. Alekseev, G.V. Solvability of inverse extremal problems for stationary heat and mass transfer equations. Sib. Math. J.; 2001; 42, pp. 811-827. [DOI: https://dx.doi.org/10.1023/A:1011940606843]

34. Korotkii, A.I.; Kovtunov, D.A. Optimal boundary control of a system describing thermal convection. Proc. Steklov Inst. Math.; 2011; 272, pp. S74-S100. [DOI: https://dx.doi.org/10.1134/S0081543811020076]

35. Mallea-Zepeda, E.; Lenes, E.; Valero, E. Boundary control problem for heat convection equations with slip boundary condition. Math. Probl. Eng.; 2018; 2018, 7959761. [DOI: https://dx.doi.org/10.1155/2018/7959761]

36. Brizitskii, R.V.; Saritskaya, Z.Y. Control problem for generalized Boussinesq model. J. Phys. Conf. Ser.; 2019; 1268, 012011. [DOI: https://dx.doi.org/10.1088/1742-6596/1268/1/012011]

37. Baranovskii, E.S. The optimal start control problem for 2D Boussinesq equations. Izv. Math.; 2022; 86, pp. 221-242. [DOI: https://dx.doi.org/10.1070/IM9099]

38. Chierici, A.; Giovacchini, V.; Manservisi, S. Analysis and computations of optimal control problems for Boussinesq equations. Fluids; 2022; 7, 203. [DOI: https://dx.doi.org/10.3390/fluids7060203]

39. ukaszewicz, G.; Krzyżanowski, P. On the heat convection equations with dissipation term in regions with moving boundaries. Math. Methods Appl. Sci.; 1997; 20, pp. 347-368. [DOI: https://dx.doi.org/10.1002/(SICI)1099-1476(19970310)20:4<347::AID-MMA864>3.0.CO;2-4]

40. Kakizawa, R. The initial value problem for motion of incompressible viscous and heat-conductive fluids in Banach spaces. Hiroshima Math. J.; 2010; 40, pp. 371-402. [DOI: https://dx.doi.org/10.32917/hmj/1291818851]

41. Amorim, C.B.; de Almeida, M.F.; Mateus, E. Global existence of solutions for Boussinesq system with energy dissipation. J. Math. Anal. Appl.; 2024; 531, 127905. [DOI: https://dx.doi.org/10.1016/j.jmaa.2023.127905]

42. Adams, R.A.; Fournier, J.J.F. Sobolev Spaces, Vol. 40 of Pure and Applied Mathematics; Elsevier: Amsterdam, The Netherlands, 2003.

43. Boyer, F.; Fabrie, P. Mathematical Tools for the Study of the Incompressible Navier–Stokes Equations and Related Models; Springer: New York, NY, USA, 2013; [DOI: https://dx.doi.org/10.1007/978-1-4614-5975-0]

44. Castillo, R.E.; Rafeiro, H. An Introductory Course in Lebesgue Spaces; Springer: Cham, Switzerland, 2016; [DOI: https://dx.doi.org/10.1007/978-3-319-30034-4]

45. Nečas, J. Direct Methods in the Theory of Elliptic Equations; Springer: Berlin/Heidelberg, Germany, 2012; [DOI: https://dx.doi.org/10.1007/978-3-642-10455-8]

46. Gaevskii, K.; Greger, K.; Zakharias, K. Nonlinear Operator Equations and Operator Differential Equations; Mir: Moscow, Russia, 1978.

47. Isac, G. Leray–Schauder Type Alternatives, Complementarity Problems and Variational Inequalities; Springer Science+Business Media: New York, NY, USA, 2006; [DOI: https://dx.doi.org/10.1007/0-387-32900-5]

48. Dinca, G.; Mawhin, J. Brouwer Degree: The Core of Nonlinear Analysis; Birkhäuser: Cham, Switzerland, 2021; [DOI: https://dx.doi.org/10.1007/978-3-030-63230-4]

49. Krasnoselskii, M.A. Topological Methods in the Theory of Nonlinear Integral Equations; Pergamon Press: New York, NY, USA, 1964.

50. Temam, R. Navier–Stokes Equations—Theory and Numerical Analysis; North-Holland Publishing Co.: Amsterdam, The Netherlands, 1977.