On some difference schemes and entropy conditions for a class of multi-species kinematic flow models with discontinuous flux

1.   Introduction

The description of a variety of natural phenomena and engineering problems deal with incompressible quasi-Newtonian flows with viscous heating and buoyancy terms (often referred to as natural convection of fluids). Mantle convection with very large viscosities, waves and currents near shorelines, heat transfer in nanoparticle fluids, creeping thermal plumes, stratified oceanic flows, chemical reactors, and many other examples can be invoked in the context of applicative problems, partial differential equations, and in the construction of numerical schemes [17,32,26,18,22,30,33,40]. In particular, finite element methods approximating the solution of these equations have been developed by the mathematical community in the last two decades. For instance, the model with constant coefficients has been addressed via primal approaches in [13,11,19,23], whereas mixed-type schemes have been employed in [20,24,25], and in particular two different formulations based on a dual-mixed approach for the momentum equation, and a primal and mixed-primal one for the energy equation, are proposed in [21].

Here we advocate to the study of well-posedness and stability of a weak formulation (as well as the construction of mixed finite element schemes) for the Boussinesq equations with thermally-dependent viscosity. In this regard, we remark that not many finite element methods that provide analysis for the case of non-constant viscosity are available in the literature (we may basically refer to [2,3,5,41,34,36,39,42,43], and some of the references therein). For instance, the authors in [36] (see [35] for the continuous analysis) propose an optimally convergent, primal finite element method for the steady-state problem applied to the flow in channels. In turn, [39] (and a stabilized version of their method, which is suggested in [42]) deal with the unsteady problem, where backward Euler discretization is used in time, and conforming finite elements in space. On the other hand, a conforming finite element method for the problem with temperature-dependent viscosity and thermal conductivity is developed in [34]. In this work, Stokes-stable elements for the velocity and pressure, as well as Lagrange elements for the temperature, and discontinuous piecewise polynomials for the normal heat flux through the boundary, are used to define the associated discrete scheme. Similar hypotheses to those in [34] are also required in primal formulations that may be computationally less expensive than a mixed method, but do not approximate directly other variables of physical interest. The user of a primal scheme would usually appeal to numerical differentiation, with the consequent loss of accuracy.

Furthermore, just a couple of the above mentioned references consider dual-mixed approaches. In particular, in the recent contribution [5] the authors construct an augmented mixed-primal finite element method for the steady Boussinesq problem without dissipation, restricting the analysis to two-dimensional bounded domains with polygonal boundary. More precisely, in the present case one seeks a velocity field ${\bf u}$, a pressure field $p$ and a temperature field ${\varphi}$ such that

where $\Omega \subset R^{n}$, $n\in\{2,3\}$, is a domain with boundary $\Gamma : = \partial\Omega$, $\mathbf{div}$ is the usual divergence operator $\mathrm{div}$ acting along the rows of a given tensor, ${\bf e}({\bf u})$ denotes the strain rate tensor (symmetric part of the velocity gradient tensor $\nabla {\bf u}$), $-{\bf g} \in {\bf L}^\infty(\Omega) : = {[L^\infty(\Omega)]^n}$ is a body force per unit mass (e.g., gravity), ${\mathbb{K}} \in {\mathbb{L}}^\infty(\Omega) : = {[L^\infty(\Omega)]^{n\times n}}$ is a uniformly positive definite tensor describing thermal conductivity and $\mu:R \to R^+$ is a temperature-dependent viscosity, which is assumed bounded and Lipschitz continuous, that is, there exist constants $\mu_2 \geq \mu_1 > 0$ and $L_\mu > 0$ such that

The viscosity can be either of Reynolds, Arrhenius, or Sutherland's type, see e.g. [32,35,33]. Note that a viscosity dissipation term ${{\mathit{\boldsymbol{\sigma}}}}:\nabla \boldsymbol{u}$ and the power of heat production should also be present in the thermal energy equation. However these contributions are assumed much smaller than heat conduction effects (see e.g. [17]). In the Boussinesq approximation it is also assumed that density fluctuations are linearly related to temperature, and so the last term on the momentum balance is the buoyancy contribution that contains a rescaling implying that, at least for the analysis of the problem, the adimensional number Ra$/$(Pr Re$^2$) is of the order of the gravity magnitude. In addition, and for sake of notational simplicity of the analysis, here we also take a zero reference temperature and assume that other constants can be absorbed by the model adimensionalization.

With respect to the boundary conditions for (1.1), we assume that ${\bf u}_D \in {\bf H}^{1/2}(\Gamma) : = {[H^{1/2}(\Gamma)]^n}$, ${\varphi}_D \in H^{1/2}(\Gamma)$, and that ${\bf u}_D$ verifies the compatibility condition

where ${{\mathit{\boldsymbol{\nu}}}}$ denotes the unit outward normal on $\Gamma$. The construction of the mixed-primal formulation in [5] begins with the introduction of the stress and vorticity tensors, respectively, defined as

where $\otimes$ stands for the tensor product operator, ${\mathbb{I}}$ is the identity matrix in $R^{n\times n}$, the symbol $^{\tt t}$ denotes matrix transpose, and

Note that if we had Neumann boundary conditions in the fluid, then $\boldsymbol\sigma {{\mathit{\boldsymbol{\nu}}}}$ would be prescribed on $\Gamma$. Therefore, after eliminating the pressure $p$, problem (1.1) is rewritten as: Find $({{\mathit{\boldsymbol{\sigma}}}},{\bf u},{{\mathit{\boldsymbol{\gamma}}}}, {\varphi})$ such that

where $\mathrm{tr}$ denotes the matrix trace, $\boldsymbol{\tau}^{\tt d} : = \boldsymbol{\tau} - \frac1n\,\rm{tr} (\boldsymbol{\tau})\,{\mathbb{I}}$ is the deviatoric tensor of a given $\boldsymbol{\tau} \in \mathbb L^2(\Omega)$, and (1.6e) constitutes a uniqueness condition for the pressure. At this point we remark that, because of the division by $\mu(\varphi)$ in (1.6a), integration by parts (after multiplication by a test function) is now possible. However, it can be seen that this leads to the usage of a continuous injection from $H^1(\Omega)$ into $L^8(\Omega)$, as required by the following estimate:

with $C(\Omega) > 0$ and valid for any ${\varphi} \in H^1(\Omega)$; ${\bf u},{\bf w} \in {\bf H}^1(\Omega)$; ${{\mathit{\boldsymbol{\tau}}}} \in {\mathbb{L}}^2(\Omega)$, and more important, for $\Omega \subset R^d$, $d\in \{1,2\}$, according to the Sobolev embedding theorem (cf., e.g., [37,Theorem 1.3.5]). This estimate is used in several ways throughout [5], at both continuous and discrete levels (see, [5,Lemmas 3.8,4.5 and 5.3]), and its main purpose is to help in the proof of Lipschitz continuity of the fixed-point operator ${\bf T}$ (respectively its discrete version ${\bf T}_h$) that consequently provides well-posedness of the continuous formulation (respectively the Galerkin scheme).

