Next Article in Journal
High-Fidelity Surrogate Based Multi-Objective Optimization Algorithm
Previous Article in Journal
Solar Photovoltaic Integration in Monopolar DC Networks via the GNDO Algorithm
Previous Article in Special Issue
Automatic Atrial Fibrillation Arrhythmia Detection Using Univariate and Multivariate Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simulation of Low-Speed Buoyant Flows with a Stabilized Compressible/Incompressible Formulation: The Full Navier–Stokes Approach versus the Boussinesq Model

by
Guillermo Hauke
1,*,† and
Jorge Lanzarote
2,†
1
Instituto de Ingeniería de Aragón—Universidad de Zaragoza, Area de Mecánica de Fluidos, Escuela de Ingeniería y Arquitectura, C/María de Luna 3, 50018 Zaragoza, Spain
2
Repsol, DE Química, C/Mendez Alvaro, 44, 28045 Madrid, Spain
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Algorithms 2022, 15(8), 278; https://doi.org/10.3390/a15080278
Submission received: 5 July 2022 / Revised: 27 July 2022 / Accepted: 27 July 2022 / Published: 5 August 2022

Abstract

:
This paper compares two strategies to compute buoyancy-driven flows using stabilized methods. Both formulations are based on a unified approach for solving compressible and incompressible flows, which solves the continuity, momentum, and total energy equations in a coupled entropy-consistent way. The first approach introduces the variable density thermodynamics of the liquid or gas without any artificial buoyancy terms, i.e., without applying any approximate models into the Navier–Stokes equations. Furthermore, this formulation holds for flows driven by high temperature differences. Further advantages of this formulation are seen in the fact that it conserves the total energy and it lacks the incompressibility inconsistencies due to volume changes induced by temperature variations. The second strategy uses the Boussinesq approximation to account for temperature-driven forces. This method models the thermal terms in the momentum equation through a temperature-dependent nonlinear source term. Computer examples show that the thermodynamic approach, which does not introduce any artificial terms into the Navier–Stokes equations, is conceptually simpler and, with the incompressible stabilization matrix, attains similar residual convergence with iteration count to methods based on the Boussinesq approximation. For the Boussinesq model, the SUPG and SGS methods are compared, displaying very similar computational behavior. Finally, the VMS a posteriori error estimator is applied to adapt the mesh, helping to achieve better accuracy for the same number of degrees of freedom.

1. Introduction

Thermally driven flows and buoyancy have important practical applications, such as ventilation problems, indoor comfort, thermally driven ocean currents, dynamics of plumes, and so on [1,2].
The dimensionless number that measures the importance of buoyancy forces due to thermal gradients compared to viscous forces is the Grashof number. Many times, the Grashof number is replaced by the Rayleigh number, which equals the Grashof number divided by the Prandtl number and represents the ratio of buoyancy forces to the viscous force and thermal diffusion [3].
From a computational mechanics point of view, flows at large Rayleigh numbers are viewed as hard. As the Rayleigh number increases, the flow evolves from steady laminar to unsteady laminar conditions. A further increase of the Rayleigh number makes the flow turbulent. The efficient computation of high Rayleigh numbers requires being capable of taking large time-steps and that the CPU time per time-step be small.
The literature on computational strategies to tackle thermally coupled low-speed flows is extremely vast. Almost every method and algorithm has been tested for this type of flow, both in segregated and the coupled formats. Regarding examples related to finite elements using the Galerkin method, one can mention [4,5]; using stabilized methods [6,7,8]; and using the Finite Increment Calculus (FIC) strategy [9]. Benitez and Bermudez present a second-order-characteristics finite-element method for thermally coupled incompressible flow problems in [10]. Examples of finite-element methods used to study the turbulent regime through the inclusion of turbulence models can be found in [2], whereas [11,12] employ quasi-static sub-scales to embrace turbulence, [13] represents the fine scales through dynamical bubbles and [14,15] use transient sub-scales to model fine-scales. Weak imposition of boundary conditions is studied in [16].
Other models used to simulate thermal convection are the anelastic model and the acoustically filtered Navier–Stokes equations (see [4,17] and references therein). In the anelastic model, the temporal derivative of density in the continuity equation is ignored, and, thermodynamically, the density is considered to be only a function of temperature, so it does not depend on the dynamical pressure. In the filtered Navier–Stokes model, by taking the zero Mach number limit of the compressible equations, a set of elliptic equations is obtained where the acoustic waves are removed from the model; thermodynamically, the density does not depend on the dynamical part of the pressure.
There are various unified approaches using stabilized finite-element methods to tackle compressible and incompressible flows, such as [18], but the work of this paper is based on [19,20]. Ref. [21] comments on advantages of this formulation in terms of treating temperature dependent flows with a unified approach based on pressure or entropy variables.
In this paper, two strategies to compute buoyant flows are compared: one based on the classical Boussinesq approximation and the other one based on considering the thermal compressibility of the liquid or gas. Furthermore, three stabilized methods based on a unified approach are applied to buoyancy flows and are compared and tested. Regarding the algorithms, a wide range of updated stabilization parameters are compared for this application, and conclusions are reached about the best choices.

2. Two Models for Low-Speed Buoyant Flows

In this paper, two models for the simulation of low-speed buoyant flows are compared: one based on the generalized liquid or gas thermodynamics and the other one based on the Boussinesq model.
In both cases, even in the presence of body forces, the flow dynamics are given by the Navier–Stokes equations, including an energy equation. For instance, using the total energy equation, in Cartesian coordinates the Navier–Stokes equations can be written as
ρ t + ( ρ u i ) , i = 0 ρ u i t + ( ρ u i u j ) , j = p , i + τ i j , j + ρ f m i ρ ( e + 1 2 u 2 ) t + ρ ( e + 1 2 u 2 ) u j , j = ( p u i ) , i + ( τ i j u i ) , j + ρ f m i u i q i , i + q ˙ v ,
where ρ is the density, u i the Cartesian velocity components, p the pressure, τ i j is the components of the viscous stress tensor, f m i is the body-force term, e is the specific internal energy, q i is the heat flux vector components, and q ˙ v is the volumetric heat source.
In vector form, the above system of equations can be written as
U , t + F i , i adv = F i , i diff + S ,
where U is the vector of conservation variables, F i adv is the i-th Euler Jacobian, F i diff is the i-th diffusive flux, and S is the source vector.

2.1. Generalized Liquid/Gas Model

