Bénard Problem for Slightly Compressible Fluids: Existence and Nonlinear Stability in 3D

'is paper shows the existence, uniqueness, and asymptotic behavior in time of regular solutions (a la Ladyzhenskaya) to the Bénard problem for a heat-conducting fluid model generalizing the classical Oberbeck–Boussinesq one.'e novelty of this model, introduced by Corli and Passerini, 2019, and Passerini and Ruggeri, 2014, consists in allowing the density of the fluid to also depend on the pressure field, which, as shown by Passerini and Ruggeri, 2014, is a necessary request from a thermodynamic viewpoint when dealing with convective problems. 'is property adds to the problem a rather interesting mathematical challenge that is not encountered in the classical model, thus requiring a new approach for its resolution.


Introduction
As is well known, one of the most accepted models in the study of convection problems in fluid mechanics is the socalled Oberbeck-Boussinesq approximation (hereafter denoted by O-B) [1][2][3][4][5]. e latter is characterized by the fact that one keeps the fluid incompressible and allows for density changes only in the buoyancy term of the linear momentum equation, by assuming a linear dependence on temperature.
However, as observed in [6] and further elaborated in [7], O-B may not be entirely satisfactory or even consistent from the physical viewpoint [8] if the density is not allowed to have a dependence also on the pressure field p. As a result, in [9], a more general theory aimed at removing this incongruence was proposed. More specifically, while still keeping the incompressibility constraint, the density in the buoyancy term was allowed to depend, linearly, also on p.
e remarkable feature of this new model is that now, unlike the classical O-B, the linear momentum equation contains a pressure term that is not in the classical gradient form. From a mathematical perspective, this feature is very challenging in that some fundamental estimates necessary to ensure stability properties and well-posedness of the relevant problems [10,11] no longer hold and, consequently, the study of such problems requires fresh ideas and different approach.
With this in mind, a systematic study of the mathematical properties of the new model was started, focusing on the classical Bénard problem. As is well known, the latter regards convective motions stemming from a motionless basic flow, r 0 , of a horizontal layer of fluid heated from below. For the new model, the state r 0 is given by [7] v ≡ 0, top of the layer, and h is the thickness of the layer. Finally, α and β are the positive material constants. Notice that (1) tends to the motionless solution for the same problem in the O-B approximation by letting β ⟶ 0. In [4,9,12], the linear and nonlinear stability of r 0 , as well as the existence and uniqueness for the associated "perturbed problem," was investigated. All this work was done in the case of the plane flow. e objective of the current paper is to investigate the same problem in the more general and appropriate threedimensional setting. e approach followed here to overcome the problem of the pressure mentioned earlier on is to introduce a new equation of the elliptic type for the pressure field, disregarding the incompressibility constraint (for steady and for nonlinear equations close to the one of non-Newtonian fluids, this approach that allows to achieve a posteriori the divergence-free equation was considered in [5]. Actually, although this work is done in a "parabolic context" (for the kinetic field), the technique in [5] in some sense has inspired the approach in the present paper). In such a way, the original system of equations is converted to three highly nonlinearly coupled equations for velocity, temperature, and pressure fields. However, for this new problem, it is shown that if the velocity field, v, is sufficiently smooth and also solenoidal at time t � 0, then v must be solenoidal in the entire interval of time where the solution exists, thus providing a solution to the original problem. In this way, it is shown that, if β < 2π, the relevant perturbed problem (to the state r 0 ) has a local, strong solution (a la Ladyzhenskaya) corresponding to initial data of arbitrary size (in a suitable Sobolev class). Such a solution becomes, in fact, global if this size is appropriately restricted and the Rayleigh number is below a certain value, see eorem 1. e plan of the paper is given in the following. After formulating the problem and collecting some preliminary results in Section 1, in Section 2, the classical Galerkin method is used (see, for instance, [4]) to construct an "approximate" solution. In Section 3, which constitutes the main core of this paper, a number of fundamental estimates on the approximate solution is proven that, again by classical techniques, allows us to conclude that the latter converges to a solution of the original problem having the properties described above.

