Versatile Mixed Methods for Non-Isothermal Incompressible Flows

The motion of a non-isothermal incompressible fluid is frequently induced by buoyancy forces, viscous forces, and pressure fields. In accordance with standard practices, we refer to the motion that is induced solely by buoyancy forces as natural or free convection, the motion that is induced solely by viscous forces and pressure fields as forced convection, and the motion that is induced by all three factors as mixed convection. In order to characterize the various types of convection, one may solve the incompressible Navier-Stokes equations for mass and momentum conservation, in conjunction with a temperature equation (usually obtained from the internal energy or enthalpy equations). In addition, one may couple the momentum and temperature equations via the approach of Oberbeck [1] and Boussinesq [2] by adding a temperature-dependent buoyancy term to the RHS of the momentum equation. The buoyancy term is assumed to be directly proportional to changes in the temperature field, and these changes are assumed to be small enough such that the density remains constant. This approximation is frequently referred to as the Boussinesq model [3], or (less commonly) the Oberbeck-Boussinesq model [4]. For practical applications, it is usually necessary to solve the Boussinesq model in the vicinity of complicated geometries, using unstructured meshes. As a result, our preference is to use finite element methods for solving the model because of their ability to operate on both structured and unstructured meshes, while simultaneously achieving high-order accuracy, stability, and robustness.

In what follows, we briefly review some previous efforts to apply finite element methods to the Boussinesq model. Some of the earliest work in this area was performed by Laskaris [5] who used a high-order continuous Galerkin (CG) method to simulate channel flows with heated walls. In addition, Young et al. [6, 7] and Tabarrok and Lin [8] used a similar approach to study natural convection in heated cavities. Next, Gartling [9] used a CG method to simulate a thin-walled tube with wall heat transfer, a rectangular heat exchanger, and a heated hexagonal cylinder in a cooled cavity. Thereafter, Marshall et al. [10] used a high-order CG method with a penalty function (for enforcement of the dilatational constraint) to simulate a heated cavity. This was the first time that a finite element method was successfully applied to natural convection problems for a wide range of Rayleigh numbers ($10^{4}-10^{7}$). Based on this work, Reddy and Satake [11] formulated an alternative CG method, and used it to simulate heated, non-convex, straight-sided cavities. It is important to note that all of the early work described above was limited to two-dimensional geometries. Fortunately, with the advent of more powerful computers and more advanced stabilization strategies, such as the Galerkin Least Squares (GLS) approach [12, 13, 14], the solutions to three-dimensional problems became possible. Some of the early work in this area was performed by Tang and Tsang [15, 16], who used least-square finite element methods to simulate three-dimensional heated cavities, and accurately reproduce the dynamics of Rayleigh-Bénard convection cells. A detailed review of the latest efforts to apply finite element methods to natural and mixed convection problems is beyond the scope of the present article. However, the interested reader may consult [17, 18] for an extensive collection of references on this topic.

Despite the many applications of finite element methods to the Boussinesq model, there have been a relatively small number of efforts to rigorously analyze the existing methods, or to develop new mixed methods which maintain inf-sup stability. Some pioneering efforts in this area were undertaken by Boland and Layton [19, 20], as they derived stability and error estimates for CG methods for steady and unsteady natural convection problems. In addition, they analyzed low-order, non-conforming discontinuous Galerkin (DG) methods. Most notably, they were among the first researchers to recognize the importance of using a skew-symmetrizing procedure to stabilize the convective operator in the temperature equation. Subsequently, their work was expanded by Dorok et al. [21] and Bernardi et al. [22], who developed stability and error estimates for mixed methods. More recently, Codina et al. [23] and Löwe and Lube [24] developed variational multiscale (VMS) methods for problems with turbulent mixed convection. Within the VMS framework, they constructed rigorous stability estimates and (in the case of [24]) error estimates for the resulting schemes. Thereafter, Dallmann and Arndt [18, 4] developed a mixed method which was stabilized using a combination of local projection stabilization [25, 26], streamline-upwind stabilization [27, 28, 29], and grad-div stabilization [30]. For this method, they rigorously derived stability and error estimates, and produced accurate numerical results for a wide range of steady and unsteady convection problems. Next, Rebollo et al. [31] developed a mixed method which they stabilized using an interpolation-based operator that acts as a low-pass filter. We note that, although the performance of this method is quite adequate from an accuracy standpoint, it is only weakly consistent. Most recently, de Frutos et al. [32] derived an optimal set of stability and error estimates for grad-div stabilized, inf-sup stable mixed methods. These methods are effectively a subset of the methods constructed by Dallmann and Arndt in [18, 4]. Lastly, we note that there are ongoing efforts to analyze mixed methods for Boussinesq models with temperature-dependent parameters (cf. [33, 34, 35, 36] for several examples).

