DynamFluid: Development and Validation of a New GUI-Based CFD Tool for the Analysis of Incompressible Non-Isothermal Flows

Abstract

A computational fluid dynamics software (DynamFluid) based on the application of the finite element method with the characteristic-based-split algorithm is presented and validated. The software is used to numerically integrate the steady and unsteady Navier–Stokes equations for both constant-density and Boussinesq non-isothermal flows. Benchmark two-dimensional computations carried out with DynamFluid show good agreement with previous results reported in the literature. Test cases used for validation include (i) the lid-driven cavity flow, (ii) mixed convection flow in a vertical channel with asymmetric wall temperatures, (iii) unsteady incompressible flow past a circular cylinder, and (iv) steady non-isothermal flow past a circular cylinder with negligible buoyancy effects. The new software is equipped with a graphical user interface that facilitates the definition of the fluid properties, the discretization of the physical domain, the definition of the boundary conditions, and the post-processing of the computed velocity, pressure and temperature fields.

1. Introduction

The science of fluid mechanics involves a broad spectrum of techniques for the study of fluid flows [1,2]. Experimental fluid mechanics plays an important role in the validation of theoretical models and the determination of their limits of application. Moreover, laboratory-scale experiments (e.g., wind tunnel model testing) combined with dimensional analysis [3,4] provides the science with an effective way of studying complex flows, saving the time and money that would otherwise be employed in full-scale experimental studies. On the other hand, computational fluid dynamics (CFD) is the art of substituting the set of partial differential equations governing the flow by a set of algebraic equations that can be solved with computers [5,6,7,8,9,10]. As a result, CFD can be regarded as the connection between theory and experiments, providing an economical alternative to experimental model testing. The ever growing importance of CFD has driven the development of a large variety of commercial software (ANSYS Fluent [11], Star-CCM+ [12], COMSOL’s CFD [13], Altair AcuSolve [14], to name a few) and open source solutions (OpenFOAM [15,16,17], FreeFEM++ [18], FEniCS [19], for example) in the last decades.

This paper presents a new CFD software tool (DynamFluid) based on the finite element method (FEM) for the analysis of incompressible flows. The finite element method has been selected in favour of other options, such as the finite difference method (FDM) or the finite volume method (FVM), because of its higher adaptability to arbitrarily shaped domains, as other methods require complex curvilinear transformation when trying to solve problems in complex geometries represented in regular Cartesian coordinates. This and other advantages are already exploited in a variety of commercial and open-source codes, such as the above mentioned COMSOL’s CFD [13], FreeFEM++ [18] or FEniCS [19]. Although the developed code has the capability of solving fully three-dimensional (3D) problems with arbitrary density variations and the general form of the gravitational term (body force), the validation process in this work will be restricted to two-dimensional (2D) flows where the fluid density undergoes tiny variations, so that the Boussinesq approximation can be employed. Further validation in 3D geometries with large density variations is left for future work. DynamFluid is implemented in C++.

The main characteristics of the DynamFluid tool can be summarized as follows:

• Windows-based software, which can be run in any Windows Operating System: both 32- and 64-bit architectures are supported. When run in a 64-bit architecture, the software leverages on the larger word-size when performing computations.
• Graphical user interface for defining the geometric domain, physical model, boundary conditions, and displaying the results obtained in any simulation.
• Custom database for storing both the project domain definition and all the information generated during the simulation. The database provides the user with ODBC (Open DataBase Connectivity) interface for interacting with any ODBC client.
• Support for NASTRAN format file importing, which allows the user to import any geometry and physical definition in DMAP (Direct Matrix Abstraction Programming) language.
• Vtk file format capabilities, to export generated simulations into ASCII text Vtk files that can be visualized using Paraview [20].
• Basic meshing capabilities to sample the geometric domain. Two methods have been implemented: (a) structured meshing (linear, logarithmic) using both quadrangular elements (QUAD) and triangular elements (TRIA) for regular geometric domains, and (b) Delaunay–Voronoi meshing for irregular geometric domains.
• Software designed to run in a computer with a motherboard that may have one or several multi-core processors. This includes parallel computation of the finite element matrices, parallel assembly of the global matrices and parallel computation of the right hand side of each step of the algorithm. The software uses the conjugate gradient stabilized algorithm provided by the Eigen library [21] to solve the linear systems, and it has been compiled with the openmp compiler flag so that the Eigen library exploits the multiple cores available in the hardware.
• Custom user language for post-processing the results (velocity, pressure and temperature), with basic algebra functions and support for different coordinate reference systems. Internal compiler for translating this user language into machine code that can be applied in parallel to every node and/or finite element.
• Support for different types of finite elements (both Lagrangian and Serendipity): (a) linear TRIA elements and (b) linear and quadratic QUAD elements.
• DynamFluid is a freeware CFD tool available at https://sites.google.com/view/dynamfluid.

Four benchmark 2D problems have been selected in this work for validation purposes, namely (i) the lid-driven cavity flow, (ii) mixed convection flow in a vertical channel with asymmetric wall temperatures, (iii) unsteady incompressible flow past a circular cylinder, and (iv) steady non-isothermal flow past a circular cylinder with negligible buoyancy effects. These problems are briefly outlined below.

(i)

The lid-driven cavity flow.This is a classical benchmark problem that has been widely used since the early days of CFD to assess and validate new techniques and methods. This test case is easy to set and simulate because its boundary conditions are particularly simple. However, the fully developed flow displays almost all fluid mechanical phenomena, with increasingly complex aspects emerging as the Reynolds number is increased, such as corner eddies, laminar to turbulence regime transition, and even turbulence at high Reynolds number. A recent comprehensive review of the literature on the subject can be found in [22], where the work of several authors is presented and discussed [23,24,25,26,27,28]. Available benchmark results have been tabulated to provide a comprehensive source of validation data [29].

(ii)

Mixed convection flow in a vertical channel with asymmetric wall temperatures. In this mixed convection heat transfer problem, an initially uniform flow develops in a slender vertical channel whose walls are at different temperatures. The cold and hot wall temperatures may also differ from the incoming flow temperature. As a result of the upward buoyancy force that appears near the hot wall, the velocity increases in the near-wall region. As the fluid accelerates downstream, the fluid near the cold wall may suffer flow reversal so as to maintain the imposed fixed flow rate. One of the most interesting features exhibited by this flow is thus the possibility of flow reversal at the cold wall as the flow develops. The occurrence (or not) of flow reversal depends on the length of the vertical channel and the buoyancy effect induced by the temperature difference between the hot and cold walls. Previous studies have shown that for the flow reversal to occur, in high Reynolds number flows the ratio of the Grashof number to the Reynolds number must be higher than a critical value that depends on the wall temperature difference ratio [30,31,32,33,34,35,36,37,38,39,40,41,42,43].

(iii)

Unsteady incompressible flow past a circular cylinder. The flow of a constant density fluid past a bluff body is another classical problem that has been widely studied in the literature. Understanding the flow regimes past bluff bodies poses a daunting challenge, so that 2D and 3D vortical structures in wakes of different bodies have been analyzed by scientists and engineers for decades. The reason for this interest is the vast range of applications of external flow past round bluff bodies: aerodynamics (planes, rockets, ground vehicles), hydrodynamics (ships, submarines) or wind energy (wind turbines), to name the few. At very low Reynolds numbers (creeping flow) the flow past a non-heated circular cylinder is symmetric in the streamwise direction [44,45,46,47]. As the Reynolds number grows to values of order unity the symmetry is lost, and when it exceeds a critical value the non-linear convective effects trigger the onset of steady flow separation, accompanied by the appearance of steady recirculation bubbles behind the cylinder. One important aspect of these flows is the so-called vortex shedding, which is an oscillating flow pattern that emerges for even larger Reynolds numbers. This regime has been thoroughly studied experimentally [48,49,50,51,52] and numerically [53,54,55,56,57]. The alternate shedding of vortices in the wake leads to the well known Kármán vortex street, which originates large fluctuating pressure forces in the direction transverse to the flow and may cause structural vibrations, acoustic noise, or resonance, which in some cases may lead to structural damage or even collapse.

(iv)

Steady non-isothermal flow past a circular cylinder with negligible buoyancy effects. This case is similar to the previous one but with the particularity that the temperature of the cylinder differs from the temperature of the incoming fluid, which causes the flow past a circular cylinder to exhibit interesting heat transfer features. This problem is of interest for the design of cylinder-shaped sensors located in fluid streams, hot-wire anemometers, tube heat exchangers, nuclear reactor fuel rods and chimneys. In this work, attention will be restricted to steady flow with negligible buoyancy effects, with the aim of characterizing the local Nusselt number at the cylinder wall for different Reynolds numbers [44,54,58,59].

2. Governing Equations

