PDF
Buoyancy Flow in Water
Introduction
This example studies the stationary state of free convection in a cavity filled with water and bounded by two vertical plates. To generate the buoyancy flow, the plates are heated at different temperatures, bringing the regime close to the transition between laminar and turbulent flow.
An important step in setting up a convection model is to assess whether the flow stays laminar or becomes turbulent. It is also important to approximate how fine should be the mesh needed to resolve velocity and temperature gradients. Both of these approximations rely on the velocity scale of the model. This makes the setup of natural convection problems challenging because the resulting velocity is part of the nonlinear solution. There are, however, tools to approximate scales for natural convection problems. These tools are demonstrated in this application using 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 for 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.
As the density changes are relatively small in this model, it is possible to use a simplified version of Navier-Stokes equation based on Boussinesq approximation:
where Tref is the reference temperature, ρref is the reference density, and αp is the isobaric compressibility coefficient evaluated at the reference pressure and temperature.
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.
Model Analysis
Before starting the simulations, it is recommended to estimate the flow regime. To this end, four indicators are presented: the Reynolds, Grashof, Rayleigh, and Prandtl numbers. They are calculated using the thermophysical properties of water listed in Table 1, found in the Water, liquid material of the Built-in material library. The thermophysical properties are given at 283 K which is in the range of the temperatures observed in the model.
Table 1: thermophysical properties for water at 283 K.
Parameter
Description
Value
 ρ
Density
1.0·103 kg/m3
 μ
Dynamic viscosity
1.3·10-3 N·s/m2
 k
Thermal conductivity
0.58 W/(m·K)
 Cp
Heat capacity at constant pressure
4.2 kJ/(kg·K)
 αp
Coefficient of thermal expansion
9·10-5 K-1
Prandtl Number
Definition
The Prandtl number is the ratio of fluid viscosity to thermal diffusivity. It is defined by:
At 283 K, Pr = 0.72 for air and Pr = 9.4 for water, as given in Ref. 2.
Boundary Layers
The Prandtl number also indicates the relative thickness of the outer boundary layer, δ, and the thermal boundary layer, δT. In the present case, it is reasonable to estimate the ratio δ ⁄ δT by the relation (7.32 in Ref. 2)
(1)
The outer boundary layer is the distance from the wall to the region where the fluid stabilizes. It is different from the momentum boundary layer, δM, which measures the distance from the wall to the velocity peak.
Application in this Tutorial
With the values given in Table 1, the Prandtl number for water at 283 K, is found to be of the order 1 or 10. According to Equation 1, δ and δT should then be of same order of magnitude.
Reynolds Number
Definition
The Reynolds number estimates the ratio of inertial forces to viscous forces. It is defined by the formula
where U denotes the typical velocity and L the typical length.
At atmospheric pressure and at 283 K, air and water have the following properties (Ref. 2).
Table 2: thermophysical properties of air and water at 283K and atmospheric pressure.
Parameter
Description
Air
Water
 ρ
Density
1.25 kg/m3
1000 kg/m3
 μ
Dynamic viscosity
1.76·10-5 N·s/m2
1.3·10-3 N·s/m2
You can thus approximate the Reynolds number by:
In these relations, U has to be in meters per second and L in meters.
The Reynolds number can be rewritten as the velocity ratio
which compares U to μ ⁄ (ρL). The latter quantity is homogeneous to a velocity and can be seen as the typical velocity due to viscous forces.
Flow Regime
The value of the Reynolds number is used to predict the flow regime. Generally, low values of Re correspond to laminar flow and high values to turbulent flow, with a critical value for the transition regime that depends on the geometry.
As an indication, Reynolds’ experiments concerning the flow along a straight and smooth pipe showed that the transition regime in this case occurs when Re is between 2000 and 104 (see chapter 1.3 in Ref. 3).
Momentum Boundary Layer
The momentum boundary layer thickness can be evaluated, using the Reynolds number, by (5.36 in Ref. 2)
(2)
Application in this Tutorial
The typical length L of the model is equal to 10 cm so the Reynolds number is evaluated as
where U is still unknown. Estimates of this typical velocity are provided later.
Grashof Number
Definition
The Grashof number gives the ratio of buoyant to viscous forces. It is defined by
where g is the gravity acceleration (equal to 9.81 m.s−2) and ΔT is the typical temperature difference.
At atmospheric pressure and at 283 K, air and water have the following properties (Ref. 2).
Table 3: thermophysical properties of air and water at 283 K and atmospheric pressure.
Parameter
Description
Air
Water
 ρ