Due to the limited number of efforts to develop mixed methods (see above), there are still opportunities to improve the robustness, accuracy, and flexibility of the methods. With this in mind, the goal of the present paper is to extend the recently developed versatile mixed methods (see [37]) to solve the Boussinesq model with constant parameters. For the sake of completeness, let us briefly describe the underlying philosophy of versatile mixed methods: i) we begin with the compressible formulation of the governing equations and then enforce the assumption of constant density, ii) we maintain the presence of dilatational terms (and similar terms) that would usually be neglected, and iii) we discretize the resulting formulation using standard, inf-sup stable, mixed methods. This approach has several advantages, as most importantly, it can be immediately applied to weakly-compressible flows, and furthermore, it ensures that the dilatational constraint is enforced in a consistent fashion in each of the governing equations. In [37], this philosophy was applied to the isothermal incompressible Navier-Stokes equations. There, we used the full compressible stress tensor (with the dilatational component) in the momentum conservation equation, and we rigorously proved the stability of the discrete velocity field. The resulting methods were successfully applied to isothermal Taylor-Green and Gresho vortex problems. In this work, we apply the same methods to non-isothermal incompressible flows.

The format of this paper is as follows. In section 2, we formally introduce the Boussinesq model equations for non-isothermal incompressible flows and we develop the notation and mathematical machinery for discretizing these equations. In sections 3 and 4, we introduce the versatile mixed methods, and prove the stability of the discrete temperature field. In section 5, we apply the methods to a set of standard benchmark problems involving natural and mixed convection. Finally, in section 6, we conclude with a summary of our work and a few final remarks.

Let us start by introducing a domain $\Omega_{t}=\left(0,t_{n}\right)\times\Omega$, where $\Omega\in\mathbb{R}^{d}$ is a spatial domain and $\left(0,t_{n}\right)\in\mathbb{R}$ is a temporal domain. In a natural fashion, we denote the spatial and temporal coordinates by $\bm{x}$ and $t$, and we denote the spatial and temporal derivatives by $\nabla\left(\cdot\right)$ and $\partial_{t}\left(\cdot\right)$, respectively. We assume $d=2$ or 3, and that the domain boundary $\partial\Omega$ is composed from straight line segments (for the case of $d=2$) and planar faces (for the case of $d=3$). Inside the domain $\Omega_{t}$, we are interested in simulating the motion of a homogeneous, non-isothermal, incompressible fluid with a constant density $\rho_{0}$, and non-constant velocity, temperature, and pressure fields $\bm{u}=\bm{u}\left(t,\bm{x}\right)$, $T=T\left(t,\bm{x}\right)$, and $p=p\left(t,\bm{x}\right)$. Since the density is constant, we find it convenient to divide the governing equations by $\rho_{0}$, and then introduce density-weighted quantities, such as $\widetilde{p}=p/\rho_{0}$ (the kinematic pressure). We introduce the tilde symbol to avoid abuses of notation which can result from ignoring differences between density scaled and unscaled quantities. Now, having established the necessary background, we present the Boussinesq model for non-isothermal flows

These equations are subject to the following boundary and initial conditions

Furthermore, in order to close the equations, we define $\widetilde{\bm{\tau}}$ as the stress tensor

$\bm{g}$ as the gravitational acceleration (where $g_{i}=-g\delta_{id}$ with $g=const$), $\widetilde{\bm{f}}_{\bm{u}}$ as a source term for the linear momentum, $\widetilde{f}_{T}$ as a source term for the temperature, $C_{v}$ as the specific heat at constant volume, $C_{p}$ as the specific heat at constant pressure, $\gamma=C_{p}/C_{v}$ as the ratio of specific heats, $\alpha=\kappa/\left(C_{p}\,\rho_{0}\right)$ as the thermal diffusivity coefficient, $\beta$ as the thermal expansion coefficient, $\kappa$ as the thermal conductivity coefficient, $\nu=\mu/\rho_{0}$ as the kinematic viscosity coefficient, and $\mu$ as the dynamic viscosity coefficient.