Let $\mathbf{x}={(x,y,z)}^T\in {\Omega }_t\subset R^3$ be the spatial domain at time $t\in (0,T)$, with $x_i$ denoting the i-th cartesian coordinate. Using the Einstein summation convention, the governing equations for three-dimensional unsteady flow with variable density (due, e.g., to temperature variations) can be written as follows

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_j\right)}{\partial x_j}& =0\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \left[\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_j{u}_i\right)}{\partial x_j}\right]& =-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+\rho g_i,i=\{1,2,3\}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \left[\frac{\partial \left(\rho e_{\mathrm{T}}\right)}{\partial t}+\frac{\partial (\rho u_j{e}_{\mathrm{T})}}{\partial x_j}\right]& =-\frac{\partial \left(p{u}_j\right)}{\partial x_j}+\frac{\partial \left({\tau }_{ij}{u}_i\right)}{\partial x_j}+\rho g_j{u}_j+\frac{\partial }{\partial x_j}\left(k\frac{\partial T}{\partial x_j}\right),\hfill \end{array}$$

(3)

where $\rho$ is the fluid density, $u_i$ the i-th component of the velocity vector $\mathbf{u}={(u,v,w)}^T$, $e_{\mathrm{T}}=e+u_i{u}_i/2$ the total energy per unit mass (with $e=c_{v}T$ the specific internal energy per unit mass, assuming a calorically perfect fluid), p the pressure, T the absolute temperature, $g_i$ the i-th component of the acceleration of gravity, and k the thermal conductivity. In the above equations ${\tau }_{ij}$ represent the deviatoric stress components

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_k}{\partial x_k}\right),$$

(4)

where $\mu$ is the dynamic viscosity and ${\delta }_{ij}$ the Kronecker delta.

In non-isothermal flows where the density variations are small but the flow is driven by buoyancy forces, such as in natural or mixed convection problems, one may simplify the equations by means of the Boussinesq approximation. This approximation assumes that variations in density have no effect on the flow field other than to give rise to buoyancy forces. In this case, the continuity and momentum equations take the simplified form

$$\begin{array}{cc}\hfill \frac{\partial u_j}{\partial x_j}& =0\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\rho }_{\infty }\left[\frac{\partial u_i}{\partial t}+\frac{\partial \left(u_j{u}_i\right)}{\partial x_j}\right]& =-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+{\rho }_{\infty }(1-\beta (T-T_{\infty }))g_i,i=\{1,2,3\}\hfill \end{array}$$

(6)

where ${\rho }_{\infty }$ is the reference fluid density, $T_{\infty }$ the reference temperature, and $\beta =-{\rho }^{-1}{(\partial \rho /\partial T)}_p$ the thermal expansion coefficient.

In the Boussinesq approximation, the energy equation is often simplified by also ignoring the density variations. However, in this work, the energy equation will be kept in its general form, without any further simplification. The reason for this approach is to minimize the approximation error when computing the temperature field in a fluid that undergoes tiny variations of density, and to give a unified treatment to the energy equation independently of the way density is computed. In the range of Reynolds number considered here, the different treatments of the energy equation, and the assumption of constant specific heats, do not make a significant difference. Leaving the general form allows validating the implementation of the energy equation when the Boussinesq approximation is used in the Navier–Stokes equations but not in the energy equation.

Let $L_{\infty }$ be the characteristic length, $U_{\infty }$ the characteristic velocity, $L_{\infty }/U_{\infty }$ the characteristic residence time of the problem, ${\rho }_{\infty }$ the reference density, ${\mu }_{\infty }$ the reference dynamic viscosity, and $k_{\infty }$ the reference thermal conductivity. With these scales, the following dimensionless variables can be introduced

$$\begin{array}{cc}{x}_i^{*}={\displaystyle \frac{x_i}{L_{\infty }}},t^{*}={\displaystyle \frac{U_{\infty }t}{L_{\infty }}},\hfill & {\rho }^{*}={\displaystyle \frac{\rho }{{\rho }_{\infty }}},u_i^{*}={\displaystyle \frac{u_i}{U_{\infty }}},p^{*}={\displaystyle \frac{p-p_{\infty }}{{\rho }_{\infty }{U}_{\infty }^2}},{\tau }_{ij}^{*}={\displaystyle \frac{{\tau }_{ij}{L}_{\infty }}{{\mu }_{\infty }{U}_{\infty }}},\hfill \\ \hfill & g_i^{*}={\displaystyle \frac{g_i}{g}},{\mu }^{*}={\displaystyle \frac{\mu }{{\mu }_{\infty }}},e_{\mathrm{T}}^{*}={\displaystyle \frac{e_{\mathrm{T}}}{U_{\infty }^2}},T^{*}={\displaystyle \frac{c_{p}T}{U_{\infty }^2}},k^{*}={\displaystyle \frac{k}{k_{\infty }}}\hfill \end{array}$$

(7)

where $p_{\infty }$ represents a convenient pressure datum, and $g\simeq 9.81\mathrm{m}/{\mathrm{s}}^2$ is the acceleration of gravity.

The asterisk will be dropped in the following for simplicity, so that all variables will be assumed to be non-dimensional. Thus, the general governing Equations (1)–(3) can be rewritten as follows

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_j\right)}{\partial x_j}& =0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \left[\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_j{u}_i\right)}{\partial x_j}\right]& =-\frac{\partial p}{\partial x_i}+\frac{1}{\mathrm{Re}}\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{1}{\mathrm{Fr}}\rho g_i,i=\{1,2,3\}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \left[\frac{\partial \left(\rho e_{\mathrm{T}}\right)}{\partial t}+\frac{\partial \left(\rho u_j{e}_{\mathrm{T}}\right)}{\partial x_j}\right]& =-\frac{\partial \left(p{u}_j\right)}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{\partial \left({\tau }_{ij}{u}_j\right)}{\partial x_i}+\frac{1}{\mathrm{Fr}}\rho g_j{u}_j+\frac{1}{\mathrm{Re}\mathrm{Pr}}\frac{\partial }{\partial x_j}\left(k\frac{\partial T}{\partial x_j}\right),\hfill \end{array}$$

(10)

where $\mathrm{Re}={\rho }_{\infty }{U}_{\infty }{L}_{\infty }/{\mu }_{\infty }$ is the Reynolds number, $\mathrm{Fr}=U_{\infty }^2/\left(L_{\infty }g\right)$ is the Froude number, and $\mathrm{Pr}={\mu }_{\infty }{c}_p/k_{\infty }$ is the Prandtl number, defined in terms of the specific heat at constant pressure $c_p$, assumed here to be constant. Note that for liquids both specific heats are equal $c_p=c_v$, whereas for gases $c_p=c_v+R_g$, where $R_g=R^{\circ }/W$ is the gas constant, defined as the ratio between the universal gas constant $R^{\circ }=8.314\mathrm{J}/(\mathrm{mol}·\mathrm{K})$ and the molecular mass of the gas W that is constant for gases of uniform composition.

Under the Boussinesq approximation, it is convenient to define the normalized dimensionless temperature

$$\mathsf{\Theta }=\frac{T-T_{\infty }}{T_w-T_{\infty }},$$

(11)

in terms of the reference temperature $T_{\infty }$ and a characteristic temperature $T_w$ (e.g., that of a hot or cold wall, thereby the subscript w) that determines the characteristic temperature difference of the problem, $T_w-T_{\infty }$. In this case, the momentum conservation Equation (6) and the energy conservation Equation (3) can be rewritten as follows, where the asterisks are dropped again so that all variables are non-dimensional

$$\begin{array}{cc}\hfill \frac{\partial u_i}{\partial t}+\frac{\partial \left(u_j{u}_i\right)}{\partial x_j}& =-\frac{\partial p}{\partial x_i}+\frac{1}{\mathrm{Re}}\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{1}{\mathrm{Fr}}{g}_i-\frac{\mathrm{Gr}}{{\mathrm{Re}}^2}{g}_i\mathsf{\Theta },i=\{1,2,3\}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \rho \left[\frac{\partial e_{\mathrm{T}}}{\partial t}+\frac{\partial \left(u_j{e}_{\mathrm{T}}\right)}{\partial x_j}\right]& =-\frac{\partial \left(p{u}_j\right)}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{\partial \left({\tau }_{ij}{u}_j\right)}{\partial x_i}+\frac{1}{\mathrm{Fr}}\rho g_j{u}_j+\frac{1}{\mathrm{Re}\mathrm{Pr}}\frac{\partial }{\partial x_j}\left(k\frac{\partial T}{\partial x_j}\right),\hfill \end{array}$$

(13)

while the continuity Equation (5) remains unchanged, implying the divergence-free nature of the velocity field. The simplified form of the continuity equation has been used, in particular, to rewrite the convective term in the energy Equation (13) including the velocity $u_j$ into the spatial derivative. Note that, following standard practice, in the momentum conservation Equation (12) the dimensionless density has been assumed to be constant ($\rho =1$) except in the buoyancy term, which is written in terms of the Grashof number

$$\mathrm{Gr}=\frac{\beta g(T_w-T_{\infty })L_{\infty }^3}{{\nu }_{\infty }^2},$$

(14)

and the normalized dimensionless temperature $\mathsf{\Theta }$, given by (11) in terms of the dimensionless temperature T computed from (13).

Finally, to be well posed, the mathematical problem also requires appropriate boundary conditions and a set of initial conditions for all the variables to be solved.

3. Numerical Method

3.1. Temporal Discretization: The Characteristics-Based-Split Algorithm

