Rayleigh–Benard Convection

PIC

Rayleigh–Bxnard Convection

Department of Aeronautics, Imperial College London, UKScientific Computing and Imaging Institute, University of Utah, USA November 9, 2025

Tutorials

1 Introduction Governing equations Mesh generation 2 Direct Numerical Simulation Nusselt number 3 Linear Stability Analysis 3.1 Estimate the value of critical Rayleigh number

1 Introduction

[next] [prev] [prev-tail] [tail] [up] Chapterx1Introduction This tutorial provides an overview of utilising the spectral/hp element framework Nektar++ for the numerical solution of a classical fluid dynamics problem: Rayleigh–Bxnard Convection (RBC). The tutorial will delve into the direct numerical simulation and linear stability analysis of a two-dimensional RBC. Information on how to install the libraries, solvers, and utilities on your own computer is available on the webpage www.nektar.info. This tutorial assumes the reader has already completed the previous tutorial in the Flow Stability series on the channel and already has the necessary software installed. Task: 1.1 Prepare for the tutorial. Make sure that you have:

Downloaded the tutorial files: rbc.tar.gz

Unpack it using unzip rbc.tar.gz to produce a directory rbc with subdirectories called tutorial and complete. We will refer to the tutorial directory as

     $NEKTUTORIAL.

Governing equations

The dynamical equations that describe Rayleigh–Bxnard Convection (RBC) under Boussinesq approximation in the non-dimensional form are

\begin{aligned}\frac{∂u}{∂t} + u · ∇u & = −∇p + RaP rT e_{y} + P r∇^{2}u, & (1.1) \\ \frac{∂T}{∂t} + u · ∇T & = ∇^{2}T, & (1.2) \\ ∇ · u & = 0, & (1.3) \end{aligned}

where u and T are the velocity and temperature fields, respectively, and e_{y} is the buoyancy direction. Here p is the modified pressure including the buoyancy term accounting for the reference temperature. Flow is controlled by two non-dimensional numbers, the Rayleigh number Ra and the Prandtl number P r.

By decomposing the velocity, temperature and pressure in the equations (1.2-1.3) as a summation of a base flow (U, T, P) and perturbation (u′, T′, p′), such that

\begin{aligned}u & = U + u′ & (1.4) \\ T & = T + T′ & (1.5) \\ p & = P + p′, & (1.6) \end{aligned}

and retaining first-order terms only yields the linearised equations governing the evolution of infinitesimal perturbations:

\frac{∂u′}{∂t} + u′· ∇U + (U · ∇)u′  = −∇p′+ RaP rT′e_{y} + P r∇^{2}u′,  (1.7)
\frac{∂T′}{∂t}+ u′· ∇T + U · ∇T′  = ∇^{2}T′,  (1.8)
∇ · u′  = 0.  (1.9)

In Chapterx2, we will validate a case of two-dimensional RBC in a domain of aspect ratio 2.02. In Chapterx3, we will discuss the linear stability analysis for a no-slip boundary condition at top and bottom plates, and will estimate the critical value of Rayleigh number. The files contained in the $NEKTUTORIAL directory are as follows:

     Folder geometry

rbc.geo - Gmshxfile that contains the geometry of the problem.

rbc.msh - Gmshxgenerated mesh data listing mesh vertices and elements.

     Folder DNS/Ra_5e3_Pr_0p71

rbc-DNS.xml - Nektar++ session file, generated with the $NEK/NekMesh utility, for computing the flow for Ra = 5000 and P r = 0.71.

     Folder LSA/Ra_1900

rbc-LSA.xml - Nektar++ session file, generated with $NEK/NekMesh for the linear stability analysis for Ra = 1900 and P r = 0.71.

Mesh generation