Notations and Preliminary Results
e flow domain is an infinite, horizontal layer bounded in the vertical z-direction with thickness h and unbounded and invariant in the x and y directions. It is assumed that the generic perturbed field is periodic in x and y of a given period L x and L y , respectively, an assumption that is justified by the fact that the convective roll pattern has a periodic structure. e relevant spatial region thus becomes a "periodicity cell" defined (in nondimensional form) as with a � L x /h and b � L y /h. Just for simplicity, set, herein, a � b � 1. e mean value of a function f is indicated as and a partial derivative with a subscript z ξ ϕ � (zϕ/zξ) � ϕ ξ with ξ ∈ x, y, z, t , and similarly for higher-order derivatives. As customary, the inner product in L 2 (Ω 0 ) is denoted by where u and v denote the scalar vector or tensor fields. e components of a vector u are denoted by u ≡ (u x , u y , u z ).
Let W m,p (Ω 0 ) be the usual Sobolev space and ‖u(t)‖ m,p be the associated norm. e subspace of W m,p (Ω 0 ) constituted by those u with 〈u〉 � 0 is denoted by W m,p (Ω 0 ). By L ∞ ((0, T); W m,p (Ω 0 )), we denote the space of functions u such that Likewise, for q ∈ [1, ∞), the Bochner space L q ((0, T); W m,p (Ω 0 )) is the space defined by the norm A decisive role in our stability analysis is played by the Poincaré inequality in Ω 0 : As is known (see [13]), set Γ � Ω 0 ∩ z � 0, 1 { }; this inequality is true if at least one of the following cases is satisfied: (a) functions with a zero mean value, (b) vector functions with a vanishing normal component on Γ, and (c) scalar functions vanishing on Γ. roughout, the same symbol C is used to denote the involved constant, which depends only on the domain. Since boundary conditions are used for which the following integration by parts is valid: if the Poincaré inequality holds true for u, then by the Schwarz inequality, it also follows As a matter of fact, in Ω 0 , there is an equivalence between the norm of the Laplacian and that of the full set of the second derivatives, as shown by Ladyzhenskaya (see [14], equations (17) and (18) at page 18]): Further tools, herein suitable, are some typical 3D immersion inequalities: the Sobolev-Gagliardo-Nirenberg inequality: which in case the Poincaré inequality also holds true becomes 2 International Journal of Differential Equations In what follows, Morrey's inequality is also used (still simplified by the Poincaré inequality): A further preliminary result consists in the estimates for the Laplace problem (e.g., [15,16]).

Lemma 1. Suppose u is a weak solution to
, and the following inequality holds: with c independent of u and f. e same result holds in case In the present framework, once the perturbation fields τ � T − T r and P � p − p r are defined, the following nondimensional system of equations are addressed: where Ra and Pr are the Rayleigh and Prandtl numbers, respectively, and k � (0, 0, 1) is directed along the upward vertical direction. Notice that, given p 0 , the perturbation P has to be zero at the reference state. e boundary conditions are homogeneous for τ, while for v, one imposes the impermeability condition v z � 0 and the stress-free condition v x z � v y z � 0, so that one could compare the stability results with the classic exact ones obtained under the same conditions. Let us set, just for simplicity, ∇ · v � ϕ.
For "regular" solutions, i.e., solutions, such that corresponding to v 0 ∈ W 2,2 (Ω 0 ) at the initial time, one can show that system (17) is equivalent to with Robin's conditions for P at z � 0, 1: e extra nonlinear term lying at the left-hand side of (19) 2 is inserted to make the stress-free condition compatible with the equation, by matching the mean value of the horizontal components. As a matter of fact, v · ∇v + ϕv � ∇ · (v⊗v) gives rise, once integrated, to a vanishing boundary integral, as any other terms in the horizontal projection of equation (19) 2 . Actually, in one direction, the implication is obvious: it is sufficient to take the divergence of the second equation in (17) to obtain the first of (19). e boundary conditions for P are established by evaluating at z � 0, 1 the normal component of (17) 2 : the trace of v z t and, by regularity, the trace of v z zz are defined and vanished, leading to (20). As a matter of fact, by deriving the divergence-free condition with respect to z, one gets v x zx + v y zy + v z zz � 0, and deriving the stress-free condition along the boundary (since one can exchange the order of derivation), one gets v x zx � v x xz � 0 and v y zy � v y yz � 0, which can be replaced in the previous one. Conversely, let us assume that (19) holds true and the solutions are regular; in particular, assume that P ∈ L 2 (0, T; W 2,2 (Ω 0 )). Let us replace (19) 1 in the divergence of (19) 2 , then one obtains 1 Pr Let us remark that 〈ϕ〉 � 0. Indeed, if v x x or v y y is different from zero, then their mean value is zero by periodicity. So, one deduces 〈ϕ〉 � 〈v z z 〉 � 0 because after integration in z ∈ (0, 1), one obtains v z . As a consequence, one can use the Poincaré inequality for ϕ.
Next, one can perform the L 2 -inner product of the equation with ϕ: by taking into account the periodicity conditions, the right-hand side of (21) becomes from the regularity hypotheses, one can deduce that the boundary terms vanish in the trace sense because of impermeability and stress-free conditions, which implies v x zx � v x xz � 0 and v y zy � v y yz � 0. On the other hand, Robin's condition for the pressure implies Δv z � 0 in the momentum balance at the boundary. Hence, since one can derive the impermeability condition with respect to x and y, it follows also v z zz � 0. us, by summing up, one sees that ϕ z vanishes at the boundaries.
In conclusion, one can write 1 Pr Finally, either the divergence of the velocity is zero at t � 0 and then for all t > 0, or it decays exponentially to zero.
International Journal of Differential Equations is means that the divergence-free condition can be fulfilled and is, moreover, physically observable.
Hereafter, a solution with the required regularity is shown to exist, provided the initial data for (19) 2 and (19) 3 are in W 2,2 (Ω 0 ).
Concerning boundary conditions, the change of variable Π � Pe βz turns Robin's ones into the simpler Neumann conditions: and simplifies the system under study, which becomes like the one proposed in [9], page. 7, and studied in [4]: e new variable Π is defined up to an arbitrary constant, unlike P which appears in the buoyancy force and depends on the boundary value p b which is in the basic solution.
us, Π can be prescribed with a mean value zero and verify the Poincaré inequality. e trivial solution (v � c 1 i + c 2 j, τ � 0, P � 0) corresponding to the Galilean invariance has to be considered equivalent to the rest state, so for this reason, solutions are looked in Banach subspaces defined by the following condition: Next, one assigns in Ω 0 the periodicity conditions in the horizontal directions and the stress-free conditions at the boundaries and homogeneous conditions for the deviatory temperature τ: At the initial time, data are given concordant with the previous ones: (v(x, y, z, 0), τ(x, y, z, 0)) � v 0 (x, y, z), τ 0 (x, y, z) .
(28) e plan is studying (17) before proving the equivalence with (25). To this end, one looks for solutions such that a.e. in t, and for k ∈ N 0 , Π belongs to W k,2 (Ω 0 ), which is defined as the closure in the Sobolev space norm of C ∞ (Ω 0 ), whose functions are periodic in x and y, verifies the homogeneous Neumann condition at z � 0, 1, and has a zero mean value (this last feature is necessary to prove the existence).
Analogously, without the isochoric conditions, v ≡ (v x , v y , v z ) has free components such that v x and v y belong to W k,2 (Ω 0 ) (with the mean value zero, as previously noticed), while v z and τ belong to W k,2 0 (Ω 0 ), which is the closure of C ∞ 0 (Ω 0 ).

