Buoyancy-Driven Flow in a Square Cavity


Buoyancy-driven convection denotes a type of the heat transfer in a fluid, in which the fluid flow is driven solely by a density difference due to a temperature gradient.

Consider a two-dimensional flow within a square cavity of solid walls, where gravity is acting in the direction. The left and right boundaries have different temperature values at and , respectively. The top and bottom boundaries, however, are assumed to be thermally insulated.

As the fluid near the left boundary absorbs heat, it becomes less dense and rises due to thermal expansion. The surrounding cooler fluid then moves in to replace it, forming a convection current within the domain. This type of flow is known as the natural/free convection.

The following simulation models the temperature, pressure and velocity fields within the square cavity. Since a stronger convection flow results in a higher rate of the heat transfer, we can measure the level of the natural convection by calculating the heat flux at the boundary.

The symbols and corresponding units used here are summarized in the Nomenclature section.

Please refer to the information provided in "Heat Transfer" for a more general theoretical background for heat transfer analysis.

Load the finite element package.

Multiphysics Model

Since this problem considers more than one kind of physics, a multiphysics model is to be constructed. The heat equation is coupled to the NavierStokes equation for modeling heat transfer with the fluid flow.

Heat Transfer Model

For a steady-state heat transfer model without sources, the temperature distribution is described by the time-independent source-free heat equation (1). The heat convection by the fluid flow is modeled with the convective term. As such the velocity field coming for the fluid flow is set in the convective term of the heat equation:

Set up a steady-state 2D heat transfer model.

Fluid Dynamics Model

The flow field is described by the steady-state Navier-Stokes equation (2). To account for the temperature effects on the fluid, a Boussinesq approximation [3] is used. This approximation neglects the variations in fluid properties except for the density differences induced by the temperature. In other words, the Boussinesq approximation assumes the buoyancy to be the only driving force acting on the fluid flow.

The buoyancy force appears in the source term in the Navier-Stokes equation:

Construct the left hand side of the 2D fluid dynamics model.


The simulation domain is constructed as a unit square. The left, right, bottom and top boundaries are denoted by , , and , respectively.

Define the 2D domain .

For a natural convection, the buoyancy is the driving force for the fluid motion. The fluid viscosity, on the other hand, is the dominant force that resisted this convection flow.

The dimensionless Rayleigh number [4], is then defined to represents the relationship between the buoyancy force and the viscosity force within the fluid:


is the gravity ,
is the thermal expansion coefficient ,
is the thermal diffusivity ,
is the kinematic viscosity ,
is the free stream temperature .

In the following simulation we will solve and compare the convection flow field under the Rayleigh numbers of and . The remaining flow properties are set as constants with a value of .

Set up the model parameters.

Boundary Conditions

Since this multiphysics simulation contains both the heat transfer model and the fluid dynamics model, the boundary conditions for each physics mode are set up respectively.

Heat Transfer Model

There are two types of the boundary conditions involved in the heat transfer model. At the left and right boundaries , the temperature are held at and , respectively.

Set up temperature surface boundary conditions at the left and right boundaries.

The top and bottom boundaries , are assumed to be perfectly insulated. The default thermally insulated boundary conditions are implicitly applied.

Fluid Dynamics Model

In the fluid dynamics model, the only boundary condition involved is the wall/no-slip boundary condition. At all four boundaries the flow velocity is set to zero to model the solid walls where the fluid has no-slip.

Set up the wall/no-slip boundary condition for all four boundaries.

To solve for the flow pressure within the domain , a reference pressure is required. Here the reference point is chosen at the bottom-left point with the reference pressure .

Set up a reference pressure point.
Construct the boundary conditions for the multiphysics PDE model.

Solve the PDE Model

Recall that the NavierStokes equation is coupled to the heat equation by a Boussinesq approximation [5]. Since the buoyancy force is alined with the axis, the source term in (6) is set up as .

Specify the multiphysics PDE with the model parameters.

A stable solution can be found if the velocities are interpolated with a higher order than the pressure. NDSolve allows an interpolation order for each dependent variable to be specified. Here the velocities and and the temperature are set to be interpolated with second order and the pressure with first order.

Solve the coupled PDE model and monitor the total time and memory used.

Post-processing and Visualization

To inspect the heat transfer within the fluid flow, the following visualization combines the resulting velocity streamlines with the temperature distribution.

Set up a legend bar and ContourPlot options for the visualization.
Visualize the temperature field and the velocity streamlines.

To measure the amount of natural convection, the heat flux is calculated at the boundary . Note that a higher implies a stronger convection flow:

Because the top and bottom boundaries are thermally insulated, the heat flux across the left boundary must equal to the one across the right boundary based on the energy conservation. Therefore, the heat flux can be calculated at either or .

Calculate the heat flux at the left boundary .

In this case the heat flux across the cavity is found as . This value matches with the reference value presented in [7].

Next, we increase the Rayleigh number from to and solve for the convection flow field again.

Update the model parameters for the increasing Rayleigh number .
Set the multiphysics PDE with the modified model parameters.

To efficiently solve the updated heat transfer model, the previous result can be used as an initial guess for the new solution. Since the initial seeding based on the previous solution is already close to the final solution the nonlinear solver will converge in fewer steps, resulting in less time elapsed.

Solve the coupled PDE model.
Compare the flow fields of different Rayleigh number .

When the Rayleigh number increases, the driving force (buoyancy) of the convective flow increases compared to the resisting force (viscosity) of the fluid motion. A stronger convection flow is then formed.

Calculate the heat flux for the updated model.

The heat flux across the cavity has been increased from to , which implies a stronger convective heat transfer within the domain. This updated heat flux value also matches with the reference value presented in [8].



1.  D. de Vahl Davis and I. P. Jones. Natural Convection of Air in a Square Cavity - A Comparison Exercise, Int. J. Numer. Meth. Fluids. vol. 3 (1983).