The first step is to generate a mesh that is readable by Nektar++. The files necessary in this section can be found in $NEKTUTORIAL/geometry/. To achieve this task, we use Gmshxin conjunction with the Nektar++ pre-processing utility called $NEK/NekMesh. Specifically, we first generate the mesh in figure 1.1 using Gmshxand successively, we convert it into a suitable Nektar++ format using $NEK/NekMesh. For all calculations in this tutorial, we fix L_{x} = 2.02 and L_{y} = 1.0 to maintain the aspect ratio of 2.02. Task: 1.2 Convert the Gmsh geometry provided into the XML Nektar++ format and with two periodic boundaries rbc.msh can be generated using Gmsh by running the following command:

 gmsh -2 rbc.geo

rbc.xml can be generated using the module peralign available with the pre-processing tool $NEK/NekMesh:

 $NEK/NekMesh -m peralign:surf1=2:surf2=3:dir=x rbc.mesh rbc.xml

where surf1 and surf2 correspond to the periodic physical surface IDs specified in Gmsh (in our case the inflow has a physical ID=2 while the outflow has a physical ID=3) and dir is the periodicity direction (i.e. the direction normal to the inflow and outflow boundaries - in our case x).

pict

Figurex1.1 100 quadrilaterals mesh

[next] [prev] [prev-tail] [front] [up] 2 Direct Numerical Simulation

[prev] [prev-tail] [tail] [up]

Chapterx2Direct Numerical Simulation

This chapter will consider two-dimensional Rayleigh–Bxnard convection in a box of aspect ratio 2.02. At the top and bottom walls, we employ a no-slip boundary condition for the velocity field and a conducting boundary condition for the temperature field. At the sidewalls, we set periodic boundary conditions for both fields. We will compare our results with Clever and Busse (J. Fluid Mech., 1974; 65:625-645) for Ra = 5000, P r = 0.71, and L_{x} = 2π/3.117 = 2.02.

In the folder $NEKTUTORIAL/DNS/Ra_5e3_Pr_0p71 you will find the file rbc-DNS.xml which

contains the geometry along with the necessary parameters to solve the problem. The GEOMETRY section is responsible for defining the mesh of the problem, and it is automatically generated, as demonstrated in the preceding task. Within the EXPANSIONS section, the expansion type and order are specified. An expansion basis is applied to a geometric composite, where composite refers to a collection of mesh entities. The $NEK/NekMesh utility always includes a default entry, and in this case, the composite C[0] represents the set of all elements. The FIELDS attribute indicates the fields for which this expansion is applicable, and the TYPE attribute designates the type of polynomial basis functions used in the expansion. For example,

   1 <EXPANSIONS> 2 xxx<ECOMPOSITE="C[0]"NUMMODES="10"FIELDS="u,v,T,p" 3 xxxTYPE="GLL_LAGRANGE_SEM"/> 4 </EXPANSIONS>.

To finalise the problem definition and facilitate solution generation, we define a section named CONDITIONS within the session file. Specifically, the CONDITIONS section encompasses the following entries:

Solver information (SOLVERINFO) such as the equation, the projection type, the evolution operator, and the analysis driver to use, along with other properties. The solver properties are specified as quoted attributes and have the form

        1 <SOLVERINFO> 2 xxx<IPROPERTY="[STRING]"VALUE="[STRING]"/> 3 xxx... 4 </SOLVERINFO>

Task: 2.1 In the SOLVERINFO section of rbc-DNS.xml:

Note: The bits to be completed are identified by xin this file. set the EvolutionOperator to Nonlinear in order to select the equations (1.2)–(1.3),

set the Projection property to Continuous in order to select the continuous Galerkin approach,

set the Driver to Standard in order to perform standard time-integration.

The parameters (PARAMETERS) are specified as name-value pairs:

  1 <PARAMETERS> 2 xxx<P>[KEY]=[VALUE]</P> 3 xxx... 4 </PARAMETERS>

Task: 2.2 Declare the two parameters Ra, that represents the Rayleigh number, and Pr, which represents the Prandtl number. Now set the Rayleigh number to 5000 and the Prandtl number to 0.71 - i.e.

 <P> Ra = 5000 </P>