Following Zienkiewicz’s et al. [60,61,62,63,64], the Navier-Stokes equations can be sampled in time using a characteristic method, the so-called characteristic-based-split (CBS) scheme. To obtain the numerical solution, the time interval will be divided into $N_t$ subintervals $I_n:=(t_n,t_{n+1}]$ with constant time step size $\mathsf{\Delta }t=t_{n+1}-t_n$. Then, a variable $\varphi$ at time $n+\theta$ can be approximated as ${\varphi }^{n+\theta }=(1-\theta ){\varphi }^n+\theta {\varphi }^{n+1}$, where the superscript denotes the time at which the variable is evaluated. As a result, it can be written that

$${\varphi }^{n+\theta }=(1-\theta ){\varphi }^n+\theta {\varphi }^{n+1}={\varphi }^n+\theta \mathsf{\Delta }\varphi ,$$

(15)

with $\mathsf{\Delta }\varphi ={\varphi }^{n+1}-{\varphi }^n$.

Consequently, the CBS scheme, proposed by Zienkiewicz and Codina [61], applied to the governing Equations (1), (12) and (13) can be sampled in time as follows

$$\begin{array}{cc}\hfill \frac{u_i^{n+1}-u_i^n}{\mathsf{\Delta }t}=\frac{\mathsf{\Delta }{u}_i}{\mathsf{\Delta }t}& ={\left[-\frac{\partial \left(u_j{u}_i\right)}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{1}{\mathrm{Fr}}{g}_i-\frac{\mathrm{Gr}}{{\mathrm{Re}}^2}{g}_i\mathsf{\Theta }\right]}^n-\frac{\partial p^{n+{\theta }_2}}{\partial x_i},i=\{1,2,3\}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \frac{{\rho }^{n+1}-{\rho }^n}{\mathsf{\Delta }t}=\frac{\mathsf{\Delta }\rho }{\mathsf{\Delta }t}& ={\left(\frac{1}{a^2}\right)}^n\frac{\mathsf{\Delta }p}{\mathsf{\Delta }t}=-\frac{\partial \left(\rho u_j^{n+{\theta }_1}\right)}{\partial x_j}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\rho }^n\frac{e_T^{n+1}-e_T^n}{\mathsf{\Delta }t}={\rho }^n\frac{\mathsf{\Delta }{e}_T}{\mathsf{\Delta }t}& ={\left[-\rho \frac{\partial \left(u_j{e}_T\right)}{\partial x_j}-\frac{\partial \left(p{u}_j\right)}{\partial x_j}+\frac{1}{\mathrm{Re}}\frac{\partial \left({\tau }_{ij}{u}_j\right)}{\partial x_i}+\frac{1}{\mathrm{Fr}}\rho g_j{u}_j+\frac{1}{\mathrm{Re}\mathrm{Pr}}\frac{\partial }{\partial x_j}\left(k\frac{\partial T}{\partial x_j}\right)\right]}^n\hfill \end{array}$$

(18)

where the continuity equation for incompressible fluids has been substituted by an equation of conservation of mass (17) that includes an artificial compressibility [65]. In that equation a represents the speed of sound in the fluid, which in the incompressible limit tends to infinity. The principal asset of this split algorithm is that it does not only apply to compressible flows, but also to incompressible flows, which makes it suitable for a wide variety of applications. In the above equations, ${\theta }_1$ and ${\theta }_2$ represent, respectively, the velocity and pressure relaxation factors. The explicit form of the algorithm is ${\theta }_2=0$, whereas in the semi-implicit form ${\theta }_2\in [0.5,1]$. In both cases the value of ${\theta }_1\in [0.5,1]$. In the numerical simulations presented in this paper we have considered the semi-implicit algorithm, setting the velocity relaxation factor to ${\theta }_1=1.0$, and the pressure relaxation factor to ${\theta }_2=1.0$. For this parameter choice, second-order accuracy in time is expected [61]. Moreover, the algorithm imposes a restriction to the time step size $\mathsf{\Delta }t$ and it is conditionally stable, as reported by [61]. In its semi-implicit form [66], the time step size must be $\mathsf{\Delta }t\le min\{h/|{\mathbf{u}}_i|,h^2/\left(2\nu \right)\}$, whereas in the explicit form, $\mathsf{\Delta }t\le h/\left|\mathbf{a}\right|$, leading to a more stringent time step in incompressible fluids. That is the reason why, for incompressible fluids, the semi-implicit form is cheaper in terms of computational cost, providing the required results in less time.

Equation (16) is coupled to Equation (17) and for this reason it cannot be used explicitly. To overcome this limitation, an auxiliary variable $\mathsf{\Delta }{u}^{**}$ is introduced, defined by

$$\begin{array}{cc}\mathsf{\Delta }{u}_i^{**}=u_i^{**}-u_i^n=\mathsf{\Delta }t[-{\displaystyle \frac{\partial \left(u_j{u}_i\right)}{\partial x_j}}+{\displaystyle \frac{1}{\mathrm{Re}}}{\displaystyle \frac{\partial {\tau }_{ij}}{\partial x_j}}+\hfill & {\displaystyle \frac{1}{\mathrm{Fr}}}{g}_i-{\displaystyle \frac{\mathrm{Gr}}{{\mathrm{Re}}^2}}{g}_i\mathsf{\Theta }\hfill \\ \hfill & +{\displaystyle \frac{\mathsf{\Delta }t}{2}}{u}_k{\displaystyle \frac{\partial }{\partial x_k}}\left({\displaystyle \frac{\partial \left(u_j{u}_i\right)}{\partial x_j}}-{\displaystyle \frac{1}{\mathrm{Fr}}}{g}_i+{\displaystyle \frac{\mathrm{Gr}}{{\mathrm{Re}}^2}}{g}_i\mathsf{\Theta }\right){]}^n.\hfill \end{array}$$

(19)

This equation does not contain the pressure term and the variable $u_i^{**}$ can thus be computed explicitly. When the pressure term is available, the following correction is applied

$$\mathsf{\Delta }{u}_i=u_i^{n+1}-u_i^n=\mathsf{\Delta }{u}_i^{**}-\mathsf{\Delta }t\frac{\partial p^{n+{\theta }_2}}{\partial x_i}-\frac{\mathsf{\Delta }{t}^2}{2}{u}_k\frac{\partial }{\partial x_k}\left(-\frac{\partial p^{n+{\theta }_2}}{\partial x_i}\right).$$

(20)

Expressing $\mathsf{\Delta }{u}_i^{n+1}$ as a function of $\mathsf{\Delta }{u}_i^{**}$ using (20) and substituting the result in (17) yields

$$\mathsf{\Delta }\rho ={\left(\frac{1}{a^2}\right)}^n\mathsf{\Delta }p=-\mathsf{\Delta }t\left[\frac{\partial u_i^n}{\partial x_i}+{\theta }_1\frac{\partial \mathsf{\Delta }{u}_i^{**}}{\partial x_i}-\mathsf{\Delta }t{\theta }_1\left(\frac{{\partial }^2{p}^n}{\partial x_i\partial x_i}+{\theta }_2\frac{{\partial }^2\mathsf{\Delta }p}{\partial x_i\partial x_i}\right)\right].$$

(21)

Finally, expressing $\mathsf{\Delta }{e}_T$ as a function of ${e_T}^n$ gives

$$\begin{array}{cc}\rho \mathsf{\Delta }{e}_T=\mathsf{\Delta }t[-{\displaystyle \frac{\partial \left(u_j(\rho e_T+p)\right)}{\partial x_j}}+{\displaystyle \frac{1}{\mathrm{Re}}}{\displaystyle \frac{\partial \left({\tau }_{ij}{u}_j\right)}{\partial x_i}}+{\displaystyle \frac{1}{\mathrm{Fr}}}\rho g_j{u}_j\hfill & +{\displaystyle \frac{1}{\mathrm{Re}\mathrm{Pr}}}{\displaystyle \frac{\partial }{\partial x_j}}\left(k{\displaystyle \frac{\partial T}{\partial x_j}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\mathsf{\Delta }t}{2}}{u}_j{\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{\partial (-u_i(\rho e_T+p))}{\partial x_i}}\right){]}^n.\hfill \end{array}$$

(22)

The way the variables of the problem are solved is summarized in the following step sequence:

• Calculate $\mathsf{\Delta }{u}_i^{**}$ from Equation (19).
• Calculate $\mathsf{\Delta }\rho$ or $\mathsf{\Delta }p$ from Equation (21).
• Calculate $\mathsf{\Delta }{u}_i$ from Equation (20), which yields the velocity at time $t_{n+1}$.
• Calculate $\mathsf{\Delta }e$ from Equation (22), which can be calculated in parallel with the other steps since the right hand side of Equation (22) does not depend on the variables at time $t_{n+1}$. Step 4 allows obtaining the value of the energy at time $t_{n+1}$.

The reason for choosing the CBS algorithm instead of others is to achieve second order convergence in both time and space using linear finite elements (TRIA and/or QUAD). Additionally, since the right hand side of the different equations of the algorithm are treated explicitly, they are easy to parallelize, as it has been done in the implementation of the algorithm.

