{'fraction_non_alphanumeric': 0.04410002396740433, 'fraction_numerical': 0.015259247423504034, 'mean_word_length': 4.335677749360614, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Versatile Mixed Methods for Non-Isothermal Incompressible Flows

The purpose of this paper is to extend the versatile mixed methods originally developed by Chen and Williams for isothermal flows in"Versatile Mixed Methods for the Incompressible Navier-Stokes Equations,"Computers&Mathematics with Applications, 2020, (under review), to simulate non-isothermal incompressible flows. These new mixed methods are particularly interesting, as with only minor modifications they can be applied to a much broader range of flows, including non-isothermal weakly-compressible flows, and fully-compressible flows. In the main body of this paper, we carefully develop these mixed methods for solving the Boussinesq model equations. Thereafter, we prove the L2-stability of the discrete temperature field, and assess the practical behavior of the methods by applying them to a set of well-known convection problems.


Introduction
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 twodimensional 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.

Preliminaries
Let us start by introducing a domain Ω t = (0, t n ) × Ω, where Ω ∈ R d is a spatial domain and (0, t n ) ∈ R is a temporal domain. In a natural fashion, we denote the spatial and temporal coordinates by x and t, and we denote the spatial and temporal derivatives by ∇ (·) and ∂ t (·), respectively. We assume d = 2 or 3, and that the domain boundary ∂Ω is composed from straight line segments (for the case of d = 2) and planar faces (for the case of d = 3). Inside the domain Ω t , we are interested in simulating the motion of a homogeneous, non-isothermal, incompressible fluid with a constant density ρ 0 , and non-constant velocity, temperature, and pressure fields u = u (t, x), T = T (t, x), and p = p (t, x). Since the density is constant, we find it convenient to divide the governing equations by ρ 0 , and then introduce density-weighted quantities, such as p = p/ρ 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 in Ω t (2.1) These equations are subject to the following boundary and initial conditions in Ω, (2.6) Furthermore, in order to close the equations, we define τ as the stress tensor g as the gravitational acceleration (where g i = −gδ id with g = const), f u as a source term for the linear momentum, 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, γ = C p /C v as the ratio of specific heats, α = κ/ (C p ρ 0 ) as the thermal diffusivity coefficient, β as the thermal expansion coefficient, κ as the thermal conductivity coefficient, ν = µ/ρ 0 as the kinematic viscosity coefficient, and µ 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 what follows, we will introduce the necessary machinery for discretizing these equations.
In accordance with standard practices, we tessellate the spatial domain Ω with a mesh 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 nonoverlapping, and the mesh is required to be devoid of hanging nodes. The boundary of each element K is denoted by ∂K and the outward-pointing unit normal vector on this boundary is denoted by 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 n + , and naturally n − = −n + . The total collection of faces in the mesh is denoted by F h , and the faces of a single element K are denoted by Finally, for a given face F , we can define a normal vector n F which points from the positive to the negative side of the face.

Next, one may define jump [[·]] and average {{·}} operators for an interior
where φ is a generic scalar function, and v is a generic vector function. Similarly, for all boundary faces F ∈ F ∂ h , 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 w and generic tensors T and U . Note: here, we assume that v, w, T , U , and φ 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 for- 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 Next, we introduce the Hilbert spaces where P k (K) is the space of polynomials of degree ≤ k, and L 2 * (Ω) 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 Lastly, we can introduce W BDM h , the Brezzi-Douglas-Marini space (see [38] for an explicit definition of this space).

Versatile Mixed Methods
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: Here, we note that the quantities with hats (for example σ 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.
Next, we must define the operators a h , b h , c h , a h , c h , Ξ, Φ, and Ψ. In order to setup these definitions, we introduce functions q h ∈ Q h , r h and θ h ∈ R h , and v h , w h and ξ h ∈ 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 (3.14) 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.

Analysis of Versatile Mixed Methods
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 a h (Eq. (3.12)) and the semicoercivity of the trilinear form 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.

Norm Definitions
Definition 4.1 (Gradient Norm). Consider the scalar-valued function r ∈ W 1,p (T h ). Then, , is a norm on Ω. In a similar fashion, for the vector-valued function w ∈ W 1,p (T h ), we have is a norm on Ω. where C I = ε − C 2 tr,2 N ∂ / (1 + ε) is a positive constant independent of h.

Analysis of Bilinear and Trilinear Forms
Proof. The proof appears in [39]
Proof. One may begin by substituting θ h = r h into Eq. (3.11) as follows Next, we note that the following identity holds Upon substituting this identity into Eq. (4.6), one obtains Finally, on substituting the definition of the seminorm into Eq. (4.7), we obtain the desired result (Eq. (4.4)). where p h (t) ∈ L 2 t 0 , t n ; L 4 (Ω) , and a discrete velocity field u h ∈ W h where u h (t) ∈ L 2 t 0 , t n ; W 1,4 (T h ) . Subject to these assumptions, the discrete temperature T h ∈ R h is governed by the following equation at time t n ≥ t 0

Analysis of Discrete Stability
, (4.8) where C I and C II are constants that are independent of h, and are seminorms and norms on (t 0 , t n ) × Ω.
Proof. We start by setting r h = T h in Eq. (3.6) as follows Next, we invoke the coercivity of a h (Lemma 4.3) and the semi-coercivity of c h (Lemma 4.4) as follows Based on this equation, we observe that

and equivalently, by the Cauchy-Schwarz and Triangle inequalities
.
Upon expanding the RHS of this expression, and using the Triangle and Cauchy-Schwarz inequalities again, we obtain . (4.14) Note: on the last line we have used the broken norm inequalities from Lemma A.1 in the Appendix. Next, we integrate Eq. (4.14) from t = t 0 to t = t n as follows , or equivalently, after applying Holder's inequality .

(4.15)
We will utilize this result shortly. For now, we turn our attention back to Eq. (4.13). On integrating this equation from t = t 0 to t = t n , we find that ds.

(4.16)
We can rewrite the last term on the RHS of Eq. (4.16) as follows .
Here, we have used the Cauchy-Schwarz inequality, the Triangle inequality, and Eq. (4.15). Next, we bound the remaining terms in the integrand above (employing the same techniques that we used to derive Eq. (4.15)), and we obtain , (4.17) where C II = 4d 2 1 +  suppose we impose a forcing function f T where f T (t) ∈ L 1 t 0 , t n ; L 2 (Ω) and a discrete velocity field u h ∈ W h where u h (t) ∈ L 2 t 0 , t n ; W 1,4 (T h ) . Subject to these assumptions, the discrete temperature T h ∈ R h is governed by the following equation at time t n ≥ t 0 .
Subject to these assumptions, the velocity field is governed by the following equation at time t n ≥ t 0 where C III is a constant independent of h, and is a norm on (t 0 , t n ) × Ω.
Proof. The important aspects of the proof are standard, and follow the arguments in Lemma 3.1 of [40] and Theorem 7.1 of [37].
Remark 4.1. The stability of the discrete temperature field depends on the stability of the discrete velocity field, as shown by Theorem 4.5 and Corollary 4.6.
Conversely, the stability of the discrete velocity field depends on the stability of the discrete temperature field, as shown by Theorem 4.7. Therefore, it is difficult to establish independent stability of either field, and (theoretically speaking) this may result in undesirable interference between the two fields.
Fortunately, in most cases the coupling between fields is weak as one of the following assumptions holds true: • The buoyancy term (with coefficient 5β 2 h g 2 ) on the right hand side of Eq. (4.19) is small.
• The viscous dissipation term (with coefficient 7C II ν 2 h /C 2 v ) on the right hand side of Eq. (4.8) or Eq. (4.18) is small.

Numerical Experiments
In this section, the results of several numerical simulations are presented to demonstrate the performance of the proposed methods. The following simulations were all performed using Taylor-Hood elements with polynomials of degree k, k+1, and k+1 for the pressure, temperature, and velocity spaces respectively; i.e. for cases with k = 1 the polynomials for each space were degree 1, 2, and 2 respectively. In addition, we imposed a zero integral mean condition for the pressure via a Lagrange multiplier. The convective numerical fluxes were computed using upwind biased fluxes with ζ = δ = 0.5, and the viscous numerical fluxes were computed using η = ε = 3(k + 1)(k + 2). In each case, a high-order

Order of Accuracy Test
For the first test case, we compared solutions from our method to an exact solution in order to check the convergence rate. To this end, we considered the traveling temperature wave proposed by [18], which can be defined on the rectangular domain Ω = [−0.5, 1.5] × [0, 1] as follows u = (100, 0) , p = 1, for t ∈ [0, 0.005]. We also defined the gravitational acceleration and forcing functions as follows Here, we considered a dimensionless formulation with α = β = γ = ν = ρ = 1.
In subsequent cases, a dimensional formulation was considered.
For this case, a temperature peak was initially located at x = 1 2 at t = 0 and moved to x = 1 at the final time t = 0.005. We compared our results to the exact solution at the final time. The time step for all polynomial degrees considered was ∆t = 1 × 10 −5 . Periodic boundary conditions were applied at the domain boundaries. The meshes were uniform and consisted of N × N 2 elements. On these meshes, we utilized Taylor-Hood spaces with degrees k = 1, 2, and 3. For the case of k = 1, mesh resolutions of N = 16, 32, 64, and 128 were used, and for k = 2 and 3 mesh resolutions of N = 4, 8, 16, and 32 were used. We expected the discrete temperature to converge at a rate of k + 2 since the associated polynomial space was degree k + 1. We see from table 1 that we recovered the expected orders of accuracy.

Heated Cavity Test
The second test case was a heated cavity as described by [18]  Specifically, we decided to vary the Rayleigh number by varying the parameter β, using the following formulas where L is the width of the cavity. We were interested in computing the average steady state Nusselt number N u based on the horizontal heat flux as follows where A is the domain area, and u x is the velocity in the x direction. The Nusselt number was calculated at Rayleigh numbers of Ra = 10 4 , 10 5 , 10 6 , and converged for this mesh. The mesh used for each simulation was biased towards the walls using the mapping proposed by [42].
Note: in order to generate our meshes, we used the average Nusselt numbers reported in [18]. On each mesh, Taylor-Hood elements of degree k = 2 and 3 were used.
At the lowest Rayleigh number, the flow was dominated by a large central vortex seen in figure 1. As the Rayleigh number was increased, this vortex disappeared and thin boundary layers developed on the left and right walls as seen in figures 2 and 3. This is the same behavior observed by [18].
We also saw very close agreement with the average Nusselt number for all mesh resolutions and polynomial degrees as seen in table 2. This leads us to believe that our method is able to accurately capture purely buoyancy-driven flows. The 64 × 64 mesh with k = 3 was used to generate these results. The 64 × 64 mesh with k = 3 was used to generate these results. The 64 × 64 mesh with k = 3 was used to generate these results.

Heated Cavity with Moving Wall Test
The final test case was a heated cavity with one moving wall, i.e. a mixed convection case. Here, the top wall moved at constant velocity V lid and was heated, while the bottom wall was cooled as proposed by [43]. This case was run at a fixed Grashof number Gr = 10 4 along with varied Richardson numbers We used the same fluid properties prescribed in the previous heated cavity case.  Therefore at these higher lid speeds, the velocity gradients in the flow became non-negligible and have a more pronounced effect on the temperature field (via the dissipation term). This deviation in Nusselt number was absent from the heated cavity case because the flow in that case was subject to natural convection only, and thus possessed smaller velocity gradients as compared to the mixed convection case.

Conclusion
In the present study, the mixed methods first put forward by Chen and Williams [37] are extended to non-isothermal incompressible flows. The primary advantages that these new methods possess are their generality and flexibility, as they utilize the full compressible formulation of the stress tensor and the expanded formulation of the temperature equation (which retains the dilatational pressure work and viscous dissipation terms). In this paper, the new versatile methods are constructed for weakly divergence-free Taylor-Hood elements, and pointwise divergence-free BDM and RT elements. Next, we rigorously derive a new condition that governs the L2-stability of the discrete temperature fields for these methods. Finally, the accuracy of the Taylor-Hood method is tested using three well-known cases from the literature; these tests are used to confirm the formal order of accuracy of the method, and demonstrate its performance on problems with natural and mixed convection. The analysis and numerical experiments in this article serve as a stepping stone towards the application of these methods to weakly-compressible and fully-compressible flows.

Declaration of Interests
None.

Funding Sources
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
This set of equations can be rewritten compactly as follows b h (u h , q h ) = 0, (A.10) The operators b h , Ξ, and a h were previously defined in Eqs. (3.7), (3.10), and (3.12). The remaining operators c h , a h , c h , and Φ can be written as follows

A.3 Supporting Results
Lemma A.1 (Broken Norm Inequalities). Suppose that w ∈ W 1,p (T h ) and p ≥ 2. Then, the following inequalities hold Proof. Let us begin by noting that Here, we have used the power mean inequality. Upon combining this result with the definition for the norm · grad,p , we obtain the first result (Eq. (A.17)). Then, by the power mean inequality, we have Upon combining this result with Eq. (A. 20), and the definition for the norm · grad,p , we obtain the inequality in Eq. (A.18).