# A Fully-Mixed Finite Element Method for the n-Dimensional Boussinesq Problem with Temperature-Dependent Parameters

Javier A. Almonacid and Gabriel N. Gatica

## Abstract

In this paper, we introduce and analyze a high-order, fully-mixed finite element method for the free convection of n-dimensional fluids, n{2,3}, with temperature-dependent viscosity and thermal conductivity. The mathematical model is given by the coupling of the equations of continuity, momentum (Navier–Stokes) and energy by means of the Boussinesq approximation, as well as mixed thermal boundary conditions and a Dirichlet condition on the velocity. Because of the dependence on the temperature of the fluid properties, several additional variables are defined, thus resulting in an augmented formulation that seeks the rate of strain, pseudostress and vorticity tensors, velocity, temperature gradient and pseudoheat vectors, and temperature of the fluid. Using a fixed-point approach, smallness-of-data assumptions and a slight higher-regularity assumption for the exact solution provide the necessary well-posedness results at both continuous and discrete levels. In addition, and as a result of the augmentation, no discrete inf-sup conditions are needed for the well-posedness of the Galerkin scheme, which provides freedom of choice with respect to the finite element spaces. In particular, we suggest a combination based on Raviart–Thomas, Lagrange and discontinuous elements for which we derive optimal a priori error estimates. Finally, several numerical examples illustrating the performance of the method and confirming the theoretical rates of convergence are reported.

## 1 Introduction

Free convective flows can be found in a wide amount of settings throughout nature and industry, for instance, in mantle convection, stratified oceanic flows and the cooling of electronic devices, to name a few. Many of these processes can be modeled by coupling the equations of continuity, momentum (Navier–Stokes) and energy using the Boussinesq approximation, which (in this context) assumes the density of the fluid to be constant in all terms of the equations, except in the buoyancy term of the momentum equation, where a linear dependence is considered. Nevertheless, other properties may also vary with temperature, as is the case of, for example, viscosity and thermal conductivity in oils and nanofluids, which poses a significant effect on the fluid flow. In this regard, several finite element methods to approximate the solution to this and related problems, both with constant and temperature-dependent properties have been proposed (see [2, 3, 4, 5, 6, 8, 12, 13, 14, 22, 23, 25, 26, 27] and the references therein).

In particular, the authors in [22, 23] propose finite element methods based on primal formulations of the Boussinesq system. The first one deals with the problem in its primitive variables, while the second one introduces the normal heat flux through the boundary as an additional variable to achieve conformity of the scheme. Nonetheless, both methods are proved to be optimally convergent whenever the exact solution is smooth enough, and the data and the W1,-norm of the velocity and temperature are small enough. More recently in [2], while the authors still consider a primal formulation of the energy equation with a space-dependent thermal conductivity, they also consider a mixed formulation of the momentum and continuity equations but with a temperature-dependent viscosity in a pseudostress-velocity-vorticity formulation. Hence, using fixed-point strategies from [5, 12], and using Raviart–Thomas elements to approximate the pseudostress, Lagrange elements for the velocity and temperature, and discontinuous elements for the vorticity and normal heat flux, they are able to construct an optimally-convergent method whenever the exact solution is smooth enough, and the data is sufficiently small. However, the presence of a variable viscosity leads to restrict the analysis to the two-dimensional case, as it becomes necessary to use Sobolev embeddings into smaller Lp spaces. To overcome this drawback, recent work [4] has shown that, by defining the rate of strain tensor as a new variable (in addition to the pseudostress, velocity and vorticity), the analysis is now valid for two and three-dimensional domains.

The purpose of this work is to extend the analysis and results of [4] by deriving now an augmented fully-mixed finite element method for the Boussinesq problem, but considering this time both the viscosity and the thermal conductivity of the fluid as temperature-dependent functions, and mixed thermal boundary conditions. To this end, we consider again the mixed formulation of the Navier–Stokes equations in [4], to then, using this same approach, construct a mixed formulation for the energy equation. More precisely, we consider the temperature gradient and pseudoheat vector as additional variables, which together with the temperature, rate of strain, pseudostress, velocity and vorticity comprise the unknowns of the problem. At this point, we remark that the main difference with respect to [13] is the consideration in this work of temperature-dependent parameters, which also becomes the cause of defining the rate of strain and temperature gradient in addition to the variables defined in the aforementioned work (recall from [2] that the vorticity appears in this formulation because of the consideration of a more physical version of the Cauchy stress tensor). That being said, part of our analysis follows basically the same uncoupling and fixed-point strategies from [2, 4, 5, 12, 13], reason why we do not provide all details but make the proper references when it corresponds. At a continuous level, we prove that the uncoupled problems are well posed thanks to the Lax–Milgram theorem, and then we prove that the fixed-point operator admits a unique fixed point by means of the Banach fixed-point theorem, whenever the data is sufficiently small and assuming that the exact solution has a slightly higher regularity than the one the well-posedness results provide. Then, following these same steps, we provide a well-posedness result for the Galerkin scheme where one of the key features of this work can be appreciated: there is no need to impose inf-sup conditions on the discrete analysis, which gives us the freedom to choose any combination of finite element subspaces. In particular, we approximate the pseudostress and pseudoheat variables using Raviart–Thomas elements of order k, the velocity and temperature using Lagrange elements of order k+1, and the rate of strain, vorticity and temperature gradient using just discontinuous piecewise polynomials of degree k. When the data are sufficiently small, optimal a priori error estimates can be derived thanks to the Strang lemma, and these are later verified with numerical examples in two and three-dimensional domains.

Finally, we consider worthwhile to mention other features of this fully-mixed method. First, as Dirichlet conditions appear naturally in mixed formulations, it is not necessary to define boundary unknowns to achieve conformity in the scheme (see the Lagrange multiplier defined in [2, 4]), thus unifying the analysis of the uncoupled problems and simplifying the computational implementation of this method. Also, as not only the velocity and temperature of the fluid are part of the solution but also their gradients, many other physical variables of interest can be computed as a simple post-process without requiring numerical differentiations that could deteriorate the good quality of the results. Nevertheless, it is also fair to mention that the use of dual-mixed approaches, and the consequent increase in the number of unknowns, yields much larger linear and nonlinear systems to be solved at the discrete level, and hence, the devising of more efficient numerical methods for computing their solution arises as an unavoidable need. In this regard, we believe that the use of Hybridizable Discontinuous Galerkin Methods (HDG) specially designed for nonlinear problems (see, e.g., [18]), appears as a very attractive possibility to be explored in a separate work.

### 1.1 Outline

The rest of this paper is organized as follows. First, we end this section by introducing some notation that will be used throughout the paper. Then, in Section 2, we set the Boussinesq problem with temperature-dependent viscosity and thermal conductivity functions, and introduce the new variables that will allow us to construct a fully-mixed formulation. Next, in Section 3, we uncouple the problem using a fixed-point argument. The uncoupled problems are then analyzed by means of the Lax–Milgram theorem, and existence of a unique fixed point is proved by fulfilling the hypotheses of the Banach fixed-point theorem. Later, in Section 4 these techniques are used to prove the well-posedness of the corresponding Galerkin scheme, but this time using the Brouwer fixed-point theorem, and then we make a specific choice of finite element subspaces. Finally, in Section 5 we derive some a priori error estimates by suitably applying the usual Strang’s lemmas for linear problems, to then in Section 6 illustrate the good performance of this augmented fully-mixed finite element method and confirm the theoretical rates of convergence through several numerical examples in two and three dimensions.

### 1.2 Preliminaries

Let us denote by ΩRn, n{2,3}, a given bounded domain with polyhedral boundary Γ, and denote by 𝝂 the outward unit normal vector on Γ. Standard notation will be adopted for Lebesgue spaces Lp(Ω) and Sobolev spaces Ws,2(Ω)=:Hs(Ω) with norm s,Ω and seminorm ||s,Ω. In particular, H1/2(Γ) is the space of traces of functions in H1(Ω). By 𝐌 and 𝕄 we will denote the corresponding vectorial and tensorial counterparts of the generic scalar functional space M, and , with no subscripts, will stand for the natural norm of either an element or an operator in any product functional space. In turn, for any vector fields 𝐯=(vi)i=1,,n and 𝐰=(wi)i=1,,n, we set the gradient, divergence and tensor product operators, as

𝐯:=(vixj)i,j=1,,n,div𝐯:=j=1nvjxj,and𝐯𝐰:=(viwj)i,j=1,,n.

In addition, for any tensor fields 𝝉=(τij)i,j=1,,n and 𝜻=(ζij)i,j=1,,n, we let 𝐝𝐢𝐯𝝉 be the divergence operator div acting along the rows of 𝝉, and define the transpose, the trace, the tensor inner product, and the deviatoric tensor, respectively, as

𝝉𝚝:=(τji)i,j=1,,n,tr(𝝉):=i=1nτii,𝝉:𝜻:=i,j=1nτijζij,and𝝉𝚍:=𝝉-1ntr(𝝉)𝕀,

where 𝕀 stands for the identity tensor in Rn×n. Furthermore, we recall that

(𝐝𝐢𝐯;Ω):={𝝉𝕃2(Ω):𝐝𝐢𝐯𝝉𝐋2(Ω)},

equipped with the usual norm

𝝉𝐝𝐢𝐯;Ω2:=𝝉0,Ω2+𝐝𝐢𝐯𝝉0,Ω2,

is a standard Hilbert space in the realm of mixed problems. Finally, in what follows, || denotes the Euclidean norm in Rn. Also, we employ 𝟎 to denote a generic null vector and use C, with or without subscripts, bars, tildes or hats, to mean generic positive constants independent of the discretization parameters, which may take different values at different places.

## 2 The Model Problem

In this section, we first introduce the Boussinesq problem with its original unknowns, to then define suitable new variables that will later allow us to construct a fully-mixed formulation.

### 2.1 The Original Formulation

Let us consider the flow of a non-isothermal, incompressible, Newtonian fluid with varying viscosity and thermal conductivity within a region Ω. Then, under the Boussinesq approximation, the problem reads: Find a velocity 𝐮, a pressure p and a temperature φ such that

(2.1)

(2.1a)-𝐝𝐢𝐯(μ(φ)𝐞(𝐮))+(𝐮)𝐮+p-φ𝐠=𝐟𝚖in Ω,
(2.1b)div𝐮=0in Ω,
(2.1c)-div(k(φ)φ)+𝐮φ=f𝚎in Ω,
(2.1d)𝐮=𝟎on Γ,
(2.1e)φ=φDon ΓD,
(2.1f)k(φ)φ𝝂=0on ΓN,

where 𝐞(𝐮) is the symmetric part of the velocity gradient tensor 𝐮, the boundary Γ of Ω is split as Γ=ΓDΓN, -𝐠𝐋(Ω) is an external force per unit mass (e.g., gravity force, centrifugal force, Coriolis force), 𝐟𝚖𝐋2(Ω) and f𝚎L2(Ω) are source terms, φDH1/2(ΓD) is a prescribed temperature and μ,k:RR+ are the temperature-dependent viscosity and thermal conductivity functions, respectively, which are assumed to be bounded above and below by positive constants, that is, there exist μ2μ1>0 and k2k1>0 such that μ1μ(t)μ2 and k1k(t)k2 for all tR. We also assume that μ and k are Lipschitz continuous functions, that is, there exist Lμ,Lk>0 such that |μ(s)-μ(t)|Lμ|s-t| and |k(s)-k(t)|Lk|s-t| for all s,tR. Notice that the difference in problem (2.1) with respect to the previous work [2] relies on the introduction of a temperature-dependent thermal conductivity, the introduction of mixed boundary conditions for the energy equation, and the presence of source terms which, albeit used in the numerical examples of [2, Section 6], were not considered in the theoretical results, and therefore, we include them for clarity purposes.

### 2.2 Introduction of New Variables

Let us first consider the spaces

𝕃𝚝𝚛2(Ω):={𝐬𝕃2(Ω):𝐬=𝐬𝚝 and tr(𝐬)=0},𝕃𝚜𝚔𝚎𝚠2(Ω):={𝜼𝕃2(Ω):𝜼+𝜼𝚝=𝟎},

and, in a similar way to [4], define the following variables, known respectively as the rate of strain, pseudostress and vorticity tensors:

𝐭:=𝐞(𝐮),𝝈:=μ(φ)𝐞(𝐮)-𝐮𝐮-p𝕀,𝜸:=𝝎(𝐮),

where 𝝎(𝐮) is the skew-symmetric part of the velocity gradient tensor 𝐮. Then the momentum and continuity equations (2.1a) and (2.1b) can be rewritten as

(2.2)

(2.2a)-𝐮+𝐭+𝜸=0in Ω,
(2.2b)μ(φ)𝐭-(𝐮𝐮)𝚍-𝝈𝚍=0in Ω,
(2.2c)-𝐝𝐢𝐯𝝈-φ𝐠=𝐟𝚖in Ω,
(2.2d)𝐮=𝟎on Γ,
(2.2e)Ωtr(𝝈+𝐮𝐮)=0.