3.2. Spatial Discretization

Following the characteristics-Galerkin approximation, Equations (19)–(22) are discretized in space using the finite element method [67,68,69]. To this end, the spatial domain is divided into a regular, unstructured triangulation ${\mathbb{T}}_h$ of non-overlapping elements. Associated with that triangulation, a conforming finite element space $V_h$ is defined composed of continuous, piecewise polynomials over each mesh element. The current implementation of the software uses only linear polynomials for all fluid variables, so higher order elements are not considered. Therefore, the numerical method is expected to be second order in space for smooth solutions in the $L^2$-norm. The advantage of the finite element method, as spatial discretization method, is that it allows to work with complex geometries, and to use fine or coarse finite elements in different zones of the domain (local mesh refinement) [70]. Moreover, the implemented code supports both TRIA and QUAD finite elements.

Hence, if m is the number of mesh points (or nodal points) of the triangulation ${\mathbb{T}}_h$, the numerical solution ${\varphi }_h$ belonging to the finite element space $V_h$ can be written as

$${\varphi }_h=\sum _{k=1}^m{N}_h^k{\varphi }^k,$$

(23)

where the summation spans over the set of nodes comprising the sampled domain, being k the index of the node, for $1\le k\le m$. ${\varphi }^k$ is the unknown variable evaluated at node k, and ${\left\{N_h^k\right\}}_{k=1}^m$ is the set of basis functions of $V_h$ satisfying $N_h^k\left({\mathbf{x}}_j\right)={\delta }_{kj}$, with ${\delta }_{kj}$ the Kronecker delta. Moreover, Equation (23) may be expressed in matricial form as

$${\varphi }_h={\mathbf{N}}_h\overline{\varphi }$$

(24)

with

$$\overline{\varphi }={\left[{\varphi }^1,\cdots ,{\varphi }^k,\cdots ,{\varphi }^m\right]}^T\mathrm{and}{\mathbf{N}}_h=\left[N_h^1,\cdots ,N_h^k,\cdots ,N_h^m\right].$$

(25)

In this problem the generic variable ${\varphi }_h$ will be replaced by the unknown fluid variables: the pressure, p, the velocity components, $u_i$, the total energy, $\rho e_T$, and the temperature, T.

As a final step to apply the standard Galerkin approximation via the Finite Element method, the weak formulation of Equations (19)–(22) must be obtained. Thus, equations can be weighted and integrated over the sampled domain using the basis functions $N_{h\varphi }^k\in V_{h0}^{\varphi }$, where $V_{h0}^{\varphi }$ is the finite element space of functions that take null values in the dirichlet boundary associated with the variable ${\varphi }_h=\{u_i,p,T\}$ (which means that it does not include the basis functions associated with the dirichlet nodes, whose values for the variable ${\varphi }_h$ are known).

The weak formulation of Equation (19) (Step 1 of the algorithm to compute $\mathsf{\Delta }{u}_i^{**}$ for $i=\{1,2,3\}$) in the standard Galerkin approximation is the following

$$\begin{array}{cc}{\displaystyle {\int }_{\Omega }}\hfill & {\displaystyle N_h^k\mathsf{\Delta }{u}_i^{**}d\Omega =\mathsf{\Delta }t{\left[-{\int }_{\Omega }{N}_h^k\frac{\partial \left(u_j{u}_i\right)}{\partial x_j}d\Omega -\frac{1}{\mathrm{Re}}{\int }_{\Omega }\frac{\partial N_h^k}{\partial x_j}{\tau }_{ij}d\Omega -{\int }_{\Omega }{N}_h^k\left(\frac{\mathrm{Gr}}{{\mathrm{Re}}^2}{g}_i\mathsf{\Theta }-\frac{1}{\mathrm{Fr}}{g}_i\right)d\Omega \right]}^n}\hfill \\ & {\displaystyle +\mathsf{\Delta }t{\left[\frac{1}{\mathrm{Re}}{\int }_{\Gamma }{N}_h^k{\tau }_{ij}{n}_{j}d\Gamma \right]}^n-{\left[\frac{\mathsf{\Delta }{t}^2}{2}{\int }_{\Omega }\frac{\partial }{\partial x_j}\left(u_j{N}_h^k\right)\left(\frac{\mathrm{Gr}}{{\mathrm{Re}}^2}{g}_i\mathsf{\Theta }-\frac{1}{\mathrm{Fr}}{g}_i+\frac{\partial \left(u_j{u}_i\right)}{\partial x_j}\right)d\Omega \right]}^n,}\hfill \end{array}$$

(26)

where $n_j$ denotes the j-th component of the unit outward normal vector to the boundary surface and $N_h^k\in V_h$, since it does not impose dirichlet boundary condition over $\mathsf{\Delta }{u}_i^{**}$.

To obtain the weak formulation of Equation (21) (Step 2 of the algorithm to compute $\mathsf{\Delta }p=p^{n+1}-p^n$), the equation should be multiplied by the $N_{hp}^k$ basis functions not associated with the dirichlet boundary nodes where pressure is prescribed $p={\tilde{p}}^{n+1}$ at ${\Gamma }_p$. This gives

$$\begin{array}{cc}{\displaystyle {\int }_{\Omega }{N}_{hp}^k\mathsf{\Delta }\rho d\Omega }\hfill & {\displaystyle ={\int }_{\Omega }{N}_{hp}^k\left(\frac{1}{a^2}\right)\mathsf{\Delta }pd\Omega =\mathsf{\Delta }t{\int }_{\Omega }\frac{\partial N_{hp}^k}{\partial x_j}\left[u_j^n+{\theta }_1\mathsf{\Delta }{u}_j^{**}-{\theta }_1\mathsf{\Delta }t\frac{\partial p^n}{\partial x_j}-{\theta }_1{\theta }_2\mathsf{\Delta }t\frac{\partial \mathsf{\Delta }p}{\partial x_j}\right]d\Omega }\hfill \\ & {\displaystyle -\mathsf{\Delta }t{\int }_{{\Gamma }_p\setminus \Gamma }{N}_{hp}^k\left[u_j^n+{\theta }_1\mathsf{\Delta }{u}_j^{**}-{\theta }_1\mathsf{\Delta }t\frac{\partial p^n}{\partial x_j}-{\theta }_1{\theta }_2\mathsf{\Delta }t\frac{\partial \mathsf{\Delta }p}{\partial x_j}\right]n_{j}d\Gamma }\hfill \end{array},$$

(27)

where the boundary ${\Gamma }_p\setminus \Gamma$ is the entire boundary of the domain minus the dirichlet pressure boundary. In the computation of the boundary integral on ${\Gamma }_p\setminus \Gamma$, the Dirichlet condition for the velocity component, $u_i$, may also be specified, in which case, the integrand will match the prescribed velocity, $u_i^n+{\theta }_1\mathsf{\Delta }{u}_i^{**}-{\theta }_1\mathsf{\Delta }t\partial p^n/\partial x_i-{\theta }_1{\theta }_2\mathsf{\Delta }t\partial \mathsf{\Delta }p/\partial x_i={\tilde{u}}_i$, at ${\Gamma }_u$.

In the same way, to build the weak formulation of Equation (20) (Step 3 of the algorithm to compute $\mathsf{\Delta }{u}_i=u_i^{n+1}-u_i^n$ for $i=\{1,2,3\}$), the basis functions $N_{hu}^k$ not associated with the dirichlet boundary nodes where velocity is prescribed, $u_i^{n+1}={\tilde{u}}_i$ at ${\Gamma }_u$, must be considered. This results in

$${\int }_{\Omega }{N}_{hu}^k\mathsf{\Delta }{u}_{i}d\Omega ={\int }_{\Omega }{N}_{hu}^k\mathsf{\Delta }{u}_i^{**}d\Omega -\mathsf{\Delta }t{\int }_{\Omega }{N}_{hu}^k\frac{\partial p^{n+{\theta }_2}}{\partial x_i}d\Omega -\frac{1}{2}\mathsf{\Delta }{t}^2{\int }_{\Omega }\frac{\partial \left(u_j^n{N}_{hu}^k\right)}{\partial x_j}\frac{\partial p^n}{\partial x_i}d\Omega$$

(28)

Finally, the weak formulation of the energy conservation Equation (22) (Step 4 of the algorithm to compute $T^{n+1}$ through $\mathsf{\Delta }{e}_T=e_T^{n+1}-e_T^n$) takes the form