One of the purposes of this work is to derive a new augmented method for (1.1) that remains valid also for the 3D case. To this end, one can look at e.g. [15] where the strain rate tensor ${\bf e}({\bf u})$ is considered as a new unknown, in addition to either the stress (or pseudostress) and vorticity tensors. This provides more flexibility in the scheme, as it is no longer necessary (nor much convenient) to divide in (1.6a) by the viscosity to set the gradient free. Instead, the decomposition of the velocity gradient into its symmetric and skew-symmetric parts provides an equation to be integrated by parts. However, proceeding a bit differently than in those works, here we do not impose the symmetry of the stress by testing against the whole space of the skew symmetric tensors, but only against a proper subspace of it, thus yielding what we call an ultra-weak imposition of such constraint. As a consequence, we do not require the vorticity as a further unknown (also interpreted as an associated Lagrange multiplier), and hence the resulting mixed scheme has fewer degrees of freedom without affecting the accuracy of the approximate solutions. Moreover, the pressure and the vorticity are easily recovered by postprocessing formulae that provide the same rates of convergence of the main unknowns. In terms of the analysis, we can suitably modify the approach from [5], thus giving rise to an augmented mixed-primal finite element method using discontinuous piecewise polynomial functions of degree $\leq k$ to approximate the strain rate and normal heat flux, Raviart-Thomas elements of order $k$ for the stress, and Lagrange elements of order $k+1$ for the velocity and temperature. We stress that the already described extension to the 3D case, the new way of imposing the symmetry of the stress with the consequent decrease in the number of unknowns, and the availability of the aforementioned postprocessing formulae, constitute the main contributions of the present work. Let us also clarify that rather than proposing a cheap method for any reformulation or simplification of the governing equations, our underlying motivation is to rigorously derive properties of numerical schemes that preserve the mixed structure of the PDE system.

The rest of this work is organized as follows. In Section 2 we rewrite the Boussinesq problem (1.1) considering the strain rate and stress tensors as new variables, to then derive an augmented mixed-primal formulation whose well-posedness is proved by means of a fixed-point approach. Similarly, in Section 3 we provide the corresponding Galerkin scheme and its associated well-posedness, to then, in Section 4, proceed to derive a priori error estimates and state the respective rates of convergence, including those for the postprocessed approximations of pressure and vorticity, when a particular choice of finite element subspaces is made. Finally, to complement our theoretical results, we present in Section 5 a set of numerical examples that serve to confirm the properties of the proposed schemes.

2.   The continuous formulation

In order to avoid the division by $\mu({\varphi})$ in (1.6a), we have a look at recent work [15], where the authors develop a augmented mixed finite element method for the Navier-Stokes equations with nonlinear viscosity. This approach relies on the definition of the strain rate tensor as a new unknown

where

In addition, for each ${\bf v} \,\in\, \mathbf H^1(\Omega)$ we let ${{\mathit{\boldsymbol{\eta}}}}({\bf v})$ be the skew-symmetric part of the velocity gradient tensor $\nabla {\bf v}$, that is

which certainly lies in $\mathbb{L}^2_{\tt skew}( {\Omega} )$ (cf. (1.5)). In this way, bearing in mind (2.1), and noting from (1.4) and (2.2) that ${{\mathit{\boldsymbol{\gamma}}}} = {{\mathit{\boldsymbol{\eta}}}}({\bf u})$, the system (1.6) can be rewritten as: Find $({\bf t},{{\mathit{\boldsymbol{\sigma}}}},{\bf u},{\varphi})$ such that

Before continuing in the following section with the derivation of the mixed-primal formulation of (2.3), we find it important to remark at this point that it is also possible to account the case of a nonlinear thermal conductivity ${\mathbb{K}}$ in the forthcoming analysis. In particular, within the context of coupled flow-transport and sedimentation-consolidation models, in which the equation for the concentration ${\varphi}$ reads as the present one for the energy, we may refer to [8,eq. (2.1)] and [9,eq. (2.1)], where diffusion coefficients depending on $\nabla{\varphi}$ and ${\varphi}$, respectively, were considered for the corresponding analyses. More recently, a slight modification of the present model (1.1), in which ${\mathbb{K}}$ was allowed to depend on the temperature ${\varphi}$, namely in the form $\kappa({\varphi})$, with $\kappa : R \longrightarrow R^+$, was introduced and analyzed in [4] by using a fully-mixed variational formulation. In this case, the pseudoheat vector $\mathbf p : = \kappa({\varphi}) \nabla {\varphi} - {\varphi} \mathbf u$ and the temperature gradient are introduced as auxiliary unknowns, so that the energy equation is rewritten as $\,\boldsymbol\zeta = \nabla {\varphi}\,, \quad \mathbf p : = \kappa({\varphi}) \boldsymbol\zeta - {\varphi} \mathbf u\,, {\quad\hbox{and}\quad} \mathrm{div}(\mathbf p) = 0 {\quad\hbox{in}\quad} \Omega$. Further details on the analysis of such fully-mixed approach are available in [4].

2.1.   An augmented mixed-primal formulation

Multiplying (2.3a) by a test function ${{\mathit{\boldsymbol{\tau}}}} \in \mathbb{H}(\mathbf{div} ; {\Omega} )$, integrating by parts and using the Dirichlet condition (2.3e), we obtain

Then, we multiply (2.3b) and (2.3c) by appropriate test functions, imposing at the same time the symmetry of the stress tensor ${{\mathit{\boldsymbol{\sigma}}}}$ in the form

thus obtaining

and

We highlight that (2.4) constitutes what we call an ultra-weak imposition of the symmetry of ${{\mathit{\boldsymbol{\sigma}}}}$ since $\, \Big\{ {{\mathit{\boldsymbol{\eta}}}}({\bf v}): \ {\bf v} \in {\bf H}^1(\Omega)\Big\}$ is a proper subspace of $\mathbb{L}^2_{\tt skew}( {\Omega} )$. This idea has also been applied in [7].

The weak form of the energy conservation equation is recalled next from [5]:

where $\lambda \,: = \, \mathbb K\nabla \varphi\cdot{{\mathit{\boldsymbol{\nu}}}}\in H^{-1/2}(\Gamma)$ is a Lagrange multiplier taking care of the non-homogeneous Dirichlet boundary condition. Notice that, due to the second term in (2.5) and the right hand side of (2.6), the velocity ${\bf u}$ must live in ${\bf H}^1(\Omega)$, since appealing to the continuous injection of $H^1(\Omega)$ into $L^4(\Omega)$, there exist positive constants $c_1(\Omega)$ and $c_2(\Omega)$ such that

and

Also, obeying to the orthogonal decomposition

where

we can consider ${{\mathit{\boldsymbol{\sigma}}}}$ and ${{\mathit{\boldsymbol{\tau}}}}$ in $\mathbb{H}_0(\mathbf{div} ; {\Omega} )$ (see [5,Lemma 3.1] for a detailed justification of this change). Having in mind these considerations, at a first glance, the weak formulation reads: Find $({\bf t},{{\mathit{\boldsymbol{\sigma}}}},{\bf u},{\varphi},\lambda) \in \mathbb{L}^2_{\tt tr}( {\Omega} ) \times \mathbb{H}_0(\mathbf{div} ; {\Omega} ) \times {\bf H}^1(\Omega) \times H^1(\Omega) \times H^{-1/2}(\Gamma)$ such that

To achieve a conforming scheme, and to properly analyze (2.10), we augment this variational formulation using Galerkin terms arising from (2.3), but tested differently from (2.10), namely:

where $\kappa_j$, $j\in\{1,\dots,4\}$ are stabilization (or augmentation) constants to be specified later on. In this way, denoting by ${\mathcal{H}} : = \mathbb{L}^2_{\tt tr}( {\Omega} ) \times \mathbb{H}_0(\mathbf{div} ; {\Omega} ) \times {\bf H}^1(\Omega)$, ${\vec{\vphantom t\smash{{{\bf t}}}}} : = ({\bf t},{{\mathit{\boldsymbol{\sigma}}}},{\bf u})$ and ${\vec{\vphantom t\smash{{{\bf s}}}}} : = ({\bf s},{{\mathit{\boldsymbol{\tau}}}},{\bf v})$, we arrive at the following augmented mixed-primal formulation: Find $({\vec{\vphantom t\smash{{{\bf t}}}}},({\varphi},\lambda)) \in {\mathcal{H}} \times H^1(\Omega) \times H^{-1/2}(\Gamma)$ such that

where, given an arbitrary $({\bf w},\phi) \in {\bf H}^1(\Omega) \times H^1(\Omega)$, the forms ${\bf A}_\phi$, ${\bf B}_{\bf w}$, ${\bf a}$, ${\bf b}$, and the functionals $F_D$, $F_\phi$, $F_{{\bf w},\phi}$ and $G$ are defined as