<P> Pr = 0.71 </P>

The declaration of the variable(s) (VARIABLES).

  1 <VARIABLES> 2 xxxxxxxx<VID="0">u</V> 3 xxxxxxxx<VID="1">v</V> 4 xxxxxxxx<VID="2">T</V> 5 xxxxxxxx<VID="3">p</V> 6 </VARIABLES>

Defining boundary regions (BOUNDARYREGIONS) involves specifying composites from the GEOMETRY section and the associated conditions applied to these boundaries (BOUNDARYCONDITIONS). The format for boundary conditions imposed on a region must be in accordance with the conditions outlined for each variable in the VARIABLES section to establish a well-posed problem. The REF attribute associated with a boundary condition region ought to align with the ID="[INDEX]" of the intended boundary region outlined in the BOUNDARYREGIONS section. For example, to model a hot bottom plate corresponding to REF="0", we set the following condition:

  1 <REGIONREF="0"> 2 xxx<DVAR="u"VALUE="0"/> 3 xxx<DVAR="v"VALUE="0"/> 4 xxx<DVAR="T"VALUE="1"/> 5 xxx<NVAR="p"USERDEFINEDTYPE="H"VALUE="0"/> 6 </REGION>

Task: 2.3 Model the top-cooled plate, which corresponds to REF="1" in the session file.

The specification of (time- and) space-dependent functions (FUNCTION), involves defining functions in terms of x, y, z, and t, encompassing aspects like initial conditions and forcing functions. The VARIABLES denote the components of the particular function in a specified direction.

  1 <FUNCTIONNAME="[NAME]"> 2 xxx<EVAR="[VARIABLE_1]"VALUE="[EXPRESSION]"/> 3 xxx<EVAR="[VARIABLE_2]"VALUE="[EXPRESSION]"/> 4 xxx... 5 </FUNCTION>

Task: 2.4 Define a function called InitialConditions. The initial condition can be derived from the solution of the conduction equation for the temperature field in the absence of fluid motion. Therefore, the initial condition will be as follows:

\begin{aligned}u & = 0 & (2.1) \\ v & = 0 & (2.2) \\ T & = 1 − y & (2.3) \\ p & = 0 & (2.4) \end{aligned}

Next, we apply the buoyancy force in the vertical direction, with gravity acting downward. This is implemented by defining the following function as the body force.

   1 <FUNCTIONNAME="BodyForce"> 2 xxxxxx<EVAR="u"VALUE="0"EVARS="uvT"> 3 xxxxxx<EVAR="v"VALUE="Ra*Pr*T"EVARS="uvT"/> 4 xxxxxx<EVAR="T"VALUE="0"EVARS="uvT"/> 5 </FUNCTION>

Note that for using the body force you need the following additional tag outside the section

CONDITIONS:

   1 <FORCING> 2 xxx<FORCETYPE="Body"> 3 xxxxxx<BODYFORCE>BodyForce</BODYFORCE> 4 xxx</FORCE> 5 </FORCING>

We also wish to monitor the simulation progress. To achieve this, we plan to record flow quantities at the location (x, y) = (1.0, 0.5). For this purpose, we incorporate a filter called HistoryPoints in the session file as follows:

   1 xxx<FILTERS> 2 xxxxxx<FILTERTYPE="HistoryPoints"> 3 xxxxxxxx<PARAMNAME="OutputFile">rbc</PARAM> 4 xxxxxxxx<PARAMNAME="OutputFrequency">10</PARAM> 5 xxxxxxxx<PARAMNAME="Points"> 6 xxxxxxxxxxx1.00.50 7 xxxxxxxx</PARAM> 8 xxxxxx</FILTER> 9 xxx</FILTERS>

This completes the specification of the problem. Task: 2.5 Run the solver by typing the following command on the command line:

 $NEK/IncNavierStokesSolver rbc-DNS.xml

