This paper is part of a Special Issue to honour the contributions of Professor Graeme Hocking to the ANZIAM Society in particular and applied and industrial mathematics more generally. For decades, Graeme has been a loyal supporter of the ANZIAM annual conference, and has organized several of them. He has recently served as the chair of ANZIAM, and with Andrew Bassom has been a joint chief editor of the ANZIAM Journal. He has also had a very long-standing relationship with the Mathematics in Industry Study Group, and has made significant contributions to a number of practical industrial technologies, such as the use of “air knives” to remove molten metal during galvanizing. In the past decade or two, Graeme has become increasingly involved with industrial mathematics worldwide, and has established extensive collaborative links with colleagues in Ireland, Europe, the UK and South Africa.
Both authors of this article have been privileged to know Graeme and work with him. Larry Forbes has been particularly fortunate to have known Graeme as both a colleague and a close friend for the past 35 years, visiting each other regularly and having written about 40 papers together. Graeme’s professional life has been characterized by the two properties of great generosity and unimpeachable integrity. He has often gone to extraordinary lengths to help students, particularly in cases where their circumstances have become difficult. His research is uncompromisingly honest and accurate; in addition, his commitment to honesty and to student welfare has seen him act in the role of whistle-blower, to his own considerable cost. It is therefore an honour to have been invited to contribute this article to this Special Issue.
1 Introduction
For many of the viscous fluids encountered in everyday experience, the Navier–Stokes equations are generally regarded as providing an accurate description of their flow behaviour. However, relatively few exact solutions to these equations are known, for situations of practical interest. The simplest of all these closed-form solutions describes Couette flow (see [Reference Drazin5, pp. 13–15]), in which the viscous fluid is confined between flat parallel plates. One plate is stationary and is located on the plane
$y = -H$
, and the second plate is at
$y = H$
and moves with some speed
$U_p$
in the x-direction. Then the Navier–Stokes equations admit an incompressible solution in which the flow occurs only in the single direction of the x-axis, with speed
$u(y) = U_p (y+H) / (2H)$
in the channel
$-H < y < H$
. Since the fluid is viscous, it must adhere to the two plates and therefore have the same velocity as each plate, and this is the reason why
$u(y)$
adopts this profile that simply varies linearly with y. Nevertheless, in spite of its simplicity, this Couette flow is of industrial interest, and occurs in lubrication theory (see Batchelor [Reference Batchelor2]).
Because so few exact solutions to the Navier–Stokes equations are known, researchers in the first half of the twentieth century became particularly adroit at making full use of them to solve closely related problems of interest; this point is argued eloquently in Van Dyke’s text on perturbation methods [Reference Van Dyke19, p. 9]. Stewartson [Reference Stewartson17] extended Prandtl boundary-layer theory to consider viscous flow near a horizontal plate that is heated. He considered steady-state boundary-layer flow, but evidently made an error in the sign of a heat-flux term, leading to mistaken conclusions about the well-posedness of the problem, as was pointed out later by Gill et al. [Reference Gill, Zeh and del Casal8]. These solutions were then generalized by Chao and Cheema [Reference Chao and Cheema4] to allow for steady flow of the fluid but an unsteady temperature response at a semi-infinite lower plate. A similar type of analysis had been carried out by Riley [Reference Riley16], who likewise assumed that the fluid flow would not be affected significantly by changes in temperature at the plate, and so could be regarded as being in a steady state, but that the temperature imposed on a semi-infinite lower plate would jump impulsively to some new constant value along the plate. He carried out a similarity-solution analysis in the boundary layer close to the plate, and solved the resulting equations both for very small and very large values of the similarity variable. A brief review of the literature on convective cooling from a hot plate in a fluid-filled channel is presented by Nair et al. [Reference Nair, Karuppasamy, Benziger and Balakrishnan15], and these authors indicate the importance of such heat-transfer and flow problems in a variety of industrial applications, and in some geophysical flow situations. There are also several works on convective heat transfer from plates oriented at some angle other than purely horizontal, and an example of these is the analysis by Merkin [Reference Merkin14] of convection due to a heated vertical plate. He carries out an asymptotic analysis, and shows that there are a number of different regions in the boundary layer, characterized by different flow and heat-transfer behaviour. Umemura and Law [Reference Umemura and Law18] allowed heat transfer to occur from a semi-infinite plate at arbitrary angle to the horizontal and used similarity solutions to characterize the flow near the leading edge of the plate, where the behaviour was similar to that for a pure horizontal plate, and further downstream, where the flow behaviour was more characteristic of a vertical plate. Mellado [Reference Mellado13] solved equations for the fluid velocity and a buoyancy variable using a pseudo-spectral method, and introduced random perturbations to the initial buoyancy so as to create highly mixed convective regions above the plate.
A viscous fluid flow over a semi-infinite plate, in which the entire plate is impulsively switched to a higher temperature at some initial time, possibly represents a situation of limited practical interest, although it has provided interesting insights into the structure of thermal boundary layers. The mathematical advantage of studying such problems in unbounded geometries is that they often offer the prospect of finding similarity solutions, which reduce the complexity of the calculations very considerably. Nevertheless, Lewandowski [Reference Lewandowski12] has pointed out that these models involving semi-infinite plates do not accord with experimental observations, and this discrepancy becomes worse as the plate moves from being oriented vertically to fully horizontal. The major reason is that such models assume that a boundary layer exists along the entire plate, whereas experiments often show that a plume of hot fluid forms instead, moving away from the plate and into the fluid. Lewandowski proposes an alternative simplified model with horizontal plates of finite dimensions with different regimes of flow within the fluid, and also presents photographs of plumes rising abruptly from heated plates.
An important application of viscous flows over horizontal plates, on which the heat exchange occurs over a patch of finite width, is to the creation of models of the phenomenon of urban heat islands. This is the situation in which large cities act as heat sources on a patch of the landscape, potentially generating plumes that move through the lower atmosphere, and a general overview was given by Kanda [Reference Kanda11]. A computational fluid mechanics (CFD) simulation has been used by Gagliano et al. [Reference Gagliano, Nocera and Aneli7] to study temperature rise within cities, and the effects of air recirculating in “urban canyons”. Huo et al. [Reference Huo, Chen, Geng, Tao, Liu, Zhang and Leng10] have recently given a detailed overview of the CFD techniques currently in use, including the various ways of modelling fluid turbulence near the urban heat island, and another recent review of experimental and computational approaches to mesoscale urban meteorology has been given by Wu et al. [Reference Wu, Wang, Chen, Cao, Cao, Yu and Song22], with particular focus on the Chinese experience of ventilation within urban valleys.
The present paper considers an idealized model of flow over an isolated heat island, notionally located over some interval
$-q < x <q$
on a horizontal bottom plate of bilaterally infinite extent, located on the plane
$y = -H$
. We note that while the heat island is of finite extent in the x-direction, we are considering it to extend indefinitely in the z-direction, thus allowing us to consider a two-dimensional model. The fluid motion is caused by a top plate, at
$y = H$
, moving at some constant speed
$U_p$
and so this would simply be a Couette flow in the absence of the heat island. A brief description of this model is presented in Section 2. For simplicity, the Boussinesq approximation [Reference Chandrashekhar3, p. 16] for the flow of the viscous fluid is adopted, and is described in Section 3, and a semianalytical spectral method for solving the resulting equations is discussed in detail. Section 4 then presents a linearized solution for this Boussinesq viscous problem, and discusses in detail a closed-form solution in the special case where the top plate is stationary (
$U_p = 0$
), since this provides a valuable independent check on the reliability of the numerical spectral solution. Results of the computation are presented in Section 5; to begin, a careful validation of the spectral solution method is made, by comparing temperature profiles at early times with those obtained from the linearized solution in Section 4 for small input heat flux on the heat island. Highly nonlinear solutions are then presented, for large heating events and Couette fluid flow, and they show plume formation with overturning sections at their heads. The powerful computing resources available to us allow highly accurate long-term calculations to be made, and these indicate that periodic plume shedding from the heat island can evidently occur. Some final remarks in Section 6 conclude the paper.
2 Theoretical framework
In this work, we examine flow between two plates in relative motion. There is a linear velocity profile between the plates in the undisturbed steady solution. We take the bottom plate at
$y= -H$
to be stationary, and the flow is driven by the fixed speed
$U_p$
of the top plate at
$y = H$
. We assume that the fluid is incompressible, and the task is to solve equation (2.1) which expresses the conservation of mass, combined with the incompressible Navier–Stokes equation (2.2). This results in the system
$$ \begin{align} \boldsymbol \nabla \cdot \mathbf{q} = 0 \kern100pt\end{align} $$
$$ \begin{align} \frac{\partial \mathbf{q}}{\partial t} + ( \mathbf{q} \cdot \boldsymbol \nabla )\mathbf{q} - \frac{\mu}{\rho}\nabla^2 \mathbf{q} = \mathbf{g} - \frac{\boldsymbol \nabla p}{\rho} \end{align} $$
of nonlinear partial differential equations. Here,
$\mathbf {q}$
denotes the velocity vector field,
$\mu $
is the dynamic viscosity,
$\rho $
is the density, p is the pressure and
$\mathbf {g}$
is the acceleration due to gravity.
We assume the fluid and both the upper and lower plates have initial temperature
$T_0$
throughout. Additional heat energy is then introduced in some notional region
$-q < x < q$
at the heat island on the lower plate. This is done in a continuous manner, so as to avoid a discontinuity in temperature at
$t=0$
. The bottom plate is therefore assigned temperature
$T_B$
, and perhaps the simplest such function would take the form
$$ \begin{align} T_B(x,t) = \begin{cases} T_0 + \epsilon [ 1+\cos(\pi x/q) ]\, f(t), & |x|<q , \\ T_0 , & |x|\geq q , \end{cases} \end{align} $$
for some
$q: 0<q<L$
, in which L is a convenient length scale along the x-axis (it will only be used in the numerical scheme of Section 3. The function
$f(t)$
defines the way in which heating is added to the system, over the heat island on the lower wall, and we set
$f(0)=0$
to ensure continuity. In the course of this study, we have investigated a number of different heat input functions
$T_B$
. These need not all have compact support over
$-q < x <q$
as (2.3) does, provided that
$T_B$
drops rapidly to zero outside this interval, in which case the results are not greatly affected by this choice of function.
In this model, the perturbation to the background density is assumed to vary oppositely to the temperature change, that is,
$$ \begin{align} \overline{\rho}=- ( T - T_0 ). \end{align} $$
The reference temperature
$T_0$
does not appear in any of the equations governing the fluid flow, so without loss of generality, we set it to zero.
Figure 1 gives a sketch of the system at the initial time
$t=0$
. The flow between the plates, in the channel
$-H < y < H$
has a linear profile in height and is thus simply a viscous Couette flow. The bottom plate is stationary, but the upper plate may move horizontally with speed
$u = U_p$
.
A definition sketch of the base Couette velocity profile, which varies linearly with height in the channel. The lower wall at
$y=-H$
is stationary and the upper wall at
$y=H$
moves to the right with speed
$U_p$
.
3 Boussinesq viscous model
The incompressibility condition equation (2.1) is satisfied immediately through the use of a streamfunction
$\psi $
to define the velocities. In vector form, this relationship can be written
$\mathbf {q} = \boldsymbol {\nabla } \times ( \psi \mathbf {k} )$
, in which
$\mathbf {k}$
represents the unit vector pointing out of the x–y flow plane. After nondimensionalizing using an arbitrary length scale and the acceleration of gravity to scale the time variable, the fluid vorticity
$\boldsymbol {\zeta } = \boldsymbol {\nabla } \times \mathbf {q}$
is introduced as usual, and found to have only the single component
$\boldsymbol {\zeta } = \Omega \mathbf {k}$
. In the Boussinesq approximation, the density is assumed to be constant everywhere, except in the buoyancy term where the perturbation
$\overline {\rho }$
to the base density is retained. The vector curl of the momentum equation (2.2) is taken, and yields a scalar vorticity equation
$$ \begin{align} \frac{\partial \Omega}{\partial t} = \nu\nabla^2 \Omega - \bigg( u\frac{\partial \Omega}{\partial x} + v\frac{\partial \Omega}{\partial y} \bigg) -\frac{\partial \overline{\rho}}{\partial x}. \end{align} $$
Here, the symbol
$\nu $
represents the dimensionless kinematic viscosity, and can be regarded as an inverse Reynolds number. In this investigation, the variation
$\overline {\rho }$
of density comes about through temperature variations, according to the effective equation of state (2.4). We choose here to fix the temperatures on the top and bottom plates, but note that other boundary conditions may also be considered. The temperature and velocity are free to vary at
$x=\pm L$
, apart from at
$y=\pm H$
.
The dimensionless temperature throughout the fluid is therefore represented spectrally as
$$ \begin{align} T(x,y) = T_B(x,t) \frac{H-y}{2H} + \sum_{m=1}^M \sum_{n=1}^N B_{m,n}(t) \cos[\alpha_{m}(x+L)] \sin[\alpha_{n} (y+H)]. \end{align} $$
With an appropriate choice of
$\alpha _{n}$
, this form for T ensures that temperatures are fixed at the top and bottom plates so that
$T=0$
at
$y=H$
, and
$T=T_B$
at
$y=-H$
, but are free to take any values elsewhere in the channel
$-H < y < H$
and in some computational window
$-L < x < L$
in the horizontal direction. In Boussinesq theory, the density variation
$\overline {\rho }$
is assumed to obey a transport equation, or equivalently, since density and temperature T are related through the equation of state (2.4), the temperature can be regarded as satisfying the temperature transport equation
$$ \begin{align} \frac{\partial T}{\partial t} + ( \mathbf{q} \cdot \boldsymbol{\nabla} ) T = \sigma \nabla^2 T \end{align} $$
with temperature diffusion constant
$\sigma $
. In terms of temperature variation, the vorticity equation (3.1) now becomes
$$ \begin{align} \frac{\partial \Omega}{\partial t} = \nu \nabla^2 \Omega - \bigg( u\frac{\partial \Omega}{\partial x} + v\frac{\partial \Omega}{\partial y} \bigg) + \frac{\partial T}{\partial x}. \end{align} $$
3.1 Velocity components
The constant-shear base flow,
$$ \begin{align} \mathbf{q} = U_p \frac{y+H}{2H} \mathbf{\hat{i}}, \end{align} $$
depicted in Figure 1 must be enforced on the upper and lower boundaries of the computational region,
$-L < x < L$
,
$-H < y < H$
, through the use of appropriate basis functions. We therefore choose a spectral representation of the streamfunction
$\psi $
that meets the boundary conditions on the top and bottom plate, but otherwise allows any velocity within the computational region. This requires basis functions which, along with their first derivatives in y, are zero on the top and bottom boundaries. The basis functions must not all be zero at any other point, and must not all be symmetric or all antisymmetric. We therefore choose
$$ \begin{align} \psi(x,y,t) = \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) \cos[\alpha_{m} (x+L)] F_n(y) + U_p\frac{(y+H)^2}{4H}, \end{align} $$
in which we have defined the functions
$$ \begin{align} F_n(y) & = \sin[\alpha_{n} (y+H)] - \frac{1-\cos(n \pi)}{2} \frac{n}{N+1} \sin \bigg[ \frac{N+1}{2H}\pi (y+H) \bigg] \nonumber \\ &\quad - \frac{1+\cos(n \pi)}{2} \frac{n}{N+2} \sin \bigg[ \frac{N+2}{2H}\pi (y+H) \bigg]. \end{align} $$
Here,
$U_p$
is the speed of the top plate in the positive x-direction. The sets of constants
$$ \begin{align} \alpha_{m}&=\frac{(m-1) \pi}{2 L}, \nonumber \\ \alpha_{n}&=\frac{n \pi}{2H} \end{align} $$
have been introduced for convenience of notation.
The functions
$F_n(y)$
,
$n=1,2,\ldots ,N$
, in equation (3.7) have been designed to allow the streamfunction
$\psi $
in (3.6) to satisfy the required boundary conditions on the walls
$y=\pm H$
. Since both velocity components u and v must vanish on these two walls, it follows that the functions
$F_n(y)$
in (3.6) and their first derivatives
$F_n^{\prime }(y)$
are both zero there. This is enforced here using a rearrangement of the inner series in (3.6) suggested by Forbes and Brideson [Reference Forbes and Brideson6]. Although the numerical method will ultimately seek to find only the N sets of coefficients
$A_{m,1}$
,
$A_{m,2}$
,
$\ldots , A_{m,N}$
for each value of
$m = 1, 2, \ldots , M$
, we nevertheless introduce the next two sets
$A_{m,N+1}$
and
$A_{m,N+2}$
in a standard Fourier series expression in the y-variable, so that the inner sum in (3.6) would take the form
$$ \begin{align*} \sum_{n=1}^N A_{m,n}(t) \sin[\alpha_{n} (y+H)] + A_{m,N+1}(t) \sin[\alpha_{N+1} (y+H)] + A_{m,N+2}(t) \sin[\alpha_{N+2} (y+H)]. \nonumber \end{align*} $$
The requirement that this function and its first derivative must both vanish at
$y=\pm H$
permits the additional coefficients
$A_{m,N+1}$
and
$A_{m,N+2}$
to be eliminated in favour of the first N coefficients, giving rise to a series involving the functions
$F_n(y)$
presented in (3.7).
This form (3.6) of the streamfunction allows the horizontal and vertical components u and v of the velocity vector to vary arbitrarily in the channel
$-H<y<H$
, except at the top and bottom where they are both zero. They are obtained from (3.6) by direct differentiation, giving
$$ \begin{align*} u&=\frac{\partial{\psi}}{\partial{y}} \nonumber\\ &=\sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) \cos[\alpha_{m}(x+L)] F_n^{\prime}(y) + U_p\frac{y+H}{2H}, \nonumber\\ v&=-\frac{\partial{\psi}}{\partial{x}} \nonumber\\ &=\sum_{m=1}^M \sum_{n=1}^N \alpha_{m}A_{m,n}(t) \sin[\alpha_{m}(x+L)] F_n(y), \nonumber \end{align*} $$
where the prime on
$F_n(y)$
denotes differentiation with respect to y.
3.2 Vorticity component
The single vorticity component
$\Omega $
is
$$ \begin{align} \Omega(x,y,t)&=\frac{\partial{v}}{\partial{x}}-\frac{\partial{u}}{\partial{y}} = -\nabla^2\psi \nonumber\\ &=-\frac{U_p}{2H}+\sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) \cos[\alpha_{m}(x+L)]( \alpha_{m}^2 F_n(y) - F_n^{\prime\prime}(y) ). \end{align} $$
The x and y derivatives of the vorticity and its Laplacian, required for equation (3.4), are obtained by direct differentiation of (3.9).
3.3 Initial conditions
Let
$\psi (x,y,0)=\psi _1$
be the initial value of the streamfunction. For convenience, we will write
$\psi _0 \equiv \psi _0(y) = U_p (y+H)^2 / (4H)$
. Thus,
$$ \begin{align*} \psi_1 = \psi_0 + \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(0) \cos[\alpha_{m}(x+L)] F_n(y). \end{align*} $$
To compute the initial values
$A_{m,n}(0)$
of the coefficients, we multiply by
$\cos [\alpha _{k}(x+L)]$
and
$\sin [\alpha _{p} (y+H)]$
, in which
$\alpha _k$
and
$\alpha _p$
are obtained from the definitions (3.8), and integrate over the spatial domain. This results in
$$ \begin{align*} &\int_{-L}^L \int_{-H}^H \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(0) \cos[\alpha_{m}(x+L)] F_n(y) \cos[\alpha_{k}(x+L)] \sin[\alpha_{p} (y+H)] \, dy \, dx \nonumber\\ &\quad=\int_{-L}^L \int_{-H}^H (\psi_1-\psi_0) \cos[\alpha_{k}(x+L)] \sin[\alpha_{p} (y+H)] \, dy \, dx. \nonumber \end{align*} $$
It follows from the usual orthogonality relations for trigonometric functions that
$$ \begin{align*} A_{k,p}(0) = \frac{1}{(1+\delta_1^k)LH} \int_{-L}^L \int_{-H}^H (\psi_1-\psi_0) \cos[\alpha_{k}(x+L)] \sin[\alpha_{p} (y+H)] \, dy \, dx. \end{align*} $$
We now have the Fourier coefficients
$A_{k,p}(0)$
for the initial conditions. The Kronecker delta symbol
$\delta _1^k$
takes the value one if
$k=1$
, but zero otherwise.
3.4 Finding the Fourier coefficients
Incorporating the Fourier series representation (3.9) for
$\Omega $
, the vorticity equation (3.4) becomes
$$ \begin{align*} &\sum_{m=1}^M \sum_{n=1}^N \dot A_{m,n}(t) \cos[\alpha_{m}(x+L)] ( \alpha_{m}^2 F_n(y)-F_n^{\prime\prime}(y) ) \nonumber\\ &\quad= \nu \sum_{m=1}^M \sum_{n=1}^N A_{m,n}(t) \cos[\alpha_{m}(x+L)] ( -\alpha_m^4 F_n(y) + 2 \alpha_m^2 F_n^{\prime\prime} - F_n^{\prime\prime\prime\prime}(y) ) \nonumber\\ & \qquad - \bigg( u\frac{\partial \Omega}{\partial x} + v\frac{\partial \Omega}{\partial y} \bigg) + \frac{\partial T}{\partial x}, \end{align*} $$
where
$\dot A_{m,n}$
denotes the time derivative of
$A_{m,n}(t)$
. To calculate the
$\dot A_{m,n}$
, we multiply through by basis functions
$\cos [\alpha _{k}(x+L)] \sin [\alpha _{p} (y+H)]$
and integrate over the spatial domain. Orthogonality of the trigonometric functions reduces this to
$$ \begin{align*} \dot A_{k,p}(t) & = - \nu A_{k,p}(t) \Delta_{k,p}^2 - \frac{1}{ \Delta_{k,p}^2 (1+\delta_k^1) L H} \\ &\quad \times \int_{-L}^L \int_{-H}^H \bigg( u\frac{\partial \Omega}{\partial x} + v\frac{\partial \Omega}{\partial y} - \frac{\partial T}{\partial x} \bigg) \cos[\alpha_{k}(x+L)] \sin[\alpha_{p} (y+H)] \, dx \, dy. \nonumber \end{align*} $$
In this expression, it is convenient to introduce the extra constants
$$ \begin{align} \Delta_{k,p}^2 = \alpha_k^2 + \alpha_p^2 \end{align} $$
and utilize parameters defined previously in (3.8).
The Fourier coefficients
$B_{m,n}(t)$
for the temperature are obtained in the same way, by analysing the temperature transport equation (3.3). The orthogonality conditions similarly generate the system of differential equations
$$ \begin{align*} \dot B_{k,p}(t) & = - \sigma B_{k,p}(t) \Delta_{k,p}^2 - \frac{1}{(1+\delta_k^1) L H} \\ &\quad \times \int_{-L}^L \int_{-H}^H \bigg( u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} - \sigma \nabla^2 T \bigg) \cos[\alpha_{k}(x+L)] \sin[\alpha_{p} (y+H)] \, dx \, dy, \nonumber \end{align*} $$
in which
$\dot B_{k,p} (t)$
denotes the time derivative of these coefficients, and the constants
$\Delta _{k,p}^2$
are as defined in (3.10). This system of differential equations is now integrated forward in time, from the initial coefficients calculated in Section 3.3.
4 Linearized model
Here, a linearized approximation to the full Boussinesq theory of Section 3 is discussed. Since variations from the underlying Couette flow (3.5) are caused by the heating on the bottom plate, this will be assumed to be a small effect, governed by a small parameter
$\epsilon $
. On the bottom wall, condition (2.3) will be generalized in this section to take the form
$$ \begin{align} T_B(x,t)= \epsilon g(x) f(t), \end{align} $$
in which the shape function
$g(x)$
is left arbitrary for now. The velocity components, vorticity, temperature and streamfunction are next expanded in the series
$$ \begin{align} u(x,y,t) & = U_p \frac{y+H}{2H} + \epsilon u_1(x,y,t) + \mathcal{O} ( \epsilon^2 ),\nonumber \\ v(x,y,t) & = 0 + \epsilon v_1(x,y,t) + \mathcal{O} ( \epsilon^2 ),\nonumber \\ \Omega(x,y,t) & = - \frac{U_p}{2H} + \epsilon \Omega_1(x,y,t) + \mathcal{O} ( \epsilon^2 ), \\ T(x,y,t) & = 0 + \epsilon T_1(x,y,t) + \mathcal{O} ( \epsilon^2 ),\nonumber \\ \psi (x,y,t) & = U_p \frac{(y+H)^2}{4H} + \epsilon \Psi_1(x,y,t) + \mathcal{O} ( \epsilon^2 ),\nonumber \end{align} $$
and substituted into the full governing equations. Only terms of first order in
$\epsilon $
are retained. In terms of the perturbed velocity components, the linearized vorticity becomes
$$ \begin{align} \Omega_1 = \frac{\partial v_1}{\partial x} - \frac{\partial u_1}{\partial y} = - \nabla^2 \Psi_1, \end{align} $$
and (3.4) shows that it satisfies the linearized momentum equation
$$ \begin{align} \frac{\partial \Omega_1}{\partial t} + U_p \frac{y+H}{2H} \frac{\partial \Omega_1}{\partial x} = \frac{\partial T_1}{\partial x} + \nu\nabla^2 \Omega_1. \end{align} $$
The linearized temperature transport equation (3.3) becomes
$$ \begin{align} \frac{\partial T_1}{\partial t} + U_p \frac{y+H}{2H} \frac{\partial T_1}{\partial x} = \sigma\nabla^2 T_1. \end{align} $$
The boundary conditions are that the two perturbation velocity components
$u_1$
and
$v_1$
must both be zero at the top and bottom plates
$y = \pm H$
. The perturbation temperature
$T_1$
must also be zero at the top plate
$y = H$
, but on the bottom plate
$y = - H$
it takes the value
$T_1(x, -H, t) = g(x)f(t)$
so as to satisfy (4.1).
Although this system of partial differential equations (4.4)–(4.5) is now linear, it is nevertheless not straightforward to solve. When the top plate is in motion,
$U_p \ne 0$
, Fourier transforms in x followed by Laplace transforms in t yield an Airy equation for the transformed linearized temperature function
$T_1$
(which can equivalently be represented in terms of Bessel functions of order
$1/3$
). However, the result is of little practical value as the transforms are not readily invertible. Since the aim of this section is to provide a linearized solution that can be compared with the fully nonlinear numerical results, we have therefore opted to consider the simpler case
$U_p = 0$
in this section.
The linearized temperature equation (4.5) with a stationary upper plate becomes a classical heat equation
$$ \begin{align} \frac{\partial T_1}{\partial t} = \sigma \bigg( \frac{\partial^2 T_1}{\partial x^2} + \frac{\partial^2 T_1}{\partial y^2} \bigg), \end{align} $$
and can be solved separately for
$T_1$
. We allow the plates to extend indefinitely in the positive and negative x-directions and define the Fourier transform
$$ \begin{align} \widehat{T}_1 (\omega; y,t) & \equiv \mathfrak{F} \{ T_1(x,y,t) \} \nonumber \\ & = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} T_1(x,y,t) {e}^{- i\omega x} \, dx. \end{align} $$
The Fourier transform of the shape function
$g(x)$
for the heat input (4.1) on the bottom plate is
$$ \begin{align} \widehat{G} (\omega ) = \mathfrak{F} \{ g(x) \}. \end{align} $$
The Laplace transform of the quantity
$\widehat {T}_1$
in (4.7) is
$$ \begin{align*} \widetilde{T}_1(\omega; y; s) & \equiv \mathfrak{L} \{ \widehat{T}_1 (\omega; y,t) \} \nonumber \\ & = \int_0^{\infty} \widehat{T}_1 (\omega; y,t) {e}^{- s t} \, dt, \end{align*} $$
and similarly it is convenient to define a Laplace transform
$$ \begin{align*} \widetilde{F}(s) = \mathfrak{L} \{ f(t) \} \end{align*} $$
of the function
$f(t)$
in (4.1) that describes how the heat input on the bottom plate is activated. It is assumed here that
$f(0)=0$
so that, initially, the temperature in the channel is everywhere constant.
The Fourier–Laplace transform of the heat equation (4.6) now yields the ordinary differential equation
$$ \begin{align*} ( s + \sigma \omega^2 ) \widetilde{T}_1 = \sigma \frac{{d}^2 \widetilde{T}_1}{{d} y^2}, \end{align*} $$
subject to the two boundary conditions
$$ \begin{align*} \widetilde{T}_1 & = 0 \quad \textrm{on } y=H,\\ \widetilde{T}_1 & = \widehat{G} (\omega ) \widetilde{F}(s) \quad \textrm{on } y=H. \end{align*} $$
This boundary-value problem can be solved at once to give
$$ \begin{align} \widetilde{T}_1(\omega; y; s) = \widehat{G} (\omega ) \widetilde{F}(s) \widetilde{\mathcal{H}} (\omega; y; s ), \end{align} $$
in which we have defined Green’s function
$$ \begin{align} \widetilde{\mathcal{H}} (\omega; y; s ) = \frac{\sinh [ (H-y) \sqrt{\omega^2 + s/\sigma}\, ]}{\sinh [ 2H \sqrt{\omega^2 + s/\sigma}\, ]}. \end{align} $$
The presence of the square-root radical in (4.10) suggests that branch cuts in the complex s-plane will be needed, to confine the function
$\widetilde {\mathcal {H}}$
to some domain in which it remains analytic. However, this turns out not to be correct, since expanding the
$\sinh $
functions in both the numerator and denominator in Taylor series shows that the square-root radical cancels from the numerator and denominator. As a result, Green’s function in (4.10) is meromorphic, possessing only simple poles at values of s for which the denominator vanishes. There is a countably infinite set of such poles distributed along the negative real s-axis, at points
$$ \begin{align} s_k = - \sigma \bigg[ \omega^2 + \bigg( \frac{k\pi}{2H} \bigg)^2 \bigg], \quad k = 0 , 1 , 2 , \ldots. \end{align} $$
Furthermore, the simple pole at
$s_0$
for which
$k=0$
in (4.11) can be discounted, since it gives only a removable singularity in the expression (4.10).
By the convolution theorem for Laplace transforms, the inverse Laplace transform of the function in (4.9) gives
$$ \begin{align} \widehat{T}_1 (\omega; y,t) = \widehat{G} (\omega ) \int_0^t f(t-\tau ) \widehat{h} (\omega; y,\tau ) \, d\tau , \end{align} $$
thus returning to the Fourier-transformed
$\omega $
-plane of the solution in the expression (4.7). Here,
$\widehat {h}$
is the inverse Laplace transform of Green’s function
$\widetilde {\mathcal {H}}$
defined in (4.10), and for positive t it can be shown to be the sum of the residues at the simple poles in (4.11). After some algebra, we obtain
$$ \begin{align} \widehat{h} (\omega; y,t) = - \frac{\sigma}{H} \sum_{k=1}^{\infty} \gamma_k \cos(k\pi ) \sin[ \gamma_k (H-y) ] \exp [ - \sigma t ( \omega^2 + \gamma_k^2 ) ], \end{align} $$
in which we have defined constants
$$ \begin{align} \gamma_k = \frac{k\pi }{2H} , \quad k = 1 , 2 , 3 , \ldots. \end{align} $$
The function in (4.13) is again an entity that exists in the Fourier-transformed
$\omega $
-plane, and a further inversion of transforms is required to obtain the final solution for the temperature perturbation
$T_1$
in physical variables.
Finally, the two expressions (4.12) and (4.13) are combined to yield
$$ \begin{align} \widehat{T}_1 (\omega; y,t) & = - \frac{\sigma}{H} \sum_{k=1}^{\infty} \gamma_k \cos(k\pi ) \sin[ \gamma_k (H-y) ] \nonumber \\ & \quad \times \int_0^t f(t-\tau ) \exp[ - \sigma \tau \gamma_k^2 ] \widehat{G} (\omega ) {e}^{-\sigma \tau \omega^2} \, d\tau. \end{align} $$
It follows from the convolution theorem for Fourier transforms that the two remaining terms in (4.15) that involve the Fourier wavenumber variable
$\omega $
may be written
$$ \begin{align*} \widehat{G} (\omega ) {e}^{-\sigma \tau \omega^2} = \mathfrak{F} \bigg\{ \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} g(\xi ) \frac{1}{\sqrt{2\sigma \tau}} \exp \bigg[ - \frac{(x- \xi )^2}{4\sigma \tau} \bigg] \, d\xi \bigg\}. \end{align*} $$
Consequently, the inverse Fourier transform applied to expression (4.15) gives the temperature perturbation function
$$ \begin{align} T_1 (x,y,t) & = - \frac{\sigma}{H} \sum_{k=1}^{\infty} \gamma_k \cos(k\pi ) \sin[ \gamma_k (H-y) ] \nonumber \\ & \quad \times \int_0^t \frac{f(t-\tau )}{2\sqrt{\pi \sigma \tau}} \exp[ - \sigma \tau \gamma_k^2 ] \int_{-\infty}^{\infty} g(\xi ) \exp \bigg[ - \frac{(x- \xi )^2}{4\sigma \tau} \bigg] \, d\xi \, d\tau. \end{align} $$
This is the general linearized solution for a stationary upper wall, and with arbitrary heat shape function
$g(x)$
and switching function
$f(t)$
and an initial condition
$T_1 = 0$
throughout the channel. It is strictly valid in the domain
$-H < y \le H$
, and does not formally hold on the lower wall
$y = - H$
due to the presence of removable singularities there.
4.1 A test example
A particularly simple choice for the heat shape function
$g(x)$
on the lower wall is the Gaussian profile
$$ \begin{align} g(x) = \exp \bigg( - \frac{x^2}{\beta^2} \bigg). \end{align} $$
This has the advantage that the inner integration in (4.16) can be carried out in closed form. The result is
$$ \begin{align} T_1 (x,y,t) & = - \frac{\sigma \beta}{H} \sum_{k=1}^{\infty} \gamma_k \cos(k\pi ) \sin[ \gamma_k (H-y) ] \nonumber \\ & \quad \times \int_0^t \frac{f(t-\tau )}{\sqrt{4\sigma \tau + \beta^2}} \exp[ - \sigma \tau \gamma_k^2 ] \exp \bigg[ - \frac{x^2}{4\sigma \tau + \beta^2} \bigg] \, d\tau. \end{align} $$
Likewise, a simple choice is also made here for the function
$f(t)$
in the bottom boundary condition (4.1). An example is the function
$$ \begin{align} f(t) = 1 - e^{- \lambda t} \end{align} $$
which satisfies the condition
$f(0) = 0$
and moves monotonically towards its limit
$1$
as
$t \rightarrow \infty $
. This function (4.19) can be substituted into (4.18) and it turns out that the integration can be carried out in closed form, and involves complementary error functions. However, the result is extremely difficult to evaluate numerically, because of large errors due to loss of significance caused by very large terms multiplying very small ones. Instead, the integrals in (4.18) are evaluated far more accurately by direct numerical quadrature. We make the change of variable
$\tau = t w$
and incorporate the special choice (4.19) to give
$$ \begin{align} T_1 (x,y,t) & = - \frac{\beta \sigma t}{H} \sum_{k=1}^{\infty} \gamma_k \cos(k\pi ) \sin[ \gamma_k (H-y) ] \nonumber \\ & \quad \times \int_0^1 [ 1 - {e}^{- \lambda t (1-w)} ] \frac{1}{D(w;t)} \exp[ - \sigma t \gamma_k^2 w ] \exp \bigg[ - \frac{x^2}{D^2(w;t)} \bigg] \, d w , \end{align} $$
in which
$$ \begin{align*} D(w;t) = \sqrt{4\sigma t w + \beta^2}. \end{align*} $$
The integrals in equation (4.20) are evaluated readily, and to very high precision, using up to 20 000 points and Gauss–Legendre quadrature as made available in the MATLAB code lgwt written by von Winckel [Reference von Winckel20].
4.2 Linearized vorticity and streamfunction
Calculating the linearized vorticity
$\Omega _1$
and streamfunction
$\Psi _1$
using Fourier and Laplace transforms is no longer a straightforward affair, even in the case of a stationary top plate, and in all likelihood advanced numerical techniques would be needed to invert these transforms. Nevertheless, there are some interesting analytical features of the flow that are revealed by this analytical approach, and these are now investigated briefly here.
When
$U_p = 0$
, the linearized vorticity equation (4.4) becomes simply
$$ \begin{align*} \frac{\partial \Omega_1}{\partial t} = \frac{\partial T_1}{\partial x} + \nu \bigg( \frac{\partial^2 \Omega_1}{\partial x^2} + \frac{\partial^2 \Omega_1}{\partial y^2} \bigg). \end{align*} $$
This shows that, at least in the linearized approximation (4.2), the flow is driven by the temperature profile (4.16). Applying the Fourier transform followed by the Laplace transform results in the ordinary differential equation
$$ \begin{align} \nu \frac{{d}^2 \widetilde{\Omega}_1}{{d} y^2} - ( s + \omega^2 \nu ) \widetilde{\Omega}_1 = - i\omega \widetilde{T}_1 \end{align} $$
for the doubly transformed linearized vorticity function
$\widetilde {\Omega }_1( \omega; y; s )$
. The transformed temperature perturbation function
$\widetilde {T}_1$
is as determined previously in (4.9)–(4.10). Similarly applying transforms to (4.3) yields a further differential equation
$$ \begin{align} \nu \frac{{d}^2 \widetilde{\Psi}_1}{{d} y^2} - \omega^2 \widetilde{\Psi}_1 = - \widetilde{\Omega}_1 \end{align} $$
for the transformed streamfunction.
It is straightforward to see that the differential equation (4.21) has a general solution which is convenient to express in the form
$$ \begin{align} \widetilde{\Omega}_1 (\omega; y; s) & = - \frac{\nu}{s} ( \tilde{\alpha}^2 + \omega^2 ) [ \tilde{K}(\omega ,s) \cosh ( \tilde{\alpha} (y-H) ) + \tilde{L}(\omega ,s) \sinh ( \tilde{\alpha} (y-H) ) ] \nonumber \\ & \quad - \widetilde{\Gamma}(\omega ,s) \sinh [ \tilde{\beta} (H-y) ] , \end{align} $$
in which
$$ \begin{align*} \tilde{\alpha} = \sqrt{ \omega^2 + s/ \nu}, \quad \tilde{\beta} = \sqrt{ \omega^2 + s/ \sigma} \end{align*} $$
are constants in the space of the transform variables, and the coefficient of the particular integral term is
$$ \begin{align} \widetilde{\Gamma}(\omega ,s) = \frac{i\sigma \omega \widehat{G}(\omega ) \widetilde{F}(s)}{s (\nu - \sigma )\sinh (2H\tilde{\beta} )}. \end{align} $$
The quantities
$\tilde {K}(\omega ,s)$
and
$\tilde {L}(\omega ,s)$
are constants of integration in the transform plane. The differential equation (4.22) for the transformed streamfunction can now also be solved without difficulty, and gives
$$ \begin{align} \widetilde{\Psi}_1 (\omega; y; s) & = \tilde{M}(\omega ,s) \cosh ( \omega (y-H) ) + \tilde{N}(\omega ,s) \sinh ( \omega (y-H) ) \nonumber \\ & \quad - \frac{\nu}{s} [ \tilde{K}(\omega ,s) \cosh ( \tilde{\alpha} (y-H) ) + \tilde{L}(\omega ,s) \sinh ( \tilde{\alpha} (y-H) ) ] \nonumber \\ & \quad + \frac{\sigma}{s} \widetilde{\Gamma}(\omega ,s) \sinh [ \tilde{\beta} (H-y) ]. \end{align} $$
This expression involves a further two arbitrary integration constants
$\tilde {M}(\omega ,s)$
and
$\tilde {N}(\omega ,s)$
.
The four arbitrary constants
$\tilde {K}$
,
$\tilde {L}$
,
$\tilde {M}$
and
$\tilde {N}$
in these transformed solutions (4.23) and (4.25) are finally determined by making the two Fourier–Laplace transformed velocity components u and v both zero on the two walls at
$y = \pm H$
. This is necessary to satisfy the no-slip boundary conditions there. The process is straightforward, but the expressions for
$\tilde {K}$
, and so on, are very complicated. This appears to prevent any possibility of inverting the Laplace transforms in analytical form, and so we do not pursue this closed-form solution any further. Nevertheless, these solutions have revealed that there is a resonance in the linearized solution when
$\nu = \sigma $
. This is at once evident from the fact that the denominator in the constant
$\widetilde {\Gamma }(\omega ,s)$
in equation (4.24) possesses a zero divisor at
$\nu = \sigma $
. This is also true of the constants
$\tilde {K}$
, and so on, since these all involve
$\widetilde {\Gamma }$
too.
4.3 Steady-state linearized temperature profile
When the series (4.16) for the temperature perturbation function
$T_1 (x,y,t)$
is evaluated on the computer, it is evident that it approaches a steady form as
$t\rightarrow \infty $
. However, it is not straightforward to take this formal limit in the solution (4.16). For completeness, we briefly derive the steady-state temperature profile here, working directly from the time-independent form of the heat equation (4.6).
The steady form of (4.6) reduces to Laplace’s equation for the perturbed temperature. This must be solved subject to boundary conditions that
$T_1$
is zero at the top plate
$y = H$
, and on the bottom plate
$y = - H$
it must take the value
$T_1(x, -H , t) = g(x)$
so as to satisfy (4.1) in the steady limit. The Fourier transforms (4.7)–(4.8) of this steady problem are taken and the resulting boundary-value problem solved in the Fourier space. The convolution theorem is used to invert the transform, yielding the solution
$$ \begin{align} T_1 (x,y) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} g(x-\xi ) h(\xi , y) \, d\xi , \end{align} $$
in which Green’s function in the Fourier space is
$$ \begin{align*} \widehat{\mathcal{H}} (\omega; y ) = \frac{\sinh ( \omega (H-y) )}{\sinh ( 2\omega H )} \end{align*} $$
and has inverse Fourier transform
$$ \begin{align} h(x,y) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \frac{\sinh ( \omega (H-y) )}{\sinh ( 2\omega H )} {e}^{i\omega x} \, d\omega. \end{align} $$
The integral in this expression can be evaluated using complex residue theory, since the integrand has simple pole singularities on the imaginary
$\omega $
-axis, at the infinite set of points
$\omega _k = \pm k\pi i / (2H)$
,
$k=0, \pm 1, \pm 2, \pm 3, \ldots .$
(The singularity at
$\omega _0 = 0$
is removable, and so can be ignored.) The integral in (4.27) along the real
$\omega $
-axis is closed with a semicircle either in the upper half plane if
$x>0$
or in the lower half-plane if
$x<0$
, and after some calculation, the function
$h(x,y)$
in (4.27) becomes
$$ \begin{align*} h(x,y) = - \frac{1}{H} \sqrt{\frac{\pi}{2}} \sum_{k=1}^{\infty} \cos(k\pi ) \sin[ \gamma_k (H-y) ] \exp [ - \gamma_k |x| ], \end{align*} $$
where the constants
$\gamma _k$
are as defined in (4.14). It follows from (4.26) that
$$ \begin{align} T_1 (x,y) & = - \frac{1}{2H} \sum_{k=1}^{\infty} \cos(k\pi ) \sin[ \gamma_k (H-y) ] \nonumber \\ & \quad \times \int_{-\infty}^{\infty} g(x-\xi ) \exp [ - \gamma_k |\xi | ] \, d\xi. \end{align} $$
With the example heat energy function
$g(x)$
given in (4.17), the integral term in the solution (4.28) can be evaluated in closed form, and gives exponentials multiplied by complementary error functions. However, these are extremely difficult to evaluate, since as the index k increases, they represent a product of extremely large numbers with very small ones, and as a result, loss-of-precision errors render the expressions useless. It is therefore better simply to evaluate the integrals in (4.28) by direct numerical quadrature, since this can be done easily and with very high precision.
5 Numerical results
We begin this presentation of numerical results by comparing the predictions of the numerical spectral method in Section 3 with the linearized solution in Section 4. In order to compare exactly, the bottom heating function
$T_B(x,t)=\epsilon g(x) f(t)$
makes direct use of the same functions
$g(x)$
and
$f(t)$
given in (4.17) and (4.19), respectively. The diffusion coefficient
$\sigma $
is set to
$0.0001$
in the numerical algorithm, unless stated otherwise. This value is somewhat arbitrary, but having a small value keeps the “interface” from becoming overly blurred. The dimensionless dynamic viscosity coefficient
$\nu $
is taken to be
$0.001$
. Again, this choice is arbitrary. We note that different values of this coefficient will have a significant effect on the flow.
5.1 Stationary upper plate
The first case considered here is for a stationary top plate,
$U_p=0$
. This allows direct comparison with the linearized solution (4.20) from Section 4. In this linearized theory, it turns out that the temperature equation decouples from the others, so that the temperature function
$T_1$
could be found without reference to the fluid velocity and streamfunction. For this reason, the linearized temperature is not affected by buoyancy, since this only appears in the momentum equations. In Figure 2, we have therefore compared the linearized solution (4.20) with the results of the nonlinear solution scheme in which the gravitational term has been removed. Here, the channel-height parameter is
$H=2$
, so that the fluid lies within the interval
$-2 < y < 2$
, and the horizontal computational window
$-L<x<L$
has taken
$L=5$
. At these early dimensionless times, the plume has not yet developed greatly, and so only the bottom twentieth
$-2 < y < -1.8$
of the channel is shown, for ease of viewing. The top panel of results shows the nonlinear solution at the four times
$t = 7.5$
,
$15$
,
$22.5$
and
$30$
and the lower panel shows the linearized solution at these same four times. The agreement between these two sets of results is excellent, which confirms the reliability of the nonlinear spectral solution scheme. This is as expected, particularly since in the absence of any gravitational force, there is no fluid flow and the temperature transport equation (3.3) reduces simply to the heat equation (4.6).
Comparison of the linearized and nonlinear solutions at early times, in the absence of gravity. The nonlinear solution is shown in the top panels and the linear solution is shown below. The half-channel height is
$H=2$
and results are shown for the four early (dimensionless) times
$t = 7.5$
,
$15$
,
$22.5$
and
$30$
. For clarity, only the bottom twentieth of the computational region is shown. At all times shown there is good agreement between the two methods.
The effects of gravity are reintroduced in Figure 3, by reinstating the buoyancy term
$\partial T / \partial x$
in equation (3.4). There is again close agreement between the linearized and nonlinear solutions at the first two times
$t = 7.5$
and
$t = 15$
in Figure 3, just as in the buoyancy-free case shown in Figure 2. However, the nonlinear temperature profile at the later two times
$t = 22.5$
and
$t = 30$
in the top panel is beginning to be affected by the fluid flow (which is driven by the temperature profile). Such effects are absent in the linearized results in the lower panel of profiles in Figure 3. By time
$t=30$
, early plume formation is clearly visible, as the temperature develops a sharply peaked profile due to the effects of fluid convection. This coupling of the temperature with the fluid flow is absent from the linearized theory, but will drive formation of more complex structures at later times in the full nonlinear model.
Comparison of the linearized and nonlinear solutions at early times. Gravity effects are now included. The nonlinear solution is shown in the top panels and the linear solution is shown beneath. The channel-height parameter is
$H = 2$
, but for clarity, only the bottom twentieth of the computational region is shown. At very early times there is good agreement between the two methods. At later times, however, nonlinear effects may be seen in the upper panels with the beginnings of a plume visible by
$t=30$
.
The solutions shown in Figure 3 over the early (dimensionless) time interval
$0 < t < 30$
have been extended to considerably later times in Figure 4, where we now present temperature profiles at the four times
$t = 30$
,
$60$
,
$90$
and
$120$
. The top panel of results again corresponds to nonlinear theory, and in the bottom panel are the results at the same four times obtained from the linearized theory of Section 4. For the first column on the left of Figure 4, these results for
$t=30$
are identical to those on the rightmost column in Figure 3; the only difference now is in the scale on the vertical axes, since here we show the complete channel height
$-H < y < H$
for the nonlinear results in the upper panel, with the same channel half-height
$H=2$
. The pointed temperature profile that was so evident at time
$t=30$
in Figure 3 is still visible at
$t=30$
in the leftmost upper panel of Figure 4, and it is clear that the plume which began to be visible near
$t=30$
now develops rapidly into a fast moving stream, accelerating up to the top wall by
$t=120$
as shown in the top panel of Figure 4. The formation of an overturning mushroom-shaped structure is evident in the frames at
$t=60$
and
$t=90$
. The narrow connecting stem of hotter (lighter-density) fluid between the heat island at the bottom and the broad overturning head at the top is typical of a lazy plume, and closely resembles the nonlinear plume structures computed by Allwright et al. [Reference Allwright, Forbes and Walters1]. A “lazy plume” is one in which the flow is driven only by buoyancy terms (for this categorization of plumes, see Hunt and Kaye [Reference Hunt and Kaye9]). By the last time
$t = 120$
illustrated in Figure 4 the rising plume has struck the top wall
$y = 2$
of the channel and has begun to spread laterally across this wall. By contrast, the linearized solution at these same four later times, depicted in the lower panel of results, merely predicts a slowly rising temperature profile, since the only physical phenomenon accounted for in that approximation is the slow diffusion of heat energy through a stationary fluid. By
$t=120$
the plume in the nonlinear numerical solution has reached the upper wall, while the bubble of hot fluid in the linearized solution has not traversed one tenth of that distance. Thus, Figures 3 and 4 have confirmed the agreement between linearized and nonlinear theory at early times, validating the numerical solution scheme, but also show the importance of nonlinearity at later times.
Comparison of the linearized and nonlinear solutions at later times, for the same case illustrated in Figure 3. The nonlinear solution is shown in the top panels, and exhibits a rapidly rising plume, starting at about
$t=30$
, rolling over by
$t=60$
and striking the top plate by
$t=120$
. In contrast, the linear solution only continues slow diffusion. Note that for clarity the y-axes differ between the upper (nonlinear) and lower (linearized) panels. The top panels show the entire channel,
$-2 < y < 2$
, but the lower panels show only the lower eighth (
$-2 < y < -1.5$
) of the channel height.
In order to see pure plume evolution, free from the effects both of background fluid movement and interference from the upper wall, we have computed a solution with half-height parameter
$H = 8$
. The horizontal computational domain is set as
$-L < x <~L$
with
$L = 4$
, and temperature profiles are presented in Figure 5 at the three later times
$t=120$
,
$208$
and
$240$
. This solution required considerable computational resources in order to maintain accuracy, and was obtained here using
$(M,N) = (384,768)$
coefficients in the spectral representation (3.2), with five times that number of points in x and y. Even with these formidable computing resources, however, a small amount of numerical error related to Gibbs’s phenomenon in the Fourier series is still present in the solution, and this may be visible as small-amplitude oscillations in a couple of the temperature contours in the solution at (dimensionless) time
$t=208$
.
Nonlinear plume development in a higher channel, with half-height
$H = 8$
. The computational region horizontally is
$-L < x < L$
with
$L = 4$
. Temperature profiles are shown at the three times
$t = 120$
,
$208$
and
$240$
. The scale on the horizontal and vertical axes is the same, so that the plumes are shown as they would actually appear.
At the first time
$t = 120$
shown in Figure 5, a vertical plume has formed, and is similar to those shown in Figure 4. There is a bulbous overturning head, connected to the heat island on the bottom plate by a long thin vertical neck. The plume accelerates upward against gravity, forming a large nearly circular overturning region at its head. The neck region becomes narrower at the second time
$t=208$
shown, so that the head almost detaches from the heated region on the bottom plate. This is again reminiscent of the lazy plumes calculated by Allwright et al. [Reference Allwright, Forbes and Walters1]. Eventually, the top of the plume strikes the upper plate at
$y = H = 8$
and then begins to spread across this upper wall, as is evident in the temperature profile at the last time
$t = 240$
. Once the plume has reached the upper boundary, the choice of boundary conditions for both temperature and flow will play a significant role in determining the flow. In the present study, we have chosen an upper boundary that absorbs all temperature variation, thereby cooling the fluid. Presumably, this will encourage the fluid to flow downwards to a greater extent that would happen with an insulating upper boundary. The no-slip condition will have a drag-like effect on the motion, reducing its horizontal speed near the upper boundary.
Alternatively, if the upper boundary were moved to a much higher location, the plume would not spread as in Figure 5, but continue to rise and overturn, perhaps eventually pinching off and detaching.
The nonlinear results shown in Figures 4 and 5 highlight the way in which nonlinear effects are responsible for the development of plumes. By contrast, the linearized theory of Section 4 causes the heat behaviour to decouple from the effects of fluid motion, so that nonlinear convection cannot occur. As a result, linearized theory predicts that a heated region of fluid will develop near the heat island on the bottom plate. Eventually, it is to be expected that, in linearized theory, heat energy put into the system at the heat island would diffuse isotropically into the surrounding fluid, so that a steady-state temperature profile would then be formed. This is difficult to prove directly from the time-dependent linearized solution (4.16). However, linearized theory does permit a steady-state temperature profile to be calculated directly, as discussed in Section 4.3, and the result is given in equation (4.28).
It is interesting to calculate and compare the steady solution (4.28) with the unsteady linearized solution (4.16) evaluated at a very large time, and this is done in Figure 6. Here, contours of the temperature have been plotted for the nonsteady linearized solution at a very late time (
$t=12\,000\,000$
), and overlaid on these are contours computed from the steady solution (4.28). These linearized formulae were evaluated in a window
$-4<x<4, -2<y<2$
using
$200$
points in x and
$300$
in y. The nonsteady solution (4.16) is not trivial to calculate at very large times, and here it required
$32\,000$
modes and
$19\,200$
integration points for convergence. The steady solution (4.28) required
$2000$
modes and
$20\,000$
integration points over an interval
$-10<\xi <10$
to evaluate the inner quadratures. This does indeed demonstrate that the linearized solution achieves a steady-state temperature profile as
$t \rightarrow \infty $
, unlike its nonlinear counterpart that also experiences the effects of convection. It also shows the level of precision required to evaluate even this linearized solution accurately.
Comparison of the linearized steady solution at
$t=\infty $
with the nonsteady solution evaluated at time
$t=12\,000\,000$
. The two sets of temperature profiles are indistinguishable.
5.2 Moving top plate
To study the effects of background fluid motion, the top plate is now given dimensionless velocity
$U_p = 0.1$
to the right. The horizontal computational domain
$-L < x < L$
is expanded to the value
$L = 10$
, in order to see the evolution of the flow. For definiteness, the channel half-height parameter will be set to the value
$H = 2$
.
Figure 7 illustrates eight different temperature profiles, computed at the eight different times
$t = 40, 60, \ldots , 180$
indicated on the figure. These solutions again required substantial computer resources, and we used
$(M,N)=(720,144)$
Fourier modes in the spectral solution, with five times that number of points in x and y. Incipient plume formation is evident in the first image at time
$t = 40$
, and an asymmetric plume is clearly visible by time
$t=60$
. As time evolves, the plume accelerates upward, similarly to the results in Figure 4, except that now the plume is also distorted to the right as a result of the moving upper plate and the background Couette flow that it establishes. The plume continues to rise and has hit the upper wall by about time
$t=140$
. By time
$t=160$
, the plume detaches altogether from the heat island on the bottom plate, and a second plume has begun to form at the downwind edge of the hot zone. It also continues to grow and move to the right, and another plume structure also begins to form at time
$t=180$
.
Evolution of the plume for top plate moving with dimensionless speed
$U_p = 0.1$
. Nonlinear temperature profiles are shown for eight different times as indicated. The initial plume forms over the centre
$x=0$
of the heated region on the bottom plate. The plume detaches, and a new plume forms at the downwind end of the hot region.
The same flow as in Figure 7 was modelled out to time
$t=960$
. Eight sample solutions are shown in Figure 8, starting at time
$t = 820$
and with a new solution shown every
$20$
time units until
$t=960$
. These diagrams illustrate the process of plumes forming and detaching, without any suggestion of a steady state ever being reached. Instead, there appears to be a repeating cycle of plumes forming above the heat island, moving away to the right under the influence of the background Couette flow and eventually detaching as a new plume forms. From the numerical results, it appears that the process repeats approximately every
$40$
time units.
Temperature profiles for the same case as shown in Figure 7, but at eight significantly later times indicated on the diagrams. The figure shows the ongoing development of plumes from the heated region, and their subsequent detachment due to the movement of the top plate. This process appears to repeat approximately every
$40$
units of dimensionless time.
For interest, we show in Figure 9 a solution obtained with the diffusion coefficient
$\sigma =0.001$
, ten times as large as used for previous solutions, and with all other parameters unchanged. Four different nonlinear results are shown, at the four times
$t = 60$
,
$100$
,
$300$
and
$500$
. At this higher diffusion, plumes continue to develop, but the initial plume shown at time
$t=40$
maintains its shape for considerably longer, before also detaching and being swept downstream by the underlying Couette flow. Further numerical results (not shown) indicate that the cyclic pattern of plume development and detachment continues.
Later times showing the development of plumes from the heat island, with the higher diffusion coefficient
$\sigma = 0.001$
.
6 Conclusions
This paper has considered the behaviour of viscous fluid forced to move over a heat island on the bottom wall of a channel in two-dimensional (planar) flow. A simple model has been adopted, in which the fluid density decreases linearly in proportion to the amount by which the temperature increases. As a result, the heated region on the fluid bottom causes a patch of low-density fluid to form, and this then rises buoyantly against gravity, effectively forming an unstable interface, similar to that in Rayleigh–Taylor flow, for example. Nonlinear convective effects cause plumes to develop above the heat island, and these accelerate upwards. When the top plate is in motion and a background Couette flow is established, these plumes are also swept downstream. Ultimately, they detach from the heated region on the bottom plate, and a new plume forms and the process repeats.
A linearized theory was also studied here, and while it agrees very well with the nonlinear results at early times, it ultimately ceases to be of relevance, because it cannot account for convection, which is crucial to the formation of plumes. The linearized theory is also available in the case when a Couette flow is established by the movement of the upper plate, but has not been discussed here because the algebra required appears to be overwhelmingly complicated, involving Airy functions and their roots. Our principal interest in the linearized theory was as an important check on the reliability of the algorithm for the nonlinear solution, but it may prove possible to extend the linearized theory to account for fluid motion.
The model used here involves Boussinesq theory, in which the fluid density is assumed to vary only very slightly from some constant background value. Here, the temperature served as a proxy for density, but in any event, variations in the density or temperature only appear in the buoyancy term in the fluid momentum equation, consistently with the Boussinesq approximation. In a recent paper, however, Walters et al. [Reference Walters, Turner and Forbes21] developed what they referred to as a “completed Boussinesq theory”, which retains the exact form of the Navier–Stokes equation and does not make this approximation concerning the fluid density. They found that Boussinesq theory could lead to exaggerated predictions about the degree to which the head of the plume overturns and rolls up, whereas their completed Boussinesq theory did not share that defect, instead agreeing closely with more accurate models. It is possible that Boussinesq theory may have exaggerated the extent of plume overturning here too, although that remains for future investigation to determine.
Open access
