COMSOL 6.4 - Buoyancy Flow in Air

This example studies the stationary state of free convection in a cavity filled with air and bounded by two vertical plates. To generate the buoyancy flow, the plates are heated at different temperatures, in a range where the flow remains in a laminar regime.

This model is similar to the model Buoyancy Flow in Water, except that water is replaced by air. The detailed model analysis for the estimation of the flow regime in the Model Analysis section is still applicable here by replacing the water properties by the air properties. Only the values of the indicators used to estimate the flow regime are given in this document, in the Indicators of the Flow Regime section.

The simulation is run in simple 2D and 3D geometries.

A first 2D model, representing a square cavity (see Figure 1), focuses on the convective flow.

Figure 1: Domain geometry and boundary conditions for the 2D model (square cavity).

The 3D model (see Figure 2) extends the geometry to a cube. Compared to the 2D model, the front and back sides are additional boundaries that may influence the fluid behavior.

Figure 2: Domain geometry and boundary conditions for the 3D model (cubic cavity).

Both models calculate and compare the velocity field and the temperature field. The predefined Nonisothermal Flow interface available in the Heat Transfer Module provides appropriate tools to fully couple the heat transfer and the fluid dynamics.

Model Definition

2D model

Figure 1 illustrates the 2D model geometry. The fluid fills a square cavity with impermeable walls, so the fluid flows freely within the cavity but cannot leak out. The right and left edges of the cavity are maintained at high and low temperatures, respectively. The upper and lower boundaries are insulated. The temperature differential produces the density variation that drives the buoyant flow.

The compressible Navier–Stokes equation contains a buoyancy term on the right-hand side to account for the lifting force due to thermal expansion that causes the density variations:

In this expression, the dependent variables corresponding to the flow are the fluid velocity vector, u, and the pressure, p. The constant g denotes the gravitational acceleration, ρ gives the temperature-dependent density, and μ is the temperature-dependent dynamic viscosity.

Because the model only contains information about the pressure gradient, it estimates the pressure field up to a constant. To define this constant, you arbitrarily fix the pressure at a point. No slip boundary conditions apply on all boundaries. The no slip condition results in zero velocity at the wall but does not set any constraint on p.

At steady-state, the heat balance for a fluid reduces to the following equation:

Here T represents the temperature, k denotes the thermal conductivity, and Cp is the specific heat capacity of the fluid.

The boundary conditions for the heat transfer interface are the fixed high and low temperatures on the vertical walls with insulation conditions elsewhere, as shown in Figure 1.

3D model

Figure 2 shows the geometry and boundary conditions of the 3D model. The cavity is now a cube with high and low temperatures applied respectively at the right and left surfaces. The remaining boundaries (top, bottom, front, and back) are thermally insulated.

Indicators of the Flow Regime

Before starting the simulations, it is recommended to estimate the flow regime. The Grashof number, Rayleigh number, and Prandtl number are then computed and summarized in Table 1 below.

Table 1: Indicators of the Flow Regime.
Indicator	Description	Value
Gr	Grashof number	1.6·106
Ra	Rayleigh number	1.1·106
Pr	Prandtl number	0.7

See the model analysis in Buoyancy Flow in Water for more details about the definition of these indicators.

The Grashof and Rayleigh numbers are about 106, and thus significantly below the critical value of 109 for the flow to become turbulent between vertical plates (see Ref. 1). A laminar regime is therefore expected.

In this regime, the typical velocity due to buoyancy can be estimated by

where αp is the coefficient of thermal expansion and L the typical length. This estimate gives U0 = 0.2 m/s when ΔT = 10 K.

The Prandtl number is lower than 1, which means that the viscous forces should not have any significant limiting effect on buoyancy in this model (see Ref. 2). The estimate of the typical velocity due to buoyancy, taking into account the limiting effect of the viscous forces on buoyancy, gives U1 = 0.1 m/s, which is close to U0.

Since the Prandtl number is still close to 1, the thermal boundary layer thickness, δT, and the shear layer thickness, δs, are of the same order of magnitude. They are expressed as (Ref. 2):

For ΔT = 10 K, these estimates give δT and δs of about 3 mm, with δs slightly smaller than δT. You can use these estimates to validate the model after computing the simulation.

The thermophysical properties of air used to compute these indicators are taken from the Material Library. They are listed in Table 2.

Table 2: thermophysical properties for air at 288.15 K.
Parameter	Description	Value
ρ	Density	1.225 kg/m3
μ	Dynamic viscosity	1.789·10-5 N·s/m2
k	Thermal conductivity	0.02537 W/(m·K)
Cp	Heat capacity at constant pressure	1.005 kJ/(kg·K)
αp	Coefficient of thermal expansion	3.47·10-3 K-1

These properties are given at 288.15 K which corresponds to the average of the hot and cold plates temperatures. In addition, the coefficient of thermal expansion αp is computed using the ideal gas approximation:

By using the average of the hot and cold plate temperatures in this approximation, you get αp = 3.47·10−3 K−1.

Results and Discussion

2D model

Figure 3 shows the velocity distribution in the square cavity.

Figure 3: Velocity magnitude for the 2D model.

The regions with faster velocities are located at the lateral boundaries. The maximum velocity is 0.05 m/s which is a bit lower than the estimated typical velocity U0 = 0.2 m/s, but still in the same order of magnitude. According to Figure 4, the shear layer thickness is about 10 mm.

Figure 4: Velocity profile at the left boundary.

Figure 5 shows the temperature field (surface) and velocity field (arrows) of the 2D model.

Figure 5: Temperature field (surface plot) and velocity (arrows) for the 2D model.

A large convective cell occupies the whole square. The fluid flow follows the boundaries. As seen in Figure 3, it is faster at the vertical plates where the highest variations of temperature are located. The thermal boundary layer is of the order 10 mm according to Figure 6, which is in agreement with the order of the estimated value of 3 mm.

Figure 6: Temperature profile at the left boundary.

3D Model

Figure 7 illustrates the velocity plot parallel to the heated plates.

Figure 7: Velocity magnitude field for the 3D model, slices parallel to the heated plates.

A second velocity magnitude field is shown in Figure 8. The plot is close to what was obtained in 2D in Figure 3.

Figure 8: Velocity magnitude field for the 3D model, slices perpendicular to the heated plates.

In Figure 9, velocity arrows are plotted on temperature surface at the middle vertical plane parallel to the plates.

Figure 9: Temperature (surface plot) and velocity (arrows) fields in the cubic cavity, for a temperature difference of 10 K between the vertical plates.

Notes About the COMSOL Implementation

The material properties for air are available in the Material Library.

At high Gr values, using a good initial condition becomes important in order to achieve convergence. Moreover, a well-tuned mesh is needed to capture the solution, especially the temperature and velocity changes near the walls. Use the Stationary study step’s continuation option with ΔT as the continuation parameter to get a solver sequence that uses previous solutions to estimate the initial condition. For this tutorial, it is appropriate to ramp up ΔT from 10−2 K to 10 K, which corresponds to a Grashof number range of 103–106. At Gr = 103, the model is easy to solve. The regime is dominated by conduction. At Gr = 106, the model becomes more difficult to solve. The regime is more influenced by convection and buoyancy.

To get a well-tuned mesh when Gr reaches 106, the element size near the prescribed temperature boundaries has to be smaller than the shear layer and thermal boundary layer thicknesses, which are of order 3 mm. It is recommended to have three to five elements across the layers when using P1 elements (the default setting for fluid flows).

References

1. F.P. Incropera, D.P. DeWitt, T.L. Bergman, and A.S. Lavine, Fundamentals of Heat and Mass Transfer, 6th ed., John Wiley & Sons, 2006.

2. A. Bejan, Heat Transfer, John Wiley & Sons, 1993.

Heat_Transfer_Module/Tutorials,_Forced_and_Natural_Convection/buoyancy_air

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
L	10[cm]	0.1 m	Square side length
W	L/2	0.05 m	Volume thickness
DeltaT	10[K]	10 K	Temperature difference
Tc	283.15[K]	283.15 K	Low temperature
Th	Tc+DeltaT	293.15 K	High temperature
T_avg	(Tc+Th)/2	288.15 K	Average temperature approximation
p_ref	1[atm]	1.0133E5 Pa	Reference pressure
p0	0[Pa]	0 Pa	Relative pressure constraint
dt_range_start	-2	-2	Range start point of the parametric sweep on DeltaT (log10 of the value in K)
dt_range_stop	1	1	Range end point of the parametric sweep on DeltaT (log10 of the value in K)
rho	1.225[kg/m^3]	1.225 kg/m³	Density at T_avg and 1 atm
mu	1.789e-5[N*s/m^2]	1.789E-5 N·s/(m·m)	Dynamic viscosity at T_avg
k	0.02537[W/(m*K)]	0.02537 W/(m·K)	Thermal conductivity at T_avg
Cp	1.005[kJ/(kg*K)]	1005 J/(kg·K)	Heat capacity at T_avg
alpha	1/T_avg	0.0034704 1/K	Coefficient of thermal expansion at T_avg
U0	sqrt(g_constalphaDeltaT*L)	0.18448 m/s	Typical velocity due to buoyancy
U1	U0/sqrt(Pr)	0.21914 m/s	Typical velocity due to buoyancy
Pr	mu*Cp/k	0.70869	Prandtl number
Gr	(U0rhoL/mu)^2	1.5957E6	Grashof number
Ra	Pr*Gr	1.1309E6	Rayleigh number
eps_t	L/(Pr*Ra)^0.25	0.0033422 m	Thermal boundary layer thickness
eps_s	eps_t*sqrt(Pr)	0.0028136 m	Shear layer thickness