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

Javier A. Almonacid; Gabriel N. Gatica

doi:10.1515/cmam-2018-0187

Article Publicly Available

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

Javier A. Almonacid and Gabriel N. Gatica

Published/Copyright: April 9, 2019

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 20 Issue 2

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.

Keywords: Boussinesq Equations; Augmented Fully-Mixed Formulation; Fixed-point Theory; Finite Element Methods; A Priori Error Analysis

MSC 2010: 65N30; 65N12; 65N15; 35Q79; 80A20; 76R05; 76D07

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

∇⁡𝐯:=(∂⁡vi∂⁡xj)i,j=1,…,n,div⁡𝐯:=∑j=1n∂⁡vj∂⁡xj,and 𝐯⊗𝐰:=(vi⁢wj)i,j=1,…,n.

In addition, for any tensor fields 𝝉=(τi⁢j)i,j=1,…,n and 𝜻=(ζi⁢j)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

𝝉𝚝:=(τj⁢i)i,j=1,…,n,tr⁢(𝝉):=∑i=1nτi⁢i,𝝉:𝜻:=∑i,j=1nτi⁢j⁢ζi⁢j,and 𝝉𝚍:=𝝉-1n⁢tr⁢(𝝉)⁢𝕀,

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, φD∈H1/2⁢(ΓD) is a prescribed temperature and μ,k:R→R+ 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 k2≥k1>0 such that μ1≤μ⁢(t)≤μ2 and k1≤k⁢(t)≤k2 for all t∈R. 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,t∈R. 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)∫Ω𝐭:𝝉𝚍+∫Ω𝜸:𝝉+∫Ω𝐮⋅𝐝𝐢𝐯⁡𝝉=0 ⁢for all ⁢𝝉∈ℍ0⁢(𝐝𝐢𝐯;Ω),

(3.2)∫Ωμ⁢(φ)⁢𝐭:𝐬-∫Ω(𝐮⊗𝐮)𝚍:𝐬-∫Ω𝝈𝚍:𝐬=0 ⁢for 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〉ΓD for 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𝚎⁢{ψ-κ6⁢div⁡𝐪}.

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⁢(φ)⁢𝐩~-∇⁡φ=0 ⁢in ⁢Ω, -div⁡𝐩~+𝐮⋅∇⁡φ=f𝚎⁢ ⁢in ⁢Ω

φ=φD ⁢on ⁢ΓD, ⁢𝐩~⋅𝝂=0⁢ ⁢on ⁢Γ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 ∥w∥1,Ω≤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⁢δ4⁢k1k2),κ6,κ8∈(0,∞),κ7∈(0,2⁢δ5⁢(k1-κ5⁢k22⁢δ4)).