$$\begin{array}{cc}{\displaystyle {\int }_{\Omega }{N}_{hT}^k\rho \mathsf{\Delta }{e}_T}\hfill & {\displaystyle d\Omega =\mathsf{\Delta }t\left[-{\int }_{\Omega }{N}_{hT}^k\frac{\partial }{\partial x_j}\left(u_j\left(\rho e_T+p\right)\right)d\Omega -\frac{1}{\mathrm{Re}}{\int }_{\Omega }\frac{\partial N_{hT}^k}{\partial x_j}\left({\tau }_{ij}{u}_j+\frac{1}{\mathrm{Pr}}k\frac{\partial T}{\partial x_j}\right)d\Omega \right.}\hfill \\ \hfill & {\displaystyle {\left.+{\int }_{\Omega }{N}_{hT}^k\frac{1}{\mathrm{Fr}}\rho g_j{u}_{j}d\Omega \right]}^n+\frac{\mathsf{\Delta }{t}^2}{2}{\left[{\int }_{\Omega }\frac{\partial }{\partial x_j}\left(u_j{N}_{hT}^k\right)\frac{\partial }{\partial x_j}\left(-u_j\left(\rho e_T+p\right)\right)d\Omega \right]}^n}\hfill \\ \hfill & {\displaystyle +\mathsf{\Delta }t\frac{1}{\mathrm{Re}}{\left[{\int }_{{\Gamma }_T\setminus \Gamma }{N}_{hT}^k\left({\tau }_{ij}{u}_j+\frac{1}{\mathrm{Pr}}k\frac{\partial h}{\partial x_i}\right)n_{i}d\Gamma \right]}^n,}\hfill \end{array}$$

(29)

where $N_{hT}^k$ are the basis functions not associated with the dirichlet boundary nodes where temperature (energy) is prescribed, $T^{n+1}=\tilde{T}$ at ${\Gamma }_T$. Moreover, in the computation of the boundary integral on ${\Gamma }_T\setminus \Gamma$, the value of the temperature gradient (heat flux) can be specified on some part of it (neumann boundary) as boundary condition.

4. Software Validation

In order to validate the implementation of the CBS algorithm described above, four test cases have been selected for study. Although the implementation of DynamFluid allows the integration of the Navier–Stokes equations in their general form (1)–(3), it can also be applied to fluid dynamics problems within the Boussinesq approximation (5) and (6). This is the approach chosen for the validation with the test cases presented in this section. All simulations are two-dimensional, hence we shall restrict our attention to cases with $x_i={(x,y)}^T$, $u_i={(u,v)}^T$, and $g_i={(g_x,g_y)}^T$. As discussed in the numerical section, in all test cases the velocity and pressure relaxation factors ${\theta }_1$ and ${\theta }_2$ have been set to unity

.

4.1. Lid-Driven Cavity Flow

The first step in the validation of a new CFD tool is to check if it is able to predict the flow behaviour in confined domains with simple boundary conditions. The flow in a cavity with an upper mobile lid, shown schematically in Figure 1, is considered as a standard software validation case for steady incompressible flows in confined geometries. In the lid-driven cavity, the fluid is confined in a 2D square cavity of side length $L(\equiv L_{\infty })$, with a zero velocity enforced by the no-slip condition in all the walls of the square cavity except in the top wall, where a uniform slip velocity $U(\equiv U_{\infty })$ is imposed. Additional parameters of the problem include the density, $\rho (\equiv {\rho }_{\infty })$, and viscosity, $\mu (\equiv {\mu }_{\infty })$, of the fluid, assumed here to be constant. Under isothermal conditions, the only governing parameter is the Reynolds number $\mathrm{Re}=UL/\nu$, based on the lid velocity, the side length of the squared cavity, and the kinematic viscosity of the fluid, $\nu =\mu /\rho (\equiv {\nu }_{\infty })$.

4.1.1. Literature Review

The lid-driven cavity flow has been extensively addressed by many authors [23,24,25,26,27,28]. Erturk et al. [23] performed numerical calculations of the two-dimensional steady incompressible lid-driven cavity flow for Reynolds number up to 21,000 with an extremely small residual error of the steady solution, finding a fourth vortex at the bottom left corner and a third vortex at the top left corner as the Reynolds number was increased, and obtaining a solution free of spurious oscillations despite the high Reynolds number.

Erturk [24] studied in detail the two-dimensional lid-driven cavity flow considering all physical, mathematical and numerical aspects, concluding that physically this flow is not two-dimensional but three-dimensional above a certain Reynolds number when the two-dimensional solution becomes unstable to small three-dimensional perturbations. Thus, while numerical solutions for the planar lid-driven cavity flow for high Reynolds numbers can be obtained, this flow can be considered fictitious. One important point for obtaining a solution at high Reynolds number is that a sufficiently fine mesh is needed so as to rule out any spurious oscillatory solution provided by the numerical method that prevents reaching the steady-state solution.

Yapici et al. [25] managed to perform several simulations up to Reynolds number 65,000, obtaining a steady solution with two new vortices in the bottom left and right corners of the square cavity for Reynolds number larger than 25,000.

Erturk et al. [29] also studied the benchmark problem of a driven skewed cavity flow for skew angles of 30${}^{\circ }$ and 45${}^{\circ }$. Using a very fine grid and highly accurate numerical solvers, they obtained numerical results for Reynolds number varying from 100 to 1000.

4.1.2. Boundary and Initial Conditions

In dimensionless form, a unit horizontal velocity is imposed at the upper wall, $u-1=v=0$, whereas the non-slip condition is imposed at the left, bottom and right walls, $u=v=0$. As an initial condition, both components of the velocity are set to 0 in all the points of the domain, except at the upper wall (lid) where the horizontal component is set to 1. The pressure is arbitrarily set to zero, $p=0$, at the lower-left corner of the cavity.

4.1.3. Convergence Analysis