for all ${\vec{\vphantom t\smash{{{\bf t}}}}},{\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$;

for all ${\varphi},\psi \in H^1(\Omega)$;

for all $(\psi,\xi) \in H^1(\Omega) \times H^{-1/2}(\Gamma)$;

for all ${\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$;

for all $\psi \in H^1(\Omega)$; and

for all $\xi \in H^{-1/2}(\Gamma)$.

2.2.   The fixed-point argument

A crucial tool in [5] to prove the well-posedness of the continuous and discrete formulations is a technique that decouples the problem into the mixed formulation of the momentum equation and the primal formulation of the energy equation, which further enables us to rewrite the formulation as a fixed-point problem. Hence, we denote ${\bf H} : = {\bf H}^1(\Omega) \times H^1(\Omega)$ and consider in what follows the operator ${\bf S}:{\bf H} \to {\mathcal{H}}$ defined by

where ${\vec{\vphantom t\smash{{{\bf t}}}}}$ is the solution of the problem: Find ${\vec{\vphantom t\smash{{{\bf t}}}}} \in {\mathcal{H}}$ such that

In addition, let $\widetilde{{\bf S}}:{\bf H} \to H^1(\Omega)$ be the operator defined by

where ${\varphi}$ is the first component of a solution to: Find $({\varphi},\lambda) \in H^1(\Omega) \times H^{-1/2}(\Gamma)$ such that

In this way, by introducing the operator ${\bf T}:{\bf H} \to {\bf H}$ as

we realize that (2.11) can be rewritten as the fixed-point problem: Find $({\bf u},{\varphi}) \in {\bf H}$ such that

As in [5], the objective is to use the Banach fixed-point theorem to prove existence and uniqueness of (2.23). We recall that the key difference in the present work with respect to [5] is in the problem that defines the operator ${\bf S}$, and therefore, those results associated to the operator $\widetilde{{\bf S}}$ and the primal formulation of the energy equation will be considered here as well, but we only cite them.

2.3.   Well-posedness of the uncoupled problems

In what follows, we consider the norms

We first recall some results that will be used for ellipticity purposes.

Lemma 2.1. There exists $c_3(\Omega)> 0$ such that

Proof. See [14,Proposition 3.1], [28,Lemma 2.3].

Lemma 2.2. There exists $\vartheta_0(\Omega) > 0$ such that

Proof. See [27,Lemma 3.1].

The following results establish sufficient conditions for the operators ${\bf S}$ and $\widetilde{{\bf S}}$ being well-defined, equivalently, (2.20) and (2.21) being well-posed.

Lemma 2.3. Assume that for $\delta_1 \in \left(0,\frac{2}{\mu_2}\right)$ and $\delta_2 \in (0,2)$ we choose

Then, there exists $r_0 > 0$ such that for each $r \in (0,r_0)$, the problem (2.20) has a unique solution ${\vec{\vphantom t\smash{{{\bf t}}}}} : = {\bf S}({\bf w},\phi) \in {\mathcal{H}}$ for each $({\bf w},\phi) \in {\bf H}$ such that ${\left\| \, {{\bf w}} \, \right\|}_{1,\Omega} \leq r$. Moreover, there exists a constant $C_{{\bf S}} > 0$, independent of $({\bf w},\phi)$ and $r$, such that there holds

Proof. Given $({\bf w},\phi) \in {\bf H}$, we notice from (2.12) that ${\bf A}_\phi$ is bilinear. Then, by using the upper bound of the viscosity function, the Cauchy-Schwarz inequality and the trace theorem with constant $c_0(\Omega)$, we see that for any ${\vec{\vphantom t\smash{{{\bf t}}}}},{\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$,

and therefore, there exists a constant $C_{\bf A} > 0$ depending only on $\mu_2$, $c_0(\Omega)$ and the stabilization parameters $\kappa_j$, such that

On the other hand, using (2.13) and (2.8), we obtain that for any ${\vec{\vphantom t\smash{{{\bf t}}}}}: = ({\bf t},{{\mathit{\boldsymbol{\sigma}}}},{\bf u}), \, {\vec{\vphantom t\smash{{{\bf s}}}}} : = ({\bf s},{{\mathit{\boldsymbol{\tau}}}},{\bf v}) \in {\mathcal{H}}$ there holds

which together with (2.25) implies the existence of a positive constant denoted by ${\left\| \, {{\bf A} + {\bf B}} \, \right\|}$, independent of $({\bf w},\phi)$ such that

To prove that ${\bf A}_\phi + {\bf B}_{\bf w}$ is elliptic, we first prove that ${\bf A}_\phi$ is elliptic. Indeed, for any ${\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$ we have

and then, using the bounds for the viscosity and the Cauchy-Schwarz and Young inequalities, we obtain for any $\delta_1,\delta_2 > 0$ and any ${\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$ that

Then, defining the following constants:

and using Lemmas 2.1 and 2.2, it is possible to find a constant $\alpha(\Omega) : = \min\{\alpha_1,\alpha_4, \alpha_5\}$, independent of $({\bf w},\phi)$, such that

Combining the foregoing inequality with (2.26), we get that, for any ${\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$, there holds

Therefore, we easily see that

provided that

that is,

thus proving ellipticity for ${\bf A}_\phi+{\bf B}_{\bf w}$ under the requirement (2.29). Finally, the linearity of the functionals $F_D$ and $F_\phi$ is clear, and from (2.16), (2.17), using the Cauchy-Schwarz inequality, the trace estimates in $\mathbb{H}(\mathbf{div} ; {\Omega} )$ and ${\bf H}^1(\Omega)$, with constants $1$ and $c_0(\Omega)$, and the continuous injection from $H^{1/2}(\Gamma)$ into $L^2(\Gamma)$ with constant $C_{1/2}$ we have

for all ${\vec{\vphantom t\smash{{{\bf s}}}}} \in {\mathcal{H}}$. Thus, there exists a constant $M_{\bf S} : = \max\Big\{(1+\kappa_2^2)^{1/2},1+ { \kappa_4} c_0(\Omega) C_{1/2} \Big\}$ such that

and by the Lax-Milgram theorem (see, e.g. [28,Theorem 1.1]), there exists a unique ${\vec{\vphantom t\smash{{{\bf t}}}}} \in {\mathcal{H}}$ solution of (2.20), and (2.24) is satisfied with $C_{\bf S} : = \frac{2M_{\bf S}}{\alpha(\Omega)}$, a constant clearly independent of $({\bf w},\phi)$. In addition, the fact that $\alpha(\Omega)$ and $M_{\bf S}$ are both independent of $r$ guarantees that $C_{\bf S}$ does not depend on $r$ either.

We stress here that, while the viscosity $\mu$ could be assumed to be a monotone decreasing function of the temperature $\phi$, this fact by itself would not help to simplify the solvability analysis of (2.20) (which defines the operator ${\bf S}$), since $\phi$ does not form part of the corresponding unknowns, but it actually constitutes a datum for that problem. Indeed, all what one needs from $\mu$ for the proof of Lemma 2.3 are its upper and lower bounds. A similar remark is valid for the Lipschitz-continuity of ${\bf S}$ (cf. Lemma 2.6 below) in which the homonimous property of the viscosity plays a key role. As it will become clear later on, the smallness assumptions on the data for the well-posedness of our problem arises from the main properties of the resulting fixed point operator, and particularly from those inherited from ${\bf S}$, so that the eventual monotonicity of $\mu$ does not help to remove those conditions either.

Lemma 2.4. For each $({\bf w},\phi) \in {\bf H}$ there exists a unique pair $({\varphi},\lambda) \in H^1(\Omega) \times H^{-1/2}(\Gamma)$ solution of the problem (2.21). Moreover, there exists $C_{\widetilde{{\bf S}}} > 0$, independent of $({\varphi},\lambda)$ and $r$, such that

Proof. The results comes from a direct application of the Babuška-Brezzi theory (see [5,Lemma 3.6] or more precisely [20,Lemma 3.4]). In particular, the right-hand side of (2.31) is obtained after bounding the functionals $F_{{\bf w},\phi}$ (cf. (2.18)) and $G$ (cf. (2.19)), respectively, whose resulting bounds are independent of $r$ (cf. [20,eq. (3.39)] for $F_{{\bf w},\phi}$). In this way, since the constants yielding the inf-sup condition of $\mathbf b$ and the ellipticity of $\mathbf a$ in the kernel of $\mathbf b$, are both independent of $r$, the corresponding continuous dependence result for (2.21) confirms that $C_{\widetilde{{\bf S}}}$ does not depend on $r$ either.

From the previous two lemmas, it is clear that ${\bf T}$ is well-defined for any pair $({\bf w},\phi) \in {\bf W}$, where

which is nothing but the closed ball in ${\bf H}$ with center $({\bf 0},0)$ and radius $r\in (0,r_0)$, with $r_0$ defined in (2.29). We also notice that a particular choice of stabilization parameters is necessary. We begin by selecting the middle points of the ranges for $\delta_1$, $\delta_2$, $\kappa_1$ and $\kappa_3$, that is,

to then pick $\kappa_2$ and $\kappa_4$ so that $\alpha_2$ and $\alpha_3$ can attain the largest value possible, that is

We stress here that the foregoing feasible choices of $\kappa_j$, $j\in\{1,2,3,4\}$, are explicitly computable since they do not depend on the usually unknown constants $c_3(\Omega)$ and $\vartheta_0(\Omega)$ from Lemmas 2.1 and 2.2, respectively, but only on the known constants $\mu_1$ and $\mu_2$ bounding the viscosity (cf. (1.2)).

2.4.   Further-regularity assumption

Although the problem that defines ${\bf S}$ (i.e. (2.20)) is well-posed, a small further-regularity assumption is needed in order to continue with the analysis. Inspired by [8], we assume that ${\bf u}_D \in {\bf H}^{1/2+\varepsilon}(\Gamma)$ for some $\varepsilon \in (0,1)$ when $n = 2$, or $\varepsilon \in \left[ \frac{1}{2}, 1 \right)$ when $n = 3$, and that for each $({\bf z},\psi) \in {\bf H}$ with ${\left\| \, {{\bf z}} \, \right\|}_{1,\Omega}\leq r$, $r>0$ given, there hold $({\bf q},{{\mathit{\boldsymbol{\zeta}}}},{\bf v}) : = {\bf S}({\bf z},\psi) \in \mathbb{L}^2_{\tt tr}( {\Omega} ) \cap {\mathbb{H}}^\varepsilon(\Omega) \times \mathbb{H}_0(\mathbf{div} ; {\Omega} ) \cap {\mathbb{H}}^\varepsilon(\Omega) \times {\bf H}^{1+\varepsilon}(\Omega)$ and

with $\widetilde{C}_{{\bf S}}(r) > 0$ independent of ${\bf z}$ but depending on the upper bound $r$ of its ${\bf H}^1$-norm.

2.5.   Solvability analysis of the fixed-point equation

Throughout the rest of the paper we assume that the regularity hypotheses of the previous section are valid. We now proceed to directly fulfill the hypotheses of the Banach fixed-point theorem. The following result shows that ${\bf T}$ can map the ball ${\bf W}$ into itself.

Lemma 2.5. Consider the closed ball ${\bf W}$ defined in (2.32) with $r\in (0,r_0)$ and $r_0$ as given in (2.29). Suppose the data satisfy

where

and $C_{{\bf S}}$, $C_{\widetilde{{\bf S}}}$ are given in Lemmas $\rm2.3$ and $\rm2.4$, respectively. Then, there holds ${\bf T}({\bf W}) \subseteq {\bf W}$.

Proof. It follows as a consequence of the continuous dependence results (2.24) and (2.31), much in an identical way as in [5,Lemma 3.7].

Next, we prove some results that will help us to arrive at the Lipschitz continuity of ${\bf T}$.

Lemma 2.6. Let $r \in (0,r_0)$ with $r_0$ as given in (2.29). Then, there exists a positive constant $\widehat{C}_{{\bf S}}$, independent of $r$, such that

for all $({\bf w}_1,\phi_1),({\bf w}_2,\phi_2) \in {\bf H}$ such that ${\left\| \, {{\bf w}_1} \, \right\|}_{1,\Omega},{\left\| \, {{\bf w}_2} \, \right\|}_{1,\Omega} \leq r$.

Proof. Let $({\bf w}_1,\phi_1),({\bf w}_2,\phi_2) \in {\bf H}$ as indicated and let ${\vec{\vphantom t\smash{{{\bf t}}}}}_j : = ({\bf t}_j,{{\mathit{\boldsymbol{\sigma}}}}_j,{\bf u}_j) = {\bf S}({\bf w}_j,\phi_j) \in \mathcal H$, $j \in \{1,2\}$ be the corresponding solutions of (2.20). Then, adding and subtracting the equality

we find that

Thus, using the ellipticity of ${\bf A}_{\phi_2} + {\bf B}_{{\bf w}_2}$ (cf. (2.28)) and the foregoing expression, we obtain

First, we bound the last two terms in the same way as Lemma 2.3 (that is, using (2.26) and (2.30)):

and

Next, for the first term, we use the Lipschitz continuity of $\mu$ and the Cauchy-Schwarz and Hölder inequalities to show in a similar way to [5,Eq. (3.56)] that

with $p,q \in [1,+\infty)$ such that $\frac{1}{p}+\frac{1}{q} = 1$. We then take into consideration the further-regularity assumed in Section 2.4, and recall from the Sobolev embedding theorem that $H^\varepsilon(\Omega)$ is continuously embedded into $L^{2p}(\Omega)$, with

(cf. [1,Theorem 4.12], [37,Theorem 1.3.4]) meaning that there exists $C_\varepsilon > 0$ such that

In this way,

and (2.38) now yields

Therefore, putting (2.36), (2.37) and (2.39) together into (2.35), we obtain

with

and since ${\bf t}_1 = {\bf S}_1({\bf w}_1,\phi_1)$ and ${\bf u}_1 = {\bf S}_3({\bf w}_1,\phi_1)$, the last inequality is exactly the estimate (2.34). We end the proof by noting that the foregoing expression defining $\widehat{C}_{{\bf S}}$ confirms that this constant is independent of $r$.

We emphasize once again that, differently from [5], the present approach allows to handle both the 2D and 3D cases for the Boussinesq model with temperature-dependent viscosity. Indeed, notice how the foregoing Lemma shows the main difference with respect to [5]: the tensor-product term in (2.36) is no longer multiplied by an $H^1$-term, thus avoiding the use of the injection $H^1(\Omega) \hookrightarrow L^8(\Omega)$ (not ensured for $\Omega \subseteq R^3$) when splitting them as in [5,Eq. (3.54)], yielding a more robust formulation for the two and three-dimensional cases. On the other hand, the analogous result for $\widetilde{{\bf S}}$ remains intact.

Lemma 2.7. There exists a positive constant $\widehat{C}_{\widetilde{{\bf S}}}$, independent of $r$, such that

for all $({\bf w}_1,\phi_1), \, ({\bf w}_2,\phi_2) \in \mathbf H$.

Proof. We refer to [20,Lemma 3.7] for full details. We just remark here that the resulting constant $\widehat{C}_{\widetilde{{\bf S}}}$ depends on the ellipticity constant of $\mathbf a$ in the kernel of $\mathbf b$, and on the constant $c(\Omega)$ arising from the boundedness estimates of the functionals $F_{{\bf w}_1,\phi_1-\phi_2}$ and $F_{{\bf w}_1 - {\bf w}_2,\phi_2}$, namely (cf. [20,eq. (3.39)]):

from which it is clear that $\widehat{C}_{\widetilde{{\bf S}}}$ does not depend on $r$.

As a consequence of the previous Lemmas, ${\bf T}$ is a Lipschitz continuous operator, as shown next.

Lemma 2.8. Let $r \in (0,r_0)$, with $r_0$ given as in (2.29). Then, there exists a constant $C_{\bf T} > 0$ depending on $r$, such that

for all $({\bf w}_1,\phi_1), \, ({\bf w}_2,\phi_2) \in \mathbf W$.

Proof. The result comes from the definition of ${\bf T}$ (cf. (2.22)) and the estimates obtained in the previous two lemmas (cf. (2.34) and (2.40)) in an identical way to [20,Lemma 3.8] (see also [5,Lemma 3.10]). We omit further details.

In summary, from Lemmas 2.3 and 2.4, ${\bf T}:{\bf W} \to {\bf W}$ is well-defined and does map the ball into itself (thanks to Lemma 2.5). Then, Lemma 2.8 shows that the operator ${\bf T}$ is Lipschitz continuous with a data-depending constant given by $C_{\bf T}\,\Big\{ {\left\| \, {{\bf g}} \, \right\|}_{\infty,\Omega} + {\left\| \, {{\bf u}_D} \, \right\|}_{1/2+\varepsilon,\Gamma} \Big\}$. When this Lipschitz constant is less than $1$, the existence and uniqueness of a fixed point follows thanks to the Banach fixed-point theorem. In this regard, we stress that the dependence of $C_{\bf T}$ on a given $r$ does not affect the imposition of the aforementioned restriction since it is clearly achieved by sufficiently small data. By what has been explained in Section 2.2, this fact is equivalent to the well-posedness of the augmented mixed-primal formulation (2.11), and the result is stated as follows.

Theorem 2.9. Assume that for $\delta_1 \in \left(0,\frac{2}{\mu_2}\right)$ and $\delta_2 \in (0,2)$ we choose

and consider the ball ${\bf W}$ (cf. (2.32)) with radius $r \in (0,r_0)$ and $r_0$ as in (2.29). In addition, assume that the data satisfy

with $c(r)$ as in Lemma $\rm2.5$, and

Then, the problem (2.11) has a unique solution $({\vec{\vphantom t\smash{{{\bf t}}}}},({\varphi},\lambda)) \in {\mathcal{H}} \times H^1(\Omega) \times H^{-1/2}(\Gamma)$, with $({\bf u},{\varphi}) \in {\bf W}$. Moreover, there hold

and

with $C_{\bf S}$ and $C_{\widetilde{{\bf S}}}$ as in Lemmas $\rm2.3$ and $\rm2.4$, respectively.

3.   The Galerkin scheme

In this section we derive a Galerkin method for the augmented mixed-primal formulation (2.11). Let us consider ${\mathcal{T}}_h$, a regular triangulation of $\bar{\Omega}$ by triangles $T$ (when $n = 2$) or tetrahedra $T$ (when $n = 3$) of diameter $h_T$ and define the meshsize $h: = \max\{h_T : T \in {\mathcal{T}}_h\}$. Then, consider arbitrary finite-dimensional subspaces ${\mathbb{H}_h^{{\bf t}}} \subset \mathbb{L}^2_{\tt tr}( {\Omega} )$, ${\mathbb{H}_h^{{\mathit{\boldsymbol{\sigma}}}}} \subset \mathbb{H}_0(\mathbf{div} ; {\Omega} )$, ${{\bf H}_h^{{\bf u}}} \subset {\bf H}^1(\Omega)$, ${H_h^{\varphi}} \subset H^1(\Omega)$, ${H_h^\lambda} \subset H^{-1/2}(\Gamma)$ and denote by ${\mathcal{H}}_h : = {\mathbb{H}_h^{{\bf t}}} \times {\mathbb{H}_h^{{\mathit{\boldsymbol{\sigma}}}}} \times {{\bf H}_h^{{\bf u}}}$, ${\vec{\vphantom t\smash{{{\bf t}}}}}_h : = ({\bf t}_h,{{\mathit{\boldsymbol{\sigma}}}}_h,{\bf u}_h)$ and ${\vec{\vphantom t\smash{{{\bf s}}}}}_h : = ({\bf s}_h,{{\mathit{\boldsymbol{\tau}}}}_h,{\bf v}_h)$. Hence, according to (2.11), the Galerkin scheme reads: Find $({\vec{\vphantom t\smash{{{\bf t}}}}}_h,({\varphi}_h,\lambda_h)) \in {\mathcal{H}}_h \times {H_h^{\varphi}} \times {H_h^\lambda}$ such that

where the forms ${\bf A}_{{\varphi}_h}$, ${\bf B}_{{\bf u}_h}$, ${\bf a}$ and ${\bf b}$; and the functionals $F_{{\varphi}_h}$, $F_D$, $F_{{\bf u}_h,{\varphi}_h}$ are defined by (2.12)-(2.19). Since the proof of well-posedness follows the steps of the previous section, and moreover, analogously to [5,Section 4], we only state the requirements to be imposed over the finite-dimensional subspaces and the main result, which is analogous to [5,Theorem 4.8].

3.1.   Well-posedness of the Galerkin scheme

It can be seen that no restrictions have to be added to ${\mathbb{H}_h^{{\bf t}}}$, ${\mathbb{H}_h^{{\mathit{\boldsymbol{\sigma}}}}}$ and ${{\bf H}_h^{{\bf u}}}$ other than being finite-dimensional subspaces of the described spaces. However, for ellipticity purposes of ${\bf a}$ in the discrete kernel of the operator induced by ${\bf b}$ (according the the Babuška-Brezzi theory), the following two inf-sup conditions must be met:

$\bf (H.1)$ There exists a constant $\widehat{\alpha} > 0$, independent of $h$ such that

where

$\bf (H.2)$ There exists a constant $\widehat{\beta} > 0$, independent of $h$ such that

Then, denoting by ${\bf W}_h$ the closed ball in ${\bf H}_h : = {{\bf H}_h^{{\bf u}}} \times {H_h^{\varphi}}$ of radius $r$ and center $({\bf 0},0)$, that is

the main result of this section reads as follows.

Theorem 3.1. Assume that for $\delta_1 \in \left(0,\frac{2}{\mu_2}\right)$ and $\delta_2 \in (0,2)$ we choose

and consider the ball ${\bf W}_h$ with radius $r \in (0,r_0)$, $r_0$ as in (2.29). Then, there exist positive constants $C_{\bf S}$, $\widetilde{C}_{\widetilde{{\bf S}}}$ and $\widetilde{c}(r)$ such that, if the data satisfy

then the problem (3.1) has at least one solution $({\vec{\vphantom t\smash{{{\bf t}}}}}_h,({\varphi}_h,\lambda_h)) \in {\mathcal{H}}_h \times {H_h^{\varphi}} \times {H_h^\lambda}$ with $({\bf u}_h,{\varphi}_h) \in {\bf W}_h$. Moreover, there hold

and

Proof. We only mention that (3.1) is transformed into a fixed-point problem that is analyzed by means of the Brouwer fixed-point theorem in the convex and compact set ${\bf W}_h$ (see [5,Theorem 4.8]).

3.2.   Specific choice of finite element subspaces

Given a set $S \subset {\bf R}: = R^n$ and an integer $k\geq 0$, we define $P_k(S)$ as the space of polynomial functions on $S$ of degree $\leq k$, and for each $T \in {\mathcal{T}}_h$, we define the local Raviart-Thomas spaces of order $k$ as

where ${\bf x}$ is a generic vector in ${\bf R}$. Hence, the strain rate, stress, velocity and temperature variables can be approximated using the following finite element subspaces:

whereas for the normal component of the heat flux, we let $\{\widetilde{\Gamma}_1, \widetilde{\Gamma}_2, \dots, \widetilde{\Gamma}_m\}$ be an independent triangulation of $\Gamma$ (made of straight segments in $R^2$ or triangles in $R^3$), and define $\widetilde{h} : = \max_{j \in \{1,\dots,m\}} |\widetilde{\Gamma}_j|$. Then, with the same integer $k \geq 0$ used in definitions (3.5)-(3.8), we approximate $\lambda$ by piecewise polynomials of degree $\leq k$ over this new mesh, that is

It can be seen that ${H_h^{\varphi}}$ and $H_{\widetilde{h}}^\lambda$ satisfy the conditions (3.2) and (3.3) as long as $h \leq C_0 \widetilde{h}$ for some $C_0 > 0$ (cf. [5,Section 4.3]). Other problems and corresponding references where some of or all the above finite element subspaces have been employed include, among others, coupled flow-transport [8,9], natural convection with phase-change [7], Navier-Stokes [15,16], and Stokes-Darcy (see, e.g. [29], where the above restriction between the mesh sizes $h$ and $\widetilde h$ was discussed).

In turn, the approximation properties of the subspaces defined in (3.5) - (3.8) and (3.9) are (cf. [14,28])

$(\mathbf{AP}_h^{{\bf t}})$ There exists $C>0$, independent of $h$, such that for each $s \in (0,k+1]$, and for each ${\bf t} \in {\mathbb{H}}^s(\Omega) \cap \mathbb{L}^2_{\tt tr}( {\Omega} )$, there holds

$(\mathbf{AP}_h^{{\mathit{\boldsymbol{\sigma}}}})$ There exists $C>0$, independent of $h$, such that for each $s \in (0,k+1]$, and for each ${{\mathit{\boldsymbol{\sigma}}}} \in \mathbb{H}^s(\Omega) \cap \mathbb{H}_0(\mathbf{div} ; {\Omega} )$ with ${\mathit{\boldsymbol{\mathrm{div}}} \, {{{\mathit{\boldsymbol{\sigma}}}}}} \in {\bf H}^s(\Omega)$, there holds

$(\mathbf{AP}_h^{{\bf u}})$ there exists $C>0$, independent of $h$, such that for each $s \in (0,k+1]$, and for each ${\bf u} \in {\bf H}^{s+1}(\Omega)$, there holds

$(\mathbf{AP}_h^{\varphi})$ there exists $C>0$, independent of $h$, such that for each $s \in (0,k+1]$, and for each ${\varphi} \in H^{s+1}(\Omega)$, there holds

$(\mathbf{AP}_{\widetilde{h}}^\lambda)$ there exists $C>0$, independent of $\widetilde{h}$, such that for each $s \in (0,k+1]$, and for each $\lambda \in H^{-1/2+s}(\Gamma)$, there holds

4.   A priori error analysis

Let $({\vec{\vphantom t\smash{{{\bf t}}}}},({\varphi},\lambda)) \in {\mathcal{H}} \times H^1(\Omega) \times H^{-1/2}(\Gamma)$ with $({\bf u},{\varphi}) \in {\bf W}$ be the solution of (2.11), and $({\vec{\vphantom t\smash{{{\bf t}}}}}_h,({\varphi}_h,\lambda_h)) \in {\mathcal{H}}_h \times {H_h^{\varphi}} \times {H_h^\lambda}$ with $({\bf u}_h,{\varphi}_h) \in {\bf W}_h$ be a solution of the discrete problem (3.1), that is,

and

In what follows, we denote as usual

and

First, the error estimate related to the variables of the momentum equation is obtained by means of the Strang Lemma, applied to the pair (4.1). We recall the Lemma, and its consequent result next.

Lemma 4.1 (Strang). Let $V$ be a Hilbert space, $F \in V'$, and $A:V \times V \to R$ be a bounded and $V$-elliptic bilinear form. In addition, let $\{V_h\}_{h>0}$ be a sequence of finite-dimensional subspaces of $V$, and for each $h>0$, consider a bounded bilinear form $A_h:V_h \times V_h \to R$ and a functional $F_h \in V_h'$. Assume that the family $\{A_h\}_{h>0}$ is uniformly elliptic in $V_h$, that is, there exists a constant $\widetilde{\alpha} > 0$, independent of $h$, such that

In turn, let $u \in V$ and $u_h \in V_h$ such that

Then, for each $h > 0$, there holds

where $C_{ST} : = \widetilde{\alpha}^{-1} \max\{1,{\left\| \, {A} \, \right\|}\}$.

Proof. See [38,Theorem 11.1].

Lemma 4.2. Let $C_{ST} : = \frac{2}{\alpha(\Omega)} \max\{ 1,{\left\| \, {{\bf A}_{\varphi}+{\bf B}_{\bf u}} \, \right\|} \}$, where $\frac{\alpha(\Omega)}{2}$ is the ellipticity constant of ${\bf A}_{\varphi}+{\bf B}_{\bf u}$ (cf. (2.28)). Then, there holds

Proof. From Lemma 2.3, we see that ${\bf A}_{\varphi}+{\bf B}_{\bf u}$ and ${\bf A}_{{\varphi}_h}+{\bf B}_{{\bf u}_h}$ are bilinear, bounded (both with constant ${\left\| \, {{\bf A}_{\varphi}+{\bf B}_{\bf u}} \, \right\|}$, w.l.o.g. since it is independent of $({\bf u},{\varphi})$) and uniformly elliptic (both with constant $\frac{\alpha(\Omega)}{2}$). Also, $F_{\varphi}+F_D$ and $F_{{\varphi}_h}+F_D$ are linear bounded functionals in ${\mathcal{H}}$ and ${\mathcal{H}}_h$, respectively. Hence, a straightforward application of Lemma 4.1 to the pair (4.1) yields

where $C_{ST} : = \frac{2}{\alpha(\Omega)} \max\big\{1,{\left\| \, {{\bf A}_{\varphi}+{\bf B}_{\bf u}} \, \right\|}\big\}$. First, we notice that

Then, in order to estimate the last supremum in (4.4), we add and subtract suitable terms to write

and so, using the boundedness of the bilinear forms ${\bf A}_{\varphi}+{\bf B}_{\bf u}$ and ${\bf A}_{{\varphi}_h}+{\bf B}_{{\bf u}_h}$ (cf. (2.27)), the estimate (2.8), the continuous embedding $H^1(\Omega) \hookrightarrow L^{n/\varepsilon}(\Omega)$ with constant $\widetilde{C}_\varepsilon$ and the further-regularity assumption in a similar way to (2.39), we obtain

This inequality, together with (4.5), back into (4.4), results in (4.3), therefore concluding the proof.

Next, we recall from [20,Lemma 5.4] (see also [5,Lemma 5.4]) the error estimate of the variables in the energy equation.

Lemma 4.3. There exists a positive constant $\widehat{C}_{ST}$, depending only on ${\left\| \, {{\bf a}} \, \right\|}$, ${\left\| \, {{\bf b}} \, \right\|}$, $\widehat{\alpha}$ and $\widehat{\beta}$ (cf. (3.2), (3.3)), such that

Proof. It suffices to apply the generalized Strang-type estimate for saddle point problems provided in [20,Lemma 5.2] (see also [5,Lemma 5.2]), which constitutes a particular case of [38,Theorem 11.12], to the context given by (4.2).

Hence, adding the estimates obtained in the previous two lemmas, we have a preliminary estimate for the total error:

where

Then, bounding the terms ${\left\| \, {{\bf u}} \, \right\|}_{1,\Omega}$, ${\left| \, {{\varphi}} \, \right|}_{1,\Omega}$ and ${\left\| \, {{\bf u}_h} \, \right\|}_{1,\Omega}$ using the continuous dependence results (2.41), (2.42) and (3.4), and the further-regularity assumption (2.33) to bound ${\left\| \, {{\bf t}} \, \right\|}_{\varepsilon,\Omega}$, (4.6) becomes

Therefore, using the continuous injection $H^{1/2+\varepsilon}(\Gamma) \hookrightarrow H^{1/2}(\Gamma)$ with constant $\widehat{C}_i$ and denoting by

we see from (4.7) that

which leads us the main result of this section.

Theorem 4.4. Assume that

Then, there exists $C > 0$ depending only on parameters, data and other constants, all of them independent of $h$, such that

Proof. The assumption (4.9) allows us to subtract the total error term in the right-hand side of (4.8), thus verifying the Céa's estimate with $C = 2 \max \Big\{ C_{ST}(1+2{\left\| \, {{\bf A}_{\varphi}+{\bf B}_{\bf u}} \, \right\|}),\widehat{C}_{ST} \Big\}$.

Finally, we state the rates of convergence of the Galerkin scheme (3.1) when the finite element subspaces (3.5)-(3.9) are used.

Theorem 4.5. In addition to the hypotheses of Theorems $\rm2.9$, $\rm3.1$ and $\rm4.4$, assume that there exists $s>0$ such that ${\bf t} \in \mathbb{H}^s(\Omega)$, ${{\mathit{\boldsymbol{\sigma}}}} \in \mathbb{H}^s(\Omega)$, ${\mathit{\boldsymbol{\mathrm{div}}} \, {{{\mathit{\boldsymbol{\sigma}}}}}} \in {\bf H}^s(\Omega)$, ${\bf u} \in {\bf H}^{s+1}(\Omega)$, ${\varphi} \in H^{s+1}(\Omega)$ and $\lambda \in H^{-1/2+s}(\Gamma)$. Then, there exists $\widehat{C} > 0$, independent of $h$ and $\widetilde{h}$ such that for all $h \leq C_0 \widetilde{h}$ there holds

Proof. It follows from Céa's estimate (4.10) and the approximation properties $(\mathbf{AP}_h^{{\bf t}})$, $(\mathbf{AP}_h^{{{\mathit{\boldsymbol{\sigma}}}}})$, $(\mathbf{AP}_h^{{\bf u}})$, $(\mathbf{AP}_h^{{\varphi}})$, and $(\mathbf{AP}_{\widetilde{h}}^{\lambda})$ described in Section 3.2.

5.   Numerical results

We now present three numerical examples showing the performance of the augmented mixed-primal method (3.1) with the subspaces specified in Section 3.2. The computational implementation uses the finite element library FEniCS (cf. [6]), in particular the module multiphenics¹ that allows a straightforward description of the Lagrange multiplier and a block structure manipulation of the linearized systems arising at each Picard step. The solution of these linear systems is done employing the unsymmetric multi-frontal direct solver MUMPS (cf. [10]).

http://mathlab.sissa.it/multiphenics

The implemented iterative method mimics the fixed-point strategy from Section 2.2: it begins with an initial state at rest $({\bf u},{\varphi}) = ({\bf 0},0)$ and it stops whenever the relative error between two consecutive iterations of the solution vector measured in the $\ell^2$-norm (the sum of squared entries of a vector) is sufficiently small, i.e.,

where $\mathtt{tol} = 10^{-6}$ is a specified tolerance. We also recall that the pressure is post-processed as

which, as explained in [5,Section 5.2], converges to the exact pressure $p$ at the same rate as the other variables (cf. Theorem 4.5). Similarly, as suggested by the second expression in (1.4), the vorticity can be postprocessed as

which easily yields

thus proving that it also converges at the same rate provided by Theorem 4.5. In this way, we define the error per variable

as well as their corresponding rates of convergence

where $h$ and $h'$ (respectively $\widetilde{h}$ and $\widetilde{h}'$) denote two consecutive mesh sizes with errors $e$ and $e'$.

5.1.   Example 1

We first consider $\Omega : = (-1,1)^2$, viscosity, thermal conductivity and body force given by $\mu({\varphi}) = \exp(-0.25 \, {\varphi})$, ${\mathbb{K}} = {\mathbb{I}}$, ${\bf g} = (0,1)^{{\tt t}}$, and boundary conditions such that the exact solution to (1.1) is

Non-homogeneous terms will then appear in the right-hand sides of the momentum and energy equations. Nevertheless, well-posedness is still ensured, since the smoothness of the exact solution makes these terms immediately belong to $L^2(\Omega)$, thus requiring only a minor modification in the functionals of the variational formulation. Concerning the stabilization parameters, these are taken as pointed out in Section 2.3, where the viscosity bounds are estimated as $\mu_1 = 0.5$, $\mu_2 = 1.25$, thus resulting in $\kappa_1 = \kappa_2 = 0.32$, $\kappa_3 = 0.25$ and $\kappa_4 = 0.125$. We also remark that the trace condition on the stress is enforced through penalization, here and also in the upcoming examples.

In Figure 5.1 we show part of the obtained numerical solution with 214,788 DOF and a first order approximation, whereas in Table 5.1 we show the convergence history given the specified data and the finite element spaces from Section 3.2 with successive quasi-uniform mesh refinements. In both cases, it can be seen that the rates of convergence are the expected ones according to Theorem 4.5, that is, $\mathcal{O}(h)$ for the first case, and $\mathcal{O}(h^2)$ for the second one. We also observe that around eight iterations are needed to reach convergence of the Picard algorithm.

5.2.   Example 2

We next consider $\Omega : = (0,1)^2$; viscosity, thermal conductivity and body force as in Example 1, and boundary conditions and source terms such that the exact solution is

where $r_1$ and $r_2$ are positive parameters, and

It is expected to find a counter-clockwise rotating vortex with center $(\hat{x},\hat{y})$, where

In particular, taking $r_1 = r_2 = 4.5$, the center of the vortex is expected to appear at the top-right corner of the cavity (that is, $(\hat{x},\hat{y}) = (0.829,0.829)$). Then, considering the stabilization parameters as in Section 2.3, estimating the viscosity bounds in $\mu_1 = 0.74$, $\mu_2 = 1.35$, we obtain the values $\kappa_1 = \kappa_2 = 0.406$, $\kappa_3 = 0.37$ and $\kappa_4 = 0.185$.

In Table 5.2 we show the corresponding convergence history. As expected, when using the finite-element subspaces of Section 3.2 with $k = 0$ and $k = 1$, the rates of convergence are near the optimal linear and quadratic orders, respectively. Part of the solution is shown in Figure 5.2, where a second-order approximation has been used with 724,448 DOF. The second-order method is more convenient than the first-order scheme, in terms of CPU time used to reach the same precision. Notice that a relatively high refinement was required to capture the small features of the solution (e.g. pressure and pseudostress).

5.3.   Example 3

The implementation of the numerical scheme and the accuracy for the three-dimensional case are assessed with this next computational test. The domain is the unit cube $\Omega = (0,1)^3$ and we consider the following closed-form solutions to the governing equations (1.1)

with $\mathbb{K} = \mathbb{I}$, $\mu({\varphi}) = \exp(-0.25 \, {\varphi})$, and ${\bf g} = (0,0,1)^{{\tt t}}$. The manufactured velocity is divergence free, it satisfies the compatibility condition (1.3) and it is used as Dirichlet datum on $\Gamma$. The exact temperature is uniformly bounded and it is also exploited as Dirichlet datum. In this configuration the viscosity bounds can be set as $\mu_1 = 0.5$, $\mu_2 = 1$ and the augmentation constants take the values $\kappa_1 = \kappa_2 = 0.5$, $\kappa_3 = 0.25$, and $\kappa_4 = 0.125$. The error history, associated with the schemes of order one and two, are performed using six steps of uniform mesh refinement applied to an initial structured tetrahedral mesh. On each level we compute approximate solutions, as well as errors and convergence rates defined as above. The boundary partition is considered conforming with the interior mesh, for sake of convenience and simplicity of the 3D mesh generation. Our findings are collected in Table 5.3, where errors and Picard iteration count are tabulated by number of degrees of freedom and meshsize. As in the 2D case, optimal error decay is observed for all individual errors (including the post-processed variables), and we also note that the errors $e({{\mathit{\boldsymbol{\sigma}}}}),e({\bf u})$ are dominant. One can also see that (perhaps assisted by the conformity between the interior and boundary meshes) the error associated with the boundary heat flux has a convergence systematically better than the optimal predicted by Theorem 4.5. Finally, we portray in Figure 5.3 a sample of the approximate solutions generated by the lowest-order mixed method on a relatively coarse mesh.

5.4.   Example 4

To close this section we conduct a benchmark test of a differentially heated, two-dimensional cavity of width 2 and height 1 (dimensionless units). The nonlinear viscosity used here is also explicitly dependent on depth, as in the models used in the context of mantle convection [12,40],

and the boundary treatment proceeds as follows. For the thermal energy conservation, the boundary is split into $\Gamma_D$ (top and bottom edges of the box) and $\Gamma_N$ (vertical walls) where temperature and heat flux are prescribed, respectively. The boundary temperature on $\Gamma_D$ is set to 0 on the top and 1 on the bottom edges, and it suffices with adequately defining $\varphi_D$ in the formulation. On $\Gamma_N$ we simply set zero thermal flux (the vertical walls are considered insulated). For the Navier-Stokes equations, and following [40], we impose free slip conditions defined as ${\bf u}\cdot{{\mathit{\boldsymbol{\nu}}}} = {{\mathit{\boldsymbol{\sigma}}}}{{\mathit{\boldsymbol{\nu}}}}\cdot{\bf m} = 0$, where ${\bf m}$ denotes the tangential vector. These conditions imply that the Galerkin term for boundary velocity of Section 2 is now modified as

and we also require the additional condition

where $\kappa_5>0$ can be taken, e.g., equal to $\kappa_4$. The coupled Boussinesq equations can be scaled so that the thermal conductivity is inverse to the adimensional Rayleigh number Ra = $10^4$, and the buoyancy term is simply ${\varphi}{\bf g}$ with ${\bf g} = (0,1)^{{\tt t}}$. In order to produce convection regimes we require extending the formulation to the transient case, adding to the thermal energy equation the term $\partial_t \varphi$ and to the momentum equation the acceleration $\partial_t {\bf u}$ (which for mantle convection could eventually be dropped as it scales with the inverse of a very large Prandtl number). These time derivatives are approximated with the backward Euler scheme, using a constant timestep of $\Delta t = 0.1$ and running the coupled Boussinesq equations until the final time $t = 12$. The initial velocity is zero and the initial temperature is ${\varphi}(0) = 1-y^2+0.01\cos(4\pi x)\sin(\pi y)$. We employ a relatively coarse structured mesh with 32K triangles (but graded to be refined in the $y-$direction near the bottom and top boundaries), which results in a mixed-primal method (for the lowest-order case) having 208,986 DOF. For this test we observe that an average of five fixed-point iterations are required to reach the fixed tolerance of $10^{-7}$. The Rayleigh number used here was sufficiently large to generate thermal blob-shaped forms (see the plumes in the temperature plot of Figure 5.4). The obtained convection modes are consistent with the flow patterns observed in [40] (see also [31,33]). They indicate that the velocity is driven by the temperature difference.