Generalized Boussinesq System with Energy Dissipation: Existence of Stationary Solutions

: 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.


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: Here, • Ω is the flow domain, Ω ⊂ R 3 ; • ∂Ω denotes the boundary of Ω; • u is the flow velocity; • π is the pressure; • D(u) denotes the rate-of-deformation tensor with the components • µ(|D(u)|) > 0 is the effective viscosity, dependent on the Euclidean norm of the rate-ofdeformation tensor; • θ is the temperature; • f (x, θ) is the body force acting on the fluid; • k stands for the thermal conductivity coefficient, k > 0; • α is the Darcy (permeability of porous medium) coefficient, α > 0; represents the Rayleigh dissipation function with the spatially averaged rate-of-deformation tensor D ρ (u) (see [18]) that is defined as follows: where and ρ : R 3 → R is a smooth function with compact support such that R 3 ρ(y)dy = 1 and ρ(x) = ρ(y) for x and y with the same Euclidean norm; • c p signifies the specific heat capacity of the fluid, c p > 0; • β stands for the heat exchange coefficient at the walls of the vessel Ω, β > 0; • g(x, θ) is the heat source intensity; • ν denotes the unit outward normal vector to the surface ∂Ω.
It is worth noting that, as the diameter of supp(ρ) tends to zero, the averaged tensor D ρ (u) converges to the original tensor D(u) in the 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 λ ∈ [0, 1].A priori bounds for solutions of this system are derived, which do not depend on λ.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 (Ω).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.

Preliminaries
In this section, we collect notions and statements that will be needed for establishing our main results.
Let E and F be normed linear spaces.By L (E, F), we denote the space of all continuous linear operators from E into F.
Let Ω be a bounded domain in R 3 with the Lipschitz-continuous boundary ∂Ω.We will use the Lebesgue spaces L p (Ω) with p ≥ 1 and the Sobolev space H 1 (Ω) := W 1,2 (Ω).The corresponding classes of vector-valued functions w : Ω → R 3 are denoted by the same symbols but in bold font, that is, L p (Ω) := L p (Ω) 3 and H 1 (Ω) := H 1 (Ω) 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 H 1 (Ω), it is convenient to use the norm , which is equivalent to the standard H 1 -norm, computed as follows: .
We now provide two function spaces developed in the theory of Navier-Stokes equations: We equip V (Ω) with the scalar product defined as follows: for arbitrary vector functions v, w ∈ V (Ω).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 ∥v∥ V (Ω) := (v, v) 1/2 V (Ω) is equivalent to the norm naturally induced from the space H 1 (Ω).
As the main space in which we will find solutions to problem (1), we consider the Cartesian product V (Ω) × H 1 (Ω) equipped with the max-norm:

Some Properties of Spatially Averaged Rate-of-Deformation Tensor
Lemma 1. Suppose Ω is a bounded Lipschitz-continuous domain of R 3 , the tensor D ρ (u) is defined by (2), and Then the operator D ρ : V (Ω) → C ∞ (Ω) has the following properties.
(i) For any vector function u ∈ V (Ω), the estimate is valid.
Proof.Let u ∈ V (Ω).Applying integration by parts, we obtain and, by using the Cauchy-Schwarz inequality, we derive inequality (3).
Taking into account relations (3) and ( 6), we arrive at the following statement.
Lemma 2. Under the conditions of Lemma 1, the operator D ρ : L 2 (Ω) → C(Ω) is well defined and continuous.

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 E be a reflexive separable Banach space with the dual space E * .By ⟨•, •⟩ E * ×E denote a duality pairing between E * and E. Definition 1.An operator A : E → E * is called d-monotone if there exists a strictly increasing function κ : R + → R such that It is easy to see that this condition is less restrictive than the requirement of the strong monotonicity: Proposition 1.Let A : E → E * be a continuous d-monotone operator with κ(s) ≡ c 0 s, where c 0 is a positive constant.Then: Proof.First, we show that the operator A is strictly monotone.Clearly, we have for any u, v ∈ E. Taking into account this inequality and c 0 > 0, it is easy to see that and hence u + v 2 In the Hilbert space E, the last relations imply the equality u = v, which can be verified, for instance, by using the parallelogram law.
Next, we observe that the mapping A is coercive.Specifically, by setting v = 0 in (9), we obtain Due to inequality (9), relation (10) implies that Since E is a Hilbert space, the weak convergence together with (11) imply the strong convergence Thus, we deduce that the operator A satisfies the so-called (S)-property (see [46], Chapter III, § 1).
Taking into account the properties of the operator A established above, proving Proposition 1 is sufficient by applying Theorem 2.2 from [46], Chapter III, § 2.

The Leray-Schauder Alternative for Completely Continuous Mappings
Proposition 2. Let us suppose that E is a real Banach space, B[0, r] := v ∈ E : ∥v∥ E ≤ r , and For more details on fixed-point results and their applications in nonlinear analysis, we refer to the monographs [47,48].