Now, after completing the calculation, we wish to verify if the solution has reached a steady state. To do this, we plot the time history of v(t) at (x, y) = (1.0, 0.5), displayed in figurex2.1. As evident from the plot, the solution has reached a steady state. It is worth noting that you can use any plotting software for this analysis. pict

Figurex2.1 Time history of v(t) velocity component measured at (x, y) = (1.0, 0.5).

Next, we want to visualise the steady solution. Specifically, we need to convert the .fld file into a format readable by a visualisation post-processing tool. In this tutorial we decided to convert the .fld file into a VTK format and to use the open-source visualisation package called Visit. Task: 2.6 Convert the file:

 $NEK/FieldConvert rbc-DNS.xml rbc-DNS.fld rbc-DNS.vtu

Now open Visit and use File ->Open, to select the VTK file. Next, you can plot different fields.Figurex2.2 displays the density plots of the temperature field and the velocity vector field in the x − y plane.

PIC

Figurex2.2 Converged solution for Ra = 5000 and P r = 0.71.

Nusselt number

To validate our results, we will calculate the Nusselt number (Nu), representing the ratio of the total heat flux (convective plus conductive) to the conductive heat flux at the plate. The Nusselt number is expressed as

follows1 1Gary A. Glatzmaier, Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification, Rotation, Princeton University Press, 2013. :

Nu = −\frac{1}{L_{x}} ∫ \frac{∂T}{∂y}dx.  (2.5)

For the calculation of the Nusselt number, we will use the FieldCovert and NekMesh utility. We will follow the following steps:

Extract the wall boundary surface in a xml file using the NekMesh utility as:

NekMesh-mextract:surf=1rbc-DNS.xmlrbc_surf.xml

Tip: Remember to change the NUMODES and the polynomial basis functions according to the session file in rbc-DNS.xml.

Compute \frac{∂T}{∂y}, and remove all unnecessary fields.

FieldConvert-mgradientrbc-DNS.xmlrbc-DNS.fldtemp_1.fldFieldConvert-mremovefield:fieldname="u,v,T,p,u_x,u_y,\\v_x,v_y,T_x,p_x,p_y"rbc-DNS.xmltemp_1.flddT_dy.fld

Next, extract \frac{∂T}{∂y} at the top boundary.

FieldConvert-mextract:bnd=0rbc-DNS.xmldT_dy.flddT_dy_bnd.fld

Finally, we will compute \frac{1}{L_{x}} ∫ \frac{∂T}{∂y}dx using mean utility.

FieldConvert-mmeanrbc_surf.xmldT_dy_bnd_b0.fldstdout

Above command gives the following output:

Domainlength:2.02Integral(variableT_y):-4.26191Mean(variableT_y):-2.10986

Now since

. \frac{1}{L_{x}} ∫ \frac{∂T}{∂y}dx = −2.10986 ⟹ Nu = −\frac{1}{L_{x}} ∫ \frac{∂T}{∂y}dx = 2.10986.

