Combined Newton-Raphson and Streamlines-Upwind Petrov-Galerkin iterations for nano-particles transport in buoyancy driven flow
CCombined Newton-Raphson and Streamlines-Upwind Petrov-Galerkin iterationsfor nano-particles transport in buoyancy driven flow
M. K. RIAHI a,d, ∗ , M. Ali c,d , Y. Addad c,d , E. Abu-Nada b a Department of Applied Mathematics, Khalifa University, PO Box 127788, Abu Dhabi, UAE b Department of Mechanical Engineering, Khalifa University, PO Box 127788, Abu Dhabi, UAE c Department of Nuclear Engineering, Khalifa University, PO Box 127788, Abu Dhabi, UAE d Emirates Nuclear Technology Center, Khalifa University, PO Box 127788, Abu Dhabi, UAE
Abstract
The present study deals with the finite element discretization of nanofluid convective transport in an enclosure withvariable properties. We study the Buongiorno model, which couples the Navier-Stokes equations for the base fluid,an advective-di ff usion equation for the heat transfer, and an advection dominated nanoparticle fraction concentrationsubject to thermophoresis and Brownian motion forces. We develop an iterative numerical scheme that combinesNewton’s method (dedicated to the resolution of the momentum and energy equations) with the transport equationthat governs the nanoparticles concentration in the enclosure. We show that Stream Upwind Petrov-Galerkin regular-ization approach is required to solve properly the ill-posed Buongiorno transport model being tackled as a variationalproblem under mean value constraint. Non-trivial numerical computations are reported to show the e ff ectiveness ofour proposed numerical approach in its ability to provide reasonably good agreement with the experimental resultsavailable in the literature. The numerical experiments demonstrate that by accounting for only the thermophoresis andBrownian motion forces in the concentration transport equation, the model is not able to reproduce the heat transferimpairment due to the presence of suspended nanoparticles in the base fluid. It reveals, however, the significant rolethat these two terms play in the vicinity of the hot and cold walls. Keywords:
Nanofluid, Navier-Stokes equation, Newton-Raphson method, Advection dominated equation, Finiteelement method, Strem-Upwind Petrov-Galerkin.
1. Introduction
Natural convection, or natural circulation, is a phenomenon in which fluid recirculates due to applied temperaturedi ff erence, where hot fluid (light) tends to rise up, while colder fluid (heavy) tends to fall down. This phenomenonoccurs in several engineering applications such as chips cooling, large compartment ventilation, passive cooling inheat exchangers, and further ocean dynamics and weather applications.The suspension of nano-sized particles in base fluids, a mixture referred to as nanofluid, represents an attractivemethod in heat transfer engineering problems during the last two decades [1, 2, 3]. Several engineering applications inheat transfer were investigated including natural convection, combined convection, heat transfer in electronic cooling,and renewable energy [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Due to the high number of publications in using nanofluidsto enhance the heat transfer rate in thermal engineering systems, several reviews were conducted, such as Jabbariet al. [13] and Khodadadi et al. [14], Fan and Wang [15], Kakac¸ and Pramuanjaroenkij [16], Buongiorno et al.[17], Sheikholeslami and Ganji [18], and Manca et al. [2]. The role of nanofluids in augmenting the heat transferrate in forced convection is accepted in the research community where the nanofluids are found to be very usefulin enhancing the performance of forced convective flows. Conversely, the role of nanofluids is still controversial innatural convection. For example, most theoretical studies reported enhancement in heat transfer due to the dispersion ∗ Corresponding author.
Email address: [email protected] (M. K. RIAHI)
Preprint submitted to Journal January 01,2021 a r X i v : . [ phy s i c s . f l u - dyn ] J a n f the nanoparticles in base fluids. However, several experiments reported a deterioration in heat transfer due tothe addition of nanoparticles to base fluids, Wen and Ding [19], Li and Peterson [20], Ho et al. [21], and Putra etal. [22]. In some experimental measurements, a significant increase in the e ff ective thermal conductivity for themixture of certain types of small solid particles with a diluted water has been noticed (see [21] and [23]). Due to thisclaimed significant improvement of transport properties of the nanofluids, nanofluid technology has attracted manyapplications including the cooling of electronic chips, cooling of smaller internal combustion engines, and biomedicaltechnology.Furthermore, along its development and use, the nanofluid technology has been controversial on some of its as-pects. In fact, increasing the nanoparticle concentration does not necessarily mean an improvement of the heat transferrate. This is indeed, a debate that until today researchers still does not converge to a common conclusion on the theheat transfer enhancement using nanoparticles. Nonetheless, the experimental results of Ho et al. [21] provide a clearevidence that heat transfer enhancement using nanofluid is actually possible with a low nanoparticles concentration(i.e. for volume fractions below 2%). For high Rayleigh number in particular, the dispersed nanoparticles in thebase fluid are able to contribute to heat transfer rate enhancement of around 18% when compared to pure water. Forconcentrations equal to or greater than 2%, on the other hand, it was observed that the addition of nanoparticles wouldhave a negative impact on the heat transfer rate.One of the mostly used mathematical model used to study the heat transfer of nanofluid is the the Buongiornomodel [5]. In this model, the flow around the nanoparticles is regarded as a continuum. It assumes that the onlymechanisms causing the nanoparticles to develop a slip velocity with respect to the base fluid are; i) the Browniandi ff usion resulting from continuous collisions between the nanoparticles and the molecules of the base fluid and ii)the thermophoresis representing the nanoparticles di ff usion due to temperature gradient in the domain. Other slipmechanisms considered by Buongiorno in his analysis, namely, the di ff usiophoresis, the lift force, the fluid drainage,and gravity settling are considered negligible. As a result, the model translates the fact that the Brownian di ff usionand thermophoresis, are the only forces responsible for the di ff usivity of nanoparticle concentrations, hence, togetherwith the main stream they update the density of the nanofluid via a convection dominated partial di ff erential equations(PDE).In this paper, we use the model Buongiorno model to investigate the heat transfer enhancement using nanofluid.We develop and analyze a new numerical iterative scheme based on a finite element method (FEM) of the steadystate PDEs. The nanofluid model at hand, consists then of four equations i) the continuity equation ii) the momentumequation, which is the classical Navier-Stokes equation subject to the buoyancy force, iii) the energy equation whichis the convected heat equation and iv) the nanoparticle transport equation. The latter, is an advection dominated PDEs,as per the P´eclet number (ratio of the convection rate and the di ff usion rate) is very high, which leads to numericalspurious while using finite di ff erence method or FEM [24]. The convected dominated nanoparticle transport equationis, indeed, ill-posed and has to be regularized in order to stabilize the calculation, see for instance [25] and also [26, 27]for an adaptive procedure based on a posteriori error estimate.The momentum and energy equations are coupled through the convection and buoyancy terms, while the nanopar-ticle transport equation plays the role of fluid density regulator that has a crucial role in the variation of the fluidviscosity, hence on the shear stresses of the fluid. Particles migration is also impacted by the thermophoresis forcesand other subs-scale forces modeled through a Brownian motion. A theoretical mathematical functional analysis ofthe Buongiorno model is studied in [28] where a mollified regularizing problem is set up and is shown to be weaklyconvergent toward the nanofluid Buongiorno model. A similar nanofluid model has been derived and studied in [29].The resulting coupled system represents the dynamics of the nanofluid inside a di ff erentially heated cavity. Wesupplement the problem at hand with the following appropriate boundary conditions; non-slip velocity, adiabatichorizontal walls, specifically heated and cooled walls (Dirichlet non homogeneous conditions), and the Neumannhomogeneous boundary condition for the nanoparticle transport equation describing simply the fact that none of theparticles is allowed to exit the enclosure (particles-flux is null).For the numerical discretization of such Oberbeck-Boussinesq equations we refer to [30, 31, 32]. We also refer tothe book [33] for a complete finite element analysis for the Navier-Stokes equations. In the present work we considerthe pair P P P P P Nomenclature table ( x , y ) coordinate system (m) n outward normal vector L width of the cavity g gravitational acceleration ( ms − ) u (cid:63) dimensional velocity ( ms − ) u horizontal velocity component p (cid:63) dimensional pressure term p dimensionless pressure term θ (cid:63) temperature (K) θ temperature φ (cid:63) nanoparticles’ volume concentration (%) φ nanoparticles’ volume concentration φ b nanoparticles’ bulk volume concentration (%) C denoting Correlation c p specific heat ( J kg − K − ) Nu Nusselt number, Nu = − ( L / ∆ θ )( ∂θ/∂ x ) w Ra Rayleigh numberPr Prantl numberPe Peclet numberLe Lewis numberSc Schmit number ρ density ( kg m − ) k thermal di ff usivity ( m s − ) µ dynamic viscosity ( N kg − s − ) α kinematic viscosity ( m s − ) β volumetric expansion coe ffi cient of the fluid D (cid:63)ω dimensional Brownian di ff usion coe ffi cient D (cid:63)θ dimensional thermophoretic di ff usion coe ffi cient D ω Brownian di ff usion coe ffi cient D θ thermophoretic di ff usion coe ffi cient π m variables for the momentum equation π e variables for the energy equation π p variables for the np transport equation (cid:63) dimension superscriptp particle subscriptnf nanofluid subscriptbf base fluid subscriptnp nanoparticle subscriptC cold wallH hot wall
2. Equations settings
Following Buongiorno model for incompressible nanofluid flow and using Boussinesq approximation for density,we describe the natural convection in a deferentially heated squared cavity. Thus the domain of computation Ω (Figure 1) is a simply connected domain with Lipschitz boundary ∂ Ω . It is assumed that the fluid has reached astatistically time invariant state, a reason for which we shall study the steady state of the problem. The description ofthe dimensional and non-dimensional problems is reported in the sequel subsections. Our focus will be on the long time statistically invariant state known as steady state. The flow is therefore consid-ered time-invariant. The heat transfer equations governing the physics consist of a coupling between i) the momentumequation associated with the continuity equation, ii) the heat equation and iii) the nanoparticle transport equation and3 (cid:63) y (cid:63) u (cid:63) = , ∂θ (cid:63) ∂ n = , ∂φ (cid:63) ∂ n = u (cid:63) = , ∂θ (cid:63) ∂ n = , ∂φ (cid:63) ∂ n = u (cid:63) = , θ (cid:63) = θ (cid:63) h , ∂ φ (cid:63) ∂ n = u (cid:63) = , θ (cid:63) = θ (cid:63) c , ∂ φ (cid:63) ∂ n = g (cid:63) L Figure 1: Geometric sketch of the considered di ff erential cavity and the boundary conditions for the velocity, temperature and nanoparticle concen-tration. write as follows: ∇ (cid:63) · u (cid:63) = Ω (cid:0) u (cid:63) · ∇ (cid:63) (cid:1) u (cid:63) = − ρ nf ∇ (cid:63) p (cid:63) + ρ nf ∇ (cid:63) · (cid:18) µ (cid:63) nf (cid:0) ∇ (cid:63) u (cid:63) + ( ∇ (cid:63) u (cid:63) ) t (cid:1) (cid:19) + g (cid:63) ( ρ ∞ − ρ nf ) ρ nf on Ω (cid:0) u (cid:63) · ∇ (cid:63) θ (cid:63) (cid:1) = ρ c p ) nf ∇ (cid:63) · (cid:16) k (cid:63) nf ( θ (cid:63) , φ (cid:63) ) ∇ (cid:63) θ (cid:63) (cid:17) + (cid:32) D (cid:63)ω ∇ (cid:63) φ (cid:63) · ∇ (cid:63) θ (cid:63) + D (cid:63)θ (cid:63) θ (cid:63) C ∇ (cid:63) θ (cid:63) · ∇ (cid:63) θ (cid:63) (cid:33) on Ω (cid:0) u (cid:63) · ∇ (cid:63) φ (cid:63) (cid:1) = ∇ (cid:63) · (cid:32) D (cid:63)ω ∇ (cid:63) φ (cid:63) + D (cid:63)θ (cid:63) θ (cid:63) C ∇ (cid:63) θ (cid:63) (cid:33) on Ω ρ ∞ − ρ nf = ( ρβ ) nf ( θ (cid:63) − θ (cid:63) c ) (1)Supplemented by the non-slip boundary condition for the fluid velocity, Dirichlet non-homogeneous for the tempera-ture di ff erential in two opposite sides of the domain, Neumann homogeneous for the adiabatic boundaries, and all overNeumann homogeneous for the nanoparticle concentration. Figure 1 showcases the considered geometric domain forthe numerical simulation. In Eq.(1) the e ff ective density, specific heat, and thermal expansion coe ffi cients of nanofluidare given by ρ nf = ρ bf (1 − φ (cid:63) ) + ρ np φ (cid:63) , (2)( ρ c p ) nf = ( ρ c p ) bf (1 − φ (cid:63) ) + ( ρ c p ) np φ (cid:63) , (3)( ρβ ) nf = ( ρβ ) bf (1 − φ (cid:63) ) + ( ρβ ) np φ (cid:63) . (4)The Brownian di ff usion coe ffi cient D (cid:63)ω is defined using the Einstein-Stokes equation in the following form D (cid:63)ω ( θ (cid:63) ) = k b πµ bf ( θ (cid:63) ) d p θ (cid:63) , (5)where k b = . · − J / K stands for the Boltzmann’s constant. The thermophoretic di ff usion coe ffi cient D (cid:63)θ (cid:63) isdefined as D (cid:63)θ ( θ (cid:63) , φ (cid:63) ) = . k bf k bf + k p µ bf ( θ (cid:63) ) ρ bf φ (cid:63) . (6)4he dynamic viscosity for water is defined as follows µ bf ( θ (cid:63) ) = . · − · . / ( θ (cid:63) − . (7)For the thermal conductivity and viscosity of nanofluid many correlations have been derived [36, 37, 23, 38].These correlations are generally nonlinear with respect to their variables. For this reasons we shall keep the correlationas a general function as such for the e ff ective viscosity for nanofluid: µ nf µ bf = C µ nf ( φ, θ ) and for the e ff ective thermalconductivity k nf k bf = C k nf ( φ, θ ) both as a potential correlations of φ and θ . We start the non-dimensionalization of the equations by defining the following non-dimensional variables: u = ( L /α ) u (cid:63) , x = x (cid:63) / Ly = y (cid:63) / L , z = z (cid:63) / L φ = φ (cid:63) /φ , θ = ( θ (cid:63) − θ (cid:63) C ) / ( θ (cid:63) H − θ (cid:63) C ) µ nf ( θ, φ ) = µ (cid:63) nf ( θ (cid:63) , φ (cid:63) ) /µ bf ( θ (cid:63) C ) , k nf ( θ, φ ) = k (cid:63) nf ( θ (cid:63) , φ (cid:63) ) / k bf D ω = D (cid:63)ω ( θ (cid:63) ) / D (cid:63)ω ( θ (cid:63) C ) , D θ = D (cid:63)θ ( θ (cid:63) , φ (cid:63) ) / D (cid:63)θ ( θ (cid:63) C , φ ) , (8)Hence the corresponding equation writes ∇ · u = u · ∇ ) u = − π m1 ( φ ) ∇ p + π m2 ( φ ) ∇ · (cid:18) µ nf ( θ, φ ) (cid:0) ∇ u + ( ∇ u ) t (cid:1) (cid:19) + π m3 ( φ ) θ ( u · ∇ θ ) = π e1 ( φ ) ∇ · ( k nf ( θ, φ ) ∇ θ ) + (cid:16) π e2 ( φ ) D ω ∇ φ · ∇ θ + π e3 ( φ ) D θ ∇ θ · ∇ θ (cid:17) ( u · ∇ φ ) = π p1 ∇ · ( D ω ∇ φ ) + π p2 ∇ · ( D θ ∇ θ ) , (9)where the variable properties are defined as follows: π m1 ( φ ) = (1 − φ + φ ρ np ρ bf ) − π m2 ( φ ) = Pr π m1 π m3 ( φ ) = PrRa nf (1 − φ (cid:63) )(1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf + φ (cid:63) (1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf β np β bf π e1 ( φ ) = (1 − φ + φ ρ np φ (cid:63) np ρ bf φ (cid:63) bf ) − π e2 ( φ ) = PrSc φ b (cid:32) (1 − φ ) ρ bf φ (cid:63) bf ρ np φ (cid:63) np + φ (cid:33) π e3 ( φ ) = St PrSc (cid:32) θ (cid:63) H − θ (cid:63) C θ (cid:63) H (cid:33) (cid:32) (1 − φ ) ρ bf φ (cid:63) bf ρ np φ (cid:63) np + φ (cid:33) − π p1 = D ω α = π p2 = St PrSc (cid:32) θ (cid:63) H − θ (cid:63) C θ (cid:63) C (cid:33) φ b . and the non-dimensional numbers are given by:Rayleigh number Ra = g ( ρβ ) bf ( θ (cid:63) H − θ (cid:63) C ) L / ( µ bf ( θ (cid:63) C ) α ) , Prantle number Pr = µ bf ( θ (cid:63) C ) / ( ρ bf α ) , Lewis number Le = α/ D (cid:63)ω ( θ (cid:63) C ) , Peclet number Pe np = (cid:107) u (cid:107) / (cid:16) π p D ω (cid:17) .
5e proceed with the numerical discretization of the governing equations. We propose to decouple the whole systeminto two parts; The first part consists in solving the momentum equations along with the energy equation, the secondpart consists in solving the stabilized nanoparticle equation with an additional constrain to assess that the bulk amountof nanoparticle in the fluid is constant during the calculation. The main motivation behind this decoupling is to reducethe non-linearity of the equations (momentum and energy) in term of the nanoparticle volume fraction and give roomto the numerical scheme to stabilize itself iteratively for the transport equation. Finally, we re-iterate through thesetwo parts in order to update the coe ffi cients and the density of the fluid. More details are described in section 3 below.
3. Variational formulations and approximation tools settings
The numerical scheme we propose in this work is based on the finite element discretization of the dimensionlessequations. As the steady problem at hand presents non-linearity in term of the velocity and temperature, a direct linearalgebra solver, is, therefore, impractical to solve the coupled system. One has to relax the non-linearity and set upan iterative numerical scheme that converges to the desired solution. Newton’s method plays a key role here and hasbeen shown to be very e ffi cient in term of convergence [33].We consider homogeneous Dirichlet boundary conditions for the velocity, i.e. u = ∂ Ω . Let us consider thefollowing Hilbert spaces for the temperature, velocity and pressure as follows: T = H ( Ω ) , U = H ( Ω ) × H ( Ω ) , P = { p ∈ L ( Ω ) | (cid:90) Ω p δω = } . where L ( Ω ) = { u | (cid:90) Ω | u | δω < ∞} , and H ( Ω ) = (cid:110) u ∈ L ( Ω ) | ∇ u ∈ L ( Ω ) (cid:111) , and H ( Ω ) = { u ∈ H ( Ω ) | u | ∂ Ω = } . The Sobolev space L ( Ω ) is endowed with the usual inner product denoted by the duality pairing ( · , · ), that generatesthe L -norm (cid:107) · (cid:107) . In order to present a general algorithm, we shall separate nanoparticles equation from the momentum and energyequations while we are processing. The generality here is mainly targeting the use of any correlations (as there is avariety available in literature). This gives our approach a flexibility to treat di ff erent type of nanofluid and considerdi ff erent (possibly highly nonlinear) thermal and viscosity correlations. The idea behind the following notation andcalculation of the tangent equations is that we formulate the problem as a multivariable function F and use the classicalNewton-Raphson iterations that reads u k + p k + θ k + = u k p k θ k − (cid:16) DF (cid:16) ( u k , p k , θ k ) (cid:17)(cid:17) − F (cid:16) ( u k , p k , θ k ) (cid:17) , where DF ( X k ) stands for the Jacobian (isomerism in U × P × T ).In practice, we define the coupled momentum-energy variational form as such for every trial ( v , q , ζ ) ∈ U × P × T we have F ( v , q ,ζ ) : U × P × T −→ R ( u , p , θ ) (cid:55)−→ F ( v , q ,ζ ) ( u , p , θ ) : = (( u · ∇ u ) , v ) + π m (cid:0) µ nf ( θ, φ ) (cid:0) ∇ u + ( ∇ u ) t (cid:1) , ∇ v (cid:1) − ( π m ∇ p , v ) − (cid:16) π m u , ∇ q (cid:17) − (cid:16) π m θ ( u · e z ) , v (cid:17) + ε ( p , q ) + (( u · ∇ θ ) , ζ ) + ( π e1 ( φ ) ( k nf ( θ, φ ) ∇ θ ) · ∇ ζ ) + ( π e2 ( φ ) D ω ( ∇ φ, ∇ θ ) , ζ ) + (cid:16) π e3 ( φ ) D θ ( ∇ θ · ∇ θ ) , ζ (cid:17) . (10)6hen we define the vector field F ( u , p , θ ) ∈ U × P × T as the vector that satisfies the following inner product formulafor every ( v , q , ζ ) ∈ U × P × T ( F ( u , p , θ ) , v q ζ ) = F ( v , q ,ζ ) ( u , p , θ ) . Besides, we define the tangent coupled momentum-energy variational trilinear form as such for every trial ( v , q , ζ ) ∈ U × P × T we have DF ( v , q ,ζ ) ( u , p , θ ) : U × P × T −→ R DF ( v , q ,ζ ) ( u , p , θ )( δ u , δ p , δθ ) : = (( δ u · ∇ u ) , v ) + (( u · ∇ δ u ) , v ) + (cid:16) π m µ nf ( θ, φ ) (cid:16) ∇ δ u + ( ∇ δ u ) t (cid:17) , ∇ v (cid:17) − (cid:16) π m ∇ δ p , v (cid:17) − (cid:16) π m δ u , ∇ q (cid:17) − (cid:16) π m δθ ( u · e z ) , v (cid:17) − (cid:16) π m θ ( δ u · e z ) , v (cid:17) + ε ( δ p , q ) + (( δ u · ∇ θ ) , ζ ) + (( u · ∇ δθ ) , ζ ) + (cid:16) π e1 ( φ ) ( k nf ( θ, φ ) ∇ δθ ) , ∇ ζ (cid:17) + (cid:16) π e2 ( φ ) D ω ( ∇ φ · ∇ δθ ) , ζ (cid:17) + (cid:16) π e3 ( φ )2 D θ ( ∇ θ · ∇ δθ ) , ζ (cid:17) , (11)and we define the linear operator DF ( u , p , θ ) that satisfies the following equation( DF ( u , p , θ )( δ u , δ p , δθ ) , v q ζ ) = DF ( v , q ,ζ ) ( u , p , θ )( δ u , δ p , δθ ) . ∀ ( v , q , ζ ) ∈ U × P × T . The variational formulation for the nanoparticle concentration transport equation reads: For avery ϕ ∈ X : = H ( Ω )find φ ∈ H ( Ω ) such that a ( φ, ϕ ) : = (( u · ∇ φ ) , ϕ ) + (cid:16) π p1 ( D ω ∇ φ ) , ∇ ϕ (cid:17) − (cid:32) π p1 D ω ∂φ∂ n (cid:33) Γ − (cid:32) π p2 D θ ∂θ∂ n (cid:33) Γ = . ( b ( θ ) , ϕ ) = π p2 ( D θ ∇ θ, ∇ ϕ ) . (12)After applying the boundary conditions. This problem is not well posed. Indeed, it does not satisfy the inf-supcondition. We shall discuss this concern in the sequel where we introduce a regularization parameter ε acting as areaction term. In the numerical code this parameter is taken very small ( ε ≈ . e −
10) to ensure satisfaction of theinf-sup condition and not to harm the model, as the nanoparticles concentration does not exhibit a reaction within themixture.
4. Finite element algorithm
For the steady states under consideration we use i) standard Taylor-Hood finite element approximation [34] forthe space discretization of the Navier-Stokes equations, approximating the velocity field with P2 finite elements, andthe pressure with the P1 finite element. We assume we have the triangulation T h of the computational domain Ω , suchthat Ω = ∪ neln = e Ω e we seek for approximated solution over the finite dimensional vector spaces U h × P h × T h ⊂ U × P × T , where h denotes the discretization parameter. We denote by ( u h , p h , θ h ) (respectively ( v h , q h , ζ h )) the discrete FE solutionapproximating the continuous solution ( u , p , θ ) (respectively ( v , q , ζ )). We also denote by φ h the FE approximation of φ . 7 .1. Matrix assembly for Newton’s method Hereafter, we associate to the linear operator DF ( u , p , θ )( · , · , · ) a matrix representation as follows: DF u , u h DF u , ph DF u ,θ h DF p , u h DF p , ph DF θ, u h DF θ,θ h δ u h δ p h δθ h , (13)where the blocks in this matrix are linear operators associated to the following bilinear forms as follows( DF u , u h δ u h , v h ) = (( δ u h , ∇ u h ) , v h ) + (( u h , ∇ δ u h ) , v h ) + ( π m µ nf ( θ h , φ h ) (cid:16) ∇ δ u h + ( ∇ δ u h ) t (cid:17) , ∇ v h )( DF u , ph δ p h , v h ) = ( π m ∇ δ p , v h )( DF u ,θ h δθ, v h ) = − ( π m δθ ( u h , e z ) , v h )( DF p , u h δ u h , q h ) = − (cid:16) π m δ u h , ∇ q h (cid:17) ( DF p , ph p h , q h ) = ε ( δ p h , q h )( DF θ, u h δθ h , ζ h ) = (( u h , ∇ δθ h ) , ζ h )( DF θ,θ h δθ h , ζ h ) = ( π e1 ( φ )( k nf ( θ h , φ h ) ∇ δθ h ) , ∇ ζ h ) + ( π e2 ( φ h ) D ω ( ∇ φ h , ∇ δθ h ) ζ h ) + ( π e3 ( φ )2 D θ h ( ∇ θ h , ∇ δθ h ) , ζ h ) . Newton’s iterations updates the FE solution ( u h , p h , θ h ) as follows u k + h p k + h θ k + h = u kh p kh θ kh − δ u k δ p k δθ k (14)where ( δ u kh , δ p kh , δθ kh ) is solution to DF u , u h DF u , ph DF u ,θ h DF p , u h DF p , ph DF θ, u h DF θ,θ h δ u kh δ p kh δθ kh = F uu h F u ph F u θ h F p u h F pph F θ u h F θθ h u kh p kh θ kh (15)with right-hand side (presented as a matrix-vector product) contrains blocks of linear operators associated to thefollowing bilinear forms as such( F uu h u h , v h ) = (( u h · ∇ u h ) , v h ) + (cid:16) π m µ nf ( θ h , φ h ) (cid:16) ∇ u h + ( ∇ u h ) t (cid:17) , ∇ v h (cid:17) ( F u ph p h , v h ) e = (cid:16) π m ∇ p h , v h (cid:17) ( F u θ h θ h , v h ) = − (cid:16) π m θ h ( u h , e z ) , v h (cid:17) ( F p u h u h , q h ) = − (cid:16) π m u h , ∇ q h (cid:17) ( F pph p h , q h ) = ε ( p h , q h )( F θ u h u h , ζ h ) = (( u h · ∇ θ h ) , ζ h )( F θθ h θ h , ζ h ) = (cid:16) π e1 ( φ ) ( k nf ( θ h , φ h ) ∇ θ h ) , ∇ ζ h (cid:17) + (cid:16) π e2 ( φ h ) D ω ( ∇ φ h , ∇ θ h ) , ζ h (cid:17) + (cid:16) π e3 ( φ h ) D θ h ( ∇ θ h , ∇ θ h ) , ζ h (cid:17) . We develop hereafter a numerical scheme that approximates the solution of the transport dominated advection-di ff usion nanoparticles concentration equation. Despite the fact that finite element solution is very adaptive for thiskind of governing equations, such method and other Galerkin based approaches may potentially fail where theirrelative discrete solution will be market by spurious oscillations in space. This happens if their element P´eclet number8oes beyond certain critical value. Other discretization methods su ff er, actually, from this issue, for instance the finitedi ff erence technique. In the later, the spatial oscillation could be reduced by upwinding scheme. For the finite elementmethod, we have equivalent methods to the later upwinding scheme, such as Petrov-Galerkin and Streamline-UpwindPetrov-Galerkin (SUPG) [39, 40, 26]. Here, the shape function is modified in order to mimic the upwinding e ff ect andtherefore reduce, or eliminate, the spatial oscillations. An other interesting stabilization method is the Galerkin / Least-Square (GLS) see for instance [27, 41] and references therein. All the above and other techniques [42, 43, 44] moreadaptable for the time-dependent equations, use artificial viscosity in the direction of the streamlines.Let us define the space X h = (cid:110) ϕ h ∈ C ( Ω ); ∀ Ω e ∈ T h , ϕ h | Ω e ∈ P (cid:111) ⊂ X , Over which, the bilinear form related to the variational formulation of the advection-di ff usion nanoparticle concentra-tion writes as follows: for any trial function ϕ h ∈ X h find φ h ∈ X h a ε ( φ h , ϕ h ) = b ( θ h ) (16)where a ε ( φ h , ϕ h ) = (cid:16)(cid:16) u · ∇ φ h (cid:17) , ϕ h (cid:17) + π p1 (cid:16) D ω ∇ φ h , ∇ ϕ h (cid:17) + ε ( φ h , ϕ h ) . ( b ( θ h ) , ϕ h ) = π p2 (cid:0) D θ h ∇ θ h , ∇ ϕ h (cid:1) We endow X with the following norm (cid:107) ϕ (cid:107) ε : = (cid:113) π p (cid:107)∇ ϕ (cid:107) + ε (cid:107) ϕ (cid:107) , which we will use in the sequel. Proposition 1. (Well-posedness) The bilinear form a ε ( · , · ) is coercive such thata ε ( ϕ h , ϕ h ) ≥ (cid:107) ϕ h (cid:107) ε , a ε ( ϕ h , ϕ h ) ≥ α p (cid:107) ϕ h (cid:107) . Proof.
We have a ε ( ϕ h , ϕ h ) = (( u h · ∇ ϕ h ) , ϕ h ) + π p1 (cid:0) D ω ∇ ϕ h , ∇ ϕ h (cid:1) + ε ( ϕ h , ϕ h )Note that we have (( u h · ∇ ϕ h ) ϕ h ) = (cid:16) ∇ · (cid:16) u h ϕ h (cid:17) , (cid:17) = − (cid:16) ( u h · n ) ϕ h (cid:17) Γ = a ε ( ϕ h , ϕ h ) ≥ π p (cid:107)∇ ϕ h (cid:107) + ε (cid:107) ϕ h (cid:107) − (cid:16) ( u h · n ) ϕ h (cid:17) Γ ≥ π p + ε c Ω c Ω (cid:107) ϕ h (cid:107) . Where c Ω stands for the Poincar´e constant. Finally, by setting α p = π p + ε c Ω c Ω , we obtain the coercivity result in X and in L ( Ω ) as stated above.In the subsequent part, we shall discuss the stabilization of the transport equation for the concentration of nanopar-ticles. Indeed, we use the so-called Streamline Upwind Petrov-Galerkin (SUPG) method that enlarges the standardtrial space via streamline trials. In practice, the SUPG formulates as follows. Given u h ∈ U h , for a test function ϕ inspan ϕ h + (cid:88) Ω e δ Ω e u h ∇ ϕ h find φ h such that a ε supg ( φ h , ϕ h ) = b supg ( φ h , ϕ h ) . (17)9here a ε supg ( φ h , ϕ h ) : = a ε ( φ h , ϕ h ) + nel (cid:88) e (cid:90) Ω e δ Ω e (cid:16) u h · ∇ φ h + εφ h − π p1 ∇ · (cid:0) D ω ∇ φ h (cid:1)(cid:17) ( u ∇ ϕ h ) δω and b supg ( φ h , ϕ ) = b ( φ h , ϕ ) + nel (cid:88) e (cid:90) Ω e δ Ω e π p2 ∇ · (cid:0) D θ h ∇ θ h (cid:1) ( u ∇ ϕ h ) δω. Here, the integral contributions stand for the SUPG formulation [45, 46, 47, 48, 49], where, δ Ω e are user-chosenweights element dependent parameters. The above convection dominated nanoparticle equation, even if it is linearwith respect to the unknown φ h sill exhibit spurious numerical oscillations due to the fact that the dominant con-tribution comes from the convective term, which is far away larger than the thermophoresis e ff ect and the di ff usiveBrownian motion. This directly reflects that the corresponding P´eclet number Pe (non-dimensional quantity), i.e.,the ratio of the convection rate and the di ff usion rate is relatively high. To overcome this handicap in the numericalsimulation, we use the SUPG method, which provides numerical stability. Indeed, the SUPG adds an artificial stream-di ff usion along the streamlines of the particles flow. In practice the SUPG weight function parameter is adjusted interm of local P´eclet number. In theory, the choice of the weighted parameter function must satisfy δ Ω e ≤ min ε , h Ω e (cid:107) D ω (cid:107) ∞ , Ω e C inv (18)a condition for which we ensure the coercivity and hence the well-posedness of the SUPG problem via the nextcoercivity result in the Banach space V h : V h = (cid:110) ϕ h ∈ X h ( Ω ) | u · ∇ ϕ h ∈ L ( Ω ) (cid:111) endowed with the SUPG-norm defined by: (cid:107) ϕ h (cid:107) supg : = (cid:113) π p1 (cid:107)∇ ϕ h (cid:107) + ε (cid:107) ϕ h (cid:107) + (cid:107) δ Ω u · ∇ ϕ h (cid:107) (19)As our numerical approach is based on Newton’s method, which increasingly varies the Rayleigh number Ra, ourSUPG weighted function is made as function of the local P´eclet number Pe np , Ω e , the nanofluid Rayleigh number Ra nf and the local element size h Ω e . Hence the numerically chosen weighted function is δ Ω e = h Ω e Ra nf Pe np , Ω e . Theorem 1.
The bilinear form a ε supg ( · , · ) satisfies the following Lax-Milgram conditionsCoercivity a ε supg ( ϕ h , ϕ h ) ≥ (cid:107) ϕ h (cid:107) supg (20) Continuity a ε supg ( ϕ h , ψ ) ≤ C (cid:107) ϕ h (cid:107) supg (cid:107) ψ (cid:107) supg (21) for which the problem (17) is well-posed and has unique solution. roof. The proof follows the classical bounds estimates for SUPG analysis. a ε supg ( ϕ h , ϕ h ) ≥ (( u · ∇ ϕ h ) , ϕ h ) + π p1 ( D ω ∇ ϕ h , ∇ ϕ h ) + ε ( ϕ h , ϕ h ) + nel (cid:88) e (cid:90) Ω e δ Ω e ( u · ∇ ϕ h ) ( u ∇ ϕ h ) δω − nel (cid:88) e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19) δω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (22) ≥ π p1 (cid:107)∇ ϕ h (cid:107) + ε (cid:107) ϕ h (cid:107) + (cid:107) δ Ω u · ∇ ϕ h (cid:107) − nel (cid:88) e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19) δω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (23) ≥ (cid:107) ϕ h (cid:107) supg − nel (cid:88) e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19) δω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (24)where we have use the fact that D ω > (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω e δ Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19) δω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) Ω e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δ Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δω ≤ (cid:18) εδ Ω e (cid:107) ϕ h (cid:107) , Ω e + π p1 δ Ω e (cid:107) D ω (cid:107) ∞ , Ω e (cid:107) ∆ ϕ h (cid:107) , Ω e (cid:19) (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (25) ≤ εδ Ω e (cid:107) ϕ h (cid:107) , Ω e + δ Ω e π p1 (cid:107) D ω (cid:107) ∞ , Ω e C inv h Ω e (cid:107)∇ ϕ h (cid:107) , Ω e (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (26) ≤ (cid:114) ε (cid:107) ϕ h (cid:107) , Ω e + (cid:115) π p1 (cid:107)∇ ϕ h (cid:107) , Ω e (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (27) ≤ ε ξ (cid:107) ϕ h (cid:107) , Ω e + π p1 (cid:107) D ω (cid:107) ∞ , Ω e ξ (cid:107)∇ ϕ h (cid:107) , Ω e + ξ (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (28)Having used Cauchy-Schwarz for the inequality (25), then the inverse inequality (see [50]) (cid:107) ∆ ϕ h (cid:107) , Ω e ≤ C inv h Ω e (cid:107)∇ ϕ h (cid:107) , Ω e in the argument for (26), we also refer to [51] and [52] for further details on inverse inequalities. Then we used (18)in (27), and young’s product inequality in (26). Thus by chosing ξ = we obtain (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω e (cid:16) εϕ h − π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ϕ h (cid:19) δω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ε (cid:107) ϕ h (cid:107) , Ω e + π p1 (cid:107) D ω (cid:107) ∞ , Ω e (cid:107)∇ ϕ h (cid:107) , Ω e + (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (29)which we combine together with (24) to obtain a ε supg ( ϕ h , ϕ h ) ≥ (cid:107) ϕ h (cid:107) supg − nel (cid:88) e (cid:18) ε (cid:107) ϕ h (cid:107) , Ω e + (cid:16) π p1 (cid:107) D ω (cid:107) ∞ , Ω e (cid:17) (cid:107)∇ ϕ h (cid:107) , Ω e + (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e (cid:19) = (cid:107) ϕ h (cid:107) supg (30)11esides, for the boundedness of the binilinear form we have a ε supg ( ϕ h , ψ h ) = (( u · ∇ ϕ h ) , ψ h ) + π p1 ( D ω ∇ ϕ h , ∇ ψ h ) + ε ( ϕ h , ψ h ) + nel (cid:88) e (cid:90) Ω e δ Ω e ( u · ∇ ϕ h ) ( u ∇ ψ h ) δω + nel (cid:88) e (cid:90) Ω e (cid:16) εϕ h ψ h + π p1 ∇ · ( D ω ∇ ϕ h ) (cid:17) (cid:18) δ Ω e u ∇ ψ h (cid:19) δω ≤ (cid:16) c Ω (cid:107) u (cid:107) + π p1 (cid:107) D ω (cid:107) ∞ + ε c Ω (cid:17) (cid:107)∇ ϕ h (cid:107) (cid:107)∇ ψ h (cid:107) + nel (cid:88) e (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) δ Ω e u ∇ ϕ h (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , Ω e (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) δ Ω e u ∇ ψ h (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , Ω e + ε c Ω (cid:107)∇ ϕ h (cid:107) , Ω e (cid:107)∇ ψ h (cid:107) , Ω e + π p1 (cid:107) D ω (cid:107) ∞ , Ω e (cid:107)∇ ϕ h (cid:107) , Ω e (cid:107)∇ ψ h (cid:107) , Ω e + (cid:107) δ Ω e u ∇ ϕ h (cid:107) , Ω e ≤ C (cid:107)∇ ϕ h (cid:107) (cid:107)∇ ψ h (cid:107) (31)where C is function of C ( c Ω , ε, (cid:107) D ω (cid:107) ∞ , (cid:107) u (cid:107) ). This completes the proof.Note here that the SUPG method demonstrates a stronger stability property in the streamline direction than thestandard Galerkin discretization.Furthermore, the nanoparticle concentration mean value has to be equal to one (as per the dimensionless formu-lation of the equations). This reflects the conservation of the bulk amount of concentration in the enclosure. Thisnanoparticle constraint hence writes as (cid:90) Ω φ ( ω ) δω = . (32)In this light, the handled problem is treated as a variational minimization problem subject to a constraint in the statevariable φ h and writes as follows min φ h ∈ H ( Ω ) , (cid:82) Ω φ h ( ω ) δω = J ( φ h , λ ) : = a ε ( φ h , φ h ) + λ (cid:90) Ω φ h ( ω ) δω where the constraint (32) is token into consideration using the Lagrange multiplier λ in the above augmented costfunctional. One can easily retrieve the equations that need to be solved for the nanoparticle concentrations as a criticalpoint of J ( · , · ), i.e., by deriving with respect to φ h and λ .In practice, we assemble the finite element matrix system for the nanoparticle transport equation as follows: N zz T (cid:34) φ h λ (cid:35) = (cid:34) b (cid:35) (33)12here ( N ) i , j = (cid:90) Ω (cid:16) u · ∇ ϕ ih (cid:17) ϕ jh δω + (cid:90) Ω π p1 D w ∇ ϕ ih · ∇ ϕ jh δω + (cid:90) Ω εϕ ih ϕ jh δω + nel (cid:88) e (cid:90) Ω e δ Ω e (cid:16) u · ∇ ϕ ih (cid:17) (cid:16) u ∇ ϕ jh (cid:17) δω − nel (cid:88) e (cid:90) Ω e δ Ω e π p1 ∇ · (cid:16) D ω ∇ ϕ ih (cid:17) (cid:16) u ∇ ϕ jh (cid:17) δω ( z ) j = (cid:90) Ω · ϕ jh δω ( b ) j = − (cid:90) Ω π p2 D θ (cid:63) H ∇ θ h · ∇ ϕ jh δω + nel (cid:88) e (cid:90) Ω e δ Ω e π p2 ∇ · (cid:0) D θ h ∇ θ h (cid:1) (cid:16) u ∇ ϕ jh (cid:17) δω. One may ask whether this augmented matrix (33) is invertible or not. Indeed, it is invertible as per the followingproposition
Proposition 2.
The discrete linear system (33) describing the nanoparticle transport in the enclosure (under the meanvalue constraint) is consistent.Proof.
It is clear that the finite element square matrix N : = [ c | · · · | c n ] ∈ R n × n is non-singular as per Proposition 1,which means that 0 is not an eigenvalue of N . Therefore, its set of column vectors c i , i = ..., n form a basis for thecolumn space Col( N ). For any non zero vector z ∈ R n with positive entries, there exists a sequence of coe ffi cients ( α ) j such that we have z = (cid:80) nj = α j c j where α = N − z is clearly unique.Let us assume that the column vector ˜ z : = (cid:34) z (cid:35) is also linear combination of columns of the augmented matrix (cid:34) N z T (cid:35) .So ˜ z = (cid:80) j α j ˜ c j (where ˜ c j = [ c Tj , z j ] T , in particular, if we take the bottom line of this linear combination we have (cid:80) j α j z j = z T α = z T N − z = , ∀ z ∈ R n , (cid:107) z (cid:107) (cid:44)
0. This contradicts N is non-singular. We present our method that combines Newton’s method for the resolution of nonlinear PDEs, combining theMomentum and Energy equations, together with the nanoparticle concentration advection dominated at steady state.The physical parameters dependency in the nanoparticle concentration di ff ers from one correlation to another, seefor instance [21, 38, 53] without being exhaustive. In most cases a highly nonlinear term involving the nanoparticleconcentration appear in the correlation (generally coming from non-linear fitting procedure). In order to avoid anydi ff erentiation in the Newton’s method with respect to the nanoparticle concentration we will split the resolutionmethod into two and update all variables accordingly in iterative fashion as demonstrated in following Algorithm 1.13 lgorithm 1: Combined Newton’s method and SUPG for nanofluid equation
Input: tol , u , p , θ for Ra n f in (10 · · · ) do (cid:15) = k = φ h = ( (cid:82) Ω δω ) − ; Initialization of the nanoparticle concentration ; while (cid:15) ≥ tol do /* Using φ k − h */ Solve for (cid:104) δ u k , δ p k , δθ k (cid:105) T following Eq. (15); (cid:15) m , e ← [ δ u kh , δ p kh , δθ k ] T [ δ u kh , δ p kh , δθ k ][ u kh , p kh , θ k ] T [ u kh , p kh , θ kh ] ; Update (cid:104) u k + h , p k + h , θ k + h (cid:105) T following Eq. (14); Solve for φ k + h under the mean constraint Eq. (33); /* Nanoparticle concentration with the SUPG stabilization */ k ← k + end k ← end As it is showing in Algorithm 1 a split procedure is adopted, where the tangent equation only concern the momen-tum and energy equations leading to an update, through Newton’s method, of the velocity and temperature respec-tively. Although, these later two equations are dependent upon the nanoparticle concentration through the mixturefluid density, viscosity and thermal conductivity. The proposed method assumes that the velocity and the temperatureare constant through one iteration of Newton’s method. Their update will then be ensured within the next iterationafter an exact resolution (through LU decomposition) of the nanoparticle concentration (with SUPG) based on the newvariables of the velocity and temperature coming out of the previous Newton’s iteration. This alternating combinationis applied throughout the iterations and leads to convergence for all variables involved. The split approach has leadto a simple yet e ff ective implementation of Newton’s method involving highly non-linear parameters. Therefore, itskipped the tedious calculation of the tangent equations of the whole coupled system of four PDEs. Furthermore, lessmemory storage is then deployed hence rapid calculation. The above procedure is thus repeated with a predefined setof increasing Rayleigh numbers.
5. Numerical experiments and validations
In this section, we investigate the numerical treatment of the heat transfer enhancement in a di ff erentially heatedenclosure using variable thermal conductivity and variable viscosity of the alumina-water nanofluid (Al O -water).The validation of numerical scheme is done through two processes. The first focuses on the validation of the numericalscheme relative to the resolution of the momentum and energy equations regardless of the volume fraction. In thiscase we consider the pure water heat transfer calculation in a square cavity. The second considers the comparison ofthe present numerical results with available numerical and experimental data. We present in Figure 2 the mesh sensitivity results for the base fluid (only the Newton’s solver).Figure3 shows the quadratic convergence of the proposed numerical scheme with the use of the SUPG stabilizationtechnique. In addition, as we shall explain in the sequel, we enforced the bulk of the nanoparticle concentration to beof mean value equal to 1, through minimization under constraint problem. A continuous P igure 2: Mesh sensitivity; Results present Nusselt (base fluid) number plotted along the Heated wall of the cavity.Figure 3: Order of Convergence, with respect to the mesh size, of the presented scheme for the nanoparticle concentration with SUPG stabilization.
15n Figure 4 we present the benefit of the implementation of the SUPG method in order to eliminate numericalartifacts showing up in the concentration of the nanoparticles, which is governed by an advection dominated problem.It is clearly shown in Figure 4 that in the case without SUPG stabilisation the numerical spurious are more accentuating Ra nf = +
06 Ra nf = +
07 Ra nf = + W i t h o u t S U P G W i t hS U P G Figure 4: Stabilizatoin e ff ect of the SUPG on the nanoparticle solutions by elimination of numerical artifacts. Plot of isovalues of the concentrationin the case of φ = for high Rayleigh simulations, where indeed, high Peclet number takes place. Our numerical scheme is validated using detailed comparison with the experimental data of Ho et al. in [21] usingAl O which their thermophysical properties are reported in Table 1. In their experimental investigations, the authorsstudied the heat transfer characteristics of alumina-water nanofluid enclosed in square cells. They used three di ff erentcell geometries and di ff erent heating conditions ( θ (cid:63) H − θ (cid:63) C ) to increase the Rayleigh number. Cases with 0 , O c p (Jkg − K − ) 4179 765 ρ (kg m − ) 997.1 3970 k (W m − K − ) 0.613 25 d p (nm) 0.384 47 α · − (m s − ) 1.47 82.23 β · − (K − ) 21 0.85 Table 1: Physical properties of base fluid and Al O nanoparticles igure 5: Comparison of the numerical experiments with the experimental results of Ho et al. [21] for the first cell case Our validation and comparison of the numerical results against experimental finding of nanofluid heat transferfollowed the cells presentation as in [21]. The corresponding results are reported in Figure 5,6 and 7. In each of thesefigures we plot separately the cases of base fluid (left) the nanofluid concentration of 1% (middle) and the nanofluidconcentration of 3% (left). To produce these results, we have used viscosity and thermal conductivity correlations asreported by Ho et al. in [21] as follows. C µ ( φ ) = + . · φ + . · φ C k ( φ ) = (1 + . · φ + . · φ ) , standing for the viscosity and thermal conductivity respectively. Figure 6: Comparison of the numerical experiment with the experimental results of Ho et al. [21] for the second cell case
Based on the averaged Nusselt number values of the nanofluid, the present predictions show a reasonably goodagreement with the experimental data for low and moderate Ra (Figures 5 and 6). Whereas, for high Ra and highconcentration, the numerical calculations are found to overestimate the Nusselt number. In fact, for the high Ra nf number cases (Figure 7), the experimental results show a more pronounced Nusselt number deterioration for thecases with nanofluid in comparison to the one with the base fluid. In their paper, Ho et al. [21] suspected thatthis behavior could be attributed to the transport mechanisms associated with nanoparticle-fluid interactions suchas Brownian di ff usion and thermophoresis in addition to the impact of the thermophysical properties changes. Thepresent numerical predictions, however, indicate that the Buongiorno model which is intended to specifically accountfor these two mechanisms, is not actually able to mimic the equivalent heat transfer impairment. This might be17 igure 7: Comparison of the numerical experiment with the experimental results of Ho et al. [21] for the third cell case suggesting that perhaps additional forces should be incorporated in the Buongiorno transport equations to provide abetter physical model for the nanoparticles concentration, capable of reflecting the e ff ect of the nanoparticles on theheat transfer impairment observed in the experimental data. This subsection is devoted to the comparison of the numerical results obtained by Algorithm 1 based on FEMdiscretization of the Buongiorno nanofluid transport model, against results of [57] that are based on finite volumediscretization using Fluent [58] (commercial software). Here it is worth noticing that both methods deal with samephysics of nanofluid transport, although, use di ff erent equations models. Indeed, the aforementioned results are basedon a solid-liquid mixture model which solves the momentum equations with an additional term to account for thephases drift velocity, the continuity equation, and energy equation for the mixture. The model adopts algebraicexpressions for the relative velocities which are then used to define the drift velocities (see Fluent’s documentation formore details [58]). Both numerical results are compared to experimental results of Ho et al. [21]. Figure 8: Comparison of the numerical experiments with the numerical experiment of Chen et al. [57] for the range of Ra in base fluid formulation. Figure 8 depicts both numerical results and showcase a good agreement between the two methods. The range ofdata available for this validation and comparison is 1 · to 6 · . For the case of 1% of the nanoparticle bulkconcentration, Buongiorno model seems to have better prediction of the heat transfer in term of the Nussult number ofthe nanofluid. However, this advantage becomes marginally on the side of the multi-phase model while we increasethe bulk of nanoparticle concentration to 3%. 18 a n f = . E R a n f = . E R a n f = . E R a n f = . E R a n f = . E Figure 9: (From left to right) velocity streamlines, heat distribution and nanoparticle concentration distribution. Plots, show (from top to bottom) thee ff ect of increasing the Rayleigh number on the profile of the variables listed earlier. Numerical simulation was performed through finite elementdiscretization, using Lagrange polynomials of degree two for the velocity and of degree 1 for the temperature and nanoparticle concentrationrespectively. nf . These results show that the stream-lines exhibit recirculations that getflatten with the increase of the Rayleigh number. These recirculations get localized in the middle of the cavity and nearthe hot and cold walls. Whereas, isothermal lines become more and more horizontal which lead to high temperaturegradient near the hot and cold walls. One, here, can directly see the increase of the heat transfer with the increase of theRayleigh number. Besides, both the energy and the nanoparticle concentration are advected mainly by the buoyancydriven flow. Although, the energy equation, has a considerable contribution of the thermal di ff usion compared to thenanoparticle transport equation, which has much smaller di ff usion terms. Indeed, the Brownian di ff usivity is veryweak in comparison to the thermal di ff usion and even negligible when compared to the advection term. This, infact, would drastically a ff ect the numerical approximation of this advection dominated equation. The SUPG artificialviscosity through the streamlines is, therefore, needed to stabilize the numerical scheme for high Rayleigh / Pecletnumber as shown in Figure 4. As the nanoparticle transport equation is mainly driven by the advection, one wouldexpect the distribution of the concentration to be very similar to the stream-line of the flow. Indeed, this behaviour isobserved in the third column of the plot in Figure 9.
Figure 10: Nanoparticle concentration profile along the ( x , ) varying from %0 to %3. Figure 10 displays the nanoparticles concentration profile along the line ( x , .
5) for di ff erent Rayleigh numbersand averaged concentration φ values ranging from 1% to 3%. Although, the profiles exhibit only minor variation, ofthe order of 10 − in magnitude along the horizontal line, very interesting phenomena occurring in the vicinity of thehot and cold walls can be observed. Near the cold wall, for instance, the concentration profiles decrease sharply asthe nanoparticles get closer to the cold wall. The slop of this decrease is a function of both; the averaged nanoparti-cles concentration value and the Rayleigh number. If one is to recall the thermophoresis e ff ect which tends to moveparticles from hotter to colder zones, then even in the near-the-wall zone in which only small convective velocity mag-20itude exist, this e ff ect is not dominant and the nanoparticles are still being carried away by the recirculating motion ofthe carrier fluid. On the hot wall, however, a di ff erent picture is depicted. Although it would not be a straight forwardtask to pin point the dominant term in this zone, the resulting force seems to favor a higher nanoparticles concentrationat the wall vicinity followed by a sharper increase which can be translated to the fact that the convective term regainsits dominant e ff ect further away from the wall. All these nanoparticles that are removed from the walls region, by oneof the two mechanisms described above, get accumulated in the center resulting in a relatively higher nanoparticlesconcentration. Unfortunately, as stated above, the concentration variation along this line remains marginal, hence,one would not expect it to provide a dramatically di ff erent outcomes from running the simulations with a constantconcentration value, hence, assuming a single-phase model approach.
6. Conclusion
We presented in this article a numerical technique based on Newton-Raphson iterations to solve the nanofluid heattransfer problem in a square cavity with variable properties. In addition to its generality (regardless of the correlationused for the variable properties), our technique has mainly two advantages compared to the conventional use ofNewton’s iterations: • Firstly, it avoids the di ffi culty coming from the highly non-linear dependency upon the concentration in sev-eral correlations published in the literature. Indeed, the Jacobian (tangent problem) disregards nanoparticleconcentration and only considers the velocity, pressure, and temperature, while nanoparticle concentration getsupdated iteratively. Admittedly, the momentum and energy equations are solved through Newton’s iterations, asthe dominant (Navier-Stokes) equations is quadratic for the velocity variable, while the nanoparticle transportequation gets solved right after each iteration of the momentum and energy equations. • Secondly, the proposed split leads to less memory consumption and allows the viscosity, density, and thermalconductivity of recirculating flow to be updated at each Newton’s iteration.The numerical experiments based on the Finite Element discretization of the nanofluid heat transfer problem havebeen regularized using the SUPG method, which showed to be very e ff ective in stabilizing the numerical solution bywiping the spurious oscillations without wrecking the solution. Here in particular we found that the ratio formulabetween the local Peclet number and global Raleigh number is a good combination for the regularization functionused in the SUPG. Besides, our numerical scheme has been thoroughly validated against experimental results andshowed a good agreement over a large spectrum of Rayleigh numbers ranging from 10 to 10 . The present studyalso reveals that although the Buongiorno’s four equations-based nanofluid transport model, tested herein, is able tocapture additional physical phenomena a ff ecting the nanoparticles distribution, additional forces might have to beaccounted for within this model in order to mimic heat transfer deterioration of similar amplitude as it is observed inthe experimental data. Appendices
A. Density dimensionless derivation
Following Eq.(2) we have ρ nf = (1 − φ (cid:63) ) ρ bf + φ (cid:63) ρ np note also from Eq. (4), ( ρβ ) nf = ( ρβ ) bf (1 − φ (cid:63) ) + ( ρβ ) np φ (cid:63) . ρβ ) bf = ρ bf β bf and ( ρβ ) np = ρ np β np . Hence( ρβ ) nf ρ nf = ρ bf β bf (1 − φ (cid:63) )(1 − φ (cid:63) ) ρ bf + φ (cid:63) ρ np + ρ np β np φ (cid:63) (1 − φ (cid:63) ) ρ bf + φ (cid:63) ρ np = β bf (1 − φ (cid:63) )(1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf + φ (cid:63) (1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf β np = β bf (1 − φ (cid:63) )(1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf + φ (cid:63) (1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf β np β bf Let M : = (1 − φ (cid:63) )(1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf + φ (cid:63) (1 − φ (cid:63) ) + φ (cid:63) ρ np ρ bf β np β bf B. Momentum equation (cid:16) u (cid:63) · ∇ (cid:63) (cid:17) u (cid:63) = − ρ nf ∇ (cid:63) p (cid:63) + ρ nf ∇ (cid:63) · (cid:16) µ (cid:63) nf (cid:16) ∇ (cid:63) u (cid:63) + ( ∇ (cid:63) u (cid:63) ) t (cid:17)(cid:17) + g ρ nf ( ρ ∞ − ρ c ) (34)Moving toward dimensionless variables the above equation writes α L ( u · ∇ ) u = − ρ bf α L ρ nf ∇ (cid:63) p (cid:63) + αµ bf L ρ nf ∇ (cid:63) · (cid:16) µ nf (cid:16) ∇ (cid:63) u (cid:63) + ( ∇ (cid:63) u (cid:63) ) t (cid:17)(cid:17) + g β nf ρ nf ρ ∞ ( θ h − θ c ) θ (35)Multiplying the above equation by L α we obtain( u · ∇ ) u = ρ bf ρ nf ∇ p + µ bf αρ bf (cid:32) − φ + φ ρ np ρ bf (cid:33) ∇ · (cid:16) µ nf (cid:16) ∇ u + ( ∇ u ) t (cid:17)(cid:17) + M L g β nf α ρ nf ρ ∞ ( θ h − θ c ) θ which rewrites using the non-dimensional constants as follows( u · ∇ ) u = π m ( φ ) ∇ p + π m ( φ ) ∇ · (cid:16) µ nf (cid:16) ∇ u + ( ∇ u ) t (cid:17)(cid:17) + π m ( φ ) θ, (36)where π m ( φ ) = (cid:32) − φ + φ ρ np ρ bf (cid:33) − , (37) π m ( φ ) = Pr (cid:32) − φ + φ ρ np ρ bf (cid:33) − , (38) π m ( φ ) = PrRa nf M . (39)22 . Energy equation The dimensional energy equation writes (cid:16) u (cid:63) · ∇ (cid:63) θ (cid:63) (cid:17) = c nf ρ nf ∇ (cid:63) · (cid:16) k (cid:63) nf ∇ (cid:63) θ (cid:63) (cid:17) + (cid:32) ρ np c np ρ nf c nf (cid:33) (cid:32) D (cid:63)θ θ (cid:63) C ∇ (cid:63) θ (cid:63) · ∇ (cid:63) θ (cid:63) + D (cid:63)ω ∇ (cid:63) φ · ∇ (cid:63) θ (cid:63) (cid:33) moving toward dimensionless variables the above equation writes α ( θ (cid:63) H − θ (cid:63) C ) L ( u · ∇ θ ) = k bf ( θ (cid:63) H − θ (cid:63) C ) L c nf ρ nf ∇ · ( k nf ∇ θ ) + (cid:32) ρ np c np ρ nf c nf (cid:33) D θ c ( θ (cid:63) H − θ (cid:63) C ) θ (cid:63) C L ( D θ ∇ θ · ∇ θ ) + (cid:32) ρ np c np ρ nf c nf (cid:33) φ b D ω c ( θ (cid:63) H − θ (cid:63) C ) L ( D ω ∇ φ · ∇ θ ) . Multiplying the above equation by L α ( θ (cid:63) H − θ (cid:63) C ) we obtain( u · ∇ θ ) = k bf α c nf ρ nf ∇ · ( k nf ∇ θ ) + (cid:32) ρ np c np ρ nf c nf (cid:33) D θ c ( θ (cid:63) H − θ (cid:63) C ) αθ (cid:63) C ( D θ ∇ θ · ∇ θ ) + (cid:32) ρ np c np ρ nf c nf (cid:33) φ b D ω c α ( D ω ∇ φ · ∇ θ ) , which rewrites using non-dimensional variables as follows( u · ∇ θ ) = π e ∇ · ( k nf ∇ θ ) + π e ( D θ ∇ θ · ∇ θ ) + π e ( D ω ∇ φ · ∇ θ )Where π e = k bf α c nf ρ nf = (cid:32) φ + (1 − φ ) ρ bf c bf ρ np c np (cid:33) − π e = (cid:32) ρ np c np ρ nf c nf (cid:33) φ b D ω c α = StPrSc θ (cid:63) H − θ (cid:63) C θ (cid:63) H (cid:32) φ + (1 − φ ) ρ bf c bf ρ np c np (cid:33) − π e = (cid:32) ρ np c np ρ nf c nf (cid:33) D θ c ( θ (cid:63) H − θ (cid:63) C ) αθ (cid:63) C = PrSc φ b (cid:32) φ + (1 − φ ) ρ bf c bf ρ np c np (cid:33) − . D. Nanoparticle transport equation
The particle dimensional equation writes ∇ (cid:63) · ∇ (cid:63) φ (cid:63) = ∇ (cid:63) · (cid:32) D (cid:63)ω ∇ (cid:63) φ (cid:63) + D (cid:63)θ (cid:63) θ (cid:63) C ∇ (cid:63) θ (cid:63) (cid:33) using the non-dimensional equations the above equation writes α L φ b ∇ · ∇ φ = φ b D ω c L ∇ · ( D ω ∇ φ ) + D θ (cid:63) C ( θ (cid:63) H − θ (cid:63) C ) L θ (cid:63) C ∇ · ( D θ ∇ θ ) . Multiplying the above by L αφ b we obtain 23 b ∇ · ∇ φ = D ω c α ∇ · ( D ω ∇ φ ) + D (cid:63)θ (cid:63) C ( θ (cid:63) H − θ (cid:63) C ) φ b αθ (cid:63) C ∇ · (cid:16) D θ ∇ θ (cid:63) (cid:17) . ∇ · ∇ φ = π p ∇ · ( D ω ∇ φ ) + π p ∇ · (cid:16) D θ ∇ θ (cid:63) (cid:17) . where π p = D ω c α = ,π p = D θ (cid:63) C ( θ (cid:63) H − θ (cid:63) C ) L θ (cid:63) C = StPrSc θ (cid:63) H − θ (cid:63) C θ (cid:63) C φ b . References [1] W. Minkowycz, E. M. Sparrow, J. P. Abraham, Nanoparticle heat transfer and fluid flow, volume 4, CRC press, 2012.[2] O. Manca, Y. Jaluria, D. Poulikakos, Heat transfer in nanofluids, 2010.[3] C. Kleinstreuer, Z. Xu, Mathematical modeling and computer simulations of nanofluid flow with applications to cooling and lubrication,Fluids 1 (2016) 16.[4] C. Kleinstreuer, Microfluidics and nanofluidics: theory and selected applications, John Wiley & Sons, 2013.[5] J. Buongiorno, Convective Transport in Nanofluids, Journal of Heat Transfer 128 (2005) 240–250. URL: https://doi.org/10.1115/1.2150834 . doi: .[6] M. A. Sheremet, I. Pop, O. Mahian, Natural convection in an inclined cavity with time-periodic temperature boundary conditions usingnanofluids: application in solar collectors, International Journal of Heat and Mass Transfer 116 (2018) 751–761.[7] O. Mahian, A. Kianifar, S. A. Kalogirou, I. Pop, S. Wongwises, A review of the applications of nanofluids in solar energy, InternationalJournal of Heat and Mass Transfer 57 (2013) 582–594.[8] J. Li, C. Kleinstreuer, Thermal performance of nanofluid flow in microchannels, International Journal of Heat and Fluid Flow 29 (2008)1221–1232.[9] A. Ba¨ıri, E ff ects of zno-h2o nanofluid saturated porous medium on the thermal behavior of cubical electronics contained in a tilted hemi-spherical cavity. an experimental and numerical study, Applied Thermal Engineering 138 (2018) 924–933.[10] Q. Li, J. Wang, J. Wang, J. Baleta, C. Min, B. Sund´en, E ff ects of gravity and variable thermal properties on nanofluid convective heat transferusing connected and unconnected walls, Energy conversion and management 171 (2018) 1440–1448.[11] Z. Xu, C. Kleinstreuer, Computational analysis of nanofluid cooling of high concentration photovoltaic cells, Journal of Thermal Scienceand Engineering Applications 6 (2014).[12] A. Ba¨ıri, N. Laraqi, K. Adeyeye, Thermal behavior of an active electronic dome contained in a tilted hemispherical enclosure and subjectedto nanofluidic cu-water free convection, The European Physical Journal Plus 133 (2018) 1–11.[13] F. Jabbari, A. Rajabpour, S. Saedodin, Thermal conductivity and viscosity of nanofluids: a review of recent molecular dynamics studies,Chemical Engineering Science 174 (2017) 67–81.[14] H. Khodadadi, S. Aghakhani, H. Majd, R. Kalbasi, S. Wongwises, M. Afrand, A comprehensive review on rheological behavior of mono andhybrid nanofluids: e ff ective parameters and predictive correlations, International Journal of Heat and Mass Transfer 127 (2018) 997–1012.[15] J. Fan, L. Wang, Review of heat conduction in nanofluids, Journal of heat transfer 133 (2011).[16] S. Kakac¸, A. Pramuanjaroenkij, Review of convective heat transfer enhancement with nanofluids, International journal of heat and masstransfer 52 (2009) 3187–3196.[17] J. Buongiorno, D. C. Venerus, N. Prabhat, T. McKrell, J. Townsend, R. Christianson, Y. V. Tolmachev, P. Keblinski, L.-w. Hu, J. L. Alvarado,et al., A benchmark study on the thermal conductivity of nanofluids, Journal of Applied Physics 106 (2009) 094312.[18] M. Sheikholeslami, D. Ganji, Nanofluid convective heat transfer using semi analytical and numerical approaches: a review, Journal of theTaiwan Institute of Chemical Engineers 65 (2016) 43–77.[19] D. Wen, Y. Ding, Experimental investigation into convective heat transfer of nanofluids at the entrance region under laminar flow conditions,International journal of heat and mass transfer 47 (2004) 5181–5188.[20] C. H. Li, G. Peterson, Experimental studies of natural convection heat transfer of al2o3 / di water nanoparticle suspensions (nanofluids),Advances in Mechanical engineering 2 (2010) 742739.[21] C. Ho, W. Liu, Y. Chang, C. Lin, Natural convection heat transfer of alumina-water nanofluid in vertical square enclosures: An experimentalstudy, International Journal of Thermal Sciences 49 (2010) 1345?1353. doi: .[22] N. Putra, W. Roetzel, S. K. Das, Natural convection of nano-fluids, Heat and mass transfer 39 (2003) 775–784.[23] C. H. Chon, K. D. Kihm, S. P. Lee, S. U. Choi, Empirical correlation finding the role of temperature and particle size for nanofluid (al 2 o 3)thermal conductivity enhancement, Applied Physics Letters 87 (2005) 153107.[24] A. C. Gale˜ao, E. G. D. Do Carmo, A consistent approximate upwind petrov-galerkin method for convection-dominated problems, ComputerMethods in Applied Mechanics and Engineering 68 (1988) 83–95.[25] F. Yurun, A comparative study of the discontinuous galerkin and continuous supg finite element methods for computation of viscoelasticflows, Computer Methods in Applied Mechanics and Engineering 141 (1997) 47 – 65. URL: . doi: https://doi.org/10.1016/S0045-7825(96)01102-4 .
26] C. Erath, D. Praetorius, Optimal adaptivity for the supg finite element method, Computer Methods in Applied Mechanics and Engineering353 (2019) 308 – 327. URL: . doi: https://doi.org/10.1016/j.cma.2019.05.028 .[27] M. ten Eikelder, I. Akkerman, Correct energy evolution of stabilized formulations: The relation between vms, supg and gls via dynamicorthogonal small-scales and isogeometric analysis. ii: The incompressible navier–stokes equations, Computer Methods in Applied Mechan-ics and Engineering 340 (2018) 1135 – 1154. URL: .doi: https://doi.org/10.1016/j.cma.2018.02.030 .[28] E. B¨ansch, S. Faghih-Naini, P. Morin, Convective transport in nanofluids: The stationary problem, Journal of Mathematical Analysis andApplications 489 (2020) 124151. URL: . doi: https://doi.org/10.1016/j.jmaa.2020.124151 .[29] E. B¨ansch, A thermodynamically consistent model for convective transport in nanofluids: existence of weak solutions and fem computations,Journal of Mathematical Analysis and Applications 477 (2019) 41–59.[30] B. C. Shekar, N. Kishan, Finite element analysis of natural convective heat transfer in a porous square cavity filled with nanofluids in thepresence of thermal radiation, in: Journal of Physics: Conference Series, volume 662, IOP Publishing, 2015, p. 012017.[31] C. S. Balla, K. Naikoti, Finite element analysis of magnetohydrodynamic transient free convection flow of nanofluid over a vertical conewith thermal radiation, Proceedings of the Institution of Mechanical Engineers, Part N: Journal of Nanomaterials, Nanoengineering andNanosystems 230 (2016) 161–173.[32] N. Ullah, S. Nadeem, A. U. Khan, Finite element simulations for natural convective flow of nanofluid in a rectangular cavity having corrugatedheated rods, Journal of Thermal Analysis and Calorimetry (2020) 1–13.[33] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, volume 5, Springer Science & BusinessMedia, 2012.[34] C. Taylor, P. Hood, A numerical solution of the navier-stokes equations using the finite element technique, Computers & Fluids 1 (1973)73–100.[35] T. Apel, H. M. Randrianarivony, Stability of discretizations of the stokes problem on anisotropic meshes, Mathematics and Computers inSimulation 61 (2003) 437–447.[36] C. Ho, W. Liu, Y. Chang, C. Lin, Natural convection heat transfer of alumina-water nanofluid in vertical square enclosures: an experimentalstudy, International Journal of Thermal Sciences 49 (2010) 1345–1353.[37] E. Abu-Nada, A. J. Chamkha, E ff ect of nanofluid variable properties on natural convection in enclosures filled with a cuo–eg–water nanofluid,International Journal of Thermal Sciences 49 (2010) 2339–2352.[38] K. Khanafer, K. Vafai, A critical synthesis of thermophysical characteristics of nanofluids, Nanotechnology and Energy (2017) 279?332.doi: .[39] L. P. Franca, G. Hauke, A. Masud, Stabilized finite element methods, International Center for Numerical Methods in Engineering (CIMNE),Barcelona . . . , 2004.[40] V. John, J. Novo, A robust supg norm a posteriori error estimator for stationary convection–di ff usion equations, Computer Meth-ods in Applied Mechanics and Engineering 255 (2013) 289 – 305. URL: . doi: https://doi.org/10.1016/j.cma.2012.11.019 .[41] M. ten Eikelder, I. Akkerman, Correct energy evolution of stabilized formulations: The relation between vms, supg and gls via dy-namic orthogonal small-scales and isogeometric analysis. i: The convective–di ff usive context, Computer Methods in Applied Mechan-ics and Engineering 331 (2018) 259 – 280. URL: .doi: https://doi.org/10.1016/j.cma.2017.11.020 .[42] E. Burman, Consistent supg-method for transient transport problems: Stability and convergence, Computer Methods in Applied Mechan-ics and Engineering 199 (2010) 1114 – 1123. URL: .doi: https://doi.org/10.1016/j.cma.2009.11.023 .[43] P. B. Bochev, M. D. Gunzburger, J. N. Shadid, Stability of the supg finite element method for transient advection–di ff usion problems,Computer Methods in Applied Mechanics and Engineering 193 (2004) 2301 – 2323. URL: . doi: https://doi.org/10.1016/j.cma.2004.01.026 .[44] A. Russo, Streamline-upwind petrov / galerkin method (supg) vs residual-free bubbles (rfb), Computer Methods in Applied Mechanicsand Engineering 195 (2006) 1608 – 1620. URL: .doi: https://doi.org/10.1016/j.cma.2005.05.031 , a Tribute to Thomas J.R. Hughes on the Occasion of his 60th Birthday.[45] A. N. Brooks, T. J. Hughes, Streamline upwind / petrov-galerkin formulations for convection dominated flows with particular emphasis on theincompressible navier-stokes equations, Computer methods in applied mechanics and engineering 32 (1982) 199–259.[46] L. P. Franca, S. L. Frey, T. J. Hughes, Stabilized finite element methods: I. application to the advective-di ff usive model, Computer Methodsin Applied Mechanics and Engineering 95 (1992) 253–276.[47] T. Gelhard, G. Lube, M. A. Olshanskii, J.-H. Starcke, Stabilized finite element schemes with lbb-stable elements for incompressible flows,Journal of computational and applied mathematics 177 (2005) 243–267.[48] E. Burman, G. Smith, Analysis of the space semi-discretized supg method for transient convection–di ff usion equations, Mathematical Modelsand Methods in Applied Sciences 21 (2011) 2049–2068.[49] E. Burman, Consistent supg-method for transient transport problems: Stability and convergence, Computer Methods in Applied Mechanicsand Engineering 199 (2010) 1114–1123.[50] V. JOHN, J. NOVO, Error analysis of the supg finite element discretization of evolutionary convection-di ff usion-reaction equations, SIAMJournal on Numerical Analysis 49 (2011) 1149–1176. URL: .[51] P. G. Ciarlet, The finite element method for elliptic problems, SIAM, 2002.[52] A. Ern, J.-L. Guermond, Theory and practice of finite elements, volume 159, Springer Science & Business Media, 2013.[53] M. S. Astanina, M. Kamel Riahi, E. Abu-Nada, M. A. Sheremet, Magnetohydrodynamic in partially heated square cavity with variableproperties: Discrepancy in experimental and theoretical conductivity correlations, International Journal of Heat and Mass Transfer 116 . doi: https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.050 .[54] M. Benedetto, S. Berrone, A. Borio, S. Pieraccini, S. Scial`o, Order preserving supg stabilization for the virtual element formulationof advection–di ff usion problems, Computer Methods in Applied Mechanics and Engineering 311 (2016) 18 – 40. URL: . doi: https://doi.org/10.1016/j.cma.2016.07.043 .[55] C. Wervaecke, H. Beaugendre, B. Nkonga, A fully coupled rans spalart–allmaras supg formulation for turbulent compressible flows onstretched-unstructured grids, Computer Methods in Applied Mechanics and Engineering 233-236 (2012) 109 – 122. URL: . doi: https://doi.org/10.1016/j.cma.2012.04.003 .[56] S. Alosious, S. Sarath, A. R. Nair, K. Krishnakumar, Experimental and numerical study on heat transfer enhancement of flat tube radiatorusing al 2 o 3 and cuo nanofluids, Heat and Mass Transfer 53 (2017) 3545–3563.[57] Y.-J. Chen, P.-Y. Wang, Z.-H. Liu, Numerical study of natural convection characteristics of nanofluids in an enclosure using multiphasemodel, Heat and Mass Transfer 52 (2016) 2471–2484. doi: .[58] ANSYS, Ansys fluent - cfd software — ansys, 2016. URL: .doi: b97b60a697227d0a7d5a660b242f281f ..