Buoyancy Flow of Free Fluids

This example follows the benchmark buoyancy flow posed by de Vahl Davis (Ref. 1) for free fluids. Buoyancy flow of free fluids is very important in Earth sciences with temperature and concentration affecting density in moving fluids, for example, in pipes, along shorelines, and within lakes. Here the buoyancy flow results from density that varies with a temperature change. The COMSOL Multiphysics results match those from the published study (Ref. 1). The model was provided by John Kamel of the University of Notre Dame.

This example repeatedly solves a problem of buoyant flow in a square cavity. It thereby analyzes different temperature distributions and convective flow patterns from variations in, for example, fluid properties, cavity size, and temperature drops. The iterative process is tuned for a fast, efficient solution using nondimensional parameters and a Boussinesq term for the buoyant drive with the Laminar Flow and the Heat Transfer in Fluids interfaces. As an alternative, the same problem can be solved using only the Non-Isothermal Flow interface available with the CFD Module or the Heat Transfer Module.

Figure 1: Domain geometry and boundary conditions for the heat balance in this example of buoyant flow in free fluids.

Model Definition

The previous figure illustrates the model geometry. The fluid fills a square cavity with impermeable walls so the fluid flows freely within it but does not exit from it. The right and left edges of the cavity are, respectively, the high and low temperature sources. The upper and lower boundaries are insulated. The temperature differential produces the density variation that drives the buoyant flow.

The laminar flow and heat transfer interfaces in this example are bidirectionally coupled. The Boussinesq term defines a force dependent on the temperature and the fluid velocity transports heat.

It is possible to formulate the incompressible Navier–Stokes equations with a Boussinesq buoyancy term on the right-hand side to account for the lifting force due to thermal expansion:

In these expressions, the dependent variables for flow are u, the vector of fluid velocity, and pressure, p. T represents temperature, T0 is a reference temperature, g denotes gravity acceleration, ρ0 gives the reference density, μ is the dynamic viscosity, and αp equals the coefficient of volumetric thermal expansion.

The heat balance comes from the conduction-convection equation

where k denotes the thermal conductivity, and Cp is the specific heat capacity of the fluid.

For the Navier–Stokes equations, impermeable, no-slip boundary conditions apply. The no-slip condition results in zero velocity at the wall, with pressure within the domain remaining undefined. Because the lack of information about p makes it difficult to achieve convergence, you arbitrarily fix the pressure at a point.

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.

Notes About the COMSOL Implementation

The model addresses a wide range of cavity sizes, fluid properties, and temperature drops using material properties set up with nondimensional Rayleigh and Prandtl numbers. The Rayleigh number, Ra = (Cpρ2gαpTL3) ⁄ (μk), gives the ratio of buoyant to viscous forces. Here L is the length of a side wall. The Prandtl number, Pr = (μCp ⁄ k), gives the ratio of kinematic viscosity to thermal diffusivity.

Setting the body force in the y-direction for the momentum equation to Fy = (Ra ⁄ Pr)(T − Tc), and the fluid properties to Cp = Pr, and ρ =μ = k = 1, produces a set of equations with nondimensional variables p, u, and T.

As Ra increases, viscous forces decrease in importance. You can examine a wide range of scenarios sequencing through different values of Ra using an with continuation selected for the Ra parameter in the study step. At high Ra values, starting with a good initial condition and a well-tuned mesh becomes increasingly important. Because the continuation feature in the parametric solver uses extrapolation from the previous solution as the initial condition for the next one, together with automatic step length control, it fulfills the first of these two requirements. Then getting a well-tuned mesh is straightforward: simply set up the mesh for the most difficult problem to solve — the one with the highest value of Ra. To that end, the element size near the prescribed temperature boundaries corresponds to the thickness of the boundary layer when Ra = 106.

Results and Discussion

As noted earlier, this COMSOL Multiphysics example reproduces a benchmark buoyant flow problem reported in Ref. 1. The composite images in Figure 2 summarize temperatures (surface), velocity fields (arrows), and x-velocities (contours) for four Ra values. With the given definitions for the dimensionless variable, it follows that as Ra increases so does temperature. The results in the figure indicate how the vigor and complexity of the convection increase at higher values of Ra, as expected. These results are virtually identical to the benchmark solution except the COMSOL Multiphysics plots provide higher resolution than the images in the original publication.