# 3D Incompressible Rayleigh-Benard Convection Tutorial
This tutorial documents the input deck `incompressible_3d_rayleigh_benard_convection.xml` and the associated Exodus mesh `3d_rayleigh_benard.exo`. The case models buoyancy-driven convection of an incompressible fluid confined between a heated bottom wall and a cooled top wall inside a cylindrical enclosure — the classical Rayleigh–Bénard configuration.
---
## 1. Problem description
The simulation solves for the transient, three-dimensional flow of a single-phase incompressible fluid enclosed in a right circular cylinder. The fluid is at rest initially and is subsequently driven into convective motion by a temperature difference imposed between the bottom (hot) and top (cold) end walls, with gravity $\mathbf{g}$ acting along the negative $z$-direction (axis of the cylinder). The lateral wall of the cylinder is adiabatic and no-slip, so that heat can only enter and leave the domain through the bottom and top disks.
The physics is built from the `IncompressibleNavierStokes` physics block with the temperature equation and the Boussinesq buoyancy source both enabled (`Build Temperature Equation = true`, `Build Buoyancy Source = true`), and the continuity constraint is enforced through an **Entropically Damped Artificial Compressibility (EDAC)** formulation rather than a classical projection/pressure-Poisson approach. This choice, together with the fluid properties and boundary temperatures, fixes the Prandtl and Rayleigh numbers of the problem, discussed in Section 2.
The XML snippet below shows how the physics is attached to the single element block (`block-1`) present in the mesh:
```xml
```
---
## 2. Governing equations
The fluid is modeled as incompressible with constant transport properties under the Boussinesq approximation (density variations are retained only in the buoyancy term). Three fields are solved: velocity $\mathbf{u}$, a Lagrange/kinematic pressure $p$, and temperature $T$.
### 2.1 Momentum equation
$$
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \nu \nabla^{2}\mathbf{u} + \alpha_{th}\,(T - T_{ref})\,\mathbf{g}
$$
where $\nu$ is the kinematic viscosity, $\alpha_{th}$ is the thermal expansion coefficient, $T_{ref}$ is the reference temperature about which the Boussinesq buoyancy term is linearized, and $\mathbf{g}$ is the gravitational acceleration vector. The last term is built by the `IncompressibleBuoyancySource` closure model.
### 2.2 Continuity equation (EDAC)
The `Continuity Model` parameter in the `Fluid Properties` closure model is set to `EDAC`, so the divergence-free constraint $\nabla\cdot\mathbf{u}=0$ is replaced by an evolution equation for the pressure that entropically damps pressure fluctuations while remaining consistent with the incompressible limit as the artificial compressibility coefficient becomes large:
$$
\frac{\partial p}{\partial t} + \mathbf{u}\cdot\nabla p + \beta\,\nabla\cdot\mathbf{u} = \nu\,\nabla^{2} p
$$
Here $\beta$ is the **artificial compressibility coefficient** (`Artificial compressibility` in the input file), which plays the role of a pseudo-acoustic wave speed squared, and the right-hand side diffusion term uses the fluid kinematic viscosity $\nu$ as the entropic damping diffusivity, consistent with the standard EDAC formulation. As $\beta \to \infty$ this equation degenerates to $\nabla\cdot\mathbf{u}=0$, recovering the incompressible constraint; for finite $\beta$ it acts as a stabilized pseudo-time relaxation toward that constraint. Table 1 reports the value of $\beta$ used in this case.
### 2.3 Energy equation
$$
\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \kappa\,\nabla^{2}T, \qquad \kappa = \frac{k}{\rho\,c_p}
$$
where $k$ is the thermal conductivity, $\rho$ the density, and $c_p$ the specific heat capacity, so that $\kappa$ is the thermal diffusivity.
### 2.4 Fluid properties and dimensionless groups
Table 1 lists the fluid properties as specified in the `Fluid Properties` closure model.
**Table 1.** Fluid properties used in the `fluids` closure model.
| Quantity | Symbol | Value | Units |
|---|---|---|---|
| Kinematic viscosity | $\nu$ | $0.1$ | $\mathrm{m^2/s}$ |
| Density | $\rho$ | $1.0$ | $\mathrm{kg/m^3}$ |
| Thermal conductivity | $k$ | $1.0$ | $\mathrm{W/(m \cdot K)}$ |
| Specific heat capacity | $c_p$ | $70.0$ | $\mathrm{J/(kg \cdot K)}$ |
| Artificial compressibility coefficient | $\beta$ | $1.0\times10^{5}$ | $\mathrm{m^2/s^2}$ |
| Thermal expansion coefficient | $\alpha_{th}$ | $0.014285715$ | $\mathrm{K^{-1}}$ |
| Reference temperature | $T_{ref}$ | $200.0$ | $\mathrm{K}$ |
| Heat capacity ratio | $\gamma$ | $1.0$ | — |
| Continuity model | — | EDAC | — |
From these values, the thermal diffusivity is $\kappa = k/(\rho c_p) = 1.0/(1.0\times70.0) \approx 0.01429\ \mathrm{m^2/s}$, giving a **Prandtl number**
$$
Pr = \frac{\nu}{\kappa} = \frac{0.1}{0.01429} \approx 7.0
$$
The **Rayleigh number**, based on the wall-to-wall temperature difference $\Delta T = T_{bottom}-T_{top}=50.0\ \mathrm{K}$ (Section 4) and the cylinder height $H$ (Section 3), is
$$
Ra = \frac{|\mathbf{g}|\,\alpha_{th}\,\Delta T\, H^{3}}{\nu\,\kappa} \approx 5.0\times10^{3}\, H^{3}
$$
with $H$ expressed in meters. The exact value of $Ra$ therefore depends on the cylinder height read from the mesh coordinate data (Table 2); the formula is reported here so it can be evaluated once $H$ is confirmed from the Exodus file.
The syntax used to define these properties and the continuity model is:
```xml
```
The remaining closure models assemble the discrete terms of the equations above from the physics block: `DQDT` (time derivative), `CONVECTIVE_FLUX` and `VISCOUS_FLUX` (momentum flux terms), `BUOYANCY_SOURCE` (buoyancy term), together with two auxiliary models used for time-step control, `element_length` and `local_time_step_size`, and a scalar response model `nusselt_number` used for output (Section 6):
```xml
```
No electric potential or electromagnetic equation is present in this input file (the physics block is purely hydro-thermal), so no such field appears in Section 6.
---
## 3. Geometry and mesh description
The computational domain is a right circular cylinder of radius $R$ and height $H$, with its axis aligned with the (vertical) $z$-direction along which gravity acts. Figure 1 shows a schematic of the domain together with the three named surfaces used to apply boundary conditions: the **bottom** disk (hot, heated wall), the **top** disk (cold wall), and the **wall** (adiabatic lateral surface).
**Figure 1.** Schematic of the cylindrical convection cell (not to scale). The bottom disk is held at the hot wall temperature, the top disk at the cold wall temperature, and the lateral surface is adiabatic. Gravity $\mathbf{g}$ points from the cold to the hot wall, i.e. the layer is heated from below.
The mesh is stored in the Exodus file `3d_rayleigh_benard.exo`, generated with CUBIT and referenced by the input deck as follows:
```xml
```
Inspection of the Exodus/netCDF4 metadata shows a single all-hexahedral (`HEX8`, trilinear) element block, `block-1`, extruded uniformly along $z$ from a circular footprint, and three side sets that match the three boundary-condition entries of the input file. The resulting mesh statistics are summarized in Table 2, and a schematic of the mesh topology (footprint and extrusion) is shown in Figure 2.
**Table 2.** Mesh statistics extracted from `3d_rayleigh_benard.exo`.
| Quantity | Value |
|---|---|
| Element topology | HEX8 (linear hexahedron) |
| Spatial dimension | 3 |
| Number of element blocks | 1 (`block-1`) |
| Number of elements | 1568 |
| Number of nodes | 1944 |
| Number of vertically extruded layers | 7 |
| Number of side sets | 3 (`top`, `bottom`, `wall`) |
| Faces on `top` side set | 224 |
| Faces on `bottom` side set | 224 |
| Faces on `wall` side set | 252 |
**Figure 2.** Mesh view extracted from `3d_rayleigh_benard.exo`.
---
## 4. Initial and boundary conditions
### 4.1 Boundary conditions
Three boundary conditions are applied through the `IncompressibleBoundaryFlux` strategy, one per side set, as summarized in Table 3.
**Table 3.** Boundary conditions applied on `block-1`.
| Side set | Velocity condition | Thermal condition |
|---|---|---|
| `top` | No-slip ($\mathbf{u}=\mathbf{0}$) | Fixed temperature, $T_{top}=200.0\ \mathrm{K}$ |
| `bottom` | No-slip ($\mathbf{u}=\mathbf{0}$) | Fixed temperature, $T_{bottom}=250.0\ \mathrm{K}$ |
| `wall` | No-slip ($\mathbf{u}=\mathbf{0}$) | Adiabatic, $\partial T/\partial n = 0$ |
The temperature difference driving convection is $\Delta T = T_{bottom}-T_{top} = 50.0\ \mathrm{K}$, with the bottom (hot) wall below the top (cold) wall relative to gravity, i.e. the classical unstable, heated-from-below configuration. The corresponding XML entries are:
```xml
```
### 4.2 Initial conditions
At $t=0$ the fluid is at rest and isothermal at the cold-wall temperature, with a zero Lagrange pressure field, as summarized in Table 4.
**Table 4.** Initial conditions on `block-1`.
| Field | Value | Units |
|---|---|---|
| Lagrange pressure, $p$ | $0.0$ | $\mathrm{m^2/s^2}$ |
| Velocity, $u_0$ | $0.0$ | $\mathrm{m/s}$ |
| Velocity, $u_1$ | $0.0$ | $\mathrm{m/s}$ |
| Velocity, $u_2$ | $0.0$ | $\mathrm{m/s}$ |
| Temperature, $T$ | $200.0$ | $\mathrm{K}$ |
```xml
```
Gravity is specified globally under `User Data` as a constant vector body force direction used by the buoyancy closure model:
```xml
```
Here $\mathbf{g}=(0,0,-10.0)\ \mathrm{m/s^2}$, i.e. gravity of magnitude $|\mathbf{g}|=10.0\ \mathrm{m/s^2}$ acting along the negative $z$-axis, from the cold (top) toward the hot (bottom) wall. The `CFL`, `CFL_init`, and `CFL_transition` entries control a Courant-number-based ramping of the local pseudo-time step used to reach the target CFL of $10$ (dimensionless) from an initial, more conservative value of $0.001$ (dimensionless) over the course of the run, while `Workset Size` sets the number of elements ($256$) grouped per work unit during assembly.
---
## 5. Solver configuration and convergence criteria
### 5.1 Time integration
Time integration is performed with the Backward Euler stepper managed by Tempus, using an adaptive time-step controller bounded between a minimum and maximum step size, and terminating on whichever of final time or final time index is reached first.
**Table 5.** Time-stepping parameters.
| Parameter | Value | Units |
|---|---|---|
| Stepper type | Backward Euler | — |
| Initial time | $0.0$ | $\mathrm{s}$ |
| Final time | $1.0\times10^{4}$ | $\mathrm{s}$ |
| Final time index | $2000$ | steps |
| Minimum time step | $1.0\times10^{-3}$ | $\mathrm{s}$ |
| Initial time step | $1.0\times10^{-3}$ | $\mathrm{s}$ |
| Maximum time step | $1.0$ | $\mathrm{s}$ |
```xml
```
### 5.2 Nonlinear solver
At each implicit time step, the nonlinear residual is driven to convergence with a Newton method (NOX), using a full (undamped) Newton step and a combined stopping criterion on the residual norm or the iteration count, summarized in Table 6.
**Table 6.** Nonlinear solver (NOX) settings.
| Parameter | Value |
|---|---|
| Nonlinear solver | Line Search Based |
| Direction method | Newton |
| Line search method | Full Step |
| Forcing term method | Constant |
| Convergence test 0 | $\lVert F \rVert$ (NormF) $< 1.0\times10^{-8}$ |
| Convergence test 1 | Maximum nonlinear iterations per time iteration $=10$ |
| Combination | OR (either test satisfied stops the solve) |
```xml
```
### 5.3 Linear solver and preconditioner
Each Newton step requires the solution of a linear system, handled by a Belos GMRES solver preconditioned with an additive Schwarz method (Ifpack2) using a sparse-direct subdomain solver (Amesos2/PARDISO-MKL), as summarized in Table 7.
**Table 7.** Linear solver and preconditioner settings.
| Parameter | Value |
|---|---|
| Linear solver | Belos, Pseudo Block GMRES |
| Convergence tolerance | $1.0\times10^{-8}$ |
| Maximum iterations | $1000$ |
| Preconditioner | Ifpack2, Schwarz |
| Schwarz overlap level | $2$ |
| Subdomain solver | Amesos2, PARDISO-MKL |
```xml
```
### 5.4 Output and monitoring
Solution fields are written to `incompressible_3d_rayleigh_benard_convection_solution.exo` every $2000$ time steps, including cell-averaged gradients of the Lagrange pressure, both velocity components, and temperature. A scalar Nusselt-number response, based on the temperature gradient evaluated on the `bottom` side set, is also recorded to monitor convective heat transfer:
```xml
```
The Nusselt number quantifies the ratio of the total (convective plus conductive) heat transfer to the purely conductive heat transfer across the layer, and is typically evaluated at the wall as
$$
Nu = -\frac{H}{\Delta T}\left.\frac{\partial T}{\partial z}\right|_{wall}
$$
with $H$ the cylinder height and $\Delta T$ the imposed wall-to-wall temperature difference (Section 4.1).
---
## 6. Numerical results
This section reports the expected qualitative outputs of the simulation: the velocity and temperature fields once convection has developed, extracted from the Exodus solution file `incompressible_3d_rayleigh_benard_convection_solution.exo`. The given parameters led to $Ra=5000$ and $Pr=7$, where $Ra$ is Rayleigh number defined as $\beta g(T_h-T_c)h/(\alpha \nu)$ and $Pr$ is Prandtl number defined as $\nu \rho c_p/k$. With the given parameters, simulation is conducted until the convergence and average Nusselt number is calculated for the bottom wall as:
$$
\overline{Nu}=\frac{1}{\pi \Gamma^2}\int_{0}^{2\pi}\int_{0}^{\Gamma}\frac{\partial \Omega}{\partial z}|_{z=0.5}RdRd\phi~~.
$$
The resultant Nusselt numbers compared with the reference value (1.6736) from {cite}`li20223d` is llustrated in Table 8 below.
**Table 8:** Nusselt numbers calculated from the bottom wall using different mesh resolutions.
| Number of Elements | Average Nusselt Number | Percentage Error (%) |
|-------------------:|-----------------------:|---------------------:|
| 1,568 | 1.5131 | 9.59 |
| 4,224 | 1.6058 | 4.05 |
| 9,344 | 1.6424 | 1.86 |
| 18,000 | 1.6615 | 0.72 |
### 6.1 Velocity field
Once the imposed temperature difference exceeds the critical value for the onset of convection, the initially quiescent fluid organizes into buoyancy-driven convective cells..
**Figure 3.** Vertical velocity $w$ obtained with VERTEX-CFD for $Ra=5000$ and $Pr=7$.
### 6.2 Temperature field
The temperature field evolves from the uniform initial condition of $200.0\ \mathrm{K}$ toward a statistically steady, convective distribution bounded by the wall temperatures of $200.0\ \mathrm{K}$ (top) and $250.0\ \mathrm{K}$ (bottom).
**Figure 4.** Temperature contours obtained with VERTEX-CFD for $Ra=5000$ and $Pr=7$.
```{bibliography}