In this section, we study the influence of the mesh in the numerical solution, which includes analyzing the influence of the mesh space discretization and cell size in the convergence of the solution. The method for estimating the convergence order of an algorithm requires to obtain a solution for different meshes with different element sizes. The meshes used for the convergence analysis are three uniform (i.e., equispaced) meshes using TRIA elements for the discretization of the domain with (mesh #1) $20\times 20$, (mesh #2) $50\times 50$ and (mesh #3) $100\times 100$ points, respectively. The coarser mesh is shown in Figure 2 for illustrative purposes.

The boundary conditions are applied to the mesh assuming the ramp condition: the velocity in the top of the cavity (moving lid) grows from zero in the corners of the cavity (left and right) until the non-dimensional value in the span of a cell/element. Alternatively, other researchers use the leaking lid formulation, where the velocity in the top lid is constant and equal to the imposed velocity and zero in the rest of the walls.

The Reynolds number used for showing the influence of the mesh element size in the accuracy of the solution is $\mathrm{Re}=1000$. The variable chosen for the convergence analysis is the horizontal velocity component u along the vertical line that goes through the center of the cavity. The steady-state solution is obtained by means of a transient simulation. It is assumed that the steady-state is reached when the relative value of the horizontal velocity in two successive time steps differ less than ${10}^{-9}$, that is, ${\sum }_{k=1}^N|u_k^{n+1}-u_k^n|/|u_k^n|\le {10}^{-9}$, the sum going through all N nodes in the mesh for the horizontal velocity component. Figure 2 shows that as the mesh is refined the solution tends to that obtained by Erturk et al. [23]. There is no appreciable difference between the results obtained with the $50\times 50$ point and $100\times 100$ point meshes.

For a TRIA element mesh composed by two-dimensional elements, as the one shown in Figure 3, an equivalent length for node k can be defined according to $h_k=\min (2{\mathrm{A}}_j/l_j)$, where j spans from 1 to the number of adjacent elements connected to node k, $A_j$ is the area of the j-th element adjacent to node k, and $l_i$ is the length of the opposite edge belonging to the j-th element adjacent to node k [71].

The convergence order of the algorithm has been estimated using the well-known grid convergence index calculation. Following Roache [72], the order of convergence for the algorithm is obtained solving the following equations iteratively for the three different meshes, starting from the coarse mesh #1 ($20\times 20$) to the more refined mesh #3 ($100\times 100$)

$$\begin{array}{cc}\hfill p& =\frac{\left|ln\left|E{r}_{32}/E{r}_{21}\right|+q\left(p\right)\right|}{ln{r}_{21}}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill q\left(p\right)& =ln\left(\frac{r_{21}^p-s}{r_{32}^p-s}\right)\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill s& =\mathrm{sign}(E_{32}/E_{21}),\hfill \end{array}$$

(32)

where $E{r}_{ij}=u_i-u_j$ and $r_{ij}=h_i/h_j$, with $u_i$ the horizontal velocity at mesh i and $u_j$ the horizontal velocity at mesh j, and $h_i$ the equivalent length for mesh i and $h_j$ the equivalent length for mesh j. This is translated into $E{r}_{32}$ being the difference in the estimating function between mesh #3 ($100\times 100$) and mesh #2 ($50\times 50$), and $E{r}_{21}$ between mesh #2 ($50\times 50$) and mesh #1 ($20\times 20$). Similarly, $r_{32}$ is the mesh ratio between mesh #3 and mesh #2, and $r_{21}$ between mesh #2 and mesh #1.

The iterative process starts with an initial guess for $q\left(p\right)$, e.g., $q\left(p\right)=0$ in this case. Using the value of q at a given iteration step, the value of p is calculated at the next iteration step. After a few iterations, the convergence order p is obtained.

Table 1. Mean and max error for the horizontal velocity.

Grid Comparison	${\mathit{Er}}_{{\mathit{L}}_1}$	${\mathit{Er}}_{{\mathit{L}}_{\mathit{\infty }}}$	${\mathit{r}}_{\mathit{ij}}$
$50\times 50$ vs. $20\times 20$	$0.00176$	$0.00748$	$2.0$
$100\times 100$ vs. $50\times 50$	$0.0138$	$0.0969$	$2.5$

shows the maximum norm of the error $E{r}_{L_{\infty }}=\max |u_i-u_j|$ ($p=2.71$), and the mean norm of the error $E{r}_{L_1}=(1/N){\sum }_{i=1}^N|u_i-u_j|$ ($p=1.98)$ for the implementation of the CBS algorithm discussed here applied to the lid-driven cavity flow.

4.1.4. Results for $\mathrm{Re}$ up to 10,000

Figure 4 compares numerical integrations of the lid-driven cavity flow carried out with DynamFluid with those obtained by Erturk [24] and Ghia et al. [26] for Reynolds number up to 10,000. Erturk et al. [23] reported that a $257\times 257$ mesh is needed in order to get a steady solution for Reynolds number up to 10,000, thus a refined mesh of $257\times 257$ elements was used for the simulations. As it can be seen, the agreement with previous results is excellent for $\mathrm{Re}\le$ 10,000, showing a slightly better agreement at low Reynolds number than at high Reynolds number.

4.2. Mixed Convection Flow in a Vertical Channel with Asymmetric Wall Temperatures

This is a mixed convection heat transfer problem in which an initially uniform flow develops in a slender vertical channel whose walls are at different temperatures. The flow is assumed to be two dimensional. The cold wall temperature, $T_c$, and hot wall temperature, $T_h(\equiv T_w)$, are also different from the incoming flow temperature, $T_{\infty }$, as illustrated in Figure 5. The parameters characterizing the flow are the channel width $H(\equiv L_{\infty })$ and length $L\gg H$, the cold and hot wall temperature differences, $T_c-T_{\infty }$ and $T_h-T_{\infty }$, and the uniform velocity of the fluid entering the rectangular channel, $U(\equiv U_{\infty })$. Additional parameters include the density, $\rho (\equiv {\rho }_{\infty })$, viscosity, $\mu (\equiv {\mu }_{\infty })$, and thermal conductivity, $k(\equiv k_{\infty })$, of the fluid. The non-dimensional numbers that emerge in this problem are the Reynolds number, $\mathrm{Re}=UH/\nu$, based on the inlet velocity and the channel width, the Grashof number, $\mathrm{Gr}=g\beta (T_h-T_{\infty })H^3/{\nu }^2$, defined in terms of the channel width and the hot wall temperature difference, and the wall temperature difference ratio ${\mathsf{\Theta }}_c=(T_c-T_{\infty })/(T_h-T_{\infty })$. The Prandtl number is assumed to be constant, $\mathrm{Pr}=0.72$. As shown in Figure 5, the acceleration of gravity is assumed to point opposite to the direction of the incoming flow. When the height of the channel is large compared to the channel width, the flow becomes fully developed far downstream, where the streamlines are parallel to each other and point vertically upwards.

4.2.1. Literature Review

The mixed convection flow in a vertical channel has been investigated by several authors. In a series of papers, Aung and Worku [30,31,32] studied the effects of buoyancy on the laminar vertical upward flow between parallel plates subject to two types of boundary conditions: (i) uniform wall temperatures, where the vertical walls may be at equal or different temperatures, and (ii) uniform heat fluxes, where the vertical walls may be subject to the same or different heat fluxes. They used the boundary layer approximation for their analysis when the dimensionless buoyancy parameter becomes $\mathrm{Gr}/\mathrm{Re}$ instead of $\mathrm{Gr}/{\mathrm{Re}}^2$. As a result, the controlling parameters are $\mathrm{Gr}/\mathrm{Re}$, the wall temperature difference ratio ${\mathsf{\Theta }}_c$, and the Prandtl number, $\mathrm{Pr}=0.72$, assumed to be constant. Aung and Worku [30] considered a set of values of $\mathrm{Gr}/\mathrm{Re}$ in the range $0,25,50,100$ and 250, founding that buoyancy extends the hydrodynamic entry length whereas it reduces the thermal entry length. As a general rule, no flow reversal was observed for sufficiently small values of $\mathrm{Gr}/\mathrm{Re}$, while it appeared for small values of ${\mathsf{\Theta }}_c$ and large values of $\mathrm{Gr}/\mathrm{Re}$ (the larger value of $\mathrm{Gr}/\mathrm{Re}$, the larger the downwards velocity in the flow reversal region). Aung and Worku [31] predicted that with constant heat fluxes the ratio $\mathrm{Gr}/\mathrm{Re}$ for flow reversal to happen is higher than with constant wall temperatures. In particular, for values of $\mathrm{Gr}/\mathrm{Re}$ up to 500 no flow reversal was observed when uniform heat flux boundary conditions are used. Flow reversal did not occur either for symmetrically heated walls. As part of their study, they developed an analytical theory [32] for fully developed mixed convection flow including flow reversal, obtaining analytical expressions for the velocity profile and the temperature profile that will be used below for validation purposes.

Ingham et al. [40] found that infinite duct walls also bring about flow reversal in the vicinity of the cold wall, while Ingham et al. [41] obtained numerical solutions of the problem for a steady laminar mixed convection flow in a vertical duct with parallel plates. They compared this flow with the case of pure forced convection finding that large values of Gr/Re cause reverse flow in the channel (similar to Aung and Worku [30,32]). Kim et al. [43] used an implicit finite difference scheme to solve the governing equations in the conjugate heat transfer flow established between two vertical plates subject to asymmetric wall temperatures. They found that the independent parameters are the Grashof number, the Prandtl number, the solid to fluid thermal conductivity ratio, the wall thickness to channel width ratio, the channel height to channel width ratio and the asymmetric heating parameter.

Gau et al. [38] studied experimentally the heat transfer process in a vertical channel comprised of two parallel plates (one heated uniformly and the opposite wall insulated), having a very large buoyancy parameter $\mathrm{Gr}/{\mathrm{Re}}^2$. From flow visualizations, they concluded that the reversal of the flow happened to be a V-shaped recirculating flow near the channel exit when $\mathrm{Gr}/{\mathrm{Re}}^2$ is greater than a critical value, although the reversal occurs initially downstream, but advances gradually downwards when $\mathrm{Gr}/{\mathrm{Re}}^2$ increases. El-Din [37], studied the flow development between two vertical plates with uniform heat and mass fluxes studying the effect of the thermal and mass transfer buoyancies. Cheng et al. [36] also studied this type of flow taking into account different boundary conditions (walls with constant temperature or with constant heat flux). Hamadah and Wirtz [39] investigated the effect of thermal buoyancy opposing the flow between two vertical plates with different boundary conditions (uniform asymmetric temperatures, uniform asymmetric heat fluxes, and the case of one wall at a uniform temperature and the other with uniform heat flux). Boulama and Galanis [35] provided analytical solutions for the fully-developed steady-state mixed convection flow past two vertical plates at constant temperatures and at constant heat flux (the two walls with the same type of boundary conditions or one with uniform temperature and the other with uniform heat flux).

Barletta et al. [34] studied the problem of a fully developed mixed convection flow with frictional heat generation in a vertical channel bounded by isothermal plane walls having the same temperature. Barletta [33] carried out an study of a laminar and fully developed flow with mixed convection in a rectangular and vertical duct where at least one of the two walls was isothermal, providing an analytical solution for the velocity and temperature fields. Desrayaud and Lauriat [73] investigated the flow reversal phenomena in a vertical channel with two parallel plates at symmetrically uniform heated walls when air is used as the fluid, for a laminar, mixed-convection flow with Reynolds number in the range $300\le \mathrm{Re}\le 1300$. The velocity and temperature profiles caused by the buoyancy forces were analyzed. Finally, Jeng et al. [42] investigated the mixed convection flow in a vertical channel with parallel walls at different temperatures: ${\mathsf{\Theta }}_c$ was in the range of $0\le {\mathsf{\Theta }}_c\le 1$, the Reynolds number was in the range of $1\le \mathrm{Re}\le 1000$, and the ratio $\mathrm{Gr}/\mathrm{Re}$ in the range of $0\le \mathrm{Gr}/\mathrm{Re}\le 500$. They found that when the streamwise coordinate was scaled out using the Reynolds number, the velocity and the temperature profile were independent of the Reynolds number for $\mathrm{Re}\ge 50$.

4.2.2. Boundary and Initial Conditions

In the simulations presented below, the right wall is assumed to be at a higher temperature than the left wall, $T_h>T_c$. Two cases are considered in the study. In the first one the cold wall is at the same temperature than the incoming fluid, $T_c=T_{\infty }$, or ${\mathsf{\Theta }}_c=0$, and in the second case the cold wall is at a temperature halfway between the incoming fluid and the hot wall, $T_c-T_{\infty }=0.5(T_h-T_{\infty })$, or ${\mathsf{\Theta }}_c=0.5$. The dimensionless boundary condition for the velocity, pressure and normalized temperature are

• Left, cold non-slip wall: $u=v=\mathsf{\Theta }-{\mathsf{\Theta }}_c=0$.
• Right, hot non-slip wall: $u=v=\mathsf{\Theta }-1=0$.
• Bottom side, incoming flow: $u=v-1=\mathsf{\Theta }=0$.
• Top side, outgoing flow: $p=u=\partial v/\partial y=\partial \theta /\partial y=0$.

Since the channel width is taken as length scale, the dimensionless width of the channel is one. By contrast, the dimensionless length of the channel is set to $H/L=80\gg 1$, which is long enough to ensure that the flow is fully developed in the outflow boundary (see Figure 5 for details; note that the height and width are not to scale).

4.2.3. Discussion of Results

As previously discussed, Aung and Worku [32] obtained analytical expressions for the velocity profile and the temperature profile in the fully developed mixed convection flow established far downstream the channel, namely

$$\left.\begin{array}{cc}\hfill v\left(x\right)& =\frac{\mathrm{Gr}}{\mathrm{Re}}(1-{\mathsf{\Theta }}_c)\left(-\frac{x^3}{6}+\frac{x^2}{4}-\frac{x}{12}\right)-6x^2+6x\hfill \\ \hfill \mathsf{\Theta }\left(x\right)& =x\hfill \end{array}\right\}0<x<1$$