Notice that the continuity equation is implicitly present in (2.2b), and it suggests us that rate of strain tensor 𝐭 must be sought in 𝕃𝚝𝚛2(Ω), whereas (2.2a) suggests the vorticity tensor to be sought in 𝕃𝚜𝚔𝚎𝚠2(Ω). Consequently, to characterize these tensors in the two-dimensional case, we only need to know two of the four components of the rate of strain tensor, and only one component of the vorticity tensor (similar simplifications hold for the three-dimensional case).

Next, in order to construct a mixed formulation for the energy equation, we follow the approach taken in [9, 13] and define

𝐩:=k(φ)φ-φ𝐮,

which from now on will be called “pseudoheat”. In addition, analogously to the mixed formulation for the momentum equation, we consider the temperature gradient as another new variable

𝜻:=φ.

Therefore, the energy equation (2.1c) can be rewritten in these terms as

(2.3)

(2.3a)-φ+𝜻=0in Ω,
(2.3b)k(φ)𝜻-φ𝐮-𝐩=0in Ω,
(2.3c)-div𝐩=f𝚎in Ω,
(2.3d)φ=φDon ΓD,
(2.3e)𝐩𝝂=0on ΓN,

where the Neumann condition (2.1f) has been converted to (2.3e) thanks to the no-slip condition 𝐮=𝟎 on Γ.

In this way, the Boussinesq problem (2.1) can now be seen as the set of equations (2.2) and (2.3). Then a mixed formulation of each one of them can be constructed upon integration by parts of (2.2a) and (2.3a) when multiplied by proper test functions, which is the purpose of the next section.

## 3 The Continuous Problem

We now turn to the construction and analysis of a fully-mixed formulation for the Boussinesq problem introduced in the previous section.

### 3.1 An Augmented Fully-Mixed Formulation

First, we recall from [4] the equations corresponding to the mixed formulation of the momentum equation:

(3.1)Ω𝐭:𝝉𝚍+Ω𝜸:𝝉+Ω𝐮𝐝𝐢𝐯𝝉=0for all 𝝉0(𝐝𝐢𝐯;Ω),
(3.2)Ωμ(φ)𝐭:𝐬-Ω(𝐮𝐮)𝚍:𝐬-Ω𝝈𝚍:𝐬=0for all 𝐬𝕃𝚝𝚛2(Ω),
(3.3)-Ω𝐯𝐝𝐢𝐯𝝈-Ω𝝈:𝜼=Ωφ𝐠𝐯+Ω𝐟𝚖𝐯for all (𝐯,𝜼)𝐋2(Ω)×𝕃𝚜𝚔𝚎𝚠2(Ω),

where

0(𝐝𝐢𝐯;Ω):={𝝉(𝐝𝐢𝐯;Ω):Ωtr(𝝉)=0}.

Then, to obtain a mixed formulation for the energy equation, we multiply equation (2.3a) by a test function 𝐪𝐇N(div;Ω), where

𝐇N(div;Ω):={𝐪𝐇(div;Ω):𝐪𝝂=0 on ΓN},

and integrate by parts. Using the boundary condition (2.3d), we get

(3.4)Ω𝜻𝐪+Ωφdiv𝐪=𝐪𝝂,φDΓDfor all 𝐪𝐇N(div;Ω),

where ,ΓD stands for the duality pairing between H-1/2(ΓD) and H1/2(ΓD). Next, we only multiply (2.3b) and (2.3c) by proper test functions:

(3.5)Ωk(φ)𝜻𝝌-Ωφ𝐮𝝌-Ω𝐩𝝌=0for all 𝝌𝐋2(Ω),
(3.6)-Ωψdiv𝐩=Ωf𝚎ψfor all ψL2(Ω).

Notice that, due to the second term in both (3.2) and (3.5), we require the velocity and temperature to have (weak) bounded derivatives as shown in the following inequalities, which can be obtained by using the Hölder inequality and the continuous injections i:H1(Ω)L4(Ω) and 𝐢:𝐇1(Ω)𝐋4(Ω) (cf. [1, Theorem 4.12], [24, Theorem 1.3.4]):

(3.7)|Ω(𝐮𝐰)𝚍:𝐬|c1(Ω)𝐮1,Ω𝐰1,Ω𝐬0,Ωfor all 𝐮,𝐰𝐇1(Ω),𝐬𝕃2(Ω),
(3.8)|Ωφ𝐮𝝌|c2(Ω)φ1,Ω𝐮1,Ω𝝌0,Ωfor all φH1(Ω),𝐮𝐇1(Ω),𝝌𝐋2(Ω),

where c1(Ω) and c2(Ω) are positive constants that depend solely on i and 𝐢. Therefore, at a first glance, the fully-mixed formulation for the Boussinesq problem (2.1) is composed by equations (3.1)–(3.3) and (3.4)–(3.6). Nevertheless, to properly analyze the problem, and to achieve conformity of the scheme, we augment it with redundant Galerkin-type terms that arise from the modelling equations, namely

(3.9)κ1Ω{𝝈𝚍+(𝐮𝐮)𝚍-μ(φ)𝐭}:𝝉𝚍=0for all 𝝉0(𝐝𝐢𝐯;Ω),
(3.10)κ2Ω{𝐝𝐢𝐯𝝈+φ𝐠}𝐝𝐢𝐯𝝉=-κ2Ω𝐟𝚖𝐝𝐢𝐯𝝉for all 𝝉0(𝐝𝐢𝐯;Ω),
(3.11)κ3Ω{𝐞(𝐮)-𝐭}:𝐞(𝐯)=0for all 𝐯𝐇01(Ω),
(3.12)κ4Ω{𝜸-𝝎(𝐮)}:𝜼=0for all 𝜼𝕃𝚜𝚔𝚎𝚠2(Ω),

for the momentum and continuity equations, and

(3.13)κ5Ω{𝐩+φ𝐮-k(φ)𝜻}𝐪=0for all 𝐪𝐇N(div;Ω),
(3.14)κ6Ωdiv𝐩div𝐪=-κ6Ωf𝚎div𝐪for all 𝐪𝐇N(div;Ω),
(3.15)κ7Ω{φ-𝜻}ψ=0for all ψH1(Ω),
(3.16)κ8ΓDφψ=κ8ΓDφDψfor all ψH1(Ω),

for the energy equations, where κj, j{1,,8}, are positive constants to be specified later on.

### Remark 3.1.

We now stress that the incorporation of equations (3.9)–(3.12), which are extracted directly from [4], allows us to establish the ellipticity of the resulting bilinear form (cf. proof of Lemma 3.3 below). The same reason applies to the introduction of (3.13)–(3.16) (cf. proof of Lemma 3.4 below). In turn, the concept “redundant” employed here refers to the fact that the equations yielding these augmented terms have already been considered before for the derivation of the weak formulation, but tested differently to the way they are tested now in (3.9)–(3.12) and (3.13)–(3.16).

In this way, denoting by

:=𝕃𝚝𝚛2(Ω)×0(𝐝𝐢𝐯;Ω)×𝐇01(Ω)×𝕃𝚜𝚔𝚎𝚠2(Ω),𝒬:=𝐋2(Ω)×𝐇N(div;Ω)×H1(Ω),

and

𝐭:=(𝐭,𝝈,𝐮,𝜸),𝐬:=(𝐬,𝝉,𝐯,𝜼),𝜻:=(𝜻,𝐩,φ),𝝌:=(𝝌,𝐪,ψ),

the augmented fully-mixed formulation for this Boussinesq problem is: Find (𝐭,𝜻)×𝒬 such that

(3.17)

𝐀φ(𝐭,𝐬)+𝐁𝐮(𝐭,𝐬)=Fφ(𝐬)+F𝚖(𝐬)for all 𝐬,
𝐂φ(𝜻,𝝌)+𝐃𝐮(𝜻,𝝌)=GD(𝝌)+G𝚎(𝝌)for all 𝝌𝒬,

where, given (𝐰,ϕ)𝐇1(Ω)×H1(Ω), the forms 𝐀ϕ, 𝐁𝐰, 𝐂ϕ, 𝐃𝐰 and the functionals Fϕ, F𝚖, GD, G𝚎 are defined as follows: for all 𝐭,𝐬,

𝐀ϕ(𝐭,𝐬):=Ωμ(ϕ)𝐭:{𝐬-κ1𝝉𝚍}+Ω𝐭:{𝝉𝚍-κ3𝐞(𝐯)}-Ω𝝈𝚍:{𝐬-κ1𝝉𝚍}
+Ω𝐮𝐝𝐢𝐯𝝉+Ω𝜸:𝝉-Ω𝐯𝐝𝐢𝐯𝝈-Ω𝜼:𝝈-κ4Ω𝝎(𝐮):𝜼
(3.18)+κ2Ω𝐝𝐢𝐯𝝈𝐝𝐢𝐯𝝉+κ3Ω𝐞(𝐮):𝐞(𝐯)+κ4Ω𝜸:𝜼,
(3.19)𝐁𝐰(𝐭,𝐬):=-Ω(𝐮𝐰)𝚍:{𝐬-κ1𝝉𝚍};

for all 𝜻,𝝌𝒬,

𝐂ϕ(𝜻,𝝌):=Ωk(ϕ)𝜻{𝝌-κ5𝐪}+Ω𝜻{𝐪-κ7ψ}-Ω𝐩{𝝌-κ5𝐪}
(3.20)+Ωφdiv𝐪-Ωψdiv𝐩+κ6Ωdiv𝐩div𝐪+κ7Ωφψ+κ8ΓDφψ,
(3.21)𝐃𝐰(𝜻,𝝌):=-Ωφ𝐰{𝝌-κ5𝐪};

for all 𝐬,

(3.22)Fϕ(𝐬):=Ωϕ𝐠{𝐯-κ2𝐝𝐢𝐯𝝉},
(3.23)F𝚖(𝐬):=Ω𝐟𝚖{𝐯-κ2𝐝𝐢𝐯𝝉};

for all 𝝌𝒬,

(3.24)GD(𝝌):=𝐪𝝂,φDΓD+κ8ΓDφDψ,
(3.25)G𝚎(𝝌):=Ωf𝚎{ψ-κ6div𝐪}.

### Remark 3.2.

Before we continue, let us have a brief look at what the energy equation (2.1c) would have looked like if we had considered instead a heat-temperature mixed formulation. Indeed, if we define

𝐩~:=k(φ)φin Ω,

equation (2.1c) can be rewritten as

1k(φ)𝐩~-φ=0in Ω,  -div𝐩~+𝐮φ=f𝚎in Ω
φ=φDon ΓD,  𝐩~𝝂=0on ΓN.

Then, multiplying the first equation by a test function 𝐪𝐐1 and integrating by parts, and multiplying the second one by a test function ψQ2 (with 𝐐1 and Q2 suitable spaces to be determined), a mixed formulation for this part of the problem reads: Find (𝐩~,φ)𝐐1×Q2 such that 𝐩~𝝂=0 on ΓN and

Ω1k(φ)𝐩~𝐪+Ωφdiv𝐪=𝐪𝝂,φDΓD,
-Ωψdiv𝐩~=Ωf𝚎ψ-Ωψ𝐮φ

for all 𝐪𝐐1 such that 𝐪𝝂=0 on ΓN, and for all ψQ2. Thus, if we augment this formulation using the same strategy as in (3.13)–(3.16), in particular (3.14) becomes

κ~6Ωdiv𝐩~div𝐪=-κ~6Ωf𝚎div𝐪+κ~6Ω(𝐮φ)div𝐪for all 𝐪𝐐1,

where κ~6>0. In this way, if we consider div𝐪L2(Ω), the last term in the previous equation is not well defined, and therefore we need higher regularity in the test functions and more demanding spaces for the unknowns, e.g., div𝐪L4(Ω) or φW1,4(Ω). Since these changes require a substantial modification of the approach we have taken for this problem, we have left this treatment for a future work.

In the upcoming sections, we analyze problem (3.17) using fixed-point strategies from [2, 4, 5, 13]. More precisely, in Section 3.2 we rewrite (3.17) as a fixed-point problem, in Section 3.3 we prove the main results stated and employed in Section 3.2, and then in Section 3.4 we establish sufficient conditions for existence and uniqueness of this fixed point.

### 3.2 The Fixed-Point Approach

First, let us define 𝐇:=𝐇1(Ω)×H1(Ω) and consider the operator 𝐌:𝐇 defined as

𝐌(𝐰,ϕ)=(𝐌1(𝐰,ϕ),𝐌2(𝐰,ϕ),𝐌3(𝐰,ϕ),𝐌4(𝐰,ϕ)):=𝐭,

where 𝐭 is the solution to the mixed formulation of the momentum equation, that is: Find 𝐭 such that

(3.26)𝐀ϕ(𝐭,𝐬)+𝐁𝐰(𝐭,𝐬)=Fϕ(𝐬)+F𝚖(𝐬)for all 𝐬.

In turn, consider the operator 𝐄:𝐇𝒬 defined as