Before proceeding further, it is important to note that our equations for the temperature and the stress tensor (Eqs. (2.3) and (2.8)) are unconventional. In particular, it is common practice to neglect the viscous dissipation and pressure work terms on the RHS of Eq. (2.3), such that

In addition, most researchers neglect the divergence and gradient transpose terms on the RHS of Eq. (2.8), as follows

However, we prefer to use Eqs. (2.3) and (2.8) due to their superior physical accuracy, flexibility, and discrete consistency. We refer the interested reader to [37] for a detailed discussion of our motivation for using the full stress tensor (Eq. (2.8)). In what follows, we will only discuss our motivation for using the full temperature equation (Eq. (2.3)).

1. 1.

   The formulation in Eq. (2.3) contains the viscous dissipation term, and thereby successfully captures the physical conversion of kinetic energy into internal energy (heat). Of course, the viscous dissipation term will be small in most incompressible flows, however it will rarely completely vanish. Therefore, by neglecting this term, we introduce a small but unnecessary amount of error into the final solution. Furthermore, this error is difficult to control, as it does not vanish in the asymptotic limit as the element size goes to zero, or the polynomial order goes to infinity.

2. 2.

   The formulation in Eq. (2.3) is more suitable for adaptation to compressible flows, as it retains the pressure work and viscous dissipation terms which become increasingly important in these types of flows. Retaining these terms helps facilitate flexibility of the resulting methods, and encourages code-reuse between incompressible and compressible CFD codes.

3. 3.

   The formulation in Eq. (2.3) is more ‘consistent’, as it enables consistent enforcement of the dilatational constraint. In order to see this, we begin by noting that Eq. (2.3) retains the pressure work term, which is guaranteed to vanish at the continuous level (by Eq. (2.1)), but which may or may not vanish at the discrete level. Evidently, for pointwise divergence-free methods, the pressure term vanishes in both cases, but for more general methods, the dilatation term typically only vanishes in the weak sense, and the pressure term is non-zero. Therefore, neglecting the pressure term a priori is inconsistent, as this effectively forces the dilatation contribution to vanish pointwise in the temperature equation, even though it may only vanish weakly in the mass conservation equation. Naturally, we prefer to use Eq. (2.3), as it avoids this inconsistency.

In summary, we have introduced a ‘versatile’ approach in which we solve Eqs. (2.1)–(2.3) in conjunction with the stress tensor in Eq. (2.8). In what follows, we will introduce the necessary machinery for discretizing these equations.

In accordance with standard practices, we tessellate the spatial domain $\Omega$ with a mesh $\mathcal{T}_{h}$. The mesh is composed from straight-sided, $d$-dimensional triangular or cubic elements $K$, with characteristic size $h$. The faces of elements on the perimeter of the mesh are required to exactly conform to the domain boundaries, and the union of all the elements is required to cover the domain. In addition, for the sake of simplicity the elements are required to be non-overlapping, and the mesh is required to be devoid of hanging nodes. The boundary of each element $K$ is denoted by $\partial K$ and the outward-pointing unit normal vector on this boundary is denoted by $\bm{n}$. Elements are considered to be ‘face neighbors’ if they share a $(d-1)$-dimensional face $F$. We denote the unit normal vector that points from the positive side to the negative side of the shared face as $\bm{n}_{+}$, and naturally $\bm{n}_{-}=-\bm{n}_{+}$. The total collection of faces in the mesh is denoted by $\mathcal{F}_{h}$, and the faces of a single element $K$ are denoted by $\mathcal{F}_{K}=\left\{F\in\mathcal{F}_{h}:F\subset\partial K\right\}$. The set of interior faces is denoted by $\mathcal{F}_{h}^{i}=\{F\in\mathcal{F}_{h}:F\cap\partial\Omega=\emptyset\}$ and the set of boundary faces by $\mathcal{F}_{h}^{\partial}=\{F\in\mathcal{F}_{h}:F\cap\partial\Omega\neq\emptyset\}$. Finally, for a given face $F$, we can define a normal vector $\bm{n}_{F}$ which points from the positive to the negative side of the face.

Next, one may define jump $\left[\!\left[\cdot\right]\!\right]$ and average $\left\{\!\left\{\cdot\right\}\!\right\}$ operators for an interior face $F\in\mathcal{F}_{h}^{i}$ as follows