(33)

the later corresponding to a linear temperature profile (i.e., a purely conductive heat flux) between the hot and cold walls. It is interesting to note that the analytical problem solved by Aung and Worku [30,31,32] included a simplified form of the energy equation that accounted only for thermal energy convection and heat conduction. As a consequence, the numerical results obtained with DynamFluid, which incorporates the general form of the energy Equation (10) or (13), can not be expected to match exactly the fully developed profiles given above.

Figure 6 shows the comparison between the fully developed profiles (33) with the numerical profiles predicted by DynamFluid in the outflow boundary corresponding to $\mathrm{Re}=100$, $\mathrm{Gr}=$ 25,000, and two values of ${\mathsf{\Theta }}_c=\{0,0.5\}$. In both cases the flow reversed near the cold wall ($x=0$) with higher downward velocities for ${\mathsf{\Theta }}_c=0$, and the temperature profiles are very close to linear. The minor differences observed between the numerical and analytical predictions can be attributed to the differences in the treatment of the energy equation. In summary, the figure shows an excellent agreement between the numerical profiles and the analytical solution for both values of ${\mathsf{\Theta }}_c$. This indicates that DynamFluid is able to predict accurately the velocity and temperature profiles in mixed convection problems under different conditions.

4.3. Isothermal Flow Past a Circular Cylinder

The flow past a circular cylinder is the prototypical benchmark case for the validation of external flows: either steady (the flow does not vary with time) or unsteady (the flow varies with time exhibiting either a transient or oscillatory behavior). Figure 7 shows an schematic representation of the flow: a circular cylinder of diameter $D(\equiv L_{\infty })$ surrounded by an unbounded fluid of density $\rho (\equiv {\rho }_{\infty })$ and viscosity $\mu (\equiv {\mu }_{\infty })$ with an uniform incoming velocity $U(\equiv U_{\infty })$. In the non-isothermal flow past a circular cylinder, to be considered in the next section, the cylinder is assumed to be at a uniform temperature $T_{\mathrm{w}}$ different from that of the incoming stream $T_{\infty }$. In both cases, the only non-dimensional parameter that appears in the limit of non-buoyant flows ($\mathrm{Ri}=\mathrm{Gr}/{\mathrm{Re}}^2\to 0$) to be considered here is the Reynolds number, $\mathrm{Re}=UD/\nu$, based on the incoming flow velocity, the diameter of the cylinder, and the kinematic viscosity of the fluid.

4.3.1. Literature Review

The flow past a non-heated circular cylinder at moderate Reynolds numbers exhibits flow separation and vortex shedding and has been studied numerically by many authors. De Sampaio et al. [56] obtained the solution of the incompressible Navier-Stokes equations with a Petrov-Galerkin method using an adaptive remeshing strategy applied to transient viscous flows. Ding et al. [54] used a mesh-free least square-based finite difference method to study the steady and unsteady viscous flow past a circular cylinder up to $\mathrm{Re}=200$. Park et al. [55] computed several flow properties such as the Strouhal number, the drag and lift coefficients, the pressure and vorticity distributions, the separation angle and the bubble separation length as a function of the Reynolds number up to $\mathrm{Re}=160$.

Nithiarasu [74] used the CBS algorithm in its fully explicit form and with artificial compressibility to predict the flow past a circular cylinder for a wide range of Reynolds numbers. Massarotti et al. [75] performed a comparison between the explicit and semi-implicit form of the CBS algorithm using benchmark test cases for both steady and unsteady flows, founding only slightly differences between both schemes for transient flows and identical results for steady flows. They concluded that the fully explicit version of the algorithm was an interesting option. Selvam [76] used an implicit CBS scheme with large eddy simulation in a 2D domain to compute the Strouhal number for high Reynolds numbers $\mathrm{Re}=\{{10}^4$, ${10}^5$, $5\times {10}^5$, ${10}^6$} and compared his results with the available experimental and numerical data. He reported that the measured reduction in the drag coefficient with the Reynolds number was not appropriately captured numerically and needed to be double-checked using a 3D model.

Qu et al. [77] performed several simulations using a FVM-based code for a wide range of Reynolds numbers ($\mathrm{Re}$ = 50–200) and studied the effect of the blockage ratio and the grid density in the vicinity of the cylinder wall, concluding that at lower Reynolds numbers the simulation requires a more uniform grid whereas at higher Reynolds number a finer grid near the vicinity of the cylinder wall is needed to resolve the thin viscous boundary layer that develops around the cylinder and eventually separates from it. Subhankar et al. [57] studied the critical Reynolds number that first causes flow separation from the cylinder wall, founding that for non-confined flows the value is $\mathrm{Re}=6.29$. Other authors that have addresses this problem include Sahin and Owens [78], Posdziech and Grundmann [79], Mittal and Raghuvanshi [80], Mittal and Kumar [81], Kieft et al. [82], and Jordan and Ragab [83].

4.3.2. Computational Domain, Boundary and Initial Conditions and Mesh Generation

The domain used for the simulations is shown in Figure 7. The cylinder is located a distance $L_{\varphi }$ from the inlet. The total length in the streamwise direction, $L>L_{\varphi }$, must be large enough so that the perturbations introduced by the cylinder become sufficiently small at the right limit of the domain for an outflow boundary condition to be used. In the simulations, the domain had a dimensionless size of $L/D=64$ with $L_{\varphi }/D=16$. Behr et al. [84] studied the influence of the distance of the lateral boundaries on the computation of two-dimensional unsteady incompressible flow past a circular cylinder. They concluded that the distance between the center of the cylinder and the lateral boundaries has a significant effect on the Strouhal number and other flow properties. In particular, they found that the minimum distance at which this influence vanishes for $\mathrm{Re}=100$ is 32 cylinder diameters. For this reason, in our simulations we chose also $H/D=64$ in order to rule out any influence of the boundary conditions in the solution obtained. The Reynolds numbers chosen for the validation campaign are 100 [56,85,86] and 200 [54,77], so that the lateral boundaries are expected not to perturb the results in either case.

The boundary conditions are schematized in dimensional form in Figure 7 for non-isothermal flow past a circular cylinder at uniform temperature, to be considered in the next section. For the isothermal flow considered here, the dimensionless boundary condition for the velocity and pressure are

• Left boundary, uniform incoming flow: $u-1=v=0$.
• Right boundary, outflow boundary condition: $p=\partial u/\partial x=v=0$.
• Top and bottom boundaries, symmetry boundary condition: $\partial u/\partial y=v=0$.
• Cylinder wall, non-slip condition: $u=v=0$.

As initial condition, the components of the velocity are set to 0 in all the points of the domain, except at the inflow boundary, where u is set to 1. The non-dimensional vertical velocity component v is set to 0 everywhere.

The strategy followed for the generation of the mesh was to use small elements in the vicinity of the cylinder and in the cylinder wake, whereas a coarser mesh was used in the rest of the domain, with the size of the elements varying gradually from the finest to the coarser regions as shown schematically in Figure 8. Several grids have been tested, from coarse to fine grids. The final grid used for the computations did not provide significantly different results than the previous coarser grid, but it was selected in order to get better estimates. In the final grid, the number of nodes was 77,011 and the number of TRIA elements was 153 318. Being D the cylinder diameter, the mesh size in the vicinity of the cylinder up to a concentric cylinder of diameter $1.1D$ was $0.01D$. The mesh size in the vicinity of the previous inner cylinder up to a concentric cylinder of radius $2D$ was less than $0.05D$. In a rectangular region with height $2D$, continuous to the outer concentric cylinder, the mesh size was less than $0.05D$. In the rest of the domain, the mesh size was less than $0.5D$.

4.3.3. Discussion of Results

The flow past a circular cylinder is periodic for $\mathrm{Re}=100$. Figure 9 shows the long-term variation of the vertical velocity v in a point located at the wake of the cylinder just in the middle of the outflow boundary ($x=64$, $y=32$). As can be seen, after an initial transient the flow becomes periodic with a given frequency f. As can be seen in Figure 10, counter rotating eddies are shed periodically with period $T=1/f$, giving rise to the well-known Karman vortex street in the cylinder wake.

