Natural Convection in a Square Cavity

This multiphysics model illustrates natural convection effects in a unit square domain using the Boussinesq approximation. The model involves a Navier-Stokes equations physics mode, representing the fluid flow with solid wall or no-slip boundary conditions everywhere. In addition a heat transfer physics mode is added to model the temperature field. The top and bottom boundaries are perfectly insulated while the left boundary is prescribed a temperature of 1 and the right zero.

The physics modes are two way coupled through the vertical source term in the Navier-Stokes equations, Pr·Ra·T, and the velocities transporting the temperature coming directly from the fluid flow. First, the Prandtl and Rayleigh numbers are set to Pr = 0.71 and Ra = 10³, respectively, after which the Ra number will be increased to 10⁴. The references contain benchmark reference and comparison results for a number of quantities such as maximum velocities and the Nusselt number [1,2].

Tutorial

This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Multiphysics > Natural Convection in a Square Cavity from the File menu. Or alternatively, follow the step-by-step instructions below. Note that the CFDTool interface differ slightly from the FEATool Multiphysics instructions described in the following.

1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.

2. Select the Navier-Stokes Equations physics mode from the Select Physics drop-down menu. (Note that for CFDTool the physics selection is done in the Equation settings dialog box.)

   
3. Press OK to finish the physics mode selection.
4. Select Rectangle from the Geometry menu.
5. Press OK to finish and close the dialog box.
6. Switch to Grid mode by clicking on the corresponding Mode Toolbar button.

7. Enter 0.03 into the Grid Size edit field, and press the Generate button to call the grid generation algorithm.

   
8. Switch to Equation mode by clicking on the corresponding Mode Toolbar button.

The non-dimensionalized form of the Boussinesq source term Pr·Ra·T couples the temperature to the y-direction source term of the Navier-Stokes equations. (Note that the Equation Settings dialog box may look different for CFDTool.)

1. Enter Pr*Ra*T into the Volume force in y-direction edit field.

   
2. Switch to the + tab.

Add the heat transfer physics mode and coupled the velocities u and v to the convective transport terms.

1. Select the Heat Transfer physics mode from the Select Physics drop-down menu.

   
2. Press the Add Physics >>> button.
3. Enter u into the Convection velocity in x-direction edit field.

4. Enter v into the Convection velocity in y-direction edit field.

   
5. Press OK to finish the equation and subdomain settings specification.
6. Press the Constants Toolbar button, or select the corresponding entry from the Equation menu, and enter the following values for the Prandtl and Rayleigh numbers in the Model Constants and Expressions dialog box.

Name	Expression
Pr	0.71
Ra	1e3

1. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.

Switch to the ns tab, which corresponds to the boundary conditions for the Navier-Stokes equations physics mode. Then select the Wall/no-slip for all four boundaries.

Click on the ht tab to change to specifying boundary conditions for the heat transfer physics mode. Select Thermal insulation/symmetry conditions for the top and bottom boundaries.

1. Switch to the ht tab.
2. Select 1 and 3 in the Boundaries list box.

Select a Temperature boundary conditions for the left and right boundaries, and set a fixed temperature equal to 1 at the left side.

1. Select Thermal insulation/symmetry from the Heat Transfer drop-down menu.
2. Select 4 in the Boundaries list box.
3. Select Temperature from the Heat Transfer drop-down menu.

4. Enter 1 into the Temperature edit field.

   
5. Select 2 in the Boundaries list box.

6. Select Temperature from the Heat Transfer drop-down menu.

   
7. Press OK to finish the boundary condition specification.
8. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.
9. Press the = Toolbar button to call the solver. After the problem has been solved FEATool will automatically switch to postprocessing mode and plot the computed solution.

The solution shows how the temperature difference is causing a vortex in the flow field, which in turn causes an offset in the temperature field.

Change the plot to visualize the temperature field as surface and contour plots, and the velocity field as an arrow plot.

1. Press the Plot Options Toolbar button.
2. Select the Arrow Plot check box.
3. Select Temperature, T from the Predefined surface plot expressions drop-down menu.
4. Select the Contour Plot check box.

5. Select Temperature, T from the Predefined contour plot expressions drop-down menu.

   

6. Press OK to plot and visualize the selected postprocessing options.

   

Use the boundary integration postprocessing tool to calculate the average Nusselt number for the vertical boundaries, and compare it to the reference value of 1.118.

1. Select Boundary Integration... from the Post menu.
2. Select 2 and 4 in the Boundaries list box.
3. Enter abs(Tx)/2 into the Integration Expression edit field.

4. Press OK to finish and close the dialog box.

   

The simulation will now be repeated but with an increased Rayleigh number Ra = 10⁴. Instead of starting over from the beginning the existing solution will be used as a starting guess which helps with non-linear convergence.

1. Select Model Constants and Expressions... from the Equation menu.

2. Enter 1e4 into the value field for Ra.

   
3. Press OK to finish and close the dialog box.
4. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.

To assist with convergence, as well as selecting the old solution as initial value, also increase the non-linear relaxation and the maximum number of non-linear iterations.

1. Press the Settings Toolbar button.
2. Select the Computed solution radio button.
3. Enter 100 into the Maximum number of non-linear iterations edit field.

4. Enter 0.8 into the Non-linear relaxation parameter (ratio of new to old solution to use) edit field.

   
5. Press the Solve button.

It is evident that the increase in Rayleigh number causes a significantly stronger rotation and stretching of the temperature field.

Use boundary integration tool again to calculate the mean Nusselt number for the vertical boundaries, and compare it to the corresponding reference value of 2.243 for Ra = 10⁴.

1. Select Boundary Integration... from the Post menu.
2. Select 2 and 4 in the Boundaries list box.

3. Enter abs(Tx)/2 into the Integration Expression edit field.

   
4. Press OK to finish and close the dialog box.

The natural convection in a square cavity multiphysics model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB m-script text file (available as the example ex_natural_convection script file), or GUI script (.fes) file.

References

[1] de Vahl D. Natural Convection of Air in a Square Cavity: A Benchmark Solution. Int. J. Numer. Meth. Fluids, vol. 3, pp. 249-264, 1983.

[2] de Vahl D, Jones IP. Natural Convection of Air in a Square Cavity: A Comparison Exercise. Int. J. Numer. Meth. Fluids, vol. 3, pp. 227-248, 1983.