Density
1.25 kg/m3
1000 kg/m3
 μ
Dynamic viscosity
1.76·10-5 N.s.m-2
1.3·10-3 N·s/m2
 αp
Coefficient of thermal expansion
3.55·10-3 K-1
9·10-5 K-1
The value of αp for air was here obtained by the ideal gas approximation:
You can then evaluate the Grashof number by:
In these relations, ΔT is given in kelvins and L in meters.
The Grashof number can also be expressed as the velocity ratio
where U0 is defined by
(3)
This quantity can be considered as the typical velocity due to buoyancy forces.
Flow Regime
When buoyancy forces are large compared to viscous forces, the regime is turbulent; otherwise it is laminar. The transition between these two regimes is indicated by the critical order of the Grashof number which is 109 (see Figure 7.7 in Ref. 2).
Application in this Tutorial
In this tutorial, ΔT is equal to 10 K so the Grashof number is about 107 which indicates that a laminar regime is expected.
Table 4 provides the values of the quantities necessary to calculate U0. This velocity is here of order 10 mm/s.
Table 4: thermophysical properties of water at 283 K used in the Grashof number.
Parameter
Description
Value
 g
Gravitational acceleration
9.81 m/s2
 αp
Coefficient of thermal expansion
9·10-5 K-1
Rayleigh Number
Definition
The Rayleigh number is another indicator of the regime. It is defined by
so it is similar to the Grashof number except that it accounts for the thermal diffusivity, k, given by
Note: The Rayleigh number can be expressed in terms of the Prandtl and the Grashof numbers through the relation Ra = Pr Gr.
At atmospheric pressure and at 300 K, you can use the approximations of Ra below for air and water (A.4 in Ref. 1)
In these relations, ΔT is given in kelvins and L in meters.
Using Equation 3, the Rayleigh number can be rewritten as the velocity ratio
where the ratio κ ⁄ L can be seen as a typical velocity due to thermal diffusion.
Flow Regime
Like the Grashof number, a critical Rayleigh value indicates the transition between laminar and turbulent flow. For vertical plates, this limit is about 109 (9.23 in Ref. 1).
Typical Velocity
Because the viscous forces limit the effects of buoyancy, U0 may give an overestimated typical velocity. Another approach (see 7.25 in Ref. 2) is to use U1 instead, defined by
that is
(4)
or
This should be a more accurate estimate of U because the fluid’s thermal diffusivity and viscosity are used in the calculations. From now on, U1 is the preferred estimate of U.
Thermal Boundary Layer
The Rayleigh number can be used to estimate the thermal boundary layer thickness, δT. When Pr is of order 1 or greater, it is approximated by the formula (7.25c in Ref. 2)
(5)
Application in this Tutorial
Here, Ra is of order 108. The laminar regime is confirmed but the Rayleigh number found is near the transition zone. The thermal boundary layer thickness is then found to be of order 1 mm and U1 of order 10 mm/s.
Synthesis
To prepare the simulation, it is very useful to follow the steps below that give indications of what results to expect. It is important to remember that the quantities computed here are only order of magnitude estimates, which should not be considered with more than one significant digit.
First evaluate the Grashof and Rayleigh numbers. If they are significantly below the critical order of 109, the regime is laminar. In this case, Equation 3 or Equation 4 provide estimates of the typical velocity U that you can use to validate the model after performing the simulation.
According to Equation 1, the Prandtl number determines the relative thickness of the thermal boundary layer and the outer layer. Equation 2 and Equation 5 then provide orders of magnitude of the thicknesses. When defining the mesh, refinements have to be done at the boundary layers by, for instance, inserting three to five elements across the estimated thicknesses.
Here, Gr and Ra are 107 and 108, respectively, and thus below the critical value of 109 for vertical plates. A laminar regime is therefore expected but because these values are not significantly below 109, convergence is not straightforward. In this regime, the estimates U0 and U1 of the typical velocity are both of the order 10 mm/s.
For water at 283 K, Pr is about 10 so δ and δT are of same orders of magnitude. Here, δT is of the order 1 mm.
The Reynolds number calculated with U1 is about 103, which confirms that the model is close to the transition regime. Using U1 and Equation 2, the momentum boundary layer thickness δM is found to be about 1 mm.
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 4 mm/s which is in agreement with the estimated typical velocity U1 of order 10 mm/s. According to Figure 4, the momentum boundary layer thickness is of order 1 mm, as calculated previously.
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 1 mm according to Figure 6, which is in agreement with the estimate provided by Equation 5. The outer layer is slightly thicker than the boundary layer.
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.
New small convective cells appear on the vertical planes perpendicular to the plates at the four corners. They are more visible at lower Gr values, that is, far from the transition regime. In Figure 10, the temperature difference between the vertical plates is reduced to 1 K and 0.1 K to decrease the Grashof number to 105 and 106.
Observe the bigger cells at the four corners of the plane.
Figure 10: Temperature (surface plot) and velocity (arrows) fields in the cubic cavity, with, for a temperature difference of 1 K (top) and 0.1 K (bottom) between the vertical plates.
Notes About the COMSOL Implementation
The material properties for water are available in the Material Library. The density and dynamic viscosity are functions of the temperature.
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-3 K to 10 K, which corresponds to a Grashof number range of 103–107. At Gr = 103, the model is easy to solve. The regime is dominated by conduction and viscous effects. At Gr = 107, the model becomes more difficult to solve. The regime is greatly influenced by convection and buoyancy.
To get a well-tuned mesh when Gr reaches 107, the element size near the prescribed temperature boundaries has to be smaller than the momentum and thermal boundary layer thicknesses, which are of order 1 mm according to Equation 2 and Equation 5. 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, 1985.
3. P. A. Davidson, Turbulence: An Introduction for Scientists and Engineers, Oxford University Press, 2004.
Application Library path: Heat_Transfer_Module/Tutorials,_Forced_and_Natural_Convection/buoyancy_water
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Fluid Flow>Nonisothermal Flow>Laminar Flow.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Stationary.
6
Click  Done.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
L
10[cm]
0.1 m
Square side length
DeltaT
10[K]
10 K
Temperature difference
Tc
283.15[K]
283.15 K
Low temperature
Th
Tc+DeltaT
293.15 K
High temperature
rho
1000[kg/m^3]
1000 kg/m³
Density
mu
1.3e-3[N*s/m^2]
0.0013 N·s/(m·m)
Dynamic viscosity
k
0.58[W/(m*K)]
0.58 W/(m·K)
Thermal conductivity
Cp
4.2[kJ/(kg*K)]
4200 J/(kg·K)
Heat capacity
alpha
9e-5[1/K]
9E-5 1/K
Coefficient of thermal expansion
U0
sqrt(g_const*alpha*DeltaT*L)
0.029709 m/s
Typical velocity due to buoyancy
U1
U0/sqrt(Pr)
0.0096828 m/s
Typical velocity estimation
Pr
mu*Cp/k
9.4138
Prandtl number
Gr
(U0*rho*L/mu)^2
5.2225E6
Grashof number
Ra
Pr*Gr
4.9163E7
Rayleigh number
Re0
rho*U0*L/mu
2285.3
Reynolds number approximation with U0
Re1
rho*U1*L/mu
744.83
Reynolds number approximation with U1
eps_t
L/(Ra)^0.25
0.0011942 m
Thermal boundary layer thickness
eps_m
L/sqrt(Re1)
0.0036641 m
Momentum boundary layer thickness
The Grashof and Rayleigh numbers should be less than 109, indicating that a laminar regime is expected.
Geometry 1
Square 1 (sq1)
1
In the Geometry toolbar, click  Square.
2
In the Settings window for Square, locate the Size section.
3
In the Side length text field, type L.
4
In the Geometry toolbar, click  Build All.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>Water, liquid.
4
Click Add to Component in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Laminar Flow (spf)
In order to ensure mass conservation, as the volume is constant, the water density cannot depend only on the temperature. It has to be either constant or pressure and temperature dependent. Select the Incompressible flow option to define a constant density evaluated from the material properties at the reference pressure and temperature.
1
In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
2
In the Settings window for Laminar Flow, locate the Physical Model section.
3
From the Compressibility list, choose Incompressible flow.
4
Select the Include gravity check box.
Pressure Point Constraint 1
1
In the Physics toolbar, click  Points and choose Pressure Point Constraint.
2
Select Point 2 only.
Fixing the pressure at an arbitrary point is necessary to define a well-posed model.
Heat Transfer in Fluids (ht)
1
In the Model Builder window, under Component 1 (comp1) click Heat Transfer in Fluids (ht).
2
In the Settings window for Heat Transfer in Fluids, locate the Physical Model section.
3
In the Tref text field, type (Tc+Th)/2.
Define the initial temperature as the mean value between the high and low temperature values.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1)>Heat Transfer in Fluids (ht) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type (Tc+Th)/2.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 1 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type Tc.
Temperature 2
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 4 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type Th.
Now, modify the default mesh size settings to ensure that the mesh satisfies the criterion discussed in the Introduction section.
Multiphysics
Nonisothermal Flow 1 (nitf1)
1
In the Model Builder window, under Component 1 (comp1)>Multiphysics click Nonisothermal Flow 1 (nitf1).
2
In the Settings window for Nonisothermal Flow, locate the Material Properties section.
3
Select the Boussinesq approximation check box.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Extra fine.
4
Click  Build All.
Study 1
Because the Grashof number is near the critical value of around 109, the model is highly nonlinear. To achieve convergence, use continuation to ramp up the temperature difference value from 10−3 K to 10 K, which corresponds to a Grashof number from 103 to 107.
Step 1: Stationary
1
In the Model Builder window, under Study 1 click Step 1: Stationary.
2
In the Settings window for Stationary, click to expand the Study Extensions section.
3
Select the Auxiliary sweep check box.
4
Click  Add.
5
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
DeltaT (Temperature difference)
1e-3 1e-2 1e-1 1 10
K
6
Click the  Show More Options button in the Model Builder toolbar.
7
In the Show More Options dialog box, in the tree, select the check box for the node Physics>Advanced Physics Options.
8
Click OK.
The pseudo time stepping option is generally useful to help the convergence of a stationary flow model. However, a continuation approach is already used here. In this precise model, disabling the pseudo time stepping option improves the convergence. Follow the instructions below to do so.
Laminar Flow (spf)
1
In the Model Builder window, under Component 1 (comp1) click Laminar Flow (spf).
2
In the Settings window for Laminar Flow, click to expand the Advanced Settings section.
3
Find the Pseudo time stepping subsection. From the Use pseudo time stepping for stationary equation form list, choose Off.
Study 1
1
In the Home toolbar, click  Compute.
The first default plot group shows the velocity magnitude as in Figure 3. Notice the high velocities near the lateral walls due to buoyancy effects.
Results
Temperature (ht)
The third default plot shows the temperature distribution. Add arrows of the velocity field to see the correlations between velocity and temperature, as in Figure 5.
1
In the Model Builder window, click Temperature (ht).
2
In the Temperature (ht) toolbar, click  Plot.
Arrow Surface 1
1
In the Temperature (ht) toolbar, click  Arrow Surface.
2
In the Settings window for Arrow Surface, locate the Coloring and Style section.
3
From the Color list, choose Black.
4
In the Temperature (ht) toolbar, click  Plot.
In the following steps, the temperature and velocity profiles are plotted near the left boundary in order to estimate the boundary layer thicknesses of the solution.
Cut Line 2D 1
1
In the Results toolbar, click  Cut Line 2D.
2
In the Settings window for Cut Line 2D, locate the Line Data section.
3
In row Point 1, set y to 5[cm].
4
In row Point 2, set x to 1[cm], and y to 5[cm].
5
Click  Plot.
Temperature at Boundary Layer
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Temperature at Boundary Layer in the Label text field.
3
Locate the Data section. From the Dataset list, choose Cut Line 2D 1.
4
From the Parameter selection (DeltaT) list, choose Last.
Line Graph 1
1
In the Temperature at Boundary Layer toolbar, click  Line Graph.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Heat Transfer in Fluids>Temperature>T - Temperature - K.
3
In the Temperature at Boundary Layer toolbar, click  Plot.
The thermal boundary layer is around 3 mm.
Velocity at Boundary Layer
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Velocity at Boundary Layer in the Label text field.
3
Locate the Data section. From the Dataset list, choose Cut Line 2D 1.
4
From the Parameter selection (DeltaT) list, choose Last.
Line Graph 1
1
In the Velocity at Boundary Layer toolbar, click  Line Graph.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Laminar Flow>Velocity and pressure>spf.U - Velocity magnitude - m/s.
3
In the Velocity at Boundary Layer toolbar, click  Plot.
The momentum boundary layer is around 1 mm and the outer layer between 5 mm and 10 mm.
Now create the 3D version of the model.
Add Component
In the Model Builder window, right-click the root node and choose Add Component>3D.
Add Physics
1
In the Home toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Fluid Flow>Nonisothermal Flow>Laminar Flow.
4
Find the Physics interfaces in study subsection. In the table, clear the Solve check box for Study 1.
5
Click Add to Component 2 in the window toolbar.
6
In the Home toolbar, click  Add Physics to close the Add Physics window.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Stationary.
4
Find the Physics interfaces in study subsection. In the table, clear the Solve check boxes for Laminar Flow (spf) and Heat Transfer in Fluids (ht).
5
Find the Multiphysics couplings in study subsection. In the table, clear the Solve check box for Nonisothermal Flow 1 (nitf1).
6
Click Add Study in the window toolbar.
7
In the Model Builder window, click the root node.
8
In the Home toolbar, click  Add Study to close the Add Study window.
Geometry 2
In the Model Builder window, under Component 2 (comp2) click Geometry 2.
Block 1 (blk1)
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type L.
4
In the Depth text field, type L/2.
5
In the Height text field, type L.
6
In the Geometry toolbar, click  Build All.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>Water, liquid.
4
Click Add to Component in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Laminar Flow 2 (spf2)
In order to ensure mass conservation, as the volume is constant, the water density cannot depend only on the temperature. It has to be either constant or pressure and temperature dependent. Select the Incompressible flow option to define a constant density evaluated from the material properties at the reference pressure and temperature.
1
In the Model Builder window, under Component 2 (comp2) click Laminar Flow 2 (spf2).
2
In the Settings window for Laminar Flow, locate the Physical Model section.
3
From the Compressibility list, choose Incompressible flow.
4
Select the Include gravity check box.
Pressure Point Constraint 1
1
In the Physics toolbar, click  Points and choose Pressure Point Constraint.
2
Select Point 4 only.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Select Boundary 2 only.
Heat Transfer in Fluids 2 (ht2)
1
In the Model Builder window, under Component 2 (comp2) click Heat Transfer in Fluids 2 (ht2).
2
In the Settings window for Heat Transfer in Fluids, locate the Physical Model section.
3
In the Tref text field, type (Tc+Th)/2.
Initial Values 1
1
In the Model Builder window, under Component 2 (comp2)>Heat Transfer in Fluids 2 (ht2) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T2 text field, type (Tc+Th)/2.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 1 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type Tc.
Temperature 2
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 6 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type Th.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Select Boundary 2 only.
Laminar Flow 2 (spf2)
1
In the Model Builder window, under Component 2 (comp2) click Laminar Flow 2 (spf2).
2
In the Settings window for Laminar Flow, click to expand the Advanced Settings section.
3
Find the Pseudo time stepping subsection. From the Use pseudo time stepping for stationary equation form list, choose Off.
Multiphysics
Nonisothermal Flow 2 (nitf2)
1
In the Model Builder window, under Component 2 (comp2)>Multiphysics click Nonisothermal Flow 2 (nitf2).
2
In the Settings window for Nonisothermal Flow, locate the Material Properties section.
3
Select the Boussinesq approximation check box.
Mesh 2
To obtain reliable results within in a reasonable computing time, create a structured mesh by following the steps below.
Mapped 1
1
In the Mesh toolbar, click  Boundary and choose Mapped.
2
Select Boundary 2 only.
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
Select Edges 1, 3, 5, and 9 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type 16.
6
In the Element ratio text field, type 3.
7
Select the Symmetric distribution check box.
8
Click  Build Selected.
The front face mesh has smaller elements near the edges because large variations in velocity and temperature are expected there.
Now extend the front mesh to the remaining structure.
Swept 1
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
From the Distribution type list, choose Predefined.
4
In the Number of elements text field, type 8.
5
In the Element ratio text field, type 3.
6
Select the Reverse direction check box.
To resolve the boundary layers, use a Boundary Layers feature to generate smaller mesh elements near the walls. The thermal boundary layer for the temperature difference of 10 K is approximately 1 mm (see the parameter eps_t defined previously). Use this value to define the thickness of the boundary layers.
Boundary Layers 1
In the Mesh toolbar, click  Boundary Layers.
Boundary Layer Properties
1
In the Model Builder window, expand the Boundary Layers 1 node, then click Boundary Layer Properties.
2
Select Boundaries 1 and 3–6 only.
3
In the Settings window for Boundary Layer Properties, locate the Boundary Layer Properties section.
4
In the Number of boundary layers text field, type 5.
5
From the Thickness of first layer list, choose Manual.
6
In the Thickness text field, type 1[mm]/5.
7
Click  Build All.
Finer mesh is needed for accurate results with higher Rayleigh numbers.
Mesh 3
In the Mesh toolbar, click  Add Mesh.
Reference 1
1
In the Mesh toolbar, click  Modify and choose Mesh>Reference.
2
In the Settings window for Reference, locate the Reference section.
3
From the Mesh list, choose Mesh 2.
Refine 1
1
In the Mesh toolbar, click  Modify and choose Elements>Refine.
2
In the Settings window for Refine, locate the Refine Options section.
3
From the Refinement method list, choose Regular refinement.
4
Click  Build All.
First step solves for the smallest Rayleigh numbers using the first mesh.
Study 2
Step 1: Stationary
1
In the Model Builder window, under Study 2 click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Study Extensions section.
3
Select the Auxiliary sweep check box.
4
Click  Add.
5
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
DeltaT (Temperature difference)
1e-3 1e-2 1e-1
K
Step 2: Stationary 1
1
Right-click Study 2>Step 1: Stationary and choose Duplicate.
Second step solves for the largest Rayleigh numbers using the finest mesh.
2
In the Settings window for Stationary, locate the Study Extensions section.
3
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
DeltaT (Temperature difference)
1 10
K
4
Click to expand the Mesh Selection section. In the table, enter the following settings:
Geometry
Mesh
Geometry 2
Mesh 3
5
Click to expand the Values of Dependent Variables section. Find the Initial values of variables solved for subsection. From the Settings list, choose User controlled.
6
From the Method list, choose Solution.
7
From the Study list, choose Study 2, Stationary.
8
From the Selection list, choose Last.
Add a Combine Solutions study step that concatenates the two solutions and makes it possible to treat the output as a single parametric study.
Combine Solutions
In the Study toolbar, click  Combine Solutions.
Solution 2 (sol2)
In the Study toolbar, click  Show Default Solver.
Step 3: Combine Solutions
1
In the Model Builder window, click Step 3: Combine Solutions.
2
In the Settings window for Combine Solutions, locate the Combine Solutions Settings section.
3
From the First solution list, choose Study 2/Solution Store 1 (sol3).
4
From the Second solution list, choose Study 2/Solution 2 (sol2).
Solution 2 (sol2)
1
In the Model Builder window, expand the Solution 2 (sol2) node.
2
In the Model Builder window, expand the Study 2>Solver Configurations>Solution 2 (sol2)>Stationary Solver 1 node, then click Fully Coupled 1.
3
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
4
From the Nonlinear method list, choose Constant (Newton).
5
In the Model Builder window, expand the Study 2>Solver Configurations>Solution 2 (sol2)>Stationary Solver 2 node, then click Fully Coupled 1.
6
In the Settings window for Fully Coupled, locate the Method and Termination section.
7
From the Nonlinear method list, choose Constant (Newton).
8
In the Damping factor text field, type 0.8.
While the default solver solves the problem without any issue, GMG is faster for this model.
9
In the Model Builder window, expand the Study 2>Solver Configurations>Solution 2 (sol2)>Stationary Solver 2>AMG, nonisothermal flow (nitf2) (merged) node, then click Multigrid 1.
10
In the Settings window for Multigrid, locate the General section.
11
From the Solver list, choose Geometric multigrid.
12
In the Study toolbar, click  Compute.
Results
Velocity (spf2)
This default plot group shows the fluid velocity magnitude in only half of the cube. To plot the other half, proceed as follows.
Mirror 3D 1
1
In the Results toolbar, click  More Datasets and choose Mirror 3D.
2
In the Settings window for Mirror 3D, locate the Plane Data section.
3
From the Plane list, choose zx-planes.
A new dataset containing mirror values is now created. Return to the velocity plot to use this dataset.
Velocity (spf2)
1
In the Model Builder window, click Velocity (spf2).
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Dataset list, choose Mirror 3D 1.
4
In the Velocity (spf2) toolbar, click  Plot.
Temperature (ht2)
This default plot group shows the temperature distribution. The mirror dataset created previously can be reused here to plot the entire cube.
Velocity, Front Plane
1
In the Home toolbar, click  Add Plot Group and choose 3D Plot Group.
2
In the Settings window for 3D Plot Group, type Velocity, Front Plane in the Label text field.
3
Locate the Data section. From the Dataset list, choose Mirror 3D 1.
Slice 1
1
In the Velocity, Front Plane toolbar, click  Slice.
2
In the Settings window for Slice, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 2 (comp2)>Laminar Flow 2>Velocity and pressure>spf2.U - Velocity magnitude - m/s.
3
Locate the Plane Data section. From the Plane list, choose zx-planes.
4
In the Velocity, Front Plane toolbar, click  Plot.
This slice view shows the velocity magnitude in the same plane as in the 2D model (Figure 8).
Next, plot arrows of the tangential velocity field in the vertical plane parallel to the plates to reproduce Figure 9.
Temperature, 10 K Offset
1
In the Home toolbar, click  Add Plot Group and choose 3D Plot Group.
2
In the Settings window for 3D Plot Group, type Temperature, 10 K Offset in the Label text field.
3
Locate the Data section. From the Dataset list, choose Mirror 3D 1.
Slice 1
1
In the Temperature, 10 K Offset toolbar, click  Slice.
2
In the Settings window for Slice, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 2 (comp2)>Heat Transfer in Fluids 2>Temperature>T2 - Temperature - K.
3
Locate the Plane Data section. In the Planes text field, type 1.
4
Locate the Coloring and Style section. From the Color table list, choose ThermalLight.
5
In the Temperature, 10 K Offset toolbar, click  Plot.
Temperature, 10 K Offset
In the Model Builder window, click Temperature, 10 K Offset.
Arrow Volume 1
1
In the Temperature, 10 K Offset toolbar, click  Arrow Volume.
2
In the Settings window for Arrow Volume, locate the Expression section.
3
In the x component text field, type 0.
4
Locate the Arrow Positioning section. Find the x grid points subsection. In the Points text field, type 1.
5
Find the y grid points subsection. In the Points text field, type 25.
6
Find the z grid points subsection. In the Points text field, type 25.
7
Locate the Coloring and Style section. From the Color list, choose Black.
8
In the Temperature, 10 K Offset toolbar, click  Plot.
Temperature, 10 K Offset
The arrows follow convective cells at the four corners for a temperature difference of 10 K. Follow the steps below to reproduce Figure 10 and to see these cells when the temperature difference is reduced to 1 K and 0.1 K.
Temperature, 1 K Offset
1
Right-click Temperature, 10 K Offset and choose Duplicate.
2
In the Settings window for 3D Plot Group, type Temperature, 1 K Offset in the Label text field.
3
Locate the Data section. From the Parameter value (DeltaT (K)) list, choose 1.
4
In the Temperature, 1 K Offset toolbar, click  Plot.
Temperature, 0.1 K Offset
1
Right-click Temperature, 1 K Offset and choose Duplicate.
2
In the Settings window for 3D Plot Group, type Temperature, 0.1 K Offset in the Label text field.
3
Locate the Data section. From the Parameter value (DeltaT (K)) list, choose 0.1.
4
In the Temperature, 0.1 K Offset toolbar, click  Plot.