where $\phi$ is a generic scalar function, and $\bm{v}$ is a generic vector function. Similarly, for all boundary faces $F\in\mathcal{F}_{h}^{\partial}$, one may define

In addition, it is convenient to introduce some standard notation for representing inner products. With this in mind, let us introduce a generic vector $\bm{w}$ and generic tensors $\bm{T}$ and $\bm{U}$. Note: here, we assume that $\bm{v}$, $\bm{w}$, $\bm{T}$, $\bm{U}$, and $\phi$ are sufficiently smooth, such that the associated integrations are possible. Based on this assumption, we can define

Using this notation, we can introduce the well-known integration by parts formulas

which can be rewritten as

In what follows, we will conclude this section by defining the standard function spaces for mixed finite element methods. We start by introducing the broken Sobolev space

where $\bm{W}^{m,p}\left(\mathcal{T}_{h}\right)=\left(W^{m,p}\left(\mathcal{T}_{h}\right)\right)^{d}$. Next, we introduce the Hilbert spaces

where $\bm{H}^{1}\left(\Omega\right)=\left(H^{1}\left(\Omega\right)\right)^{d}$. Having established these spaces, we can define scalar-valued polynomial spaces $Q_{h}^{DC}$ and $Q_{h}^{C}$ for the pressure, and $R_{h}^{C}$ for the temperature

where $P_{k}\left(K\right)$ is the space of polynomials of degree $\leq k$, and $L^{2}_{\ast}\left(\Omega\right)$ is the space of $L^{2}$ functions with zero mean. Furthermore, we can define the vector-valued Raviart-Thomas and Taylor-Hood spaces for the velocity

where $\bm{C}^{0}\left(\Omega\right)=\left(C^{0}\left(\Omega\right)\right)^{d}$, and

Lastly, we can introduce $\bm{W}_{h}^{BDM}$, the Brezzi-Douglas-Marini space (see [38] for an explicit definition of this space).

In this section, we develop a general class of mixed methods for solving Eqs. (2.1) – (2.3). The methods can be constructed using the following steps: i) choose function spaces $Q_{h}\subset L^{2}_{\ast}\left(\Omega\right)$, $R_{h}\subset H^{1}_{0}\left(\Omega\right)$, and $\bm{W}_{h}\subset\bm{H}_{0}(\text{div};\Omega)$, ii) identify test functions $\left(q_{h},r_{h},\bm{w}_{h}\right)$ that span $Q_{h}\times R_{h}\times\bm{W}_{h}$, and iii) find unknowns $\left(\widetilde{p}_{h},T_{h},\bm{u}_{h}\right)$ in $Q_{h}\times R_{h}\times\bm{W}_{h}$ that satisfy

Here, we note that the quantities with hats (for example $\widehat{\bm{\sigma}}_{\text{inv}}$) denote numerical fluxes. An array of possible formulas for the fluxes are given in section A.1 of the Appendix. By substituting these formulas into Eqs. (3.1) – (3.3), one may rewrite the equations in standard form as follows

Next, we must define the operators $a_{h}$, $b_{h}$, $c_{h}$, $\underline{a}_{h}$, $\underline{c}_{h}$, $\Xi$, $\Phi$, and $\Psi$. In order to setup these definitions, we introduce functions $q_{h}\in Q_{h}$, $r_{h}$ and $\theta_{h}\in R_{h}$, and $\bm{v}_{h},\bm{w}_{h}$ and $\bm{\xi}_{h}\in\bm{W}_{h}$. Thereafter, we expand the operators in Eqs. (3.4) and (3.5) as follows

In addition, the operators in Eq. (3.6) can be expanded as follows

It is possible to simplify these expressions in the particular case when the method is pointwise divergence-free. One may consult section A.2 of the Appendix for details.

In this section, we rigorously analyze the stability of the versatile mixed methods which were introduced in section 3. In preparation for this analysis, we first define a special set of norms on broken Sobolev spaces. Thereafter, we establish the coercivity of the bilinear form $\underline{a}_{h}$ (Eq. (3.12)) and the semi-coercivity of the trilinear form $\underline{c}_{h}$ (Eq. (3.11)). Next, we use these results to prove the L2-stability of the discrete temperature field. Finally, we discuss the relationship between the stability properties of the discrete temperature and velocity fields.