𝐄(𝐰,ϕ)=(𝐄1(𝐰,ϕ),𝐄2(𝐰,ϕ),𝐄3(𝐰,ϕ)):=𝜻,

where 𝜻𝒬 is now the solution to the mixed formulation of the energy equation, that is: Find 𝜻𝒬 such that

(3.27)𝐂ϕ(𝜻,𝝌)+𝐃𝐰(𝜻,𝝌)=GD(𝝌)+G𝚎(𝝌)for all 𝝌𝒬

The well-definedness of the mappings 𝐌 and 𝐄 is guaranteed by the following lemmas, whose proofs are provided below in Section 3.3.

### Lemma 3.3.

Assume that for δ1(0,2μ2), δ2,δ3(0,2) we choose

κ1(0,2μ1δ1μ2),κ2(0,),κ3(0,2δ2(μ1-κ1μ22δ1)),κ4(0,2δ3κ0(1-δ22)κ3),

where κ0 is a positive constant depending solely on Ω. Then there exists r1>0 such that for each r(0,r1), and for each (w,ϕ)H satisfying w1,Ωr, problem (3.26) has a unique solution t:=M(w,ϕ)H. Moreover, there exists a constant CM>0, independent of (w,ϕ), such that there holds

(3.28)𝐌(𝐰,ϕ)=𝐭C𝐌{𝐠,Ωϕ0,Ω+𝐟𝚖0,Ω}.

### Lemma 3.4.

Assume that for δ4(0,2k2), δ5(0,2) we choose

κ5(0,2δ4k1k2),κ6,κ8(0,),κ7(0,2δ5(k1-κ5k22δ4)).

Then there exists r2>0 such that for each r(0,r2) and for each (w,ϕ)H satisfying w1,Ωr, problem (3.27) has a unique solution 𝛇:=E(w,ϕ)Q. Moreover, there exists a constant CE>0, independent of (w,ϕ), such that there holds

(3.29)𝐄(𝐰,ϕ)=𝜻C𝐄{φD1/2,ΓD+f𝚎0,Ω}.

Consequently, we can define the operator 𝐓:𝐇𝐇 as

𝐓(𝐰,ϕ):=(𝐌3(𝐰,ϕ),𝐄3(𝐌3(𝐰,ϕ),ϕ)),

and look at (3.17) as the fixed-point problem: Find (𝐮,φ)𝐇 such that

(3.30)𝐓(𝐮,φ)=(𝐮,φ).

Indeed, we begin the corresponding analysis by showing next that the operator 𝐓 is well defined in the closed ball of center 0 and radius r, with r(0,min{r1,r2}).

### Lemma 3.5.

Given r(0,r0) with r0:=min{r1,r2}, r1 as in (3.37) and r2 as in (3.40), let

(3.31)𝐖r:={(𝐰,ϕ)𝐇:(𝐰,ϕ)r}.

Assume that the stabilization parameters κj, j{1,,8}, are taken as in Lemmas 3.3 and 3.4, and that there holds

(3.32)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}<r2

with CM>0 as in (3.28). Then T:WrH is well defined and there holds

(3.33)𝐓(𝐰,ϕ)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}+C𝐄{φD1/2,ΓD+f𝚎0,Ω}.

### Proof.

Let (𝐰,ϕ)𝐖r. Since 𝐓(𝐰,ϕ):=(𝐌3(𝐰,ϕ),𝐄3(𝐌3(𝐰,ϕ),ϕ)), the norm of the first component 𝐓 is bounded above by the right-hand side of (3.28), whereas the norm of the second one is bounded above by the right-hand side of (3.29), provided (3.32) holds. Therefore, 𝐓 is well defined and (3.33) trivially holds. ∎

### 3.3 Proofs of Lemmas 3.3 and 3.4

As usual, we consider

𝐭:={𝐭0,Ω2+𝝉𝐝𝐢𝐯;Ω2+𝐯1,Ω2+𝜼0,Ω2}1/2for all 𝐭,

and

𝜻:={𝜻0,Ω2+𝐪div;Ω2+ψ1,Ω2}1/2for all 𝜻𝒬.

### Proof of Lemma 3.3.

Notice that an analogous result has been proved in [4, Lemma 2.3] for the case where a non-homogeneous velocity boundary condition is imposed, however, since we want to find 𝐮𝐇01(Ω), augmenting Galerkin-type boundary terms are not needed, thus modifying the forms and the way ellipticity is proved. Nevertheless, the other properties asked by the Lax–Milgram theorem (see, e.g., [16, Theorem 1.1]) can be easily extracted from this result, that is, given (𝐰,ϕ)𝐇, 𝐀ϕ+𝐁𝐰 is a bilinear form and there exists a positive constant denoted by 𝐀ϕ+𝐁𝐰, independent of (𝐰,ϕ), such that

(3.34)|(𝐀ϕ+𝐁𝐰)(𝐭,𝐬)|𝐀ϕ+𝐁𝐰𝐭𝐬for all 𝐭,𝐬.

We begin by analyzing the ellipticity of the bilinear form 𝐀ϕ. From (3.18), we have for any 𝐬 that

𝐀ϕ(𝐬,𝐬)=Ωμ(ϕ)𝐬:𝐬-κ1Ωμ(ϕ)𝐬:𝝉𝚍-κ3Ω𝐬:𝐞(𝐯)-κ4Ω𝝎(𝐯):𝜼
+κ1𝝉𝚍0,Ω2+κ2𝐝𝐢𝐯𝝉0,Ω2+κ3𝐞(𝐯)0,Ω2+κ4𝜼0,Ω2.

Then, using the bounds for the viscosity and the Cauchy–Schwarz and Young inequalities, we obtain for any δ1,δ2,δ3>0 and any 𝐬 that

𝐀ϕ(𝐬,𝐬)μ1𝐬0,Ω2-κ1μ22δ1𝐬0,Ω2-κ1μ2δ12𝝉𝚍0,Ω2-κ32δ2𝐬0,Ω2-κ3δ22𝐞(𝐯)0,Ω2
-κ42δ3𝝎(𝐯)0,Ω2-κ4δ32𝜼0,Ω2+κ1𝝉𝚍0,Ω2+κ2𝐝𝐢𝐯𝝉0,Ω2+κ3𝐞(𝐯)0,Ω2+κ4𝜼0,Ω2
(μ1-κ1μ22δ1-κ32δ2)𝐬0,Ω2+κ1(1-μ2δ12)𝝉𝚍0,Ω2+κ2𝐝𝐢𝐯𝝉0,Ω2
+κ3(1-δ22)𝐞(𝐯)0,Ω2-κ42δ3|𝐯|1,Ω2+κ4(1-δ32)𝜼0,Ω2.

Note here, in particular, that in order to get positive quantities multiplying 𝐞(𝐯)0,Ω2 and 𝜼0,Ω2, one needs to impose that δ2 and δ3 lie in the interval (0,2). Similar comments apply to the quantities multiplying 𝐬0,Ω2 and 𝝉𝚍0,Ω2, which explains the other ranges assumed in the hypotheses of Lemma 3.3. Then, using the classical Korn inequality (cf. [21, Theorem 10.1]), the norm of 𝐞(𝐯) can be bounded by

(3.35)𝐞(𝐯)0,Ω2κ0(Ω)|𝐯|1,Ω2for all 𝐯𝐇01(Ω),

whereas the norm of 𝝉 can be constructed thanks to the inequality (cf. [7, Proposition 3.1] or [16, Lemma 2.3])

𝝉𝚍0,Ω2+𝐝𝐢𝐯𝝉0,Ω2c3(Ω)𝝉0,Ω2for all 𝝉0(𝐝𝐢𝐯;Ω),

where κ0 and c3 are positive constants depending solely on Ω. Thus, defining the constants

α1:=μ1-κ1μ22δ1-κ32δ2,α2:=min{κ1(1-μ2δ12),κ22},α3:=κ3(1-δ22),α4:=min{α2c3(Ω),κ22},α5:=α3κ0-κ42δ3,α6:=κ4(1-δ32),

we deduce the existence of a positive constant α(Ω):=min{α1,α4,α5,α6}, independent of (𝐰,ϕ), such that

𝐀ϕ(𝐬,𝐬)α(Ω)𝐬2for all 𝐬.

The rest of the proof is identical to [4, Lemma 2.3], but we recall it for completion purposes. Thus, the foregoing inequality, the definition of 𝐁𝐰 (cf. (3.19)) and inequality (3.7) give place to

(𝐀ϕ+𝐁𝐰)(𝐬,𝐬)(α(Ω)-c1(Ω)(1+κ12)1/2𝐰1,Ω)𝐬2for all 𝐬,

and then, we easily see that

(3.36)(𝐀ϕ+𝐁𝐰)(𝐬,𝐬)α(Ω)2𝐬2for all 𝐬

provided that

α(Ω)2c1(Ω)(1+κ12)1/2𝐰1,Ω,

that is,

(3.37)𝐰1,Ωα(Ω)2c1(Ω)(1+κ12)1/2=:r1,

thus proving ellipticity for 𝐀ϕ+𝐁𝐰 under the requirement (3.37). On the other hand, it is not hard to see that Fϕ+F𝚖 is a linear functional, and that there holds

|(Fϕ+F𝚖)(𝐬)|(1+κ22)1/2(𝐠,Ωϕ0,Ω+𝐟𝚖0,Ω)𝐬for all 𝐬.

Therefore, by the Lax–Milgram theorem, there exists a unique 𝐭 solution to (3.26), and the continuous dependence result (3.28) is satisfied with C𝐌:=1α(Ω)2(1+κ22)1/2. ∎

Next, we recall the following Poincaré-type inequality that will help us to prove the ellipticity of the bilinear form defining the operator 𝐄.

### Lemma 3.6.

There exists c4(Ω)>0 such that

|ψ|1,Ω2+ψ0,ΓD2c4(Ω)ψ1,Ω2for all ψH1(Ω).

### Proof.

See [20, Theorem 5.11.2]. ∎

### Proof of Lemma 3.4.

Let (𝐰,ϕ)𝐇. It is clear from (3.20) and (3.21) that 𝐂ϕ and 𝐃𝐰 are bilinear forms, and their bounded character can be seen using the upper bound of the thermal conductivity function, the Cauchy–Schwarz inequality and the trace inequality with boundedness constant c0(Ω), that is, for the first one

|𝐂ϕ(𝜻,𝝌)|k2(1+κ52)1/2𝜻0,Ω𝝌+(1+κ72)1/2𝜻0,Ω𝝌
+(1+κ52)1/2𝐩0,Ω𝝌+φ0,Ωdiv𝐪0,Ω+div𝐩0,Ωψ0,Ω
+κ6div𝐩0,Ωdiv𝐪0,Ω+κ7|φ|1,Ω|ψ|1,Ω+κ8c0(Ω)2φ1,Ωψ1,Ω,

thus obtaining the existence of a positive constant C𝐂>0 such that

|𝐂ϕ(𝜻,𝝌)|C𝐂𝜻𝝌.

For 𝐃𝐰, in addition to the above, we use (3.8) to get

|𝐃𝐰(𝜻,𝝌)|c2(Ω)(1+κ52)1/2𝐰1,Ωφ1,Ω𝝌.

Hence, from the previous two equations, there exists a positive constant denoted by 𝐂ϕ+𝐃𝐰 independent of (𝐰,ϕ) such that

(3.38)|(𝐂ϕ+𝐃𝐰)(𝜻,𝝌)|𝐂ϕ+𝐃𝐰𝜻𝝌for all 𝜻,𝝌𝒬.

Next, to prove that 𝐂ϕ+𝐃𝐰 is elliptic, we first prove that 𝐂ϕ has this property. Indeed, for any 𝝌𝒬, we obtain from (3.20) that

𝐂ϕ(𝝌,𝝌)=Ωk(ϕ)𝝌𝝌-κ5Ωk(ϕ)𝝌𝐪-κ7Ω𝝌ψ+κ5𝐪0,Ω2+κ6div𝐪0,Ω2+κ7|ψ|0,Ω2+κ8ψ0,ΓD2.

Then, using the bounds for the thermal conductivity function, the Cauchy–Schwarz and Young inequalities, we get for any δ4,δ5>0 that

𝐂ϕ(𝝌,𝝌)k1𝝌0,Ω2-κ5k22δ4𝝌0,Ω2-κ5k2δ42𝐪0,Ω2-κ72δ5𝝌0,Ω2-κ7δ52|ψ|1,Ω2
+κ5𝐪0,Ω2+κ6div𝐪0,Ω2+κ7|ψ|0,Ω2+κ8ψ0,Γ2
=(k1-κ5k22δ4-κ72δ5)𝝌0,Ω2+κ5(1-k2δ42)𝐪0,Ω2+κ6div𝐪0,Ω2
+κ7(1-δ52)|ψ|1,Ω2+κ8ψ0,ΓD2.

Hence, applying Lemma 3.6 and defining the constants

