We study convection of an incompressible Newtonian fluid heated from below in a two-dimensional domain of height : the Bénard problem. The lower wall is maintained at a temperature and the upper wall is maintained at a temperature , where . The governing equations are the (2D) Navier–Stokes equations under the Boussinesq approximation, in which all variations in physical properties with temperature are neglected, apart from that of the density in the gravitational-body-force term in the momentum equations. This "buoyancy" term is given by
where is the variation in density and is the -th component of the gravitational body force. Under the additional assumption that variations in temperature are small, we can use the linear relationship
where is the coefficient of thermal expansion of the fluid, is the (dimensional) temperature and is the density at the reference temperature .
The equations governing the fluid motion are thus the Navier–Stokes equations with the inclusion of the additional buoyancy term. In Cartesian coordinates, we have
Here, is the -th (dimensional) velocity component and is the position in the -th coordinate direction; is the dynamic viscosity of the fluid at the reference temperature and is the dimensional time.
The equation that governs the evolution of the temperature field is the advection-diffusion equation where the "wind" is the fluid velocity. Thus,
where is the (constant) thermal diffusivity of the fluid.
We choose the height of the domain, , as the length scale and let the characteristic thermal diffusion speed over that length, , be the velocity scale, so that the Péclet number, . The fluid pressure is non-dimensionalised on the viscous scale, , and the hydrostatic pressure gradient is included explicitly, so that we work with the dimensionless excess pressure. The temperature is non-dimensionalised so that it is -0.5 at the upper (cooled) wall and 0.5 at the bottom (heated) wall and the reference temperature is then . Finally, the timescale is chosen to be the thermal diffusion timescale, . Hence
The governing equations become
The appropriate dimensionless numbers are the Prandtl number , and the Rayleigh number, ; is the acceleration due to gravity and is the kinematic viscosity of the fluid.
We consider the solution of this coupled set of equations in a two-dimensional domain , . The boundary conditions are no-slip at the top and bottom walls
constant temperature at the top and bottom walls (heated from below)
and symmetry boundary conditions at the sides:
We assume that gravity acts vertically downward so that and .
There is a trivial steady-state solution that consists of a linearly-varying temperature field balanced by a quadratic pressure field:
A linear stability analysis shows that this solution becomes unstable via an up-down, symmetry-breaking, pitchfork bifurcation at a critical Rayleigh number of with a critical wavenumber of see for example Hydrodynamic and Hydromagnetic Stability by S. Chandrasekhar OUP (1961). Thus, for there are three possible steady solutions, the (unstable) trivial steady state and two (stable) symmetry-broken states. In principle, all three states can be computed directly by solving the steady equations. However, we typically find that if the steady computation is started with a zero initial guess for the velocity and temperature, the Newton method converges to the trivial state. In order to demonstrate that this state is indeed unstable we therefore apply a time-dependent, mass-conserving perturbation to the vertical velocity at the upper wall and time-march the system while rapidly reducing the size of the perturbation. The system then evolves towards the nontrivial steady state as shown in the animation from which the plots shown above were extracted. (In the next tutorial where we discuss the adaptive solution of this problem we shall demonstrate an alternative technique for obtaining this solutions).
Note that by choosing our domain of a particular size and applying symmetry conditions at the sidewalls we are only able to realise a discrete set of wavelengths (those that exactly fit into the box). At the chosen Rayleigh number, 1800, only one of these modes is unstable; that of wavelength 2.
The problem contains three global parameters, the Péclet number, the Prandtl number and the Rayleigh number which we define in a namespace, as usual. In fact, is the natural dimensionless grouping, and so we use the inverse Prandtl number as our variable.
In the driver code we set the direction of gravity and construct our problem, using the new BuoyantQCrouzeixRaviartElement, a multi-physics element, created by combining the
QCrouzeixRaviart Navier-Stokes elements with the
QAdvectionDiffusion elements via multiple inheritance. (Details of the element's implementation are discussed in the section Creating the new BuoyantQCrouzeixRaviartElement class below.)
We assign the boundary conditions at the time and initially perform a single steady solve to obtain the trivial (and temporally unstable) trivial solution; see the section Comments for a more detailed discussion of the
The result of this calculation is the trivial symmetric base flow. We next timestep the system using the (unstable) steady solution as the initial condition. As time increases, the flow evolves to one of the stable asymmetric solutions, as shown in the animation of the results. As usual, we only perform a few timesteps when the code is used as a self-test, i.e. if any command-line parameters are passed to the driver code.
The problem class contains five non-trivial functions: the constructor, the
fix_pressure(...) function, as well as the functions
doc_solution(...), all discussed below.
We pass the element type as a template parameter to the problem constructor, which has no arguments. The constructor creates a
BFD<2> timestepper and builds a
RectangularQuadMesh of elements.
Next, the boundary constraints are imposed. We pin all velocities and the temperature on the top and bottom walls and pin only the horizontal velocity on the sidewalls. Since the domain is enclosed, the pressure is only determined up the an arbitrary constant. We resolve this ambiguity by pinning a single pressure value, using the
We complete the build of the elements by setting the pointers to the physical parameters and finally assign the equation numbers
In order to examine the stability of the symmetric state, we impose a time-dependent boundary condition that transiently perturbs the vertical velocity field at the upper boundary. Our boundary condition is
where . The perturbation is zero at , tends to zero as , and is mass conserving. This is implemented in the function below
This function is a simple wrapper to the element's
This function is used to ensure that the time-dependent boundary conditions are set to the correct value before solving the problem at the next time level.
This function writes the complete velocity, pressure and temperature fields to a file in the output directory.
The sketch below illustrates how the new multi-physics
BuoyantQCrouzeixRaviartElement is constructed by multiple inheritance from the two existing single-physics elements:
QCrouzeixRaviartElementis based on a nine-node quadrilateral geometric
QElementfamily. All of its
Nodesstore two values, the horizontal and vertical velocity, respectively. The element also stores internal
Datawhich represents the (discontinuous) pressure degrees of freedom; in the sketch this
Datais represented by the dashed box.
QAdvectionDiffusionElementis based on the same geometric
FiniteElementand stores one value (the temperature, ) at each
Both elements are fully-functional and provide their contributions to the global system of nonlinear algebraic equations that is solved by Newton's method via the two member functions
fill_in_contribution_to_residuals(...)computes the element's contribution to the global residual vector for a given "wind". The "wind" is specified by its virtual member function
get_wind_adv_diff(...)and in the single-physics advection diffusion problems studied so far, the "wind" tended to specified a priori by the user. The element's
fill_in_contribution_to_jacobian(...)computes the elemental Jacobian matrix, i.e. the derivatives of the elemental residual vector with respect to its unknown nodal values (the temperatures).
fill_in_contribution_to_residuals(...)computes the element's contribution to the global residual vector for a given body force. The body force is specified by its virtual member function
get_body_force_nst(...)and in the single-physics Navier-Stokes problems studied so far, the body force tended to specified a priori by the user. The element's member function
fill_in_contribution_to_jacobian(...)computes the elemental Jacobian matrix, i.e. the derivatives of the elemental residual vector with respect to its unknown nodal and internal values (the velocities and the pressure).
When combining the two single-physics elements to a multi-physics element, we have to take the interaction between the constituent equations into account: In the coupled problem the "wind" in the advection-diffusion equations is given by the Navier-Stokes velocities, while the body force in the Navier-Stokes equations is a function of the temperature. When implementing these interactions we wish to recycle as much of the elements' existing functionality as possible. This may be achieved by the following straightforward steps:
FiniteElement::required_nvalue(...)function to ensure that each
Nodeprovides a sufficient amount of storage for the (larger) number of nodal values required in the multi-physics problem.
fill_in_contribution_to_jacobian(...)functions and to use finite-differencing only for the off-diagonal (interaction) blocks; see the section Comments a more detailed discussion of this technique.]
That's all! Here is the implementation:
The class contains a single new physical parameter, the Rayleigh number, as usual referenced by a pointer to a double precision datum,
with suitable access functions.
The constructor calls the constructors of the component classes (
QAdvectionDiffusionElement) and initialises the value of the Rayleigh number to zero, via a static default parameter value.
We must overload the function
FiniteElement::required_nvalue() because the new element will store
DIM+1 unknowns at each node:
DIM fluid velocity components and the value of the temperature, as shown in the sketch above.
In the standard single-physics advection-diffusion elements the temperature is the only value stored at the nodes and is stored as
value(0). Similarly, in the single-physics Navier–Stokes elements, the fluid velocities are stored in the first
DIM nodal values. In our new multi-physics element, we must decide where to store the different variables and then inform the single-physics elements of our choice. As indicated in the above sketch, we choose to store the temperature after the fluid velocities, so that it is
value(DIM). The recommended mechanism for communicating the location of the variables to the single-physics elements is to use an index function. Hence, single-physics elements that are to be the components of multi-physics elements must have an index function for their variables. For instance, the function
u_index_adv_diff(...) is used in the
AdvectionDiffusionEquations class to read out the position (index) at which the advected variable (the temperature) is stored. That function is now overloaded in our multi-physics element:
We need not overload the index function for the fluid velocities because they remain stored in the first
DIM positions at the node.
The coupling between the two sets of single-physics equations is achieved by overloading the two functions
get_wind_adv_diff(), used in the advection-diffusion equations and
get_body_force_nst(), used in the Navier–Stokes equations
The elemental residual vector is composed of the residuals from the two single-physics elements and we simply call the underlying functions for each element in turn.
Finally, we compute the Jacobian matrix by finite-differencing the element's combined residual vector, using the default implementation of the
fill_in_contribution_to_jacobian(...) function in the
FiniteElement base class:
Finally, we overload the output function to print the fluid velocities, the fluid pressure and the temperature.
Problem::newton_solve()employs Newton's method to solve the system of nonlinear algebraic equations arising from the
Problem'sdiscretisation. The current
Datavalues are used as the initial guess for the Newton iteration. On return from this function, all unknown
Datavalues will have been assigned their correct values so that the solution of the problem may be plotted by calls to the elements'
outputfunctions. We tended to use this function to solve steady problems.
Problem::unsteady_newton_solve(dt,...)increments time by
dt, shifts the "history" values and then computes the solution at the advanced time, +
dt. On return from this function, all unknown
Datavalues (and the corresponding "history" values) will have been assigned their correct values so that the solution at time +
dtmay be plotted by calls to the elements'
outputfunctions. We tended to use this function for unsteady problems.
Problem::unsteady_newton_solve(...)function shows that this function is, in fact, a wrapper around
Problem::newton_solve(), and that the latter function solves the discretised equations including any terms that arise from an implicit time-discretisation. The only purpose of the wrapper function is to shift the history values before taking the next timestep. This raises the question how to compute steady solutions (i.e. solutions obtained by setting the time-derivatives in the governing equation to zero) of a
Problemthat was discretised in a form that allows for timestepping, as in the problem studied here. This is the role of the function
Problem::steady_newton_solve(): The function performs the following steps:
Problemby calling their
Problem::newton_solve()function to compute the solution of the discretised problem with all time-derivatives set to zero.
TimeSteppers(unless they were already in "steady" mode when the function was called).
Problem::assign_initial_values_impulsive()to ensure that the "history" values used by the (now re-activated)
TimeSteppersare consistent with an impulsive start from the steady solution just computed.
Datavalues (and the corresponding "history" values) will have been assigned their correct values so that the solution just computed is a steady solution to the full unsteady equations.
BuoyantQCrouzeixRaviartElementincludes such an implementation. The full finite-difference-based computation discussed above is used if the code is compiled with the compiler flag
USE_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENT. Finite-differences are used for the off-diagonal blocks only when the compiler flag
USE_OFF_DIAGONAL_FD_JACOBIAN_FOR_BUOYANT_Q_ELEMENTis passed. When comparing the two versions of the code, we found the run times for the full finite-difference-based version to be approximately 3-7% higher, depending on the spatial resolution used. The implementation of the more efficient version is still straightforward and can be found in the source code boussinesq_convection.cc.
BuoyantQCrouzeixRaviartElementincludes such an implementation and, moreover, it is the default behaviour. We found the assembly time for the analytic coupled Jacobian to be approximately 15% of the finite-difference based versions. The implementation is reasonably straightforward and can be found in the source code boussinesq_convection.cc.
QTaylorHoodElementsas the "fluid" element part of the multi-physics elements. N.B. in this case, the temperature must be stored as the first variable at the nodes because we assume that it is always stored at the same location in every node.
A pdf version of this document is available.