Even though, for an incompressible substance, the density is constant, real liquids display compressible behavior, and the real density is a function of temperature and pressure. One way of modeling such real liquids is through the model of a general divariant substance, in which the thermodynamic variables depend on two independent thermodynamic variables. In [19,20,22], the thermodynamic variables are taken as the compressibility coefficients at constant pressure α p and at constant temperature β T :
β T = 1 ρ ρ p T
α p = 1 ρ ρ T p
Remarks
(i)
This model has some non-standard features with respect to the usual approaches employed to simulate buoyancy-driven flows. Following the unified approach in [19,20], the present model is based on the conservative form of the transport equations, and it uses the total energy equation instead of the temperature equation.
(ii)
This model is valid for any temperature gradient (or temperature difference between walls), as long as the equations of state and the constitutive relationships hold. Compare to the Boussinesq model, which holds up to relative temperature variations of about 10%.
(iii)
As shown in [19,20,22,23], this model is endowed with the correct second law of thermodynamics, which stems from the entropy production of the Navier–Stokes equations. This is especially relevant, since many of the models for buoyancy-driven flows violate this principle [24,25]. Furthermore, in the present model, the flow divergence may be non-vanishing; therefore, it is not approximated to be zero, as in the Boussinesq model.
(iv)
The total energy equation may be replaced by the internal energy equation or another expression for the temperature equation. However, in order for the model to be entropy consistent, the temperature or energy equation has to retain all the terms, for instance, the dissipation function and the power expansion term (see [3] for other forms of the energy or temperature equations that are consistent). These pitfalls were shown in [24,25] for many usual models employed in the simulation of buoyant flows.
(v)
A typical equation of state for a liquid is that the density is only a function of temperature, i.e., ρ = ρ ( T ) . In this case, the equations are termed as the acoustically filtered Navier–Stokes equations [4].

2.2. Boussinesq Approximation

Let us expand the density in a Taylor series with respect to temperature variations. Neglecting higher order terms (see for instance [4]):
ρ ( T ) ρ 0 + ρ T p Δ T = ρ 0 ρ 0 α p ( T T 0 ) ,
where ρ 0 is the density at the reference temperature T 0 , and α p is the thermal expansion coefficient at constant pressure.
In this model, the fluid is assumed incompressible and, therefore, of constant density ρ 0 , and the effect of density variations is accounted for in the source term. In particular, the gravity-forcing term is substituted by the buoyant forces with respect to equilibrium:
f m = g α p ( T T 0 ) ,
where g is the gravity acceleration.
Note that, if T > T 0 , the fluid density decreases with respect to ρ 0 , so the buoyant force works against g , and the fluid particle lifts up, whereas if T < T 0 , the fluid density increases, so the buoyant force acts in favor of g , and the particle drops down. Also note that in this model p is the departure from the equilibrium hydrostatic pressure [4].
Thus, in the system Equation (2) we set
(6) ρ = ρ 0 (7) ρ f m i = ρ g i α p ( T T 0 )

Entropy Production

The entropy production 𝒫 (per unit time and volume) of the Boussinesq terms in the above system of equations can be calculated as [23]
(8) 𝒫 = u T · mom . eq . term + 1 T energy eq . term (9) = u T · ρ f m 1 T ρ f m · u .
(i)
Note that the entropy production of the buoyant terms in the Boussinesq model should cancel, because gravity is a conservative field. Indeed, this is satisfied if the same buoyancy model is applied to the momentum and energy equations. This is the case of the model presented in this paper. However, there are implementations that do not respect entropy production [24,25]. This is an important trait for physical and numerical reasons. From a numerical point of view, methods can inherit stability from the discrete second law of thermodynamics. From a physical point of view, if the system of equations does not respect the second law of thermodynamics, its solutions are physically wrong [23].
(ii)
Note, with this set of equations, the velocity divergence does not appear explicitly in the total energy equation, and, therefore, this model does not suffer the entropy production pitfall of other methods, in which the expansion power appears in the temperature or internal energy equation [24,25]. Indeed, if this is the case, on the one side, the velocity divergence of the Boussinesq model is zero, whereas the real divergence does not cancel due to the changes of density.
(iii)
When using the internal energy equation or the temperature equation, the viscous dissipation and power expansion terms should be retained for consistency. As mentioned above, there are models that do not retain these terms.

3. Stabilized Formulation