β1:=k1-κ5k22δ4-κ72δ5,β2:=min{κ5(1-k2δ42),κ6},β3:=min{κ7(1-δ52),κ8},

there exists β(Ω):=min{β1,β2,c4(Ω)β3} such that

𝐂ϕ(𝝌,𝝌)β(Ω)𝝌2,

which, together with (3.21) and inequality (3.8), allows us to write

(𝐂ϕ+𝐃𝐰)(𝝌,𝝌)(β(Ω)-c2(Ω)(1+κ52)1/2𝐰1,Ω)𝝌2.

Therefore, we easily see that

(3.39)(𝐂ϕ+𝐃𝐰)(𝝌,𝝌)β(Ω)2𝝌2for all 𝝌𝒬,

provided that

β(Ω)2c2(Ω)(1+κ52)1/2𝐰1,Ω,

that is,

(3.40)𝐰1,Ωβ(Ω)2c2(Ω)(1+κ52)1/2=:r2,

thus proving ellipticity for 𝐂ϕ+𝐃𝐰 under the requirement (3.40). Finally, it is clear from (3.24) and (3.25) that for any 𝝌 there holds

|GD(𝝌)|𝐪div,ΩφD1/2,ΓD+κ8φD0,ΓDψ0,ΓD

and

|G𝚎(𝝌)|(1+κ62)1/2f𝚎0,Ω𝝌,

and then, using the continuous injection from H1/2(ΓD) into L2(ΓD) with constant C1/2(ΓD), and the trace inequality in H1(Ω) with constant c0(Ω), we conclude the existence of a positive constant

CG:=min{1+κ8c0(Ω)C1/2,(1+κ62)1/2}

such that

|(GD+G𝚎)(𝝌)|CG{φD1/2,ΓD+f𝚎0,Ω}𝝌for all 𝝌𝒬.

Since GD+G𝚎 is also linear, existence and uniqueness of a solution to problem (3.27) is just a consequence of the Lax–Milgram theorem. Moreover, the continuous dependence result (3.29) holds with C𝐄:=2CGβ(Ω). ∎

At this point we remark that, while (3.36) and (3.39) guarantee the ellipticities of 𝐀ϕ+𝐁𝐰 and 𝐂ϕ+𝐃𝐰 for each (𝐰,ϕ)𝐇:=𝐇1(Ω)×H1(Ω) satisfying assumptions (3.37) and (3.40), respectively, this fact does not imply, however, that the problems given by each one of the rows of (3.17) could be solved separately by using precisely these properties of the foregoing bilinear forms. Indeed, it suffices to realize that the variable 𝐮 (resp. φ) forming part of the definitions of 𝐁𝐮 and 𝐃𝐮 (resp. 𝐀φ and 𝐂φ) is in turn the third component of the unknown 𝐭:=(𝐭,𝝈,𝐮,𝜸) (resp. 𝜻:=(𝜻,𝐩,φ)), which constitutes the main reason why a fixed point approach needs to be applied in order to analyze the solvability of (3.17). Certainly, as already shown by our previous analysis, the above goes through the decoupling of (3.17) into the two problems given by each row of it, thus generating the fixed-point equation (3.30). If one wanted to analyze (3.17) as a whole, the eventual tools of nonlinear functional analysis to be employed (assuming that they exist) would be much beyond any ellipticity or strong monotonicity properties.

We end this subsection by commenting that for computational purposes, a particular choice of stabilization parameters has to be made. Hence, we first consider the middle points of the intervals for δ1, δ2, δ3, δ4, δ5, κ1, κ3, κ4, κ5 and κ7; to then choose κ2, κ6 and κ8 so α2, β2 and β3 can be respectively as large as possible. This results in the following set of values:

κ1=μ1μ22,κ2=μ1μ22,κ3=μ12,κ4=κ0μ14,κ5=k1k22,κ6=k12k22,κ7=k12,κ8=k14.

Notice that κ0, the constant that appears in the Korn inequality (3.35), takes the value 12 when ΩR2.

### 3.4 Solvability Analysis of the Fixed-Point Problem

In this subsection, we pursue to comply with the hypotheses of the Banach fixed-point theorem to ensure existence and uniqueness of a fixed point. More precisely, the main results of this subsection are stated as follows.

### Lemma 3.7.

Let r(0,r0) with r0 given as in Lemma 3.5. Assume that

(3.41)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}r

with CM>0 as in (3.28) (cf. Lemma 3.3). Then T:WrH is Lipschitz-continuous and there holds

𝐓(𝐰1,ϕ1)-𝐓(𝐰2,ϕ2)L𝐓(𝐰1,ϕ1)-(𝐰2,ϕ2)for all (𝐰1,ϕ1),(𝐰2,ϕ2)𝐖r,

where

L𝐓:=(L𝐌(𝐠,Ω+𝐟𝚖ε,Ω)+1)(L𝐄(φD1/2+ε,ΓD+f𝚎0,Ω)+1)-1,

and LT and LE are the constants yielding the Lipschitz-continuity of the operators M and E, respectively (cf. Lemmas 3.10 and 3.11 below).

### Theorem 3.8.

Let r(0,r0) with r0 as in Lemma 3.5, and let Wr be the closed ball in H with center (0,0) and radius r (cf. (3.31)). Assume that the stabilization parameters κj, j{1,,8}, are taken as in Lemmas 3.3 and 3.4, and that the data g, fm, fe and φD satisfy

(3.42)

(3.42a)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}+C𝐄{φD1/2,ΓD+f𝚎0,Ω}r,
(3.42b)L𝐌(𝐠,Ω+𝐟𝚖ε,Ω)<2-1,
(3.42c)L𝐄(φD1/2+ε,ΓD+f𝚎0,Ω)<2-1

with constants CM, CE, LM and LE as in Lemmas 3.3, 3.4, 3.10 and 3.11, respectively. Then there exists a unique (t,𝛇)H×Q, with (u,φ)Wr, solution to the fully-mixed formulation (3.17) of the Boussinesq problem (2.1). Moreover, there holds

𝐭C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}

and

𝜻C𝐄{φD1/2,ΓD+f𝚎0,Ω}.

In order to prove Lemma 3.7 and Theorem 3.8, we first recall from (3.33) that 𝐓 maps the ball 𝐖r into itself whenever the data satisfy

(3.43)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}+C𝐄{φD1/2,ΓD+f𝚎0,Ω}r.

Therefore, it only remains to prove that 𝐓 is a Lipschitz-continuous and contracting mapping. As in [2, 4, 5, 9], this analysis will require a further-regularity assumption on the solution of both uncoupled problems. Hence, in what follows we assume that 𝐟𝚖𝐇ε(Ω) and φDH1/2+ε(ΓD) for some ε(0,1) when n=2, or ε[12,1) when n=3, and that for each (𝐳,ρ)𝐇, with 𝐳1,Ωr, r>0 given, there hold, on the one hand, that

(𝐬,𝝉,𝐯,𝜼):=𝐌(𝐳,ρ)𝕃𝚝𝚛2(Ω)ε(Ω)×0(𝐝𝐢𝐯;Ω)ε(Ω)×𝐇01+ε(Ω)×𝕃𝚜𝚔𝚎𝚠2(Ω)ε(Ω)

and (following (3.28))

(3.44)𝐬ε,Ω+𝝉ε,Ω+𝐯1+ε,Ω+𝜼ε,ΩC~𝐌(r){𝐠,Ωρ1,Ω+𝐟𝚖ε,Ω},

and on the other hand, that

(𝝌,𝐪,ψ):=𝐄(𝐳,ρ)𝐋2(Ω)𝐇ε(Ω)×𝐇N(div;Ω)𝐇ε(Ω)×H1+ε(Ω)

and (following (3.29))

(3.45)𝝌ε,Ω+𝐪ε,Ω+ψ1+ε,ΩC~𝐄(r){φD1/2+ε,ΓD+f𝚎0,Ω}.

Here, C~𝐌(r) and C~𝐄(r) are positive constants independent of 𝐳 but depending on the upper bound r of its 𝐇1-norm. A difference with respect to other works, namely [2, 4], is that the presence in this case of a homogeneous boundary condition for the momentum equation leads us to ask for more regularity to the source term, which for this uncoupled problem, becomes 𝐟𝚖+ρ𝐠, hence the reason why not only we consider the Hε-norm of 𝐟𝚖 in (3.44) but also the Hε-norm of ρ, though since we are already assuming that ρH1(Ω), we have used a continuous injection from H1(Ω) into Hε(Ω) to transform this Hε-norm into the H1-norm that appears. That being said, we now proceed to prove the Lipschitz-continuity of 𝐌 and 𝐄. We begin by recalling from [4] a preliminary result for 𝐌.

### Lemma 3.9.

Let r(0,r0) with r0 given as in Lemma 3.5. Then there exists a positive constant C^M, independent of r, such that

(3.46)

𝐌(𝐰1,ϕ1)-𝐌(𝐰2,ϕ2)C^𝐌{𝐌1(𝐰1,ϕ1)ε,Ωϕ1-ϕ2Ln/ε(Ω)
+𝐌3(𝐰1,ϕ1)1,Ω𝐰1-𝐰21,Ω+𝐠,Ωϕ1-ϕ20,Ω}

for all (w1,ϕ1),(w2,ϕ2)H such that w11,Ω,w21,Ωr.

### Proof.

See [4, Lemma 2.6]. ∎

We stress here that Lemma 3.9 makes use precisely of the regularity assumption (3.44) and the fact, that under the specified range of ε, H1(Ω) is continuously embedded in Ln/ε(Ω) for n{2,3}. The Lipschitz-continuity of 𝐌 is proved now.

### Lemma 3.10.

Let r(0,r0), with r0 given as in Lemma 3.5. Then there exists LM>0 such that

(3.47)𝐌(𝐰1,ϕ1)-𝐌(𝐰2,ϕ2)L𝐌{𝐠,Ω+𝐟𝚖ε,Ω}(𝐰1,ϕ1)-(𝐰2,ϕ2)

for all (w1,ϕ1),(w2,ϕ2)Wr.

### Proof.

Let (𝐰1,ϕ1),(𝐰2,ϕ2)𝐖r. Considering that H1(Ω)Ln/ε(Ω) with constant C~ε, that H1+ε(Ω)H1(Ω) with constant C^ε, and that the solution to the problem defining 𝐌 verifies the further-regularity assumption (3.44), the following can be inferred from (3.46):

𝐌(𝐰1,ϕ1)-𝐌(𝐰2,ϕ2)C^𝐌{C~𝐌(r)(C~ε2+C^ε2)1/2(r𝐠,Ω+𝐟𝚖ε,Ω)+𝐠,Ω}(𝐰1,ϕ1)-(𝐰2,ϕ2).

Then, defining the constants

C^1=C~𝐌(r)(C~ε2+C^ε2)1/2,C^2=1+C^1r,

the result (3.47) holds with L𝐌:=C^𝐌max{C^1,C^2}. ∎

In a similar way, we show next that the operator 𝐄 is Lipschitz-continuous.

### Lemma 3.11.

Let r(0,r0) with r0 given as in Lemma 3.5. Then there exists LE>0 such that

(3.48)𝐄(𝐰1,ϕ1)-𝐄(𝐰2,ϕ2)L𝐄{φD1/2+ε,ΓD+f𝚎0,Ω}(𝐰1,ϕ1)-(𝐰2,ϕ2)

for all (w1,ϕ1),(w2,ϕ2)Wr.

### Proof.

Let (𝐰1,ϕ1),(𝐰2,ϕ2)𝐇 such that 𝐰11,Ω,𝐰21,Ωr, and let 𝜻j=(𝜻j,𝐩j,φj):=𝐄(𝐰j,ϕj), for j{1,2}, be the corresponding solutions of (3.27). Then, since for any 𝝌𝒬 there holds

𝐂ϕ1(𝜻1,𝝌)+𝐃𝐰1(𝜻1,𝝌)=GD(𝝌)+G𝚎(𝝌),

it follows that

𝐂ϕ2(𝜻1-𝜻2,𝝌)+𝐃𝐰2(𝜻1-𝜻2,𝝌)=𝐂ϕ2(𝜻1,𝝌)+𝐃𝐰2(𝜻1,𝝌)-GD(𝝌)-G𝚎(𝝌)
=𝐂ϕ2(𝜻1,𝝌)-𝐂ϕ1(𝜻1,𝝌)+𝐃𝐰2(𝜻1,𝝌)-𝐃𝐰1(𝜻1,𝝌)

for any 𝝌𝒬. In this way, using the ellipticity of the bilinear form 𝐂ϕ2+𝐃𝐰2 (cf. (3.39)), we see that