Continuity of Superposition Operator in Lebesgue Spaces
Proposition 3 (Krasnoselskii's theorem).Let G be a bounded domain in the real n-space R n and let Q : G × R → R be a function such that for any y ∈ R and almost every x ∈ G ; • the function Q(•, y) : G → R is measurable for any y ∈ R;

is continuous and bounded as a mapping from L p 1 (G ) into L p 2 (G ).
A detailed proof of this proposition can be found in the monograph [49], Chapter I, § 2.

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 µ : R + → R + is continuous, and there exist constants µ min and µ max such that 0 < µ min ≤ µ(τ) ≤ µ max for all τ ∈ R + ; (H.2) with a positive constant µ ⋆ , the inequality holds for any t, s ∈ R + such that t ≥ s; (H.3) the functions f 1 (•, y), f 2 (•, y), f 3 (•, y) : Ω → R and g(•, y) : Ω → R are measurable for any y ∈ R; (H.4) the functions f 1 (x, •), f 2 (x, •), f 3 (x, •) : R → R and g(x, •) : R → R are continuous for almost every x ∈ Ω; (H.5) there exist functions f max ∈ L 2 (Ω) and g max ∈ L 2 (Ω) such that for any y ∈ R and almost every x ∈ Ω.It can easily be checked that condition (H.2) holds for any C 1 -smooth function µ : R where µ ′ := dµ/dτ.Indeed, by the Mean Value Theorem (Lagrange's theorem) and condition (H.1), we derive for arbitrary t, s ∈ R + , t > s, and some Note that the requirements (H.1) and (H.2) imposed on the viscosity function µ 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 µ ≡ µ min satisfies (H.1) and (H.2) with µ ⋆ = µ min .Now, let us provide the weak formulation of problem (1).
Definition 2. We shall say that a pair (u, θ) is a weak solution to problem (1) if the following two conditions hold: • the pair (u, θ) belongs to the space V (Ω) × H 1 (Ω); • the equalities are valid for all test functions v ∈ V (Ω) and η ∈ H 1 (Ω).
By M we denote the set of all weak solutions to problem (1).

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

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 R 3 and conditions (H.1)-(H.5)are satisfied.Then: (a) problem (1) has at least one weak solution, that is, M ̸ = ∅; (b) any pair (u, θ) ∈ M satisfies the energy equalities: k (c) the set M is compact in the space V (Ω) × H 1 (Ω) as well as in the space L p (Ω) × L p (Ω), 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 (u, θ) from the space V (Ω) × H 1 (Ω) that satisfy the following two equations: where λ is a numerical parameter belonging to the interval [0, 1].
By setting v = u in equality (14), we obtain The first term on the left-hand side of ( 16) is zero.Indeed, applying integration by parts, we obtain Of course, here we used the two identities: div u ≡ 0 in Ω and u ≡ 0 on ∂Ω.
In view of (17), equality (16) reduces to From this equality, using assumptions (H.1) and (H.5) as well as the Cauchy-Schwarz inequality, we deduce that Note that ∥u∥ Here and in the succeeding discussion, the symbol I denotes the identity operator.
Taking into account (20), we can continue (19) as follows: Next, on setting η = θ in (15), we find 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 div u ≡ 0 in Ω and u ≡ 0 on ∂Ω, we obtain Therefore, equality (22) simplifies to From ( 23), using assumption (H.5) and the Cauchy-Schwarz inequality, we derive min{k, β}∥θ∥ 2 Notice that Therefore, (24) can be continued as follows: Taking into account inequalities (3), (20), and (21), from (25) we can easily derive the following bound for θ: Step 3: Operator formulation.Now, let us proceed with an operator interpretation of system ( 14)- (15).To achieve this, we introduce the eight operators: J(u, θ, λ) := (u, θ, λ); It is easy to see that system ( 14)-( 15) is equivalent to the operator equation Step 4: Properties of operators.Let us establish some properties of the operators A and K.
First, using assumption (H.2), we obtain for any u, w ∈ V (Ω).From this, we deduce the inequality This means that the operator A 1 is d-monotone with κ(s) := µ ⋆ s, for s ∈ R + .Applying Proposition 1, we obtain that the operator A 1 is continuously invertible.The same holds true for the operator A 2 .In order to prove the continuous invertibility of this operator, we observe that and apply the classical Lax-Milgram lemma (see [45], Chapter I, § 1.3.1).
Thus, we see that the operator A is continuously invertible, and the following relation Next, note that, due to compact embedding results concerning Sobolev spaces (see [45], Chapter II, § 2.6.1),we have Therefore, the operator J is completely continuous.Moreover, by conditions (H.3)-(H.5),Proposition 3 and Lemma 2, it can be directly checked that the operator K is continuous.Thus, we see that the operator K = K • J (and consequently, the operator A −1 • K) is completely continuous.
Step 5: Solvability.Applying the operator A −1 to both sides of Equation ( 27), we obtain 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 with sufficiently large radius r.
For example, we can take r = r ⋆ , where and the pair of functions (u 0 , θ 0 ) is a weak solution to problem (1).Thus, we establish that M ̸ = ∅, and hence part (a) is proved.
Finally, let us prove statement (c).It is readily apparent from estimates ( 21) and ( 26) that set M is bounded in the Cartesian product V (Ω) × H 1 (Ω).Moreover, this set coincides with the fixed-points set of the completely continuous operator (A −1 • K)(•, •, 1).Taking into account these facts, it is easy to show the compactness of M in V (Ω) × H 1 (Ω).
Therefore, the set M belongs to the space L p (Ω) × L p (Ω) with p ∈ [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 µ = const and α = 0, the results of Theorem 1 coincide with the results from the work [17] by the first author.

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