Our calculation of the Nusselt number match very well with that of Clever and Busse (J. Fluid Mech., 1974; 65:625-645), which is 2.112. Task: 2.7 In the file rbc-DNS.xml change the values of Ra and P r, and again run the solver to compute the Nusselt number (your results should match with Tablex(2.1).

\begin{matrix}11  Ra & 11  P r & 11  Δt & 11  Nu \\ 11  10^{4} & 11 0.71 & 11 5 × 10^{−4} & 11 2.65179 \\ 11 3 × 10^{3} & 11  7.0 & 11  10^{−3} & 11 1.66236 \end{matrix}

Tablex2.1 Nusselt number as a function of Rayleigh number and Prandtl number. Above results matches very well with the Table 1 of Clever and Busse (J. Fluid Mech., 1974; 65:625-645).

This completes the tutorial.

[prev] [prev-tail] [front] [up] 3 Linear Stability Analysis

[prev] [prev-tail] [tail] [up]

Chapterx3Linear Stability Analysis

In this chapter, we will compute the leading eigenvalues and eigenvectors for Rayleigh–Bxnard Convection (RBC). We have the analytical solution for the base flow, which is simply the solution of the conduction equation for the temperature field in the absence of fluid motion (the same solution used as the initial condition in the previous chapter). Therefore, U = 0 and T = 1 − y. Our goal is to reproduce the value of the critical Rayleigh number, Ra_{c} = 1707.762, first predicated by

Chandrasekhar1 ^{1S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Clarendon, 1968.} .

First, we will perform linear stability analysis for Ra = 1900 and P r = 0.71. In the $NEKTUTORIAL/LSA/Ra_1900 folder there is a file called rbc-LSA.xml. This file is similar to rbc-DNS.xml, but contains additional instructions to perform the linear stability analysis. Task: 3.1 Configure the following additional SOLVERINFO options which are related to the stability analysis.

set the EvolutionOperator to Direct in order to activate the linearised governing equations (1.8)–(1.9).

set the Driver to ModifiedArnoldi in order to use the ModifiedArnoldi eigenvalue analysis.

Task: 3.2 Set the parameters for the Arnoldi algorithm.

kdim=8: dimension of Krylov-space, nvec=1: number of requested eigenvalues,

nits=500: number of maximum allowed iterations,

evtol=1e-7: accepted tolerance on the eigenvalues and it determines the stopping criterion of the method.

Next, we need to define the InitialConditions. We will use the random initial condition, which is implemented by using following within the CONDITIONS section.

1 <FUNCTIONNAME="InitialConditions"> 2 xxx<EVAR="u"VALUE="awgn(0.0001)"/> 3 xxx<EVAR="v"VALUE="awgn(0.0001)"/> 4 xxx<EVAR="T"VALUE="awgn(0.0001)"/> 5 xxx<EVAR="p"VALUE="0"/> 6 </FUNCTION>

Task: 3.3 Define a function called BaseFlow. Set the following function for the base flow:

\begin{aligned}U & = 0 & (3.1) \\ V & = 0 & (3.2) \\ T & = 1 − y & (3.3) \\ P & = 0 & (3.4) \end{aligned}

Task: 3.4 Run the solver to perform the analysis

 $NEK/IncNavierStokesSolver rbc-LSA.xml

The simulation should converge in 8 iterations, and the terminal screen should look similar to the one below at the end:

    EV:01.1332501.250880Writing:"rbc-LSA_eig_0.fld"(0.015323s,XML)

Task: 3.5 Plot the leading eigenvector in Paraviewxor VisIt. This should look like the solution shown in figures 3.1.

PIC

Figurex3.1 The velocity vector field u′ on top of the density plots of T′ for the leading eigenmode.

3.1 Estimate the value of critical Rayleigh number Task: 3.6 Now again run the linear stability analysis for Ra = 1600, 1700, 1800. Value of growth rates for these simulation are given in Tablex3.1.

\begin{matrix}11  Ra & 11  σ \\ 11 1600 & 11  −0.731021 \\ 11 1700 & 11 −0.0519525 \\ 11 1800 & 11  0.608177 \end{matrix}

Tablex3.1 Growth rate of the most unstable mode as a function of Rayleigh number.

Task: 3.7 Use linear interpolation to estimate the Rayleigh critical. You should get Ra_{c} = 1707.87, which makes a good agreement with the value obtained by Chandrasekhar (1961).For linear interpolation you can use following Python script:

    In[1]:Ra=np.array([1600,1700,1800,1900])In[2]:sigma=np.array([-0.731021,-0.0519525,0.608177,1.25088])In[3]:np.lnterp(0,sigma,Ra)Out[3]:1707.8700467105318

This completes the tutorial.

[prev] [prev-tail] [front] [up] [next]