Using any well-defined set of variables Y , it is possible to rewrite (2) in quasi-linear form as
A 0 Y , t + A i Y , i = ( K i j Y , j ) , i + S ,
where A 0 = U , Y , A i = F i , Y adv is the ith Euler Jacobian matrix, and K = [ K i j ] is the diffusivity matrix, where K i j Y , j = F i diff .
In this paper, special attention is paid to the choice of pressure-primitive variables,
Y = p u T ,
where p is the pressure, u i is the Cartesian velocity components, and T is the absolute temperature. The matrices and vectors for these variables can be found in [20]. This set of variables is endowed with the property that the incompressible limit is well behaved.
The quasi-linear form can be written as
𝓛 Y = S 0 ,
with the differential operators defined below
𝓛 Y = A 0 Y , t + A i Y , i ( K i j Y , j ) , i C Y 𝓛 T Y = A 0 T Y , t + A i T Y , i ( K i j T Y , i ) , j C T Y 𝓛 * Y = A 0 T Y , t + A i T Y , i + ( K i j T Y , i ) , j + C T Y 𝓛 SUPG T Y = A i T Y , i ,
and the source term written as
S = C Y + S 0
Note that the difference between 𝓛 T and 𝓛 * is the sign of the diffusive and source terms.
In this paper, the time-discontinuous space-time Streamline Upwind Petrov-Galerkin (SUPG) and Subgrid-Scale (SGS) stabilized methods [26,27] are considered. In the absence of source terms, and for linear shape functions, both methods coincide.
Consider a space-time domain, where the time interval I = ] 0 , T [ is subdivided into N intervals I n = ] t n , t n + 1 [ , n = 0 , 1 , , N 1 . We define for each time interval Q n = Ω × I n and P n = Γ × I n , where Ω is the spatial domain and Γ its boundary. Finally, the “slab” Q n is decomposed into the elements Q n e , e = 1 , 2 , , ( n el ) n .
Following [19,20], the variational formulation is defined for the set of variables Y . Within each Q n , n = 0 , 1 , , N 1 , find Y 𝓢 Y such that W 𝓥 Y :
Q n ( W , t · U ( Y ) W , i · F i adv ( Y ) + W , i · K i j Y , j W · S ) d Q
+ Ω ( W ( t n + 1 ) · U ( Y ( t n + 1 ) ) W ( t n + ) · U ( Y ( t n ) ) ) d Ω
+ e = 1 ( n el ) n Q n e ( L W ) · τ ( 𝓛 Y S ) d Q
= P n W · ( F i adv ( Y ) + F i diff ( Y ) ) n i d P
The first and last integrals constitute the Galerkin terms expressed as functions of the variables Y , again written in conservative form to ensure that the weak solution is bestowed with the correct Rankine–Hugoniot relations. The jump term is written with the help of the right and left limits
W ( t n ± ) = lim ϵ 0 ± W ( t n + ϵ )
The stabilizing term is written in terms of the differential operators 𝓛 and L , which have been defined previously. According to the choice of L , Table 1 shows the finite-element methods that can be recovered.
In the case that the source term is independent of the unknown variables and the shape functions are linear, 𝓛 * = 𝓛 SUPG T = 𝓛 T . Note that, when entropy variables are used, 𝓛 ˜ = 𝓛 ˜ T because of the symmetry of the coefficient matrices and the symmetric form is recovered.
We assume
τ = Y , V τ ˜ ,
where τ ˜ is the stabilization matrix for entropy variables [27,28]. For definitions of the τ matrix for low-speed flows based on primitive variables, see [20,29], and for compressible flows, [30,31,32,33]. Other proposals for τ for a low-speed limit can be found in [34]. As shown in Codina et al. [14], the buoyancy terms in C need not intervene into τ ; therefore, the above definitions for incompressible flows can be used directly for the Boussinesq model. Masud et al. [13,35] present tau definitions based on the condensation of bubble information within the subgrid-scales, where off-diagonal terms arise naturally. Another relevant modification to compute the continuity stabilization parameter from the momentum parameter is presented in [36,37,38,39,40,41].
However, for compressible flows [42] suggests including the acceleration of gravity in the stabilization matrix. For more complex stabilization parameter computation for compressible flows, an iterative algorithm for computing tau is proposed in [43].
Note that the discontinuity-capturing operator has been omitted, since we are dealing with low speed flows [44].

4. Numerical Example Preliminaries

In order to solve compressible and incompressible flows with the same code, in this work, the set of pressure-primitive variables is chosen:
Y = p u T
For the numerical examples presented below, the full Navier–Stokes (full NS) model uses the perfect gas equation of state. For the Boussinesq model, this equation of state is substituted by constant density. Note that the equation of state ρ = ρ ( T ) , would yield approximated acoustically filtered compressible equations [4].
In this work, the following stabilization matrices are analyzed and compared:
  • The compressible “classic” tau matrix transformed into pressure primitive variables [20,27];
  • The incompressible diagonal tau presented in [20], which is an extension of [45];
  • The incompressible non-diagonal tau [20], referred to here as HH, also an extension of [45];
  • The incompressible Polner tau, which is an extension of the above definition, with the parameter ω = 0.6 [29];
  • The compressible Polner tau, an extension of the above definition [33].
The incompressible taus defined above have been updated with the newer continuity-stabilization parameter, τ c , based on the inverse of the momentum-stabilization parameter, τ m , which was presented and proposed in these works [36,37,38,39,40,41]. See Appendix A for the structure of the above incompressible stabilization matrices.
Other stabilization matrices are possible. See, for instance, the stabilization matrix derived from residual-free bubbles [13,35], where the off-diagonal terms come up naturally. These terms will be shown later to be fundamental in buoyancy-driven flows.

5. Buoyancy-Driven Square Cavity

The purpose of this section is to compare, from the computational point of view, three strategies: the full Navier–Stokes equations (solved with the SUPG method) and the Boussinesq model solved with the SUPG and the SGS methods.
The definition of this problem is taken from [46,47]. The flow in a closed square cavity is impulsed by the temperature difference between the vertical walls (see Figure 1). The dimensionless parameters that define this problem are the temperature semi-difference, defined as
ϵ = T h T c T h + T c = Δ T 2 T 0 ,
and the Rayleigh number, defined as
Ra = 2 ϵ g L 3 ν 0 α 0 ,
which for the Boussinesq model can be taken as
Ra = g β Δ T L 3 ν 0 α 0 ,
where T 0 = ( T h + T c ) / 2 is the average vertical walls temperature, and the subindex 0 denotes the variables computed at T 0 . Here, ν = μ / ρ stands for the kinematic viscosity and α = κ / ( ρ c p ) the thermal diffusivity. The above definitions usually are written as a function of β = α p ,
Additionally, the Rayleigh number can be written as a function of the Prandtl number:
Ra = Pr   2 ϵ g L 3 ν 0 2 = Pr   g Δ T L 3 ν 0 2 T 0
The Boussinesq model, which gives anti-symmetric solutions, is valid up to ϵ < 0.1 [48]. For larger temperature differences, one can talk about large temperature differences and the full Navier–Stokes equations or a low-Mach-number model should be used [11,49,50].
Two temperature differences are considered in this section, ϵ = 0.01 and ϵ = 0.6 . The fluid is air with Pr = 0.71 . Quadrilateral bilinear elements are employed for this test case. Benchmark solutions for these cases can be found in [46,47,49,50,51,52,53,54].

5.1. Residual Convergence Analysis

The baseline test case is run in a 70 × 70 uniform mesh with Δ t = 0.1 , a Generalized Minimal Residual (GMRES) tolerance of 10 8 with a maximum of 60 iterations with 9 restarts and an ILU preconditioner. The classic tau allowed a larger time-step for ϵ = 0.6 , so a value of Δ t = 1.0 was chosen for this case. The viscosity has been assumed to be constant.
Figure 2 compares the residual convergence for various stabilization matrices for the full Navier–Stokes equations. One can conclude that the diagonal tau should be avoided due to slow convergence or even divergence for large ϵ . This result stresses the importance of off-diagonal terms in the stabilization matrix, which is in agreement with [13,35]. For any ϵ , the HH, Polner, and compressible Polner taus show a similar behavior. Regarding the classic stabilization matrix, it is observed that the residual converges more slowly than the method with the incompressible taus.
Figure 3 compares, for the Boussinesq model and the SUPG finite element method, the residual convergence of the incompressible stabilization matrices. From the results, one can conclude that the diagonal tau should be avoided, since, for the given computing parameters, it does not converge for small ϵ . Again, this result stresses the importance of off-diagonal terms [13,35]. For any ϵ , the HH and Polner taus show a similar behavior. For the incompressible model, recall that the compressible Polner tau is the same as the Polner tau.
The Boussinesq model solved with the SGS finite element method is evaluated in Figure 4, which compares the residual convergence for the incompressible stabilization matrices. Again, the diagonal tau should be avoided, since, for the given computing parameters, it shows no converge for small ϵ , stressing the importance of off-diagonal terms [13,35]. For small ϵ , the Polner tau presents a better convergence, whereas for large ϵ , both Polner and HH taus behave similarly. As before, for the incompressible model, the Polner tau is the same as the compressible Polner tau.
Figure 5 compares, for the Boussinesq model and both the SUPG and SGS methods, the effect of GMRES tolerance. The Boussinesq model allows further CPU time savings by increasing the GMRES tolerance. This is not the case for the full Navier–Stokes model, which requires a tolerance of 10 8 .
Finally, Figure 6 compares the residual convergence for various models and methods. For the Boussinesq model, both SUPG and SGS with either the HH and the Polner taus, behaved similarly; SUPG perhaps showed faster convergence than SGS for small ϵ . Regarding the full Navier–Stokes equations, both HH and the Polner tau behave identically and produce a faster convergence than the classic tau, although the classic tau improves convergence as the temperature difference is increased. Note that, for the same computing parameters and small ϵ , the Boussinesq model converges faster than the full Navier–Stokes equations. This advantage tends to disappear as ϵ is increased.

5.2. Relevance of the Body-Force Term within the Full Navier–Stokes Tangent Matrix

For the full Navier–Stokes equations, one may be tempted to neglect the contribution of the body-force term inside the tangent of the system of equations. However, this contribution depends on the density and velocity components, and, therefore, it is not zero. Thus, for the same computational parameters of the above subsection, Figure 7 shows the residual convergence for the SUPG method and the HH stabilization matrix, with and without this contribution. The improvement in residual convergence is remarkable when the tangent of the body-force term is taken into consideration for all Ra and ϵ numbers. This fact was also observed in [55] for non-Newtonian fluids, when the whole Jacobian is computed with finite differences.
The tangent of the body-force term can be found in Appendix B, both for the full Navier–Stokes equations and the Boussinesq model.

5.3. Case R a = 10 6 , ϵ = 0.01 , Constant Viscosity

In this section, the present solutions are analyzed and compared with the reference solutions [46,52,56], in particular, the Nusselt number along the hot vertical wall. Vahl Davis [46] presents extrapolated values obtained with a finite-difference method for the streamfunction-vorticity formulation applied to the Boussinesq model with up to 81 × 81 nodes. Le Quéré [56] solves the same equations with a spectral Chebysev Galerkin method (with up to 72 × 72 degrees of freedom per variable), and Vierendeels [52] uses a finite volume method to solve the full Navier–Stokes equations on a 512 × 512 stretched grid. Masud and coworkers use a stabilized finite element method applied to the Boussinesq model with stabilization based on bubble functions and 80 × 80 grids.
For the three methods and the HH stabilizing matrix, Figure 8 shows the convergence of the minimum, maximum and average Nusselt numbers along the vertical hot wall as a function of the nodes on the side n np , where n np is the number of nodal points. It can be seen that, for small ϵ , the full NS equations and the Boussinesq models give very similar solutions, as shown in many publications. The present methods for the given mesh resolution provide Nusselt numbers between those reported by Vahl Davis and Le Quéré.
Table 2 compares the present results on a uniform mesh of 140 × 140 elements with various published reference solutions. The Nusselt numbers computed with full NS and the classic tau converge from below. The Boussinesq model computed with the HH and Polner taus give very similar Nusselt numbers.
In summary, it can be concluded that both SUPG and SGS give very similar solutions, so there is not a favorite method. Regarding the stabilization matrices, in this regime, the incompressible HH and Polner taus are the best choices for all models, including the compressible Navier–Stokes equations, because they attain a faster convergence of the Nusselt numbers with respect to the number of degrees of freedom.

5.4. Case R a = 10 6 , ϵ = 0.01 , Constant Viscosity with Adaptivity

In order to adapt the mesh, the Variational Multiscale (VMS) a posteriori error estimator is applied to this problem [57,58,59]. In particular, the version for compressible-incompressible flows is found in [60]. As shown in [61], the adapted meshes attained with this method are very smooth.
Figure 9 shows the first two adapted meshes obtained with the naïve error estimator, driven by the Polner error-stabilization matrix, the L 1 norm over the new mesh and an error tolerance 1.0 × 10 3 , taking the global minimum element size among all dependent variables. The mesh size was pruned within the interval [ 1 × 10 4 , 60 ] . Indeed, the mesh-size evolution in the adapted mesh is very smooth.
Figure 10 shows the better accuracy obtained with the adapted meshes for a much smaller number of degrees of freedom. Indeed, the dimensionless Nusselt numbers obtained with the coarsest adapted mesh basically match the results obtained with the finest mesh. Thus, VMS adaptivity succeeds at improving accuracy with respect to a solution computed in a uniform mesh.

5.5. Case R a = 10 6 , ϵ = 0.6 , Sutherland Viscosity

In this section, the present solutions are analyzed and compared with other works for R a = 10 6 and ϵ = 0.6 . This benchmark is described in [47,54], i.e., where the Sutherland viscosity law is taken from.
μ ( T ) μ * = T T * 3 2 T * + S T + S ,
with T * = 273 K, S = 110.5 K, and μ * = 1.65 × 10 5 kg /(m s).
Figure 11 shows the convergence of the minimum, maximum, and average Nusselt numbers along the vertical hot wall for the full Navier–Stokes equations computed with two stabilization matrices. The black line is the reference solution taken from [53,54], which solves the full Navier–Stokes equations on a 1024 × 1024 stretched grid with a maximum aspect ratio of 80. Table 3 further compares the Nusselt numbers with the benchmarks of [49,50] that utilize low-Mach-number solvers based on stabilized FEM methods on very fine stretched meshes (with about 4.2 million and 0.86 million degrees of freedom, respectively).
It can be seen that, while the method with the classic tau converges to the Nusselt numbers from above, the method with the incompressible HH tau converges from below. In general, for this test case, the Nusselt numbers computed with the present methods are a little higher than those reported in the literature.
Figure 12 compares the dimensionless Nusselt numbers obtained with the naïve (nv) and standard (st) adaptivity approaches, showing that a higher accuracy is attained than for a uniform grid with the same number of degrees of freedom. In general, the Nusselt numbers computed with both adaptation procedures have better accuracy than those obtained for a uniform grid.

6. Wind Towers

In this section, the above models are used to analyze an ancient device to climatize indoor environments in dry areas. The device is called a wind tower (see Figure 13), and it is typical of Iran [62]. Nowadays, it is also applied to advanced indoor energy-saving climatization devices, mainly in large buildings.
Typical values of ϵ encountered in the examples below are around 0.1 . Thus, the full Navier–Stokes equations will be applied together with the HH stabilization matrix, combined with the naive error estimator to adapt the meshes.

6.1. Application: Wind Tower 1

Figure 14 shows the working principle of wind towers. The free-stream is directed into the interior of a building, refreshing the interior hot air. Furthermore, a convection heat-exchange takes place inside the building, lowering the temperature sensation.
For this simulation, in a computational domain of 125 × 40 m, the following data are used to approximate the boundary conditions (see Figure 15). Sun radiation heats the surfaces at a temperature of T0 = 323 °C. Two wind-speed values are tested, U = 1 and 5 m/s, with a free-stream temperature of T = 303 °C. This corresponds to an approximate dimensionless temperature difference of ϵ 0.1 . Other free-stream values for the simulation are ρ = 1.151 kg/m3, p = 100,570 Pa.
Figure 16 shows the initial mesh ( n np = 3410 ) and the adapted mesh ( n np = 36 , 628 ) obtained with the VMS a posteriori error estimator driven by the naïve approach with the Polner error-tau, using the L 1 norm and the objective error over the new mesh with an error tolerance of 0.01 for the temperature field [60].
Figure 17 and Figure 18 compare the flow field obtained with the full Navier–Stokes equations and the Boussinesq model. For this test case and the same boundary conditions, the flow fields present some differences. The full Navier–Stokes solution shows larger eddies with somewhat smaller maximum velocity and larger temperature spreading. Therefore, the use of the full model is advised.
Indoor comfort can be evaluated by means of the apparent temperature T a , which combines information about temperature and velocity so as to indicate the human sensation of temperature. The apparent temperature is calculated according to the US National Weather Service as
T a = 13.12 + 0.6215   T ( 11.37 + 0.3965   T ) 1.23 | u | 0.16 ,
where T and T a are measured in degrees C and u in m/s.
The field of apparent temperature is shown in Figure 19 for both models, the full Navier–Stokes equations and the Boussinesq model. Figure 20 shows the apparent temperature field for two pieces of free-stream velocity data. For the higher velocity, more dramatic temperature variations due to ventilation are obtained.

6.2. Wind Tower 2

The simulations of the above subsection are repeated with the same data and boundary conditions for the geometry shown in the zoom of Figure 21. VMS adaptivity is used again to refine the mesh, yielding the mesh shown in Figure 22 of nnp = 27,449.
Figure 23 and Figure 24 compare the velocity and apparent temperature fields obtained with the full Navier–Stokes equations and the Boussinesq model. As above, the full Navier–Stokes solution shows larger eddies with somewhat smaller maximum velocities and larger temperature spreading.

7. Conclusions

The thermodynamic method presented is consistent with entropy production and the second law of thermodynamics, which provides discrete entropy stability for the numerical method. This model is valid for any temperature difference. Combined with an incompressible stabilization matrix, such as HH or Polner, it attains reasonable residual convergence speed with respect to CPU-time when compared to the Boussinesq model. Accounting for the tangent of the source, the terms accelerate the convergence rate.
The Boussinesq model can only be applied for small temperature differences. As shown in the numerical examples, special care should be taken before this model is used in practical problems. For this model, both SUPG and SGS attain very similar numerical properties, yielding almost identical solutions. In this case, using the Boussinesq model with the SUPG method gives faster residual convergence. It is advised to combine the method with the HH or Polner stabilization matrix.
VMS adaptivity has been used to drive the adapted meshes, and it has been shown to yield very smooth meshes, which improve the solution accuracy for the same number of degrees of freedom.

Author Contributions

Conceptualization, G.H. and J.L.; methodology, G.H.; software, G.H. and J.L.; validation, G.H. and J.L.; formal analysis, G.H.; investigation, G.H. and J.L.; writing—original draft preparation, G.H and J.L.; writing—review and editing, G.H. and J.L.; visualization, J.L..; supervision, G.H.; funding acquisition, G.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially funded by the Ministerio de Economia y Competitividad under contract PID2019-106099RB-C44 (AEI/FEDER,UE), Gobierno de Aragon/FEDER-UE (Grupo de Investigacion de Referencia de Mecanica de Fluidos Computacional T32 20R).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We thank T. Lecuyer for performing the wind towers simulations.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Structure of the Stabilization Matrices

Given the stabilization parameters τ c , τ m , and τ e , the incompressible stabilization matrices for 2D can be written as follows [20,29,33].
τ diag = τ c 0 0 0 0 τ m 0 0 0 0 τ m 0 0 0 0 τ e
τ HH = τ c ρ u 1 τ m ρ u 2 τ m 0 0 τ m 0 0 0 0 τ m 0 e ¯ 1 τ e u 1 τ e u 2 τ e τ e
τ polner = τ HH + ω 0 ρ u 1 τ m ρ u 2 τ m 0 u 1 τ m 0 0 0 u 2 τ m 0 0 0 0 0 0 0
τ polner comp = τ c ρ ( ω + 1 ) f u 1 ( ω + 1 ) f u 2 ρ α p | u | 2 τ e ω f u 1 τ m 0 α p u 1 τ e ω f u 2 0 τ m α p u 2 τ e e ¯ 1 τ e ( α p T 1 ) u 1 τ e ( α p T 1 ) u 2 τ e τ e ,
where f = τ m + e ¯ 1 α p τ e and e 1 can be found in [20]. The constant
ω ( 1 2 ( 1 5 ) , 1 2 ( 1 + 5 ) )
is given in [29]. Additionally, from [20,41],
τ m = min Δ t 2 ρ , 1 ρ ( u · G u ) 1 / 2 , 4 12 μ G : G
τ e = min Δ t 2 ρ c v , 1 ρ c v ( u · G u ) 1 / 2 , 4 12 κ G : G
τ c = 1 ρ τ m a · a
G i j = k = 1 2 ξ k x i ξ k x j
a i = j = 1 2 ξ j x i .

Appendix B. Jacobians of the Body Force Term

As a function of the compressibility coefficients expressed in Equations (3) and (4), for pressure-primitive variables, the Jacobian of the body-force term can be written as follows.

Appendix B.1. Full Navier–Stokes Equations

S Y = 0 0 0 0 ρ β T f m 1 0 0 ρ α p f m 1 ρ β T f m 2 0 0 ρ α p f m 2 ρ β T ( f m 1 u 1 + f m 2 u 2 ) ρ f m 1 ρ f m 2 ρ α p ( f m 1 u 1 + f m 2 u 2 )

Appendix B.2. Boussinesq Model

S Y = 0 0 0 0 0 0 0 ρ α p f m 1 0 0 0 ρ α p f m 2 0 ρ α p f m 1 ( T 0 T ) ρ α p f m 2 ( T 0 T ) ρ α p ( f m 1 u 1 + f m 2 u 2 )

References

  1. Baker, A.; Williams, P.; Kelso, R. Development of a robust finite element CFD procedure for predicting indoor room air motion. Build. Environ. 1994, 29, 261–273. [Google Scholar] [CrossRef]
  2. Lube, G.; Knopp, T.; Rapin, G.; Gritzki, R.; Rösler, M. Stabilized finite element methods to predict ventilation efficiency and thermal comfort in buildings. Int. J. Numer. Methods Fluids 2008, 57, 1269–1290. [Google Scholar] [CrossRef]
  3. Hauke, G. An Introduction to Fluid Mechanics and Transport Phenomena; Springer Science+Business Media, B.V.: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  4. Martinez, M.; Gartling, D. A finite element method for low-speed compressible flows. Comput. Methods Appl. Mech. Eng. 2004, 193, 1959–1979. [Google Scholar] [CrossRef]
  5. Gresho, P.; Sutton, S. Application of the FIDAP code to the 8:1 thermal cavity problem. Int. J. Numer. Methods Fluids 2002, 40, 1083–1092. [Google Scholar] [CrossRef]
  6. du Toit, C. A segregated finite element solution for non-isothermal flow. Comput. Methods Appl. Mech. Eng. 2000, 182, 457–481. [Google Scholar] [CrossRef]
  7. Tezduyar, T.E.; Ramakrishnan, S.; Sathe, S. Stabilized formulations for incompressible flows with thermal coupling. Int. J. Numer. Methods Fluids 2008, 57, 1189–1209. [Google Scholar] [CrossRef]
  8. Principe, J.; Codina, R. A stabilized finite element approximation of low speed thermally coupled flows. Int. J. Numer. Methods Heat Fluid Flow 2008, 18, 835–867. [Google Scholar] [CrossRef]
  9. Ryzhakov, P.; Rossi, R.; Oñate, E. An algorithm for the simulation of thermally coupled low speed flow problems. Int. J. Numer. Methods Fluids 2012, 70, 1–19. [Google Scholar] [CrossRef]
  10. Benitez, M.; Bermudez, A. A second order characteristics finite element scheme for natural convection problems. J. Comput. Appl. Math. 2011, 235, 3270–3284. [Google Scholar] [CrossRef]
  11. Gravemeier, V.; Wall, W. Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number. Int. J. Numer. Methods Fluids 2011, 65, 1260–1278. [Google Scholar] [CrossRef]
  12. Yan, J.; Korobenko, A.; Tejada-Martinez, A.; Golshan, R.; Bazilevs, Y. A new variational multiscale formulation for stratified incompressible turbulent flows. Comput. Fluids 2017, 158, 150–156. [Google Scholar] [CrossRef]
  13. Zhu, L.; Masud, A. Residual-based closure model for density-stratified incompressible turbulent flows. Comput. Methods Appl. Mech. Eng. 2021, 386, 113931. [Google Scholar] [CrossRef]
  14. Codina, R.; Principe, J. Dynamic subscales in the finite element approximation of thermally coupled incompressible flows. Int. J. Numer. Methods Fluids 2007, 54, 707–730. [Google Scholar] [CrossRef]
  15. Codina, R.; Principe, J.; Avila, M. Finite element approximation of turbulent thermally coupled incompressible flows with numerical sub-grid scale modelling. Int. J. Numer. Methods Heat Fluid Flow 2010, 20, 492–516. [Google Scholar] [CrossRef]
  16. Xu, S.; Gao, B.; Hsu, M.C.; Ganapathysubramanian, B. A residual-based variational multiscale method with weak imposition of boundary conditions for buoyancy-driven flows. Comput. Methods Appl. Mech. Eng. 2019, 352, 345–368. [Google Scholar] [CrossRef]
  17. Principe, J.; Codina, R. Mathematical models for thermally coupled low speed flows. Adv. Theor. Appl. Mech. 2009, 2, 93–112. [Google Scholar]
  18. Mittal, S.; Tezduyar, T. A unified finite element formulation for compressible and incompressible flows using augmented conservation variables. Comput. Methods Appl. Mech. Eng. 1998, 161, 229–243. [Google Scholar] [CrossRef]
  19. Hauke, G.; Hughes, T. A unified approach to compressible and incompressible flows. Comput. Methods Appl. Mech. Eng. 1994, 113, 389–395. [Google Scholar] [CrossRef]
  20. Hauke, G.; Hughes, T. A comparative study of different sets of variables for solving compressible and incompressible flows. Comput. Methods Appl. Mech. Eng. 1998, 153, 1–44. [Google Scholar] [CrossRef]
  21. Givoli, D.; Flaherty, J.E.; Shephard, M.S. Parallel adaptive finite element analysis of viscous flows based on a combined compressible-incompressible formulation. Int. J. Numer. Methods Heat Fluid Flow 1997, 7, 880–906. [Google Scholar] [CrossRef]
  22. Chalot, F.; Hughes, T.; Shakib, F. Symmetrization of conservation laws with entropy for high-temperature hypersonic computations. Comput. Syst. Eng. 1990, 1, 459–521. [Google Scholar] [CrossRef]
  23. Hauke, G.; Landaberea, A.; Garmendia, I.; Canales, J. On the thermodynamics, stability and hierarchy of entropy functions in fluid flow. Comput. Methods Appl. Mech. Eng. 2006, 195, 4473–4489. [Google Scholar] [CrossRef]
  24. Pons, M.; Quéré, P.L. An example of entropy balance in natural convection, Part 1: The usual Boussinesq equations. Comptes Rendus Méc. 2005, 333, 127–132. [Google Scholar] [CrossRef]
  25. Pons, M.; Quéré, P.L. An example of entropy balance in natural convection, Part 2: The thermodynamic Boussinesq equations. Comptes Rendus Méc. 2005, 333, 133–138. [Google Scholar] [CrossRef]
  26. Brooks, A.; Hughes, T. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 1982, 32, 199–259. [Google Scholar] [CrossRef]
  27. Shakib, F.; Hughes, T.; Johan, Z. A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 1991, 89, 141–219. [Google Scholar] [CrossRef]
  28. Hughes, T.; Franca, L.; Hulbert, G. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput. Methods Appl. Mech. Eng. 1989, 73, 173–189. [Google Scholar] [CrossRef]
  29. Polner, M.; van der Vegt, J.; van Damme, R. Analysis of stabilization operators for Galerkin least-squares discretizations of the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 2006, 196, 982–1006. [Google Scholar] [CrossRef]
  30. Hauke, G. Simple Stabilizing Matrices for the Computation of Compressible Flows in Primitive Variables. Comput. Methods Appl. Mech. Eng. 2001, 190, 6881–6893. [Google Scholar] [CrossRef]
  31. Hauke, G.; Landaberea, A.; Garmendia, I.; Canales, J. A Segregated Method for Compressible Flow Computation. Part I.: Isothermal Compressible Flows. Int. J. Numer. Methods Fluids 2005, 47, 271–323. [Google Scholar] [CrossRef]
  32. Hauke, G.; Landaberea, A.; Garmendia, I.; Canales, J. A Segregated Method for Compressible Flow Computation. Part II.: General Divariant Compressible Flows. Int. J. Numer. Methods Fluids 2005, 49, 183–209. [Google Scholar] [CrossRef]
  33. Polner, M.; Pesch, L.; van der Vegt, J. Construction of stabilization operators for Galerkin least-squares discretizations of compressible and incompressible flows. Comput. Methods Appl. Mech. Eng. 2007, 196, 2431–2448. [Google Scholar] [CrossRef]
  34. Costa, G.K.; Lyra, P.R.M. A stabilized finite element formulation to solve high and low speed flows. Commun. Numer. Methods Eng. 2006, 22, 411–419. [Google Scholar] [CrossRef]
  35. Zhu, L.; Goraya, S.; Masud, A. A Variational Multiscale Method for Natural Convection of Nanofluids. Mech. Res. Commun. 2022. [Google Scholar]
  36. Taylor, C.; Hughes, T.; Zarins, C. Finite element modeling of blood flow in arteries. Comput. Methods Appl. Mech. Eng. 1998, 158, 155–196. [Google Scholar] [CrossRef]
  37. Whiting, C.H.; Jansen, K.E. A stabilized finite element method for the incompressible Navier–Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 2001, 35, 93–116. [Google Scholar] [CrossRef]
  38. Codina, R. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Eng. 2000, 190, 1579–1599. [Google Scholar] [CrossRef]
  39. Codina, R. A stabilized finite element method for generalized incompressible flows. Comput. Methods Appl. Mech. Eng. 2001, 190, 2681–2706. [Google Scholar] [CrossRef]
  40. Codina, R. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Comput. Methods Appl. Mech. Eng. 2002, 191, 4295–4321. [Google Scholar] [CrossRef]
  41. Bazilevs, Y.; Calo, V.; Cottrell, J.; Hughes, T.; Reali, A.; Scovazzi, G. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Eng. 2007, 197, 173–201. [Google Scholar] [CrossRef]
  42. Bayona, C.; Baiges, J.; Codina, R. Variational multi-scale finite element approximation of the compressible Navier-Stokes equations. Int. J. Numer. Methods Heat Fluid Flow 2016, 26, 1240–1271. [Google Scholar] [CrossRef]
  43. Xu, F.; Moutsanidis, G.; Kamensky, D.; Hsu, M.C.; Murugan, M.; Ghoshal, A.; Bazilevs, Y. Compressible flows on moving domains: Stabilized methods, weakly enforced essential boundary conditions, sliding interfaces, and application to gas-turbine modeling. Comput. Fluids 2017, 158, 201–220. [Google Scholar] [CrossRef]
  44. Hughes, T.; Scovazzi, G.; Tezduyar, T. Stabilized methods for compressible flows. J. Sci. Comput. 2010, 43, 343–368. [Google Scholar] [CrossRef]
  45. Franca, L.; Frey, S. Stabilized finite element methods: II. The incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 1992, 99, 209–233. [Google Scholar] [CrossRef]
  46. de Vahl Davis, G. Natural Convection of Air in a Square Cavity: A Bench Mark Numerical Solution. Int. J. Numer. Methods Fluids 1983, 3, 249–264. [Google Scholar] [CrossRef]
  47. LeQuéré, P.; Weisman, C.; Paillère, H.; JanVierendeels.; Dick, E.; Becker, R.; Braack, M.; Locke, J. Modelling of Natural Convection Flows with LargeTemperature Differences: A Benchmark Problem for Low MachNumber Solvers. Part 1. Reference Solutions. ESAIM Math. Model. Numer. Anal. 2005, 39, 609–616. [Google Scholar] [CrossRef]
  48. Mamimid, S.; Guellal, M.; Bouafia, M. Numerical study of natural convection in a square cavity under non-Boussinesq conditions. Therm. Sci. 2016, 20, 1509–1517. [Google Scholar]
  49. Becker, R.; Braak, M. Solution of a stationary benchmark problem for natural convection with large temperature difference. Int. J. Therm. Sci. 2002, 41, 428–439. [Google Scholar] [CrossRef]
  50. Heuveline, V. On higher-order mixed FEM for low Mach number flows: Application to a natural convection benchmark problem. Int. J. Numer. Methods Fluids 2003, 41. [Google Scholar] [CrossRef]
  51. de Vahl Davis, G.; Jones, I. Natural Convection in a square Cavity: A comparison exercise. Int. J. Numer. Methods Fluids 1983, 3, 227–248. [Google Scholar] [CrossRef]
  52. Vierendeels, J.; Merci, B.; Dick, E. Numerical study of natural convective heat transfer with large horizontal temperature differences. Int. J. Numer. Methods Heat Fluid Flow 2001, 11, 329–341. [Google Scholar] [CrossRef]
  53. Vierendeels, J.; Merci, B.; Dick, E. Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences. Int. J. Numer. Methods Heat Fluid Flow 2003, 13, 1057–1078. [Google Scholar] [CrossRef]
  54. Paillère, H.; Le LeQuéré, P.; Weisman, C.; Vierendeels, J.; Dick, E.; Braack, M.; Dabbene, F.; Beccantini, A.; Studer, E.; Kloczko, T.; et al. Modelling of Natural Convection Flows with LargeTemperature Differences: A Benchmark Problem for Low MachNumber Solvers. Part 2. Contributions to the 2004 June Conference. ESAIM Math. Model. Numer. Anal. 2005, 39, 617–621. [Google Scholar] [CrossRef]
  55. Wittschieber, S.; Demkowicz, L.; Behr, M. Stabilized finite element methods for a fully-implicit logarithmic reformulation of the Oldroyd-B constitutive law. J. Non-Newton. Fluid Mech. 2022, 306, 104838. [Google Scholar] [CrossRef]
  56. Quéré, P.L. Accurate solutions to the square thermally driven cavity at high Rayleigh number. Comput. Fluids 1991, 20, 29–41. [Google Scholar] [CrossRef]
  57. Hauke, G.; Doweidar, M.H.; Miana, M. Proper intrinsic scales for a-posteriori multiscale error estimation. Comput. Methods Appl. Mech. Eng. 2006, 195, 3983–4001. [Google Scholar] [CrossRef]
  58. Irisarri, D.; Hauke, G. A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method. Comput. Methods Appl. Mech. Eng. 2016, 311, 648–670. [Google Scholar] [CrossRef]
  59. Irisarri, D.; Hauke, G. A posteriori error estimation and adaptivity based on VMS for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 2021, 373, 113508. [Google Scholar] [CrossRef]
  60. Hauke, G.; Fuster, D.; Lizarraga, F. Variational multiscale a posteriori error estimation for systems: The Euler and Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. 2015, 283, 1493–1524. [Google Scholar] [CrossRef]
  61. Hauke, G.; Doweidar, M.H.; Fuentes, S. Mesh adaptivity for the transport equation led by variational multiscale error estimators. Int. J. Numer. Methods Fluids 2011, 69, 1835–1850. [Google Scholar] [CrossRef]
  62. Ortego, I. Torres de Viento. Relectura Contemporanea de Sistemas Energetivos Pasivos. Master’s Thesis, ETS Arquitectura de Madrid, Madrid, Spain, 2020. [Google Scholar]
Figure 1. Buoyancy-driven square cavity. Problem setup.
Figure 1. Buoyancy-driven square cavity. Problem setup.
Algorithms 15 00278 g001
Figure 2. Buoyancy-driven square cavity. Full Navier–Stokes equations. Residual convergence for various stabilization matrices.
Figure 2. Buoyancy-driven square cavity. Full Navier–Stokes equations. Residual convergence for various stabilization matrices.
Algorithms 15 00278 g002
Figure 3. Buoyancy-driven square cavity. Boussinesq model with SUPG. Residual convergence for various stabilization matrices.
Figure 3. Buoyancy-driven square cavity. Boussinesq model with SUPG. Residual convergence for various stabilization matrices.
Algorithms 15 00278 g003
Figure 4. Buoyancy-driven square cavity. Boussinesq model with SGS. Residual convergence for various stabilization matrices.
Figure 4. Buoyancy-driven square cavity. Boussinesq model with SGS. Residual convergence for various stabilization matrices.
Algorithms 15 00278 g004
Figure 5. Buoyancy-driven square cavity. Residual convergence for GMRES tolerances for the Boussinesq model.
Figure 5. Buoyancy-driven square cavity. Residual convergence for GMRES tolerances for the Boussinesq model.
Algorithms 15 00278 g005
Figure 6. Buoyancy-driven square cavity. Comparative study of residual convergence for various methods.
Figure 6. Buoyancy-driven square cavity. Comparative study of residual convergence for various methods.
Algorithms 15 00278 g006
Figure 7. Buoyancy-driven square cavity. Full Navier–Stokes with the HH stabilization matrix. Comparison of Jacobian.
Figure 7. Buoyancy-driven square cavity. Full Navier–Stokes with the HH stabilization matrix. Comparison of Jacobian.
Algorithms 15 00278 g007
Figure 8. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.01 and HH tau.
Figure 8. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.01 and HH tau.
Algorithms 15 00278 g008
Figure 9. Buoyancy-driven square cavity. Initial mesh and adapted meshes for ϵ = 0.01 .
Figure 9. Buoyancy-driven square cavity. Initial mesh and adapted meshes for ϵ = 0.01 .
Algorithms 15 00278 g009
Figure 10. Buoyancy-driven square cavity. Adapted solutions for ϵ = 0.01 obtained with various methods and the HH tau.
Figure 10. Buoyancy-driven square cavity. Adapted solutions for ϵ = 0.01 obtained with various methods and the HH tau.
Algorithms 15 00278 g010
Figure 11. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.6 and Sutherland viscosity.
Figure 11. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.6 and Sutherland viscosity.
Algorithms 15 00278 g011
Figure 12. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.6 and Sutherland viscosity for adapted meshes.
Figure 12. Buoyancy-driven square cavity. Nusselt numbers along hot vertical wall for ϵ = 0.6 and Sutherland viscosity for adapted meshes.
Algorithms 15 00278 g012
Figure 13. Example of Iranian wind tower.
Figure 13. Example of Iranian wind tower.
Algorithms 15 00278 g013
Figure 14. 2D model and working principle of a wind tower.
Figure 14. 2D model and working principle of a wind tower.
Algorithms 15 00278 g014
Figure 15. Wind tower problem setup.
Figure 15. Wind tower problem setup.
Algorithms 15 00278 g015
Figure 16. Wind tower 1. Zoom of initial and adapted mesh around the tower.
Figure 16. Wind tower 1. Zoom of initial and adapted mesh around the tower.
Algorithms 15 00278 g016
Figure 17. Wind tower 1. Velocity modulus obtained using the full NS and the Boussinesq model and U = 5 m/s.
Figure 17. Wind tower 1. Velocity modulus obtained using the full NS and the Boussinesq model and U = 5 m/s.
Algorithms 15 00278 g017
Figure 18. Wind tower 1. Temperature field obtained using the full NS and the Boussinesq model and U = 5 m/s.
Figure 18. Wind tower 1. Temperature field obtained using the full NS and the Boussinesq model and U = 5 m/s.
Algorithms 15 00278 g018
Figure 19. Wind tower 1. Apparent temperature field obtained by the full NS and the Boussinesq model with U = 5 m/s.
Figure 19. Wind tower 1. Apparent temperature field obtained by the full NS and the Boussinesq model with U = 5 m/s.
Algorithms 15 00278 g019
Figure 20. Wind tower 1. Comparative apparent temperature field obtained with the full NS with U = 1 m/s and U = 5 m/s.
Figure 20. Wind tower 1. Comparative apparent temperature field obtained with the full NS with U = 1 m/s and U = 5 m/s.
Algorithms 15 00278 g020
Figure 21. Wind tower 2. Zoom around dome geometry.
Figure 21. Wind tower 2. Zoom around dome geometry.
Algorithms 15 00278 g021
Figure 22. Wind tower 2. VMS T-adapted mesh.
Figure 22. Wind tower 2. VMS T-adapted mesh.
Algorithms 15 00278 g022
Figure 23. Wind tower 2. Velocity modulus obtained with the full NS and the Boussinesq model and U = 5 m/s.
Figure 23. Wind tower 2. Velocity modulus obtained with the full NS and the Boussinesq model and U = 5 m/s.
Algorithms 15 00278 g023
Figure 24. Wind tower 2. Apparent temperature field obtained by the full NS and the Boussinesq model with U = 5 m/s.
Figure 24. Wind tower 2. Apparent temperature field obtained by the full NS and the Boussinesq model with U = 5 m/s.
Algorithms 15 00278 g024
Table 1. Relationship between differential operator and type of stabilized method.
Table 1. Relationship between differential operator and type of stabilized method.
SUPGGLSSGS
L 𝓛 SUPG T 𝓛 T 𝓛 *
Table 2. Nusselt numbers along hot vertical wall. Ra = 10 6 , ϵ = 0.01 , constant μ . Present methods on a 140 × 140 uniform mesh with various taus.
Table 2. Nusselt numbers along hot vertical wall. Ra = 10 6 , ϵ = 0.01 , constant μ . Present methods on a 140 × 140 uniform mesh with various taus.
Reference Solution Nu ¯ Nu min Nu max
vahl Davis [46] 8.8170.98917.925
Le Quéré [56] 8.82520.9794617.5360
Vierendeels [52] 8.8257
Masud [35] 8.81490
Present Methodtau Nu ¯ Nu min Nu max
Full NSClassic8.08330.80313.932
Full NSHH8.85240.98417.704
Full NSPolner8.85260.98417.703
Boussinesq SUPGHH8.85400.98217.689
Boussinesq SUPGPolner8.85420.98217.688
Boussinesq SGSHH8.86320.98517.704
Boussinesq SGSPolner8.86340.98517.703
Table 3. Nusselt numbers along hot vertical wall. Ra = 10 6 , ϵ = 0.6 , Sutherland viscosity. Present methods on a 140 × 140 uniform mesh.
Table 3. Nusselt numbers along hot vertical wall. Ra = 10 6 , ϵ = 0.6 , Sutherland viscosity. Present methods on a 140 × 140 uniform mesh.
Reference Solution Nu ¯ Nu min Nu max
Vierendeels et al. [53,54] 8.68661.066720.2704
Heuveline [50,54] 8.68611.067420.3051
Becker-Braak [49] 8.6866
Present methodtau Nu ¯ Nu min Nu max
Full NSClassic8.62271.02717.940
Full NSHH9.11341.09321.629
Full NSPolner9.11421.09321.632
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hauke, G.; Lanzarote, J. Simulation of Low-Speed Buoyant Flows with a Stabilized Compressible/Incompressible Formulation: The Full Navier–Stokes Approach versus the Boussinesq Model. Algorithms 2022, 15, 278. https://doi.org/10.3390/a15080278

AMA Style

Hauke G, Lanzarote J. Simulation of Low-Speed Buoyant Flows with a Stabilized Compressible/Incompressible Formulation: The Full Navier–Stokes Approach versus the Boussinesq Model. Algorithms. 2022; 15(8):278. https://doi.org/10.3390/a15080278

Chicago/Turabian Style

Hauke, Guillermo, and Jorge Lanzarote. 2022. "Simulation of Low-Speed Buoyant Flows with a Stabilized Compressible/Incompressible Formulation: The Full Navier–Stokes Approach versus the Boussinesq Model" Algorithms 15, no. 8: 278. https://doi.org/10.3390/a15080278

APA Style

Hauke, G., & Lanzarote, J. (2022). Simulation of Low-Speed Buoyant Flows with a Stabilized Compressible/Incompressible Formulation: The Full Navier–Stokes Approach versus the Boussinesq Model. Algorithms, 15(8), 278. https://doi.org/10.3390/a15080278

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop