Brahim Ben Beya

Nader Ben Cheikh

Taieb Lili

Comptes Rendus Mécanique (CRM)

C. R. Mecanique 335 (2007) 113–118http://france.elsevier.com/direct/CRAS2B/Natural convection flow in a tall enclosure using a multigrid methodNader Ben Cheikh ∗, Brahim Ben Beya, Taieb LiliFaculté des Sciences de Tunis, Département de Physique, Campus Universitaire 2092 El-Manar II, TunisiaReceived 30 April 2006; accepted after revision 23 January 2007Available online 26 February 2007Presented by Sébastien CandelAbstractA two-dimensional numerical study of natural convection flow in an air filled enclosure of aspect ratio A = 9 is investigatedin this Note. The numerical method is based on a second order finite volume scheme and a projection method. A full multigridtechnique is used to accelerate the convergence of the Poisson pressure equation. The multigrid procedure is briefly described andthe critical Rayleigh number above which the flow becomes unsteady is determined. To cite this article: N. Ben Cheikh et al.,C. R. Mecanique 335 (2007).© 2007 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.RésuméÉcoulement de convection naturelle dans une longue cavité à l’aide d’une méthode multigrille. Nous présentons dans cepapier une étude numérique de convection naturelle bidimensionnelle dans une cavité de rapport de forme A = 9, remplie d’air. Laméthode numérique est basée sur un schéma de type volumes finis du second ordre et une méthode de projection. Une approchemultigrille est utilisée pour accélérer la convergence de l’équation de Poisson. La méthode est brièvement décrite et le nombre deRayleigh critique au-delà duquel l’écoulement devient instationnaire est déterminé. Pour citer cet article : N. Ben Cheikh et al.,C. R. Mecanique 335 (2007).© 2007 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.Keywords: Computational fluid mechanics; Natural convection flow; Multigrid techniqueMots-clés : Mécanique des fluides numérique ; Écoulement de convection naturelle ; Approche multigrille1. IntroductionAccurate solutions of the incompressible Navier–Stokes/Boussinesq equations (NSBE) are of considerable interestin many industrial flows, such as speed aerodynamics and hydrodynamics. The application of computational fluiddynamics methods in these areas shows that for sufficient accuracy, large grid sizes are needed for capturing thinboundary layer properties, detecting high heat transfer spots or resolving small eddy regions. Because of that need,procedures for solving the NSBE, in acceptable CPU times, are at the heart of large scale computations.* Corresponding author.E-mail address: nader_bc@yahoo.fr (N. Ben Cheikh).1631-0721/$ – see front matter © 2007 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.doi:10.1016/j.crme.2007.01.004114 N. Ben Cheikh et al. / C. R. Mecanique 335 (2007) 113–118So far, one of the most popular techniques for solving the unsteady NSBE may be the use of multigrid methods(MG) [1]. The applications of MG methods have indeed shown nearly optimum convergence characteristics [2–4].The computation times are directly proportional to the number of grid points, allowing very fine grids to be used andtherefore more accurate solutions to be obtained.In this Note a finite volume, full multi-grid method (FMG) applied to natural convection flows is utilised. Thesmoother used is the iterative red and black successive over relaxation scheme (RBSOR) [5]. In order to accelerate theconvergence, an acceleration parameter is implemented in the classical FMG procedure. It was indeed demonstratedin Ref. [6] that convergence could be substantially improved by just multiplying the correction field in the FMGprocedure by some suitable factor Γ > 0 (AFMG method).In the next section the numerical approach and the AFMG method are briefly described. Next, the implementationof the AFMG algorithm is validated by a time dependent benchmark problem of natural convection in a tall cavityof aspect ratio 8 [7,8]. The method is then applied to study natural convection in an enclosure of aspect ratio 9 andto determine the critical Rayleigh above which the flow becomes unsteady. In the final section the most importantfindings of this study are summarized.2. Numerical approachThe non-dimensional governing equations for an incompressible flow, corresponding to the continuity, Navier–Stokes and energy equations, under the Boussinesq approximation, are given by:∂ui∂xi= 0 (1)∂ui∂t+ ∂(ujui)∂xj= − ∂p∂xi+√PrRa∂2ui∂xj ∂xj+ θδi2 (2)and∂θ∂t+ ∂(uiθ)∂xi= 1√Ra Pr∂2θ∂xi∂xi(i = 1,2) (3)where ui = (u, v), p and θ are the velocity, the kinematic pressure, and temperature, respectively. Here δij is theKrönecker symbol. Eqs. (1)–(3) where obtained using the characteristic length W , velocity scale u0 = √gβWT ,time scale t0 = W/u0 and pressure scale p0 = ρu20. ρ is the mass density, g the gravitational acceleration, and βthe isobaric coefficient of thermal expansion. The non-dimensional temperature θ is defined in terms of the walltemperature difference and a reference temperature:θ = T − TrTh − TcTh is the temperature of the hot wall, and Tc is that of the cold wall. The non-dimensional parameters are the Rayleighnumber (Ra) and the Prandtl number (Pr) defined as Ra = gβT W 3/(να) and Pr = ν/α , where ν is the kinematicviscosity, α the thermal diffusivity, and T = Th − Tc the temperature difference between the hot and cold walls.The unsteady Navier–Stokes and energy equations are discretized by a second-order time stepping of finite differ-ence type. Non-linear terms in Eq. (2) are treated explicitly with a second-order Adams–Bashforth scheme. Convectiveterms in Eq. (3) are treated semi-implicitly and diffusion terms in both Eqs. (2) and (3) are treated implicitly.In order to get round the difficulty that resides in the strong velocity-pressure coupling, we choose to use a projec-tion method [9,10]. Using this method, the global algorithm consists of:• solving at a first step Eq. (4) related to the energy equation:3θn+1 − 4θn + θn−12t+ ∂[(2uni − un−1i )θn+1]∂xi= 1√Ra Pr∂2θn+1∂xi∂xi(4)N. Ben Cheikh et al. / C. R. Mecanique 335 (2007) 113–118 115• the temperature field θn+1 permits the solution of the system of Eqs. (5) related to an intermediate velocity fieldu∗i that may not be divergence free:3u∗i − 4uni + un−1i2t+ 2∂(ujui)n∂xj− ∂(ujui)n−1∂xj= −(∂p∂xi)n+√PrRa∂2u∗i∂xj ∂xj+ θn+1δi2 (5)• the obtained intermediate velocity field permits to solve the Poisson equation:∂2φ∂xi∂xi= ∂ui∂xi(6)where φ is an arbitrary scalar field;• once u∗i and φ are obtained, the real velocity field and the pressure are updated by the following expressions:un+1i = u∗i −∂φ∂xi(7)pn+1 = pn + 3φ2t−√PrRa(∂2φ∂xi∂xi)(8)A second order finite-volume method [11] is used to discretize the Navier–Stokes and energy equations. Thus,equations to solve are the equation of energy, the equation related to an intermediate velocity field and a Poissonequation. The two first equations are solved using the red and black successive over-relaxation method (RBSOR) [5],while the Poisson pressure correction equation is solved using an accelerated full multigrid method (AFMG) withthe RBSOR method as a smoother. The main idea of the multigrid method can be found in Refs. [1,2]. We haveimplemented our multigrid procedure in a so-called full multigrid (FMG) fashion [2]. Indeed, before starting V-cycles,the source term is calculated on the coarsest grid permitting the determination of an exact solution. This solution isprogressively interpolated from the coarsest to the finest grid, used there as a starting guess for the V-cycle procedure.In order to optimize the number of V-cycles, it is demonstrated [6] that convergence can be substantially improvedby just multiplying the correction field in the multigrid procedure by some suitable factor Γ > 0 (AFMG method).Several tests showed that for Γ = 3.75 and after a short time of integration, only one to two V-cycles are necessary tosatisfy the convergence.3. Results and discussion3.1. Code validationAs a test case we considered the time dependent solution for an 8:1 differentially heated cavity. The benchmarkdescription and the physical problem are well detailed in Ref. [7]. Thus, in this paragraph, results related to Ra =3.4 × 105, Pr = 0.71 and A = H/W = 8 are presented. Here, H and W are the height and the width of the enclosure,respectively. The number of control volumes applied on the finest grid is 96 × 480 corresponding to 6 levels in themultigrid calculation. The coarsest grid level consisted thus of 3 × 15 number of control volumes. Cartesian gridsare uniform in the vertical direction and non-uniform in the horizontal direction with smallest cells near hot and coldwalls, respectively.We ran the computations on three different grids with increasing refinement: 24 × 120 nodes (coarse), 48 × 240nodes (medium) and 96 × 480 nodes (fine). Tests on a coarse and a medium grid, respectively, 24 × 120 and 48 × 240(corresponding to 4 and 5 levels in the multigrid calculations), were first carried on to determine an optimal value offactor Γ . On both grids, a mean value of Γ = 3.75 was obtained assuming that the value of Γ is grid independent.According to the CFL criterion, time steps on grids 24×120, 48×240 and 96×480 were fixed at t = 3.2×10−2,t = 1.6×10−2 and t = 8×10−3, respectively. Computational results were obtained for each of the three grid typesand simulations were run until a non-dimensional time of t = 1500.During the flow solution, time history data at point 1, i.e., (x = 0.181; y = 7.37), in the cavity were reported ateach time step. Table 1 provides a summary of the time-history data for the three grid resolutions. The table presentsinformation (at point 1) for the velocity in the x-direction (u1), the velocity in the y-direction (v1), the temperature (θ1)and the Nusselt number (Nu). The corresponding fluctuating x-velocity, y-velocity, temperature, and Nusselt numberare u1, v1, θ1, and Nu. The period, τθ , is the period associated with the temperature oscillation at point 1.116 N. Ben Cheikh et al. / C. R. Mecanique 335 (2007) 113–118Table 1Comparison of the time-averaged data and the corresponding fluctuations at point 1Tableau 1Comparaison des grandeurs moyennes et fluctuations correspondantes au point 1Quantity Grid resolution24 × 120 48 × 240 96 × 480u1 5.01203 × 10−2 5.37072 × 10−2 5.61780 × 10−2u1 7.16083 × 10−4 3.59755 × 10−2 5.43379 × 10−2v1 4.52762 × 10−1 4.59891 × 10−1 4.61412 × 10−1v1 1.12033 × 10−3 5.28620 × 10−2 7.67392 × 10−2θ1 2.68107 × 10−1 2.66011 × 10−1 2.65582 × 10−1θ1 5.93518 × 10−4 2.88141 × 10−2 4.26336 × 10−2Nu 4.61383 4.58802 4.58179Nu 1.70366 × 10−4 4.98058 × 10−3 7.10453 × 10−3τθ 3.7440 3.4400 3.4160In order to compare the performances of contributions to the benchmark problem, a number of metrics basedon some of the quantities shown in Table 1 were defined in Ref. [7]. The values obtained by each participant werecompared against the accepted values of Xin and Le Quéré [8], and the deviations from these latter values defined themetrics. Hence, values as close to zero as possible are the ‘target’. Here, we present our performance, based on theoverall metric, E, defined by:E = ε̄θ1 + ε̄u1 + ε̄Nu1 + ετθ + ε′θ1+ ε′u1 + ε′Nu17(9)which represents the average of the errors computed in the period, mean velocity, temperature and Nusselt number,and the corresponding error in the oscillation amplitudes. This metric has been calculated on the grid 96 × 480 for1465.84  t  1500 (ten periods) and compared with some of the contributors to the MIT session (see Table 2). Asseen, our metric is ranked first, meaning that the numerical method is very efficient giving results very close to the‘true’ solution provided by Xin and Le Quéré [8].In order to reflect the computational cost required to obtain the solution, a normalized algorithm timing defined as:ηAT =(m secnode.step.stepsperiod).specFP95 (10)was proposed in Ref. [7]. This normalized performance metric is scaled in a relative sense so that the minimumvalue of ηAT is unity. Our metric measured on a Dell Dimension 4600 system, with a single 2.9 GHz processor(specFP95 = 71), is presented in Table 3 and is compared with some of the contributors to the MIT session. As it canbe seen, the value is ranked fourth, meaning that the method is quite fast.Table 2Comparison of overall metric ETableau 2Comparaison du métrique EFirst author [7] E Rank orderPresent study 0.2491 1 / (31)Johnston 0.3352 2 / (31)Davis 0.5631 3 / (31)Armfield 0.9375 7 / (31)Christon 0.9849 9 / (31)Westerberg 2.7708 15 / (31)Bruneau 9.5773 21 / (31)Chan 21.266 27 / (31)Table 3Comparison of the normalized algorithm timing metricTableau 3Comparaison du métrique du temps algorithmique normaliséFirst author [7] ηAT RankingJohnston 1.00 1 / (26)Christon II 1.31 2 / (26)Christon III 1.35 3 / (26)Pres. Study 2.11 4 / (26)Matsumoto 2.34 5 / (26)Bruneau 3.40 6 / (26)Armfield 85.14 22 / (26)Le Quéré 90.17 24 / (26)Ingberg 2010.14 26 / (26)N. Ben Cheikh et al. / C. R. Mecanique 335 (2007) 113–118 117Fig. 1. Nusselt number versus Rayleigh.Fig. 1. Nombre de Nusselt en fonction du Rayleigh.Fig. 2. Vorticity isolines for: (a) Ra = 104, (b) Ra =3.4 × 105 and (c) Ra = 106.Fig. 2. Iso-vorticités pour : (a) Ra = 104, (b) Ra =3.4 × 105 and (c) Ra = 106.Fig. 3. (Amplitude)2 versus Rayleigh.Fig. 3. (Amplitude)2 en fonction du Rayleigh.3.2. Natural convection in a 9:1 heated enclosureWe here consider the problem of a differentially heated closed cavity of aspect ratio A = 9. The numerical re-sults are obtained for Rayleigh numbers ranging from Ra = 102 to Ra = 106 and a Prandtl number of Pr = 0.71.118 N. Ben Cheikh et al. / C. R. Mecanique 335 (2007) 113–118A proportional number of control volumes according to the 8:1 heated enclosure case described above was applied,i.e., 96 × 520 nodes. During the flow solution, time history data at point 2 of coordinate (x = 0.15; y = 7.75) werereported at each time step.The results showed that for Rayleigh numbers lower than Ra = 105, the flow was still steady and for Rayleighnumbers somewhat greater than 2 × 105 the flow bifurcates to an unsteady periodic state, followed by a chaotic stateif Rayleigh number was till increased.In Fig. 1 the average Nusselt number, Nu, is plotted versus the Rayleigh number. As expected, the Nusselt numberincreases with Ra and for Raleigh numbers lower than 103, Nu is very weak and therefore diffusion is the principalmode of heat transfer. For Ra greater than 5 × 104, the intensity of convection increases and in the range 5 × 104 Ra  106, Nu can be expressed through the following correlation:Nu = 0.178Ra0.252 (11)Fig. 2 shows the instantaneous iso-vorticity features of the flow for three Rayleigh numbers related to a steady, a pe-riodic and a chaotic state, obtained for Ra = 104, Ra = 3.4 × 105 and Ra = 106, respectively. As seen, by increasingthe Rayleigh number, a more complex solution is found and an increase of vortex dynamics is observed. Note thatthe period of the flow at Ra = 3.4 × 105 has been calculated, and comparing to the 8:1 benchmark problem, a smallervalue of τθ = 3.120 was found.The determination of the critical Rayleigh number, Rac , was determined using the same method as presented inRef. [12]. In fact, the solution of a supercritical Hopf bifurcation is proportional to√Ra − Rac . Hence, in Fig. 3, wereport amplitudes squared, (Amp)2, of quantities θ2, u2 and v2 at point 2 versus the Rayleigh number. By extrapolatingresults at Amp = 0, we obtained an approximate value of Rac = 2.4 × 105. Comparing that value with proposed onesrelated to the benchmark enclosure of aspect ratio A = 8 (Rac ∼ 3.0×105), a decrease of the critical Rayleigh numbermay be observed by increasing the aspect ratio.4. ConclusionIn this Note, a numerical method for solving the Navier–Stokes/Boussinesq equations has been investigated. Thenumerical method is based on a finite volume approach, a projection method and an accelerated full-multigrid tech-nique (AFMG). The performed test calculations demonstrated the potential of the method for enabling very accuratesolutions on very fine grids in acceptable CPU times. Tests showed that a mean value for acceleration parameter Γ of3.75 was well appropriated for natural convection flows in differentially heated enclosures. Using the AFMG method,a study in an enclosure of aspect ratio A = 9 was investigated. For Rayleigh numbers ranging from 104 to 106,a correlation between Nusselt number and Rayleigh number has been reported. The supercritical solution was alsodetermined giving an approximate critical Rayleigh number of Rac = 2.4 × 105.References[1] W. Hackbusch, Multigrid Methods and Applications, Springer-Verlag, Berlin/New York, 1985.[2] M. Hortmann, M. Peric, G. Scheuerer, Finite volume multigrid prediction of laminar natural convection: bench-mark solutions, Int. J. Numer.Methods Fluids 11 (1990) 189–207.[3] E. Nobile, Simulation of time-dependent flow in cavities with the additive-correction multigrid method, Part I: Mathematical formulation,Numer. Heat Transfer B 30 (1996) 341–350.[4] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput. 31 (1977) 333–390.[5] R. Barrett, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, 1994.[6] D. Braess, Towards algebraic multigrid for elliptic problems of second order, Computing 55 (1995) 379–393.[7] M.A. Christon, P.M. Gresho, S.B. Sutton, Computational predictability of time-dependent natural convection flows in enclosures (includinga benchmark solution), Int. J. Numer. Methods Fluids 40 (2002) 953–980.[8] S. Xin, P. Le Quéré, An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity, Int. J. Numer. MethodsFluids 40 (2002) 981–998.[9] R. Peyret, T.D. Taylor, Methods for Fluid Flow, Springer-Verlag, Berlin/New York, 1983.[10] Y. Achdou, J.L. Guermond, Convergence analysis of a finite element projection/Lagrange–Galerkin method for the incompressible Navier–Stokes equations, SIAM J. Numer. Anal. 37 (2000) 799–826.[11] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980.[12] H. Wang, S. Xin, P. Le Quéré, Étude numérique du couplage de la convection naturelle avec le rayonnement de surfaces en cavité carréeremplie d’air, C. R. Mécanique 334 (2006) 48–57.

Natural convection flow in a tall enclosure using a multigrid method

https://comptes-rendus.academie-sciences.fr/mecanique/item/10.1016/j.crme.2007.01.004.pdf

Natural convection flow in a tall enclosure using a multigrid method

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mécanique (CRM)