Table 2. Strouhal number as a function of the blockage ratio for $\mathrm{Re}=100$.

$\mathit{D}/\mathit{H}$	$\mathbf{St}$
$1/16$	$0.1792$
$1/32$	$0.1703$
$1/64$	$0.1650$

shows the variation of the Strouhal number ($\mathrm{St}=fD/U_{\infty }$) associated with the vortex shedding process corresponding to $\mathrm{Re}=100$ and different blockage ratios ($D/H$). For the smallest blockage ratio, the predicted Strouhal number perfectly matches that obtained by De Sampaio [56] (0.165), while it deviates only $0.01\%$ from the value obtained by Ding et al. [54] and by Rahman et al. [85] (0.164). As the Reynolds number increases the influence of the blockage ratio becomes less important, hence the mesh with blockage ratio equal to $1/64$ has been used in the computations presented below for both $\mathrm{Re}=100$ and 200. For $\mathrm{Re}=200$ the computed Strouhal number is 0.1954, which deviates only $-0.18$% from the value obtained by Qu et al. [77] (0.1958) and $-0.3$% from the value obtained by Ding et al. [54] (0.196).

Table 3. Comparison between the Strouhal number for $\mathrm{Re}=\{100,200\}$ as predicted by DynamFluid and reported by previous authors.

$\mathbf{Re}$	[56]	[77]	[86]	[85]	[54]	Present Work
100	$0.165$	$0.1649$	$0.1569$	$0.164$	0.164	$0.165$
200	−	$0.1958$	$0.1957$	−	0.196	$0.1954$

compares the results obtained with DynamFluid with those of previous references. As can be seen, the Strouhal number obtained for the two Reynolds numbers considered in the study show very good agreement with the existing literature, which validates the performance of DynamFluid for the prediction of unsteady flows.

4.4. Flow Past a Heated Circular Cylinder with Forced Convection

The non-isothermal flow past a heated circular cylinder is similar to that past a non-heated cylinder with the particularity that the temperature of the cylinder, $T_w$, is now different from the temperature of the bulk fluid, $T_{\infty }$. In the presence of buoyancy forces, the temperature variations that cause the flow to exhibit additional features due to the effects of aiding or opposing buoyancy, but for simplicity attention will be restricted here to the case where the effects of buoyancy can be neglected. As the main difference with the mixed convection problem in a vertical channel, the Prandtl number is assumed here to be $\mathrm{Pr}=0.71$ (instead of $0.72$) in order to match the value used in previous works used here for validation.

4.4.1. Literature Review

Flow and heat transfer in forced convection past a circular cilinder has been studied by Apelt and Ledwich [87], Dennis et al. [46] and Rashid and Ahmad [88] up to $\mathrm{Re}=40$. Bitwas et al. [58] showed that in the absence of thermal buoyancy, the separation angle and the length of the recirculation bubble increase with the Reynolds number. They also showed that the average Nusselt number increases with increasing Reynolds number and the predicted results were in accordance with well-known experimental observations. Badr [89] studied numerically the influence of the flow direction, vertically upwards (parallel flow) vs. vertically downwards (opposing flow), accounting for buoyancy effects.

4.4.2. Computational Domain and Boundary Conditions

The computational domain used for the simulations is similar to that used in previous section, but in this case the length of the domain was reduced to $L/D=25$ cylinder diameters, and the height of the domain was reduced to $H/D=20$ cylinder diameters, which gives a blockage ratio of $0.05$. The center of the cylinder is located at $L_{\varphi }/D=10$ diameters from the inlet. As a result, the outflow boundary is located 15 diameters downstream the center of the cylinder. The boundary conditions imposed to the velocity and pressure fields are the same as in the previous section, with additional Dirichlet boundary conditions for the temperature ($\mathsf{\Theta }=0$) at the inlet and at the cylinder wall ($\mathsf{\Theta }=1$), zero heat flux boundary conditions at the upper and lower boundaries ($\partial \mathsf{\Theta }/\partial y=0$) as well as the outflow boundary ($\partial \mathsf{\Theta }/\partial x=0$).

4.4.3. Discussion of Results

Several simulations were carried out for increasing Reynolds numbers ranging from $\mathrm{Re}=10$ to 40 in the limit of non-buoyant flows ($\mathrm{Ri}=\mathrm{Gr}/{\mathrm{Re}}^2\to 0$). Figure 11 shows the steady-state solution reached in each simulation, including the streamlines and temperature contours obtained for the different Reynolds numbers under study. As previously reported by Williamson [90], for low Reynolds numbers ($\mathrm{Re}<49$) a vertically symmetric steady wake is formed with the structure shown schematically in Figure 12. The nomenclature used to describe the flow behind the cylinder is: A is the front stagnation point, B the rear stagnation point, C the stagnation point in the cylinder wake, E and F the upper and lower separation points, $L_s$ the eddy length of the recirculation region E-C-F-B-E and ${\varphi }_s$ the angle with respect to the horizontal line where the fluid detaches from the cylinder. The plots in the left panels of Figure 11 show a perfect symmetry with respect to the horizontal mid-line which passes through the center of the cylinder as previously reported by Williamson [90]. The temperature gradient at the front stagnation point, the eddy length, and the separation angle all increase monotonically with the Reynolds number. For reference purposes,

Table 4. Comparison between the eddy length ($L_s$) and the separation angle (${\theta }_s$) for several Reynolds numbers as predicted by DynamFluid and reported by previous authors.

Variable	$\mathbf{Re}$	[44]	[59]	[58]	[54]	[55]	DynamFluid
	10	$0.504$	$0.498$	$0.52$	$0.504$	$0.51$	$0.512$
	15	−	$1.162$	$1.189$	−	−	$1.227$
$2L_s/D$	20	$1.88$	$1.844$	$1.865$	$1.86$	$1.87$	$1.866$
	25	−	−	$2.517$	−	−	$2.548$
	40	$4.69$	$4.65$	$4.424$	$4.4$	$4.59$	$4.480$
Variable	$\mathrm{Re}$	[44]	[59]	[58]	[54]	[55]	DynamFluid
	10	$29.6$	$29.3$	$29.12$	$30.0$	−	$28.57$
	15	−	$38.6$	$38.57$	−	−	$38.57$
${\theta }_s$ (${}^{\circ }$)	20	$43.7$	$43.65$	$43.64$	$44.1$	−	$43.58$
	25	−	−	$46.89$	−	−	$47.14$
	40	$53.8$	$53.55$	$53.1$	$53.5$	−	$51.43$

shows the variation with the Reynolds number of the non-dimensional eddy length and the separation angle. Note that the eddy length is made non-dimensional with the cylinder radius ($D/2$) to be consistent with previous references, and that the angles are measured from the rear stagnation point. As can be seen, the results obtained with DynamFluid show excellent agreement with the values obtained by other authors [44,54,58,59] for all Reynolds numbers under study.

The local heat lost from the cylinder by heat conduction to the fluid is given by the local Nusselt number

$$\mathrm{Nu}=-\frac{\partial \mathsf{\Theta }}{\partial n},$$

(34)

where n denotes the outward normal to the cylinder surface. The local value of $\mathrm{Nu}$ along the cylinder surface was computed and compared with previous results taken from the literature [46,58,89], as can be seen in Figure 13. The Nusselt number computed with DynamFluid shows good agreement with the results by Biswas et al. [58] for all Reynolds number under study, showing also good agreement in the rear stagnation point with other authors [46,89]. However, it deviates slightly in the front stagnation point from the value computed by Bard [89] and Dennis et al. [46]. The results shown in Figure 13 indicate that Badr [89] overstimated the value of the Nusselt number in the range Re = 20–40, while Dennis et al. [46] overestimated the value of the Nusselt number in the range Re = 10–40. Note that the works by Dennis et al. [46] and Bard [89] were published in 1968 and 1984, when the precision of numerical computations was still far from that of today. These results indicate that DynamFluid is able to reproduce with good agreement results of forced convection heat transfer problems, particularly those published more recently, exhibiting only local small variations that can be attributed to the imprecision of the numerical results reported decades ago.

5. Conclusions

A new CFD software tool (DynamFluid) based on the finite element method and the characteristic-based-split algorithm has been presented and validated against four benchmark constant-density and Boussinesq-type non-isothermal two-dimensional flows, showing excellent agreement with previous results taken from the literature. The test cases have comprised both stationary problems, such as the lid-driven cavity flow, mixed convection in a vertical channel, or flow past a heated circular cylinder at low Reynolds numbers, as well as unsteady problems, such as isothermal flow past a circular cylinder at moderate Reynolds numbers, covering a wide range of Reynolds, Grashof and Richardson numbers. The results provide the prospect users high confidence for the application of this software to other 2D fluid dynamic problems of interest, particularly those involving isothermal and non-isothermal incompressible flows under the Boussinesq approximation. In the near future, DynamFluid is expected to cover a broader scope of flow regimes, including fully compressible three-dimensional flows, non-constant density flows, and general non-isothermal flows.