FEATool Multiphysics
v1.10 Finite Element Analysis Toolbox |

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 ^{3}*, respectively, after which the

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.

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

**Navier-Stokes Equations**physics mode from the*Select Physics*drop-down menu.- Press
**OK**to finish the physics mode selection. - Select
**Rectangle**from the*Geometry*menu. - Press
**OK**to finish and close the dialog box. - Switch to
**Grid**mode by clicking on the corresponding*Mode Toolbar*button. Enter

`0.03`

into the*Grid Size*edit field, and press the**Generate**button to call the grid generation algorithm.- 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.

Enter

`Pr*Ra*T`

into the*Volume force in y-direction*edit field.- Switch to the
**+**tab.

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

Select the

**Heat Transfer**physics mode from the*Select Physics*drop-down menu.- Press the
**Add Physics >>>**button. - Enter
`u`

into the*Convection velocity in x-direction*edit field. Enter

`v`

into the*Convection velocity in y-direction*edit field.- Press
**OK**to finish the equation and subdomain settings specification. - 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 |

- 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.

- Switch to the
**ht**tab. - 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.

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

`1`

into the*Temperature*edit field.- Select
**2**in the*Boundaries*list box. Select

**Temperature**from the*Heat Transfer*drop-down menu.- Press
**OK**to finish the boundary condition specification. - Switch to
**Solve**mode by clicking on the corresponding*Mode Toolbar*button. - 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.

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

**Temperature, T**from the*Predefined contour plot expressions*drop-down menu.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*.

- Select
**Boundary Integration...**from the*Post*menu. - Select
**2**and**4**in the*Boundaries*list box. - Enter
`abs(Tx)/2`

into the*Integration Expression*edit field. Press

**OK**to finish and close the dialog box.

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

- Select
**Model Constants and Expressions...**from the*Equation*menu. Enter

`1e4`

into the value field for*Ra*.- Press
**OK**to finish and close the dialog box. - 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.

- Press the
**Settings***Toolbar*button. - Select the
**Computed solution**radio button. - Enter
`100`

into the*Maximum number of non-linear iterations*edit field. Enter

`0.8`

into the*Non-linear relaxation parameter (ratio of new to old solution to use)*edit field.- 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 ^{4}*.

- Select
**Boundary Integration...**from the*Post*menu. - Select
**2**and**4**in the*Boundaries*list box. Enter

`abs(Tx)/2`

into the*Integration Expression*edit field.- 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, or GUI script (.fes) file.

[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.