Galerkin Approximation
Still in [9] and in [4], an existence result for Π and further estimates are given in 2D. e proof is unaffected by the dimension, and the results are summarized below.
with periodic side conditions, has a unique solution Π ∈ W 2,2 (Ω 0 ) such that 〈Π〉 � 0, and the following estimate holds true: In this lemma and hereafter, C(β) stands for a positive function of β such that lim β⟶0 C(β) � c ≥ 0. e same convention is kept for C(β, Pr, Ra).
Notice that both these estimates fall into the classic Navier-Stokes estimate: as β goes to zero. However, here the case ∇ · v ≠ 0 is, in principle, considered.
By taking the data of (29) from system (25), all hypotheses concerning the boundary conditions are verified, and by writing one sees that (31) can be used. Hence, Now, let us introduce Galerkin's "approximate" solutions of (25). Let us denote by φ j j∈N and χ j j∈N the periodic eigenfunctions of the Laplace operator, respectively, verifying Dirichlet and Neumann conditions at zΩ, such that the vector functions Ψ j � (χ j , χ j , φ j ) verify the same boundary condition as v does: 4 International Journal of Differential Equations (36) Moreover, Ψ j 's and φ j 's are complete orthogonal bases, respectively, for the vector and scalar fields in W 1,2 (Ω 0 ). Further, they belong to W 2,2 (Ω 0 ).
e Galerkin approximation solutions v N and τ N , are then defined in a standard way by projecting system (25) on the N-dimensional subspace spanned by the same eigenfunctions. e ordinary differential equation system so obtained is the first order in the unknowns C j N (t) and A j N (t). It is an autonomous homogeneous system in a normal form, whose solutions are easily shown to exist for all N ∈ N.
As a matter of fact, although the pressure terms in the momentum balance (which will be written in a compact form as e − βz ∇Π) cannot be solved a posteriori as for Navier-Stokes equations, one can nevertheless define for all N ∈ N some suitable Π N as a function of C j N (t) and A j N (t).
To this end, let us use (34) in Lemma 2 and denote by A the linear operator: where, of course, ϕ N � ∇ · v N . By Lemma 2, the operator is continuously invertible for 0 < β < 2π and where the summation over repeated indices is understood. According to Lemma 2, the linear operator A − 1 is defined in L 2 (Ω 0 ), and its image belongs to W 2,2 (Ω 0 ). Moreover, A − 1 does not depend on t.
For the sake of brevity, let us denote Next, let us compute is expression can be inserted in the Galerkin system, which can be written as It allows for solutions since the right-hand side is a Lipschitz continuous function of the variables (in fact, it is at most quadratic). For all N ∈ N, the initial conditions C j N (0) � C j 0 � (v 0 , Ψ j ) and A j N (0) � A j 0 � (τ 0 , Ψ j ) correspond to the initial conditions of (25). Since the aim is still to derive the equivalence of the differential systems (17) and (19) by showing existence with ϕ � 0, the required regularity for the initial data is W 2,2 . In particular, one has to define the space as the closure of linear combinations of the eigenfunctions χ j − 〈χ j 〉 so that v x , v y ∈ W 2,2 (Ω 0 ), while v z , τ ∈ W 2,2 0 (Ω 0 ), which is directly obtained by closing the φ j 's.
Moreover, an estimate for ‖v N t (0)‖ 2 and ‖τ N t (0)‖ 2 is hereafter needed, for which the left-hand side of (42) evaluated in t � 0 has to be increased by directly inserting (34) in estimate (31): International Journal of Differential Equations After squaring and summing the inequalities derived by (42) and evaluating at the initial time, one obtains a bound independent of N, estimated by the W 2,2 -norm of the initial data through Bessel's inequality. Since the basis is in W 2,2 (Ω 0 ), by linearity Concerning the nonlinear terms, one can compute, for instance, Next, these nonlinear terms at the right-hand side can be increased by (12) and (13), through Hölder (p � (3/2) and p ′ � 3) and interpolation (p � (4/3) and p ′ � 4) inequalities (the procedure is detailed in the proof of Lemma 3). Finally, using the Poincaré inequality and (10), one gets Here and hereafter, the L 2 -norms of the two nonlinear terms v · ∇v and ϕv are estimated exactly in the same way, but for brevity only the first one will be developed.

Evolutive Estimates and Main Results
Let us formally derive the a priori estimates for solutions of system (25): by quite a standard procedure, they are verified just by Galerkin solutions, denoted here by (v, τ), and, as seen, are bound uniformly with respect to N.

Lemma 3. Let us define the energy functions:
en, for sufficiently small ε > 0, the following inequalities hold true: where c ∈ (0, 1). Proof. One can test (25) with (v, τ), and then one multiplies the third equation by Ra and then sums: 6 International Journal of Differential Equations In order to estimate the nonlinear term arising from the coupling with the pressure equation, one uses (35): the nonlinear term at the right-hand side can be increased by (12) and (13) (61) e terms of order ϵ are still absorbed at the left-hand side, and then inequalities (50) and (51), respectively, follow by (59) and (60).
Proof. By summing side-by-side the first and second inequalities in Lemma 3, one obtains d dt Here, one has and consequently the right-hand side can be bounded by E + E 3 . Hence, then It follows that once the maximal interval corresponding to a vanishing denominator is identified, 8 International Journal of Differential Equations the energy is bounded. Next, coming back to (65) and integrating in [0, T), one deduces a bound for the norm in L ∞ (0, T; W 1,2 (Ω 0 )) ∩ L 2 (0, T; W 2,2 (Ω 0 )), as can be seen hereafter: the inequality implies in particular the integrability of ‖∇v(t)‖ 2 2 and of ‖Δv(t)‖ 2 2 .
□ Although, the Galerkin solutions exist for all t > 0, their uniform estimates exist only in the finite interval defined in (69). Global-in-time solution will be found only by asking suitable smallness of Ra and of the initial data.
Finally, (25) is derived formally with respect to t, and the same procedure is repeated as in the first estimate, but for the nonlinear convective terms. Now, in order to bound the initial data for v t , initial data in W 2,2 is chosen: (86) e proof that τ is bounded in L ∞ (0, T; W 1,2 (Ω 0 )) ∩ L 2 (0, T; W 2,2 (Ω 0 )) is similar and simpler. e proof of global existence follows from Lemmas 3 and 5 by setting In fact, if inequalities (49), (51), and (74) are added, bringing to the left-hand side the three terms with factor Ra/ε and using the Poincaré inequality to absorb them in P 0 for sufficiently small Ra (depending on ε), one gets d dt E(t) + c 1 P(t) ≤ c 2 E 3/2 (t) + E 3 (t) .
e constant c 1 is now fixed positive by suitably restricting Ra. Again, using the Poincaré inequality on the left-hand side of (89), one infers d dt us, the exponential decay follows by choosing at the initial time: In fact, by continuity, for some time, E(t) ≤ E(0). But, this reinforces the previous condition so that the positive function in square brackets cannot decrease and is bounded from below by a positive constant for all times.

Data Availability
e data that support the findings of this study are available within the article.

Conflicts of Interest
e author declares that there are no conflicts of interest.