(3.49)β(Ω)2𝜻1-𝜻22(𝐂ϕ2+𝐃𝐰2)(𝜻1-𝜻2,𝜻1-𝜻2)=(𝐂ϕ2-𝐂ϕ1)(𝜻1,𝜻1-𝜻2)+𝐃𝐰2-𝐰1(𝜻1,𝜻1-𝜻2)=Ω{k(ϕ2)-k(ϕ1)}𝜻1{(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}+Ωφ1(𝐰1-𝐰2){(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}.

The last term in the previous inequality can be easily split using (3.8), that is,

(3.50)|Ωφ1(𝐰1-𝐰2){(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}|c2(Ω)(1+κ52)1/2φ11,Ω𝐰1-𝐰21,Ω𝜻1-𝜻2,

whereas for the first term we use the Hölder inequality to show that

(3.51)|Ω{k(ϕ2)-k(ϕ1)}𝜻1{(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}|LkΩ|(ϕ2-ϕ1)𝜻1{(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}|Lk(1+κ52)1/2(ϕ2-ϕ1)𝜻10,Ω𝜻1-𝜻2Lk(1+κ52)1/2ϕ2-ϕ1L2q(Ω)𝜻1𝐋2p(Ω)𝜻1-𝜻2

with p,q[1,) such that 1p+1q=1. Having in mind the further-regularity assumption (3.45), we recall that Hε(Ω) is continuously embedded into L2p(Ω) whenever

2p={21-εif n=2,63-2εif n=3,

and hence, there exists Cε>0 such that

(3.52)𝜻𝐋2p(Ω)Cε𝜻ε,Ωfor all 𝜻𝐇ε(Ω).

In this way,

2q=2pp-1={2εif n=2,3εif n=3,=nε,

and (3.51) now yields

(3.53)|Ω{k(ϕ2)-k(ϕ1)}𝜻1{(𝜻1-𝜻2)-κ5(𝐩1-𝐩2)}|CεLk(1+κ52)1/2𝜻1ε,Ωϕ1-ϕ2Ln/ε(Ω)𝜻1-𝜻2.

Since 𝜻1=𝐄1(𝐰1,ϕ1) and φ1=𝐄3(𝐰1,ϕ1), putting together (3.50) and (3.53) into (3.49), we find that there exists

C^𝐄:=2(1+κ52)1/2β(Ω)max{CεLk,c2(Ω)}

such that

(3.54)𝐄(𝐰1,ϕ1)-𝐄(𝐰2,ϕ2)C^𝐄{𝐄1(𝐰1,ϕ1)ε,Ωϕ1-ϕ2Ln/ε(Ω)+𝐄3(𝐰1,ϕ1)1,Ω𝐰1-𝐰21,Ω}

for any (𝐰1,ϕ1),(𝐰2,ϕ2)𝐇 such that 𝐰11,Ω,𝐰21,Ωr. Considering the same injections as in the proof of Lemma 3.10 and the further-regularity assumption (3.45), we see that for any (𝐰1,ϕ1),(𝐰2,ϕ2)𝐖r,

(3.55)𝐄(𝐰1,ϕ1)-𝐄(𝐰2,ϕ2)C^𝐄C~𝐄(r)(C~ε2+C^ε2)1/2(φD1/2+ε,ΓD+f𝚎0,Ω)(𝐰1,ϕ1)-(𝐰2,ϕ2),

thus obtaining (3.48) with L𝐄:=C^𝐄C~𝐄(r)(C~ε2+C^ε2)1/2. ∎

Consequently, the proof of Lemma 3.7, that is, the Lipschitz-continuity of 𝐓:𝐖r𝐇 becomes straightforward.

### Proof of Lemma 3.7.

Let (𝐰1,ϕ1),(𝐰2,ϕ2)𝐖r. Then, according to the definition of 𝐓,

𝐓(𝐰1,ϕ1)-𝐓(𝐰2,ϕ2)𝐌3(𝐰1,ϕ1)-𝐌3(𝐰2,ϕ2)+𝐄3(𝐌3(𝐰1,ϕ1),ϕ1)-𝐄3(𝐌3(𝐰2,ϕ2),ϕ2).

The bound for the first term comes directly from (3.47), whereas for the second one, we first use the Lipschitz-continuity of 𝐄 (notice that the assumption (3.41) is necessary to be able to use this property, i.e., (3.48)) and then the Lipschitz-continuity of 𝐌, that is,

𝐄3(𝐌3(𝐰1,ϕ1),ϕ1)-𝐄3(𝐌3(𝐰2,ϕ2),ϕ2)
L𝐄(φD1/2+ε,ΓD+f𝚎0,Ω)(𝐌3(𝐰1,ϕ1),ϕ1)-(𝐌3(𝐰2,ϕ2),ϕ2)
L𝐄(φD1/2+ε,ΓD+f𝚎0,Ω)(L𝐌(𝐠,Ω+𝐟𝚖ε,Ω)+1)(𝐰1,ϕ1)-(𝐰2,ϕ2).

Summarizing, Lemma 3.5 ensures us that 𝐓:𝐇𝐇 is well defined in the closed ball 𝐖r, and maps the ball into itself whenever the data satisfy condition (3.43). Then this operator is Lipschitz-continuous according to the foregoing Lemma, and it becomes a contracting map whenever L𝐓<1. Putting together all these results, we get Theorem 3.8, which is proved as follows.

### Proof of Theorem 3.8.

Notice that (3.42a) ensures that both (3.32) and (3.43) hold, so that 𝐓:𝐖r𝐖r is indeed well defined. In addition, (3.42b) and (3.42c) imply that 𝐓 is a contracting map (that is, L𝐓<1). Therefore, the result is just a consequence of the Banach fixed-point theorem applied to the operator 𝐓 in 𝐖r. ∎

## 4 The Galerkin Scheme

We advocate in this section to present the Galerkin scheme for the continuous problem (3.17), whose well-posedness will be proved using the same steps and methods as in the previous section.

### 4.1 Preliminaries

Let us consider 𝒯h a regular triangulation of Ω¯ made by triangles T when n=2 (or tetrahedra when n=3) of diameter hT an define the meshsize h:=max{hT:T𝒯h}. Then consider arbitrary finite-dimensional subspaces h𝐭𝕃𝚝𝚛2(Ω), h𝝈0(𝐝𝐢𝐯;Ω), 𝐇h𝐮𝐇01(Ω), h𝜸𝕃𝚜𝚔𝚎𝚠2(Ω), 𝐇h𝜻𝐋2(Ω), 𝐇h𝐩𝐇N(div;Ω), HhφH1(Ω) and denote by

h:=h𝐭×h𝝈×𝐇h𝐮×h𝜸,𝒬h:=𝐇h𝜻×𝐇h𝐩×Hhφ,

and

𝐭h:=(𝐭h,𝝈h,𝐮h,𝜸h),𝐬h:=(𝐬h,𝝉h,𝐯h,𝜼h),𝜻h:=(𝜻h,𝐩h,φh),𝝌h:=(𝝌h,𝐪h,ψh).

Hence, according to the continuous formulation (3.17), the Galerkin scheme reads: Find (𝐭h,𝜻h)h×𝒬h such that

(4.1)

(4.1a)𝐀φh(𝐭h,𝐬h)+𝐁𝐮h(𝐭h,𝐬h)=Fφh(𝐬h)+F𝚖(𝐬h)for all 𝐬hh,
(4.1b)𝐂φh(𝜻h,𝝌h)+𝐃𝐮h(𝜻h,𝝌h)=GD(𝝌h)+G𝚎(𝝌h)for all 𝝌h𝒬h,

where the forms 𝐀φh, 𝐁𝐮h, 𝐂φh and 𝐃𝐮h, and the functionals Fφh, F𝚖, GD and G𝚎 are defined by (3.18)–(3.25).

We will see that it is possible to establish sufficient conditions for well-posedness of (4.1) in the same form they were established for the continuous problem (3.17). To this end, we now split the discrete formulation into the two corresponding mixed formulations. In fact, we first set 𝐇h:=𝐇h𝐮×Hhφ and consider the operator 𝐌h:𝐇hh defined for any (𝐰h,ϕh)𝐇h by

𝐌h(𝐰h,ϕh)=(𝐌1,h(𝐰h,ϕh),𝐌2,h(𝐰h,ϕh),𝐌3,h(𝐰h,ϕh),𝐌4,h(𝐰h,ϕh)):=𝐭h,

where 𝐭hh is the solution to the problem: Find 𝐭hh such that

(4.2)𝐀ϕh(𝐭h,𝐬h)+𝐁𝐰h(𝐭h,𝐬h)=Fϕh(𝐬h)+F𝚖(𝐬h)for all 𝐬hh.

In turn, we let 𝐄h:𝐇h𝒬h be the operator defined for any (𝐰h,ϕh)𝐇h by

𝐄h(𝐰h,ϕh)=(𝐄1,h(𝐰h,ϕh),𝐄2,h(𝐰h,ϕh),𝐄3,h(𝐰h,ϕh)):=𝜻h,

where 𝜻h𝒬h is the solution to the problem: Find 𝜻h𝒬h such that

(4.3)𝐂ϕh(𝜻h,𝝌h)+𝐃𝐰h(𝜻h,𝝌h)=GD(𝝌h)+G𝚎(𝝌h)for all 𝝌𝒬h.

Therefore, by introducing the operator 𝐓h:𝐇h𝐇h as

(4.4)𝐓h(𝐰h,ϕh):=(𝐌3,h(𝐰h,ϕh),𝐄3,h(𝐌3,h(𝐰h,ϕh),ϕh)),

we can rewrite (4.1) as the fixed-point problem: Find (𝐮h,φh)𝐇h such that

𝐓h(𝐮h,φh)=(𝐮h,φh).

In this case, existence of a fixed point for this problem will be proved by means of the Brouwer fixed-point theorem, which we recall next.

### Theorem 4.1 (Brouwer).

Let W be a compact and convex subset of a finite-dimensional Banach space, and let T:WW be a continuous mapping. Then T has at least one fixed-point.

### Proof.

See [10, Theorem 9.9-2]. ∎

### 4.2 Solvability Analysis

We first study under which conditions 𝐓h is well defined by looking at the well-posedness of the uncoupled problems (4.2) and (4.3). It is easy to realize that the analysis comes in a straightforward way from the one realized in Lemmas 3.3, 3.4 and 3.5, thus obtaining the following discrete versions for these results.

### Lemma 4.2.

Assume that for δ1(0,2μ2), δ2,δ3(0,2) we choose

κ1(0,2μ1δ1μ2),κ2(0,),κ3(0,2δ2(μ1-κ1μ22δ1)),κ4(0,2δ3κ0(1-δ22)κ3),

where κ0 is the constant in the Korn inequality (3.35). Then, for each r(0,r1), with r1 as in (3.37), problem (4.2) has a unique solution th:=Mh(wh,ϕh)Hh for each (wh,ϕh)Hh such that wh1,Ωr. Moreover, there holds

𝐌h(𝐰h,ϕh)=𝐭hC𝐌{𝐠,Ωϕh0,Ω+𝐟𝚖0,Ω}

with CM>0 as in (3.28).

### Proof.

It follows from a direct application of the Lax–Milgram theorem to (4.2) in the same way it was applied in Lemma 3.3. ∎

### Lemma 4.3.

Assume that for δ4(0,2k2), δ5(0,2) we choose

κ5(0,2δ4k1k2),κ6,κ8(0,),κ7(0,2δ5(k1-κ5k22δ4)).

Then, for each r(0,r2) with r2 as in (3.40), problem (4.3) has a unique solution 𝛇h:=Eh(wh,ϕh)Qh for each (wh,ϕh)Hh such that wh1,Ωr. Moreover, there holds

𝐄h(𝐰h,ϕh)=𝜻hC𝐄{φD1/2,ΓD+f𝚎0,Ω}

with CE>0 as in (3.29).

### Proof.

It also follows from an application of the Lax–Milgram theorem to (4.3) in the same way it was applied in Lemma 3.4. ∎

Analogously to the continuous case, the previous two lemmas provide the well-definedness of the operator 𝐓h.

### Lemma 4.4.

Given r(0,r0) with r0:=min{r1,r2}, r1 as in (3.37) and r2 as in (3.40), let

(4.5)𝐖h,r:={(𝐰h,ϕh)𝐇h:(𝐰h,ϕh)r}.

Assume that the stabilization parameters κj, j{1,,8}, are taken as in Lemmas 4.2 and 4.3 and that the data satisfy (3.32). Then Th:Wh,rHh is well defined and there holds

𝐓h(𝐰h,ϕh)C𝐌{r𝐠,Ω+𝐟𝚖0,Ω}+C𝐄{φD1/2,ΓD+f𝚎0,Ω}.

Moreover, Th maps the ball Wh,r into itself whenever (3.43) holds.

### Proof.

It comes straightforwardly from Lemmas 4.2 and 4.3. We omit further details. ∎

We now turn to prove the continuity of 𝐓h. Notice that we are no longer able to use the further-regularity assumptions (3.44) and (3.45), and therefore, we will only obtain results that are analogous to (3.46) and (3.54) in terms of L4 norms. Indeed, although their bounds are not known, they are finite since the solution of the Galerkin scheme will be formed by piecewise polynomial functions.

### Lemma 4.5.

Let r(0,r0) with r0 as in Lemma 4.4. Then there exists a positive constant C¯M, independent of r, such that

𝐌h(𝐰h1,ϕh1)-𝐌h(𝐰h2,ϕh2)C¯𝐌{𝐌1,h(𝐰h1,ϕh1)𝕃4(Ω)ϕh1-ϕh2L4(Ω)
(4.6)+𝐌3,h(𝐰h1,ϕh1)1,Ω𝐰h1-𝐰h21,Ω+𝐠,Ωϕh1-ϕh20,Ω}

for all (wh1,ϕh1),(wh2,ϕh2)Hh such that wh11,Ω,wh21,Ωr.

### Proof.

It comes from [4, Lemma 2.6], but employing an L4-L4-L2 argument, instead of the L2p-L2q-L2 argument taken in [4, equation (2.52)]. ∎

### Lemma 4.6.

Let r(0,r0) with r0 as in Lemma 4.4. Then there exists a positive constant C¯E such that

(4.7)𝐄h(𝐰h1,ϕh1)-𝐄h(𝐰h2,ϕh2)C¯𝐄{𝐄1,h(𝐰h1,ϕh1)𝐋4(Ω)ϕh1-ϕh2L4(Ω)+𝐄3,h(𝐰h1,ϕh1)1,Ω𝐰h1-𝐰h21,Ω}

for all (wh1,ϕh1),(wh2,ϕh2)Hh such that wh11,Ω,wh21,Ωr.

### Proof.

As in the previous lemma, the proof is based on the one for its continuous counterpart Lemma 3.11, but just taking p=q=2 in (3.51). ∎

Consequently, we have the following result for the operator 𝐓h.

### Lemma 4.7.

Given r(0,r0) with r0 as in Lemma 4.4, let Wh,r be the closed ball in Hh with center (0,0) and radius r (cf. (4.5)), and assume that (3.41) holds. Then Th:Wh,rHh is continuous and there exists a constant CTh>0 such that

(4.8)𝐓h(𝐰h1,ϕh1)-𝐓h(𝐰h2,ϕh2)C𝐓h{𝐄1,h(𝐌3,h(𝐰h1,ϕh1),ϕh1)𝐋4(Ω)+(1+𝐄3,h(𝐌3,h(𝐰h1,ϕh1),ϕh1)1,Ω)×(𝐌1,h(𝐰h1,ϕh1)𝕃4(Ω)+𝐌3,h(𝐰h1,ϕh1)1,Ω+𝐠,Ω)}×(𝐰h1,ϕh1)-(𝐰h2,ϕh2)

for all (wh1,ϕh1),(wh2,ϕh2)Wh,r.

### Proof.

Let (𝐰h1,ϕh1),(𝐰h2,ϕh2)𝐖h,r. From the definition of 𝐓h (cf. (4.4)) we see that

𝐓h(𝐰h1,ϕh1)-𝐓h(𝐰h2,ϕh2)𝐌3,h(𝐰h1,ϕh1)-𝐌3,h(𝐰h2,ϕh2)1,Ω
+𝐄3,h(𝐌3,h(𝐰h1,ϕh1),ϕh1)-𝐄3,h(𝐌3,h(𝐰h2,ϕh2),ϕh2),

and since (3.41) holds, 𝐌3,h(𝐰h1,ϕh1)1,Ωr and we can use (4.7) to bound the second term in the previous inequality. In this way we get

𝐓h(𝐰h1,ϕh1)-𝐓h(𝐰h2,ϕh2)C¯𝐄Ci𝐄1,h(𝐌3,h(𝐰h1,ϕh1),ϕh1)𝐋4(Ω)ϕh1-ϕh21,Ω
+(1+C¯𝐄𝐄3,h(𝐌3,h(𝐰h1,ϕh1),ϕh1)1,Ω)
×𝐌3,h(𝐰h1,ϕh1)-𝐌3,h(𝐰h2,ϕh2)1,Ω,

where Ci is the boundedness constant of the continuous injection of H1(Ω) into L4(Ω). Then we can use the analogous result for 𝐌h (cf. (4.6)) and, after some algebraic work, we arrive at (4.8) with

C𝐓h=max{C¯1,C¯2},C¯1=C¯𝐌max{1,C¯𝐄}max{1,Ci},C¯2=CiC¯𝐄.

Having proved that 𝐓h is continuous and that it maps the closed and convex set 𝐖h,r into itself whenever (3.43) holds, existence of at least one fixed point is ensured by the Brouwer fixed-point theorem, equivalently, existence of at least one solution for (4.1); a result that is summarized in the following statement.

### Theorem 4.8.

Let r(0,r0) with r0 as in Lemma 4.4, and let Wh,r be the closed ball of center (0,0) and radius r (cf. (4.5)). Assume the stabilization parameters κj, j{1,,8}, are taken as in Lemmas 4.2 and 4.3 and that the data satisfy (3.43). Then the Galerkin scheme (4.1) has at least one solution (th,𝛇h)Hh×Qh with (uh,φh)Wh,r. Moreover, there hold

𝐭hC𝐌{r𝐠,Ω+𝐟𝚖0,Ω},

and

𝜻hC𝐄{φD1/2,ΓD+f𝚎0,Ω},

with CM and CE as in Lemmas 4.2 and 4.3, respectively.

### 4.3 Specific Finite Element Subspaces

An interesting point to realize in this fully-mixed approach with respect to the finite-dimensional subspaces is that we have not imposed any kind of inf-sup conditions as in [2, 4], or any other requirement than being finite-dimensional, which give us the chance to freely choose these subspaces. In particular, the most natural finite element subspaces of 0(𝐝𝐢𝐯;Ω) and 𝐇N(div;Ω), that is, those preserving in the cheapest way the continuity of the respective normal components across the edges (in R2) and faces (in R3), are the ones defined by the Raviart–Thomas polynomials since they include precisely those normal components as part of their degrees of freedom. That being said, we recall now that, given an integer k0, we define for each T𝒯h the local Raviart–Thomas space of order k as

𝐑𝐓k(T):=𝐏k(T)+Pk(T)𝐱,

where according to the terminology described in Section 1, 𝐏k(T):=[Pk(T)]n and 𝐱 is a generic vector in Rn=𝐑. Similarly 𝐂(Ω¯):=[C(Ω¯)]n denotes the space of continuous functions in Ω¯. Thus, we approximate the rate of strain, the vorticity and temperature gradient by discontinuous piecewise polynomial tensors of degree k, the pseudostress and pseudoheat by Raviart–Thomas elements of order k, and the velocity and temperature by Lagrange elements of order k, that is,

(4.9)h𝐭:={𝐬h𝕃𝚝𝚛2(Ω):𝐬h|Tk(T) for all T𝒯h},
(4.10)h𝝈:={𝝉h0(𝐝𝐢𝐯;Ω):𝐜𝚝𝝉h|T𝐑𝐓k(T) for all 𝐜𝐑,T𝒯h},
(4.11)𝐇h𝐮:={𝐯h𝐂(Ω¯):𝐯h|T𝐏k+1(T) for all T𝒯h}𝐇01(Ω),
(4.12)h𝜸:={𝜼h𝕃𝚜𝚔𝚎𝚠2(Ω):𝜼h|Tk(T) for all T𝒯h},
(4.13)𝐇h𝜻:={𝝌h𝐋2(Ω):𝝌h|T𝐏k(T) for all T𝒯h},
(4.14)𝐇h𝐩:={𝐪h𝐇N(div;Ω):𝐪h|T𝐑𝐓k(T) for all T𝒯h},
(4.15)Hhφ:={ψhC(Ω¯):ψh|TPk+1(T) for all T𝒯h}.

According to [7, 16], their corresponding approximation properties are as follows.

### $(\mathbf{AP}_{h}^{\mathbf{t}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝐭s(Ω)𝕃𝚝𝚛2(Ω), there holds

(4.16)dist(𝐭,h𝐭)Chs𝐭s,Ω.

### $(\mathbf{AP}_{h}^{\boldsymbol{\sigma}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝝈s(Ω)0(𝐝𝐢𝐯;Ω) with 𝐝𝐢𝐯𝝈𝐇s(Ω), there holds

(4.17)dist(𝝈,h𝝈)Chs{𝝈s,Ω+𝐝𝐢𝐯𝝈s,Ω}.

### $(\mathbf{AP}_{h}^{\mathbf{u}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝐮𝐇s+1(Ω)𝐇01(Ω), there holds

(4.18)dist(𝐮,𝐇h𝐮)Chs𝐮s+1,Ω.

### $(\mathbf{AP}_{h}^{\boldsymbol{\gamma}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝜸s(Ω)𝕃𝚜𝚔𝚎𝚠2(Ω), there holds

(4.19)dist(𝜸,h𝜸)Chs𝜸s,Ω.

### $(\mathbf{AP}_{h}^{\boldsymbol{\zeta}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝜻s(Ω), there holds

(4.20)dist(𝜻,𝐇h𝜻)Chs𝜻s,Ω.

### $(\mathbf{AP}_{h}^{\mathbf{p}})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each 𝐩𝐇s(Ω)𝐇N(div;Ω) with div𝐩Hs(Ω), there holds

(4.21)dist(𝐩,𝐇h𝐩)Chs{𝐩s,Ω+div𝐩s,Ω}.

### $(\mathbf{AP}_{h}^{\varphi})$.

There exists a constant C>0, independent of h, such that for each s(0,k+1], and for each φHs+1(Ω), there holds

(4.22)dist(φ,𝐇h𝐩)Chsφs+1,Ω.

## 5 A Priori Error Analysis

Let (𝐭,𝜻)×𝒬 with (𝐮,φ)𝐖r0 (r0 as in Lemma 3.5) and (𝐭h,𝜻h)h×𝒬h with (𝐮h,φh)𝐖h,r0 be the solutions to the continuous and discrete problems (3.17) and (4.1), respectively, that is,

(5.1)

(5.1a)𝐀φ(𝐭,𝐬)+𝐁𝐮(𝐭,𝐬)=Fφ(𝐬)+F𝚖(𝐬)for all 𝐬,
(5.1b)𝐀φh(𝐭h,𝐬h)+𝐁𝐮h(𝐭h,𝐬h)=Fφh(𝐬)+F𝚖(𝐬)for all 𝐬hh,

and

(5.2)

(5.2a)𝐂φ(𝜻,𝝌)+𝐃𝐮(𝜻,𝝌)=GD(𝝌)+G𝚎(𝝌)for all 𝝌𝒬,
(5.2b)𝐂φh(𝜻h,𝝌h)+𝐃𝐮h(𝜻h,𝝌h)=GD(𝝌)+G𝚎(𝝌)for all 𝝌h𝒬h.

In addition, we denote as usual

dist(𝐭,h):=inf𝐬hh𝐭-𝐬h,dist(𝜻,𝒬h):=inf𝝌h𝒬h𝜻-𝝌h.

Then the main result of this section reads as follows.

## Theorem 5.1.

Assume that the data satisfy

(5.3)C0(𝐠,Ω+𝐟𝚖ε,Ω+φD1/2+ε,ΓD+f𝚎0,Ω)<12,

where C0 is a positive constant to be specified below. Then there exists a constant C>0 depending only on parameters, data, and other constants, all of them independent of h, such that the following Céa estimate holds:

(5.4)(𝐭,𝜻)-(𝐭h,𝜻h)C{dist(𝐭,h)+dist(𝜻,𝒬h)}.

Similar to [2, Lemma 5.3] and [4, Lemma 4.2], we will apply the Strang Lemma to the pair of equations (5.1) and (5.2) separately, to then join the resulting estimates to derive the Céa estimate (5.4). In this regard, we emphasize that, given arbitrary 𝐮,𝐮h𝐇1(Ω), and φ,φhH1(Ω), particularly those forming part of the solutions to the continuous and discrete problems (3.17) and (4.1), the generic forms involved in (5.1) and (5.2), namely 𝐀φ, 𝐁𝐮, 𝐀φh, 𝐁𝐮h, 𝐂φ, 𝐃𝐮, 𝐂φh, and 𝐃𝐮h, are all bilinear, and hence the aforementioned applications of the Strang Lemma are indeed possible. We begin by recalling next from [11] and [4] this lemma and the result regarding the first pair of equations, respectively.

## Lemma 5.2 (Strang).

Let V be a Hilbert space, let FV, and let A:V×VR be a bounded and V-elliptic bilinear form. In addition, let {Vh}h>0 be a sequence of finite-dimensional subspaces of V, and for each h>0, consider a bounded bilinear form Ah:Vh×VhR and a functional FhVh. Assume that the family {Ah}h>0 is uniformly elliptic in Vh, that is, there exists a constant α~>0, independent of h, such that

Ah(vh,vh)α~vhV2for all vhVh and all h>0.

In turn, let uV and uhVh such that

A(u,v)=F(v)for all vV  𝑎𝑛𝑑  Ah(uh,vh)=F(vh)for all vhVh.

Then, for each h>0, there holds

u-uhVCST{supwhVhwh0|F(wh)-Fh(wh)|whV+infvhVhvh0(u-vhV+supwhVhwh0|A(vh,wh)-Ah(vh,wh)whV)},

where CST:=α~-1max{1,A}.

## Proof.

See [11, Theorem 11.1]. ∎

## Lemma 5.3.

Let CSTM:=2α(Ω)max{1,Aφ+Bu}, where α(Ω)2 is the ellipticity constant of Aφ+Bu (cf. (3.36)) and Aφ+Bu is its norm, which is independent of (u,φ) (cf. (3.34)). Then there holds

𝐭-𝐭hCST𝐌{(1+2𝐀φ+𝐁𝐮)dist(𝐭,h)+c1(Ω)(1+κ12)1/2𝐮1,Ω𝐮-𝐮h1,Ω
+(LμCεC~ε(1+κ12)1/2𝐭ε,Ω+(1+κ22)1/2𝐠,Ω)φ-φh1,Ω}.

## Proof.

See [4, Lemma 4.2]. ∎

## Lemma 5.4.

Let CSTE:=2β(Ω)max{1,Cφ+Du}, where β(Ω)2 is the ellipticity constant of Cφ+Du (cf. (3.39)) and Cφ+Du is its norm, which is independent of (u,φ) (cf. (3.38)). Then there holds

(5.5)

𝜻-𝜻hCST𝐄{(1+2𝐂φ+𝐃𝐮)dist(𝜻,𝒬h)+c2(Ω)(1+κ52)1/2φ1,Ω𝐮-𝐮h1,Ω
+LkCεC~ε(1+κ52)1/2𝜻ε,Ωφ-φh1,Ω}.

## Proof.

It is clear from (3.29) that 𝐂φ+𝐃𝐮 and 𝐂φh+𝐃𝐮h are bounded bilinear forms (both with constant 𝐂φ+𝐃𝐮, without loss of generality since it is independent of (𝐮,φ)) and uniformly elliptic with constant β(Ω)2. Also, GD+G𝚎 is a bounded linear functional in both and h. Therefore, a straightforward application of Lemma 5.2 to the pair (5.2) leads us to

(5.6)𝜻-𝜻hCST𝐄inf𝝃h𝒬h𝝃h𝟎{𝜻-𝝃h+sup𝝌h𝒬h𝝌h𝟎|(𝐂φ+𝐃𝐮)(𝝃h,𝝌h)-(𝐂φh+𝐃𝐮h)(𝝃h,𝝌h)|𝝌h},

where

CST𝐄:=2β(Ω)max{1,𝐂φ+𝐃𝐮}.

In what follows we make use of the boundedness constants Cε and C~ε, which correspond, respectively, to the continuous embedding of 𝐇ε(Ω) in 𝐋2p(Ω) (cf. (3.52)) and H1(Ω) in Ln/ε(Ω) (cf. proof of Lemma 3.10). Then, adding and subtracting suitable terms, using the boundedness of the bilinear forms, and the same bounding technique as in (3.51), the numerator of the last term in (5.6) can be treated as follows:

|(𝐂φ+𝐃𝐮)(𝝃h,𝝌h)-(𝐂φh+𝐃𝐮h)(𝝃h,𝝌h)||(𝐂φ+𝐃𝐮)(𝝃h-𝜻,𝝌h)+(𝐂φ-𝐂φh)(𝜻,𝝌h)+(𝐃𝐮-𝐃𝐮h)(𝜻,𝝌h)-(𝐂φh+𝐃𝐮h)(𝝃h-𝜻,𝝌h)|2𝐂φ+𝐃𝐮𝜻-𝝃h𝝌h+Ω(k(φ)-k(φh))𝜻{𝝌h-κ5𝐪h}-Ωφ(𝐮-𝐮h){𝝌h-κ5𝐪h}2𝐂φ+𝐃𝐮𝜻-𝝃h𝝌h+Lk(1+κ52)1/2CεC~ε𝜻ε,Ωφ-φh1,Ω𝝌h+c2(Ω)(1+κ52)1/2φ1,Ω𝐮-𝐮h1,Ω𝝌h,

which back into (5.6) gives (5.5), concluding this way the proof. ∎

As a result of the previous two lemmas, we have a preliminary estimate for the error:

(5.7)

(𝐭,𝜻)-(𝐭h,𝜻h)C1dist(𝐭,h)+C2dist(𝜻,𝒬h)+{C3C^ε𝐮1+ε,Ω
+C4𝐭ε,Ω+C5𝐠,Ω+C6C^εφ1+ε,Ω+C7𝜻ε,Ω}(𝐭,𝜻)-(𝐭h,𝜻h),

where

C1:=CST𝐌(1+2𝐀φ+𝐁𝐮),C2:=CST𝐄(1+2𝐂φ+𝐃φ),C3:=CST𝐌c1(Ω)(1+κ12)1/2C4:=CST𝐌LμCεC~ε(1+κ12)1/2,C5:=CST𝐌(1+κ22)1/2,C6:=CST𝐄c2(Ω)(1+κ52)1/2,C7:=CST𝐄LkCεC~ε(1+κ52)1/2,

all being positive constants independent of the discretization parameters. Thus, we first bound the terms 𝐮1+ε,Ω and 𝐭ε,Ω using the further-regularity assumption (3.44), and then we bound φ1+ε,Ω and 𝜻ε,Ω in a similar way using (3.45), whence denoting by

C8:=(C3C^ε+C4)C~𝐌(r),C9:=(C6C^ε+C7)C~𝐄(r)andC0:=max{C8r+C5,C8,C9},

inequality (5.7) becomes

(5.8)(𝐭,𝜻)-(𝐭h,𝜻h)C1dist(𝐭,h)+C2dist(𝜻,𝒬h)+C0(𝐠,Ω+𝐟𝚖ε,Ω+φD1/2+ε,ΓD+f𝚎0,Ω)(𝐭,𝜻)-(𝐭h,𝜻h),

thus leading us to the main result of this section, that is, to Theorem 5.1. In fact, thanks to (5.3), the estimate (5.8) readily implies (5.4) with C:=2max{C1,C2}.

Consequently, when using the finite element subspaces (4.9)–(4.15), the following can be established regarding the rates of convergence of the method.

## Lemma 5.5.

In addition to the hypotheses of Theorems 3.8, 4.8 and 5.1, assume that there exists s>0 such that tHs(Ω), 𝛔Hs(Ω), div𝛔Hs(Ω), uHs+1(Ω), 𝛄Hs(Ω), 𝛇Hs(Ω), pHs(Ω), divpHs(Ω) and φHs+1(Ω), and that the finite element subspaces are defined by (4.9)–(4.15). Then there exists a constant C>0, independent of h, such that

(𝐭,𝜻)-(𝐭h,𝜻h)Chmin{s,k+1}{𝐭s,Ω+𝝈s,Ω+𝐝𝐢𝐯𝝈s,Ω+𝐮s+1,Ω+𝜸s,Ω
+𝜻s,Ω+𝐩s,Ω+div𝐩s,Ω+φs+1,Ω}.

## Proof.

It follows from the Céa estimate (5.4) and the approximation properties (4.16)–(4.22). ∎

## 6 Numerical Results

We now present two examples that will illustrate the performance of the augmented fully-mixed finite element method (4.1) with the subspaces indicated in (4.9)–(4.15) on a set of quasiuniform triangulations. This means that, given an integer k0, the unknowns 𝐭, 𝜸, and 𝜻 are approximated by discontinuous piecewise polynomials of degree k, whereas 𝝈 and 𝐩 are approximated by Raviart–Thomas elements of order k, and 𝐮 and φ by continuous piecewise polynomials of degree k+1. In turn, the computational implementation is based on a FreeFem++ code (cf. [19]) and the iterative method comes straightforward from the uncoupling strategy presented in Section 4.1. Then, as a stopping criteria, we finish the algorithm when the relative error between two consecutive iterations of the complete coefficient vector measured in the discrete 2 norm is sufficiently small, this is,

𝐜𝐨𝐞𝐟𝐟m+1-𝐜𝐨𝐞𝐟𝐟m2𝐜𝐨𝐞𝐟𝐟m2<𝚝𝚘𝚕,

where tol is a specified tolerance.

Let us first define the error per variable

e(𝐭):=𝐭-𝐭h0,Ω,e(𝝈):=𝝈-𝝈h𝐝𝐢𝐯;Ω,e(𝐮):=𝐮-𝐮h1,Ω,e(p):=p-ph0,Ω,
e(𝜸):=𝜸-𝜸h0,Ω,e(𝜻):=𝜻-𝜻h0,Ω,e(𝐩):=𝐩-𝐩hdiv;Ω,e(φ):=φ-φh1,Ω,

as well as their corresponding rates of convergence

r():=log(e()e())log(hh)for all {𝐭,𝝈,𝐮,p,𝜸,𝜻,𝐩,φ},

where h and h denote two consecutive meshsizes with errors e and e.

### 6.1 Example 1: Two-Dimensional Smooth Exact Solution

In our first example, we consider Ω=(-1,1)2 and viscosity, thermal conductivity and body force given by

μ(φ)=exp(-0.25φ),
k(φ)=exp(0.25φ),
𝐠=(0,1)𝚝.

Dirichlet boundary conditions will be imposed on

ΓD:=ΓD1ΓD2,

where ΓD1:=[-1,1]×{-1}, ΓD2:=[-1,1]×{1}, whereas Neumann boundary conditions will be imposed on the remainder of the border, that is, ΓN:=Γ\ΓD. In this way, we consider boundary data φD and source terms 𝐟𝚖 and f𝚎 such that the exact solution to the Boussinesq problem is given by

𝐮(x,y)=(u1(x,y),u2(x,y))𝚝

with

u1(x,y):=2ysin(πx)sin(πy)(x2-1)+πsin(πx)cos(πy)(x2-1)(y2-1),
u2(x,y):=-2xsin(πx)sin(πy)(y2-1)-πcos(πx)sin(πy)(x2-1)(y2-1),

and

p(x,y):=y2-x2,φ(x,y):=-0.6944y4+1.6944y2.

Concerning the stabilization parameters, these are taken as pointed out at the end of Section 3.3, with κ0=12 and the bounds for the viscosity and thermal conductivity functions are estimated in

μ1=0.5,μ2=1.25,k1=0.75,k2=1.3.

In Figure 1 we display part of the solution obtained with fully-mixed finite element method using a first-order approximation and 1,409,884 DOF. Notice that not only we are able to recover the original unknowns but also to compute further variables of physical interest. In turn, Tables 1 and 2 show the convergence history for a sequence of quasi-uniform mesh refinements, thus confirming the rates of convergence predicted by Lemma 5.5, that is, when using first and second-order finite elements, and considering that the exact solution is smooth enough, whence the method converges with orders 𝒪(h) and 𝒪(h2), respectively. More precisely, while some variations appear at the first refinements, we stress that the above orders of convergence are indeed attained as the meshsize h approaches 0, as Lemma 5.5 established. This remark is certainly valid as well for the other tables reported in the present Section 6.

Table 1

Convergence history for Example 1, with a uniform mesh refinement and a first-order approximation.Here, the simulation with 1,816 DOF took ten fixed-point iterations, while the rest of them took nine iterationsto achieve a tolerance of 𝚝𝚘𝚕=10-8.

 Finite Element: ℍh𝐭 – ℍh𝝈 – 𝐇h𝐮 – ℍh𝜸 – 𝐇h𝜻 – 𝐇h𝐩 – Hhφ with k=0 DOF e⁢(𝐭) e⁢(𝝈) e⁢(𝐮) e⁢(p) e⁢(𝜸) e⁢(𝜻) e⁢(𝐩) e⁢(φ) 1,816 4.5307 15.6976 7.7257 2.3619 14.3622 0.4298 1.5650 0.3698 6,972 2.1933 7.9406 3.3074 1.0874 8.8796 0.2124 0.7864 0.1758 27,052 1.0947 3.9421 1.6088 0.5396 6.0695 0.1022 0.3809 0.0831 107,142 0.5331 2.0064 0.8011 0.2678 2.5649 0.0517 0.1962 0.0426 431,398 0.2581 0.9882 0.3925 0.1216 1.5602 0.0259 0.0970 0.0211 1,707,922 0.1271 0.4937 0.1947 0.0613 0.7027 0.0127 0.0482 0.0105 h r⁢(𝐭) r⁢(𝝈) r⁢(𝐮) r⁢(p) r⁢(𝜸) r⁢(𝜻) r⁢(𝐩) r⁢(φ) 0.4129 – – – – – – – – 0.1940 0.9605 0.9023 1.1233 1.0270 0.6366 0.9331 0.9111 0.9845 0.0995 1.0402 1.0482 1.0788 1.0490 0.5695 1.0951 1.0851 1.1221 0.0527 1.1333 1.0638 1.0982 1.1032 1.3567 1.0736 1.0453 1.0521 0.0311 1.3744 1.3419 1.3519 1.4961 0.9419 1.3087 1.3350 1.3342 0.0150 0.9729 0.9535 0.9635 0.9411 1.0960 0.9778 0.9600 0.9503
Table 2

Convergence history for Example 1, with a uniform mesh refinement and a second-order approximation.Here, all the simulations took nine fixed-point iterations to achieve a tolerance of 𝚝𝚘𝚕=10-8.

 Finite Element: ℍh𝐭 – ℍh𝝈 – 𝐇h𝐮 – ℍh𝜸 – 𝐇h𝜻 – 𝐇h𝐩 – Hhφ with k=1 DOF e⁢(𝐭) e⁢(𝝈) e⁢(𝐮) e⁢(p) e⁢(𝜸) e⁢(𝜻) e⁢(𝐩) e⁢(φ) 5,812 0.9568 3.0676 1.3170 0.6622 2.2406 0.0690 0.2052 0.0434 22,564 0.2193 0.7735 0.3002 0.1932 0.6699 0.0154 0.0467 0.0098 88,036 0.0554 0.1939 0.0750 0.0472 0.2187 0.0038 0.0118 0.0023 349,660 0.0140 0.0494 0.0192 0.0119 0.0494 0.0010 0.0031 0.0006 1,409,884 0.0035 0.0121 0.0047 0.0029 0.0152 0.0002 0.0008 0.0001 h r⁢(𝐭) r⁢(𝝈) r⁢(𝐮) r⁢(p) r⁢(𝜸) r⁢(𝜻) r⁢(𝐩) r⁢(φ) 0.4129 – – – – – – – – 0.1940 1.9507 1.8241 1.9580 1.6310 1.5987 1.9811 1.9595 1.9687 0.0995 2.0591 2.0712 2.0759 2.1084 1.6752 2.1123 2.0631 2.1848 0.0527 2.1678 2.1538 2.1490 2.1666 2.3443 2.0727 2.1265 2.1000 0.0311 2.6452 2.6677 2.6506 2.6598 2.2355 2.6761 2.6587 2.6447
Figure 1

Numerical Results for Example 1. From left to right and from up to down: XX-component of the Rate of Strain tensor, true vorticity magnitude (computed as twice the YX-component of the full vorticity tensor), Y-component of the temperature gradient, XX and YY components of the pseudostress tensor, pseudoheat magnitude and its vector field, velocity magnitude and its vector field, pressure and temperature. Snapshots obtained from a simulation with 1,409,884 DOF and a second-order approximation.

### 6.2 Example 2: Three-Dimensional Smooth Exact Solution

In our second example, we consider Ω=(0,1)3 and viscosity, thermal conductivity and body force given by

μ(φ)=2.0-0.5φ2-0.5φ4,
k(φ)=-0.5+2.0μ(φ),
𝐠=(0,0,1)𝚝.

With respect to boundary conditions, we impose Dirichlet conditions on [0,1]2×{0}=:ΓD and Neumann conditions on the rest of the boundary, that is, ΓN:=Γ\ΓD. Then we consider boundary data φD and source terms 𝐟𝚖 and f𝚎 such that the exact solution to the Boussinesq problem is given by

𝐮(x,y,z)=(u1(x,y,z),u2(x,y,z),u3(x,y,z))𝚝

with

u1(x,y,z):=8x2yz(x-1)2(y-1)(z-1)(y-z),
u2(x,y,z):=-8xy2z(x-1)(y-1)2(z-1)(x-z),
u3(x,y,z):=8xyz2(x-1)(y-1)(z-1)2(x-y)

and

φ(x,y,z):=sin(πx)2sin(πy)2(z-1)2,p(x,y,z):=(x-0.5)3sin(y+z).

Concerning the stabilization parameters, we take them again as in Section 3.3, but this time with κ0=1. In turn, the bounds for the viscosity and thermal conductivity are estimated in

μ1=1.0,μ2=2.0,k1=1.5,k2=3.5.

Part of the solution is displayed in Figure 2, and a convergence history for a set of quasi-uniform mesh refinements is provided in Table 3, thus showing also that, having the problem a smooth exact solution, this fully-mixed finite element method converges optimally with order 𝒪(h) (when using a first order element).

Figure 2

Numerical Results for Example 2. From left to right and from up to down: XY and ZZ components of the rate of strain tensor (the last one postprocessed as -th,11-th,22), YY component of the pseudostress tensor, velocity vector field, flow streamlines, pressure, Z components of the temperature gradient and pseudoheat, and temperature. Snapshots obtained from a simulation with 3,903,877 DOF and a first-order approximation.

Table 3

Convergence history for Example 2, with a uniform mesh refinement and a first-order approximation.Here, to achieve a tolerance of 𝚝𝚘𝚕=10-8, a number of six fixed-point iterations were needed for the first simulation,seven iterations for the second one, and eight for the third, fourth and fifth simulations.

 Finite Element: ℍh𝐭 – ℍh𝝈 – 𝐇h𝐮 – ℍh𝜸 – 𝐇h𝜻 – 𝐇h𝐩 – Hhφ with k=0 DOF e⁢(𝐭) e⁢(𝝈) e⁢(𝐮) e⁢(p) e⁢(𝜸) e⁢(𝜻) e⁢(𝐩) e⁢(φ) 1,117 0.0219 0.2559 0.0420 0.0249 0.0262 0.7239 17.6230 1.2979 8,181 0.0128 0.1367 0.0265 0.0176 0.0196 0.4249 8.0953 0.6128 62,821 0.0079 0.0700 0.0140 0.0097 0.0132 0.2240 4.3469 0.3044 492,741 0.0042 0.0351 0.0071 0.0047 0.0077 0.1137 2.1971 0.1505 3,903,877 0.0022 0.0175 0.0035 0.0022 0.0040 0.0571 1.1017 0.0750 h r⁢(𝐭) r⁢(𝝈) r⁢(𝐮) r⁢(p) r⁢(𝜸) r⁢(𝜻) r⁢(𝐩) r⁢(φ) 0.7071 – – – – – – – – 0.3536 0.7778 0.9044 0.6636 0.4972 0.4220 0.7689 1.1223 1.0828 0.1768 0.7024 0.9661 0.9205 0.8660 0.5634 0.9235 0.8971 1.0094 0.0884 0.9094 0.9951 0.9871 1.0518 0.7899 0.9788 0.9844 1.0165 0.0442 0.9513 1.0018 1.0000 1.0742 0.9191 0.9943 0.9958 1.0044

Funding statement: This research was partially supported by CONICYT-Chile through the project AFB170001 of the PIA Program: Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal; and by Centro de Investigación en Ingeniería Matemática (CI2MA), Universidad de Concepción.

## References

[1] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, 2nd ed., Pure Appl. Math. (Amsterdam) 140, Elsevier, Amsterdam, 2003. Search in Google Scholar

[2] J. A. Almonacid, G. N. Gatica and R. Oyarzúa, A mixed–primal finite element method for the Boussinesq problem with temperature-dependent viscosity, Calcolo 55 (2018), no. 3, Article ID 36. 10.1007/s10092-018-0278-zSearch in Google Scholar

[3] J. A. Almonacid, G. N. Gatica and R. Oyarzúa, A posteriori error analysis of a mixed-primal finite element method for the Boussinesq problem with temperature-dependent viscosity, J. Sci. Comput. 78 (2019), no. 2, 887–917. 10.1007/s10915-018-0810-ySearch in Google Scholar

[4] J. A. Almonacid, G. N. Gatica, R. Oyarzúa and R. Ruiz-Baier, A new mixed finite element method for the n-dimensional Boussinesq problem with temperature-dependent viscosity, Preprint 2018-18, Centro de Investigacion en Ingenieria Matematica (CI2MA), Universidad de Concepcion, Chile, 2018. 10.1007/s10092-018-0278-zSearch in Google Scholar

[5] M. Alvarez, G. N. Gatica and R. Ruiz-Baier, An augmented mixed-primal finite element method for a coupled flow-transport problem, ESAIM Math. Model. Numer. Anal. 49 (2015), no. 5, 1399–1427. 10.1051/m2an/2015015Search in Google Scholar

[6] J. Boland and W. Layton, An analysis of the finite element method for natural convection problems, Numer. Methods Partial Differential Equations 6 (1990), no. 2, 115–126. 10.1002/num.1690060202Search in Google Scholar

[7] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Ser. Comput. Math. 15, Springer, New York, 1991. 10.1007/978-1-4612-3172-1Search in Google Scholar

[8] J. Camaño, G. N. Gatica, R. Oyarzúa and R. Ruiz-Baier, An augmented stress-based mixed finite element method for the steady state Navier–Stokes equations with nonlinear viscosity, Numer. Methods Partial Differential Equations 33 (2017), no. 5, 1692–1725. 10.1002/num.22166Search in Google Scholar

[9] S. Caucao, G. N. Gatica and R. Oyarzúa, Analysis of an augmented fully-mixed formulation for the coupling of the Stokes and heat equations, ESAIM Math. Model. Numer. Anal. 52 (2018), no. 5, 1947–1980. 10.1051/m2an/2018027Search in Google Scholar

[10] P. G. Ciarlet, Linear and Nonlinear Functional Analysis with Applications, SIAM, Philadelphia, 2013. Search in Google Scholar

[11] P. G. Ciarlet and J.-L. Lions, Handbook of Numerical Analysis. Vol. II. Finite Element Methods. Part 1, North-Holland, Amsterdam, 1991. Search in Google Scholar

[12] E. Colmenares, G. N. Gatica and R. Oyarzúa, Analysis of an augmented mixed-primal formulation for the stationary Boussinesq problem, Numer. Methods Partial Differential Equations 32 (2016), no. 2, 445–478. 10.1002/num.22001Search in Google Scholar

[13] E. Colmenares, G. N. Gatica and R. Oyarzúa, An augmented fully-mixed finite element method for the stationary Boussinesq problem, Calcolo 54 (2017), no. 1, 167–205. 10.1007/s10092-016-0182-3Search in Google Scholar

[14] M. Farhloul, S. Nicaise and L. Paquet, A mixed formulation of Boussinesq equations: Analysis of nonsingular solutions, Math. Comp. 69 (2000), no. 231, 965–986. 10.1090/S0025-5718-00-01186-8Search in Google Scholar

[15] G. N. Gatica, An augmented mixed finite element method for linear elasticity with non-homogeneous Dirichlet conditions, Electron. Trans. Numer. Anal. 26 (2007), 421–438. Search in Google Scholar

[16] G. N. Gatica, A Simple Introduction to the Mixed Finite Element Method. Theory and Applications, Springer Briefs Math., Springer, Cham, 2014. 10.1007/978-3-319-03695-3Search in Google Scholar

[17] G. N. Gatica, A. Márquez and M. A. Sánchez, A priori and a posteriori error analyses of a velocity-pseudostress formulation for a class of quasi-Newtonian Stokes flows, Comput. Methods Appl. Mech. Engrg. 200 (2011), no. 17–20, 1619–1636. 10.1016/j.cma.2011.01.010Search in Google Scholar

[18] G. N. Gatica and F. A. Sequeira, A priori and a posteriori error analyses of an augmented HDG method for a class of quasi-Newtonian Stokes flows, J. Sci. Comput. 69 (2016), no. 3, 1192–1250. 10.1007/s10915-016-0233-6Search in Google Scholar

[19] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (2012), no. 3–4, 251–265. 10.1515/jnum-2012-0013Search in Google Scholar

[20] A. Kufner, O. Jhon and S. Fučík, Function Spaces, Noordhoff International, Leyden, 1977. Search in Google Scholar

[21] W. McLean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University, Cambridge, 2000. Search in Google Scholar

[22] R. Oyarzúa, T. Qin and D. Schötzau, An exactly divergence-free finite element method for a generalized Boussinesq problem, IMA J. Numer. Anal. 34 (2014), no. 3, 1104–1135. 10.1093/imanum/drt043Search in Google Scholar

[23] R. Oyarzúa and P. Zúñiga, Analysis of a conforming finite element method for the Boussinesq problem with temperature-dependent parameters, J. Comput. Appl. Math. 323 (2017), 71–94. 10.1016/j.cam.2017.04.009Search in Google Scholar

[24] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer Ser. Comput. Math. 23, Springer, Berlin, 1994. 10.1007/978-3-540-85268-1Search in Google Scholar

[25] P. W. Schroeder and G. Lube, Stabilised dG-FEM for incompressible natural convection flows with boundary and moving interior layers on non-adapted meshes, J. Comput. Phys. 335 (2017), 760–779. 10.1016/j.jcp.2017.01.055Search in Google Scholar

[26] M. Tabata and D. Tagami, Error estimates of finite element methods for nonstationary thermal convection problems with temperature-dependent coefficients, Numer. Math. 100 (2005), no. 2, 351–372. 10.1007/s00211-005-0589-2Search in Google Scholar

[27] J. Wu, J. Shen and X. Feng, Unconditionally stable gauge-Uzawa finite element schemes for incompressible natural convection problems with variable density, J. Comput. Phys. 348 (2017), 776–789. 10.1016/j.jcp.2017.07.045Search in Google Scholar