Then there exists r2>0 such that for each r∈(0,r2) and for each (w,ϕ)∈H satisfying ∥w∥1,Ω≤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𝐄⁢{∥φD∥1/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:Wr→H is well defined and there holds

(3.33)∥𝐓⁢(𝐰,ϕ)∥≤C𝐌⁢{r⁢∥𝐠∥∞,Ω+∥𝐟𝚖∥0,Ω}+C𝐄⁢{∥φD∥1/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/2 for all ⁢𝐭→∈ℋ,

and

∥𝜻→∥:={∥𝜻∥0,Ω2+∥𝐪∥div;Ω2+∥ψ∥1,Ω2}1/2 for 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,Ω2 for 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,Ω2≥c3⁢(Ω)⁢∥𝝉∥0,Ω2 for 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⁡{α2⁢c3⁢(Ω),κ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

𝐀ϕ⁢(𝐬→,𝐬→)≥α⁢(Ω)⁢∥𝐬→∥2 for 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,Ω)⁢∥𝐬→∥2 for all ⁢𝐬→∈ℋ,

and then, we easily see that

(3.36)(𝐀ϕ+𝐁𝐰)⁢(𝐬→,𝐬→)≥α⁢(Ω)2⁢∥𝐬→∥2 for all ⁢𝐬→∈ℋ

provided that

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

that is,

(3.37)∥𝐰∥1,Ω≤α⁢(Ω)2⁢c1⁢(Ω)⁢(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,ΓD2≥c4⁢(Ω)⁢∥ψ∥1,Ω2 for 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,Ω

+κ6⁢∥div⁡𝐩∥0,Ω⁢∥div⁡𝐪∥0,Ω+κ7⁢|φ|1,Ω⁢|ψ|1,Ω+κ8⁢c0⁢(Ω)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+κ6⁢∥div⁡𝐪∥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-κ5⁢k22⁢δ4⁢∥𝝌∥0,Ω2-κ5⁢k2⁢δ42⁢∥𝐪∥0,Ω2-κ72⁢δ5⁢∥𝝌∥0,Ω2-κ7⁢δ52⁢|ψ|1,Ω2

+κ5⁢∥𝐪∥0,Ω2+κ6⁢∥div⁡𝐪∥0,Ω2+κ7⁢|ψ|0,Ω2+κ8⁢∥ψ∥0,Γ2

=(k1-κ5⁢k22⁢δ4-κ72⁢δ5)⁢∥𝝌∥0,Ω2+κ5⁢(1-k2⁢δ42)⁢∥𝐪∥0,Ω2+κ6⁢∥div⁡𝐪∥0,Ω2

+κ7⁢(1-δ52)⁢|ψ|1,Ω2+κ8⁢∥ψ∥0,ΓD2.

Hence, applying Lemma 3.6 and defining the constants

β1:=k1-κ5⁢k22⁢δ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⁢∥𝝌→∥2 for all ⁢𝝌→∈𝒬,

provided that

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

that is,

(3.40)∥𝐰∥1,Ω≤β⁢(Ω)2⁢c2⁢(Ω)⁢(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,Ω⁢∥φD∥1/2,ΓD+κ8⁢∥φD∥0,ΓD⁢∥ψ∥0,ΓD

and

|G𝚎⁢(𝝌→)|≤(1+κ62)1/2⁢∥f𝚎∥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+κ8⁢c0⁢(Ω)⁢C1/2,(1+κ62)1/2}

such that

|(GD+G𝚎)⁢(𝝌→)|≤CG⁢{∥φD∥1/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𝐄:=2⁢CGβ⁢(Ω). ∎

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=k12⁢k22,κ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:Wr→H 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𝐄⁢(∥φD∥1/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𝐄⁢{∥φD∥1/2,ΓD+∥f𝚎∥0,Ω}≤r,

(3.42b)L𝐌⁢(∥𝐠∥∞,Ω+∥𝐟𝚖∥ε,Ω)<2-1,

(3.42c)L𝐄⁢(∥φD∥1/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𝐄⁢{∥φD∥1/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𝐄⁢{∥φD∥1/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 φD∈H1/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)⁢{∥φD∥1/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-ϕ2∥Ln/ε⁢(Ω)

+∥𝐌3(𝐰1,ϕ1)∥1,Ω∥𝐰1-𝐰2∥1,Ω+∥𝐠∥∞,Ω∥ϕ1-ϕ2∥0,Ω}

for all (w1,ϕ1),(w2,ϕ2)∈H such that ∥w1∥1,Ω,∥w2∥1,Ω≤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^1⁢r,

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𝐄⁢{∥φD∥1/2+ε,ΓD+∥f𝚎∥0,Ω}⁢∥(𝐰1,ϕ1)-(𝐰2,ϕ2)∥

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

Proof.

Let (𝐰1,ϕ1),(𝐰2,ϕ2)∈𝐇 such that ∥𝐰1∥1,Ω,∥𝐰2∥1,Ω≤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-𝜻→2∥2≤(𝐂ϕ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⁢∥φ1∥1,Ω⁢∥𝐰1-𝐰2∥1,Ω⁢∥𝜻→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)⁢𝜻1∥0,Ω⁢∥𝜻→1-𝜻→2∥≤Lk⁢(1+κ52)1/2⁢∥ϕ2-ϕ1∥L2⁢q⁢(Ω)⁢∥𝜻1∥𝐋2⁢p⁢(Ω)⁢∥𝜻→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 L2⁢p⁢(Ω) whenever

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

and hence, there exists Cε>0 such that

(3.52)∥𝜻∥𝐋2⁢p⁢(Ω)≤Cε⁢∥𝜻∥ε,Ω for all ⁢𝜻∈𝐇ε⁢(Ω).

In this way,

2⁢q=2⁢pp-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-ϕ2∥Ln/ε⁢(Ω)⁢∥𝜻→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-ϕ2∥Ln/ε⁢(Ω)+∥𝐄3⁢(𝐰1,ϕ1)∥1,Ω⁢∥𝐰1-𝐰2∥1,Ω}

for any (𝐰1,ϕ1),(𝐰2,ϕ2)∈𝐇 such that ∥𝐰1∥1,Ω,∥𝐰2∥1,Ω≤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⁢(∥φD∥1/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𝐄⁢(∥φD∥1/2+ε,ΓD+∥f𝚎∥0,Ω)⁢∥(𝐌3⁢(𝐰1,ϕ1),ϕ1)-(𝐌3⁢(𝐰2,ϕ2),ϕ2)∥

≤L𝐄⁢(∥φD∥1/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 ⁢𝐬→h∈ℋh,

(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:𝐇h→ℋh 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 𝐭→h∈ℋh is the solution to the problem: Find 𝐭→h∈ℋh such that

(4.2)𝐀ϕh⁢(𝐭→h,𝐬→h)+𝐁𝐰h⁢(𝐭→h,𝐬→h)=Fϕh⁢(𝐬→h)+F𝚖⁢(𝐬→h) for all ⁢𝐬→h∈ℋh.

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:W→W 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 t→h:=Mh⁢(wh,ϕh)∈Hh for each (wh,ϕh)∈Hh such that ∥wh∥1,Ω≤r. Moreover, there holds

∥𝐌h⁢(𝐰h,ϕh)∥=∥𝐭→h∥≤C𝐌⁢{∥𝐠∥∞,Ω⁢∥ϕh∥0,Ω+∥𝐟𝚖∥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⁢δ4⁢k1k2),κ6,κ8∈(0,∞),κ7∈(0,2⁢δ5⁢(k1-κ5⁢k22⁢δ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 ∥wh∥1,Ω≤r. Moreover, there holds

∥𝐄h⁢(𝐰h,ϕh)∥=∥𝜻→h∥≤C𝐄⁢{∥φD∥1/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,r→Hh is well defined and there holds

∥𝐓h⁢(𝐰h,ϕh)∥≤C𝐌⁢{r⁢∥𝐠∥∞,Ω+∥𝐟𝚖∥0,Ω}+C𝐄⁢{∥φD∥1/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-ϕh2∥L4⁢(Ω)

(4.6)+∥𝐌3,h(𝐰h1,ϕh1)∥1,Ω∥𝐰h1-𝐰h2∥1,Ω+∥𝐠∥∞,Ω∥ϕh1-ϕh2∥0,Ω}

for all (wh1,ϕh1),(wh2,ϕh2)∈Hh such that ∥wh1∥1,Ω,∥wh2∥1,Ω≤r.

Proof.

It comes from [4, Lemma 2.6], but employing an L4-L4-L2 argument, instead of the L2⁢p-L2⁢q-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-ϕh2∥L4⁢(Ω)+∥𝐄3,h⁢(𝐰h1,ϕh1)∥1,Ω⁢∥𝐰h1-𝐰h2∥1,Ω}

for all (wh1,ϕh1),(wh2,ϕh2)∈Hh such that ∥wh1∥1,Ω,∥wh2∥1,Ω≤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,r→Hh 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-ϕh2∥1,Ω

+(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=Ci⁢C¯𝐄.∎

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 (t→h,𝛇→h)∈Hh×Qh with (uh,φh)∈Wh,r. Moreover, there hold

∥𝐭→h∥≤C𝐌⁢{r⁢∥𝐠∥∞,Ω+∥𝐟𝚖∥0,Ω},

and

∥𝜻→h∥≤C𝐄⁢{∥φD∥1/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 k≥0, 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|T∈ℙk⁢(T)⁢ for all ⁢T∈𝒯h},

(4.10)ℍh𝝈:={𝝉h∈ℍ0⁢(𝐝𝐢𝐯;Ω):𝐜𝚝⁢𝝉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|T∈ℙk⁢(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φ:={ψh∈C⁢(Ω¯):ψh|T∈Pk+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𝐭)≤C⁢hs⁢∥𝐭∥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𝝈)≤C⁢hs⁢{∥𝝈∥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𝐮)≤C⁢hs⁢∥𝐮∥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𝜸)≤C⁢hs⁢∥𝜸∥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𝜻)≤C⁢hs⁢∥𝜻∥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𝐩)≤C⁢hs⁢{∥𝐩∥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𝐩)≤C⁢hs⁢∥φ∥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 ⁢𝐬→h∈ℋh,

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𝐬→h∈ℋh⁡∥𝐭→-𝐬→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⁢(∥𝐠∥∞,Ω+∥𝐟𝚖∥ε,Ω+∥φD∥1/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 φ,φh∈H1⁢(Ω), 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 F∈V′, and let A:V×V→R 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×Vh→R and a functional Fh∈Vh′. 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)≥α~⁢∥vh∥V2 for all ⁢vh∈Vh⁢ and all ⁢h>0.

In turn, let u∈V and uh∈Vh such that

A⁢(u,v)=F⁢(v) for all ⁢v∈V 𝑎𝑛𝑑 Ah⁢(uh,vh)=F⁢(vh) for all ⁢vh∈Vh.

Then, for each h>0, there holds

∥u-uh∥V≤CST⁢{supwh∈Vhwh≠0⁡|F⁢(wh)-Fh⁢(wh)|∥wh∥V+infvh∈Vhvh≠0⁡(∥u-vh∥V+supwh∈Vhwh≠0⁡|A(vh,wh)-Ah(vh,wh)∥wh∥V)},

where CST:=α~-1⁢max⁡{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

∥𝐭→-𝐭→h∥≤CST𝐌{(1+2∥𝐀φ+𝐁𝐮∥)dist(𝐭→,ℋh)+c1(Ω)(1+κ12)1/2∥𝐮∥1,Ω∥𝐮-𝐮h∥1,Ω

+(LμCεC~ε(1+κ12)1/2∥𝐭∥ε,Ω+(1+κ22)1/2∥𝐠∥∞,Ω)∥φ-φh∥1,Ω}.

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)

∥𝜻→-𝜻→h∥≤CST𝐄{(1+2∥𝐂φ+𝐃𝐮∥)dist(𝜻→,𝒬h)+c2(Ω)(1+κ52)1/2∥φ∥1,Ω∥𝐮-𝐮h∥1,Ω

+LkCεC~ε(1+κ52)1/2∥𝜻∥ε,Ω∥φ-φh∥1,Ω}.

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)∥𝜻→-𝜻→h∥≤CST𝐄⁢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 𝐋2⁢p⁢(Ω) (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/2⁢Cε⁢C~ε⁢∥𝜻∥ε,Ω⁢∥φ-φh∥1,Ω⁢∥𝝌→h∥+c2⁢(Ω)⁢(1+κ52)1/2⁢∥φ∥1,Ω⁢∥𝐮-𝐮h∥1,Ω⁢∥𝝌→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𝐄⁢Lk⁢Cε⁢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:=(C3⁢C^ε+C4)⁢C~𝐌⁢(r),C9:=(C6⁢C^ε+C7)⁢C~𝐄⁢(r) and C0:=max⁡{C8⁢r+C5,C8,C9},

inequality (5.7) becomes

(5.8)∥(𝐭→,𝜻→)-(𝐭→h,𝜻→h)∥≤C1⁢dist⁡(𝐭→,ℋh)+C2⁢dist⁡(𝜻→,𝒬h)+C0⁢(∥𝐠∥∞,Ω+∥𝐟𝚖∥ε,Ω+∥φD∥1/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:=2⁢max⁡{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 t∈Hs⁢(Ω), 𝛔∈Hs⁢(Ω), div⁡𝛔∈Hs⁢(Ω), u∈Hs+1⁢(Ω), 𝛄∈Hs⁢(Ω), 𝛇∈Hs⁢(Ω), p∈Hs⁢(Ω), div⁡p∈Hs⁢(Ω) 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 k≥0, 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-𝐜𝐨𝐞𝐟𝐟m∥ℓ2∥𝐜𝐨𝐞𝐟𝐟m∥ℓ2<𝚝𝚘𝚕,

where tol is a specified tolerance.

Let us first define the error per variable

e⁢(𝐭):=∥𝐭-𝐭h∥0,Ω, ⁢e⁢(𝝈):=∥𝝈-𝝈h∥𝐝𝐢𝐯;Ω, ⁢e⁢(𝐮):=∥𝐮-𝐮h∥1,Ω, ⁢e⁢(p):=∥p-ph∥0,Ω,

e⁢(𝜸):=∥𝜸-𝜸h∥0,Ω, ⁢e⁢(𝜻):=∥𝜻-𝜻h∥0,Ω, ⁢e⁢(𝐩):=∥𝐩-𝐩h∥div;Ω, ⁢e⁢(φ):=∥φ-φh∥1,Ω,

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):=2⁢y⁢sin⁡(π⁢x)⁢sin⁡(π⁢y)⁢(x2-1)+π⁢sin⁡(π⁢x)⁢cos⁡(π⁢y)⁢(x2-1)⁢(y2-1),

u2⁢(x,y):=-2⁢x⁢sin⁡(π⁢x)⁢sin⁡(π⁢y)⁢(y2-1)-π⁢cos⁡(π⁢x)⁢sin⁡(π⁢y)⁢(x2-1)⁢(y2-1),

and

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

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.

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):=8⁢x2⁢y⁢z⁢(x-1)2⁢(y-1)⁢(z-1)⁢(y-z),

u2⁢(x,y,z):=-8⁢x⁢y2⁢z⁢(x-1)⁢(y-1)2⁢(z-1)⁢(x-z),

u3⁢(x,y,z):=8⁢x⁢y⁢z2⁢(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{-t_{h,11}-t_{h,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.$

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. 10.1137/1.9781611972597Search 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

Received: 2018-07-26

Revised: 2019-03-07

Accepted: 2019-03-14

Published Online: 2019-04-09

Published in Print: 2020-04-01

Articles in the same Issue

https://doi.org/10.1515/cmam-2018-0187

Keywords for this article

Boussinesq Equations; Augmented Fully-Mixed Formulation; Fixed-point Theory; Finite Element Methods; A Priori Error Analysis