%ParentFile :A00compgfdE.tex
%\documentclass{book}
%\usepackage{amsmath}
%\renewcommand{\thechapter}{A\arabic{chapter}}
%\renewcommand \theequation {\arabic{equation}}
%\newcommand{\Dt}[1]{ \frac{D #1}{D t} }
%\newcommand{\PD}[2]{ \frac{\partial #1}{\partial #2} }
%\setcounter{chapter}{2}
%\begin{document}
%-----------------------------------------------------------------------
%\chapter{Thermal Convection}

The thermal convection is driven by heterogeneous temperature field.
When a fluid layer is heated from below and/or cooled from above, its thermal expansion yields a gravitationally unstable situation.
Namely, the dense cool fluid tends to sink and the hot light fluid rises to generate a thermal convection.
In this chapter, we will formulate the thermal convection of 2-dimensional non-compressible Newtonian fluid.
The nature of the convection depends on two non-dimensional parameters, i.e., the Rayleigh number and the Prandtl number.
We observe the generation and the growth of these convection currents in the numerical experiments.
%\begin{itemize}
%\item[]\textit{Keywords: convection, advection, Rayleigh number, Prandtl number}
%\end{itemize}

%-----------------------------------------------------------------------
\section{Theoretical Background}
%----------
\subsection{Master equations}
We suppose that a non-compressible Newtonian fluid of thickness $d$
 lies between horizontal rigid walls.
The assumption of non-compressibility is called anelastic
 approximation.
Independent variables $x,z,t$ are horizontal and vertical (upward)
 coordinates and time, respectively, while dependent variables
 $\mathbf{u}=(u,v), p, \rho$ and $T$ represent velocity, pressure,
 density and temperature.
The temperature of the upper wall $(z=d)$ and the lower wall $(z=0)$
 are fixed at $T_0$ and at $T_0+\Delta T (\Delta T>0)$, respectively.

The thermal convection is described by the following three equations,
\begin{align}
  &\textrm{Navier-Stokes equation}
    &\rho\Dt{\mathbf{u}} = -\nabla p + \eta\nabla^2\mathbf{u} + \rho\mathbf{g}, 
    \label{eq:NS} \\
  &\textrm{equation of continuity}
    &\rho\nabla\cdot\mathbf{u} = 0,
    \label{eq:cont} \\
  &\textrm{equation of thermal advection and conduction}
    &\Dt{T} = \kappa\nabla^2T,
    \label{eq:therm}
\end{align}
 where constants $\mathbf{g}=(0,-g), \eta$ and $\kappa$ are the gravitational
 acceleration, the coefficient of viscosity and the thermal diffusivity.
Equation (\ref{eq:NS}) is the equation of motion for Newtonian fluid,
 while (\ref{eq:cont}) corresponds to the conservation of mass for 
 non-compressible fluid.
The viscous heating and the internal heat generation are neglected.

For simplicity, we assume that density the $\rho$ depends only on
 the temperature,
\begin{equation}
  \rho = \rho_0 \{1-\alpha(T-T_0)\},
\end{equation}
 where constants $\rho_0$ and $\alpha$ are the density at temperature $T_0$
 and the thermal expansion coefficient, respectively.
Assuming that the change of density is small, we can neglect them except in
 the body force term of (\ref{eq:NS}), $\rho\mathbf{g}$.
That is so-called Boussinesq approximation, under which the above
 master equations are expanded to be
\begin{align}
  \PD{u}{t}+u\PD{u}{x}+w\PD{u}{z}
    &= -\frac{1}{\rho_0}\PD{p}{x} + \nu\nabla^2u,
    \label{eq:NSx} \\
  \PD{w}{t}+u\PD{w}{x}+w\PD{w}{z}
    &= -\frac{1}{\rho_0}\PD{p}{z} + \nu\nabla^2w - g\{1-\alpha(T-T_0)\},
    \label{eq:NSz} \\
  \PD{u}{x}+\PD{w}{z} &= 0, 
    \label{eq:cont_ex} \\
  \PD{T}{t}+u\PD{T}{x}+w\PD{T}{z} &= \kappa\nabla^2T,
\end{align}
 where $\nu=\eta/\rho_0$ is the coefficient of kinematic viscosity.

%----------
\subsection{Introduction of vorticity and stream function}
The elimination of the pressure terms by $\PD{}{x}$ (\ref{eq:NSz})
 $-\PD{}{z}$ (\ref{eq:NSx}) gives
\begin{equation*}
  \PD{\zeta}{t}+u\PD{\zeta}{x}+w\PD{\zeta}{z}
    = \nu\nabla^2\zeta + \alpha g \PD{T}{x},
\end{equation*}
where $\zeta$ is the $y$-component of the vorticity,
\begin{equation}
  \zeta = \PD{w}{x} - \PD{u}{z}.
    \label{eq:defzeta}
\end{equation}

To satisfy (\ref{eq:cont_ex}), the stream function $\psi$ is introduced,
\begin{equation}
  u=-\PD{\psi}{z}, w=\PD{\psi}{x}.
    \label{eq:defpsi}
\end{equation}
The above two definitions show that the vorticity and the stream 
 function are related by the Poisson's equation,
\begin{equation}
  \zeta = \nabla^2 \psi. \label{eq:Poisson}
\end{equation}
Using $\zeta$ and $\psi$, we can rewrite the master equations as
\begin{align}
  \PD{\zeta}{t} - \PD{\psi}{z}\PD{\zeta}{x} + \PD{\psi}{x}\PD{\zeta}{z}
    &= \nu\nabla^2\zeta + \alpha g \PD{T}{x}
    \label{eq:zeta} \\
  \PD{T}{t} - \PD{\psi}{z}\PD{T}{x} + \PD{\psi}{x}\PD{T}{z}
    &= \kappa\nabla^2T.
    \label{eq:temp}
\end{align}

Consequently, we obtained the equations
 (\ref{eq:Poisson}), (\ref{eq:zeta}) and (\ref{eq:temp})
 for unknown variables $\psi, \zeta$ and $T$.
Note that (\ref{eq:zeta}) represents the generation and the advection 
 of $\zeta$ caused by the \textit{horizontal} temperature gradient
 ($\PD{T}{x}$).
The generated vorticity yields the velocity field according to
 (\ref{eq:Poisson}).
Finally, that velocity field gives the new temperature field 
 through (\ref{eq:temp}).
This circulation describes the generation and the growth of the 
 thermal convection.

%----------
\subsection{Non-dimensionalization}
The master equations can be non-dimensionalized by normalizing length
 by $d$, time by the time scale of thermal diffusion, $d^2/\kappa$, and
 temperature by $\Delta T$,
\begin{equation} \begin{array}{lllll}
  \displaystyle (x,z) = (d \cdot x^*, d \cdot z^*), &
  \displaystyle t=\frac{d^2}{\kappa}t^*, &
  \displaystyle \zeta=\frac{\kappa}{d^2}\zeta^*, &
  \displaystyle \psi=\kappa\psi^*, &
  \displaystyle T=\Delta T^*,
\end{array} \end{equation}
 where the superscripts $^*$ represent non-dimensionalized variables.
%Hereafter, since all variables are non-dimensionalized, the superscripts
% $^*$ are omitted.
%\begin{align}
%  \textrm{Pr}^{-1}\left(
%    \PD{\zeta}{t} - \PD{\psi}{z}\PD{\zeta}{x} + \PD{\psi}{x}\PD{\zeta}{z}
%    \right) &= \nabla^2\zeta + \textrm{Ra}\PD{T}{x},
%    \label{eq:zeta'} \\
%  \PD{T}{t} - \PD{\psi}{z}\PD{T}{x} + \PD{\psi}{x}\PD{T}{z}
%    &= \nabla^2T,
%    \label{eq:temp'}
%\end{align}
Equations (\ref{eq:zeta}) and (\ref{eq:temp}) are rewritten as
\begin{align}
  \textrm{Pr}^{-1}\left(
    \PD{\zeta^*}{t^*} - \PD{\psi^*}{z^*}\PD{\zeta^*}{x^*} + \PD{\psi^*}{x^*}\PD{\zeta^*}{z^*}
    \right) &= \nabla^2\zeta^* + \textrm{Ra}\PD{T^*}{x^*},
    \label{eq:zeta'} \\
  \PD{T^*}{t^*} - \PD{\psi^*}{z^*}\PD{T^*}{x^*} + \PD{\psi^*}{x^*}\PD{T^*}{z^*}
    &= \nabla^2T^*,
    \label{eq:temp'}
\end{align}
where two non-dimensional numbers Ra and Pr are the Rayleigh number and
the Prandtl number, respectively. They are defined as
\begin{equation} \begin{array}{ll}
  \displaystyle \textrm{Ra} = \frac{\alpha gd^3\Delta T}{\kappa\nu}, &
  \displaystyle \textrm{Pr} = \frac{\nu}{\kappa}.
\end{array} \end{equation}

Ra is a ratio between the time scale of the thermal conduction and 
 that of the thermal advection.
A large Ra means that the advection is more efficient than 
 the conduction for the heat transport, and causes 
 the thermal convection.
According to the linear stability theory, the convection occurs when 
 Ra exceeds the critical value Ra$_c$ which is order of $10^3$.
If the Ra is substantially large, the convection shifts to
 turbulent flow.

Meanwhile, Pr is a ratio between the thermal diffusivity and 
 the diffusivity of momentum, i.e., the kinematic viscosity.
In such a case of the earth's mantle convection, a large value of
 $\nu$ makes the inertia terms in (\ref{eq:zeta'}) negligible.

%-----------------------------------------------------------------------
\section{Exercises}
%----------
\subsection{Boundary and initial conditions}
Equations (\ref{eq:Poisson}), (\ref{eq:zeta'}) and (\ref{eq:temp'}) are
 solved to determine velocity, vorticity and temperature fields under
 the following boundary and initial conditions.
On the lower and upper walls $(z=0,d)$, the temperature is fixed,
\begin{equation}
  T = \begin{cases}
        T_0+\Delta T & \textrm{($z=0$)} \\
        T_0.         & \textrm{($z=d$)}
      \end{cases}
\end{equation}
The computational domain is rectangular with height (depth) $d$
 and width $l = 2d$.
On all the walls, the fluid slips freely with no viscous stress
 (stress-free condition),
\begin{equation*} \begin{array}{llll}
  w = 0, &
  \displaystyle \PD{u}{z}+\PD{w}{x} = 0. &
  (x=0,2d, & z=0,d)
  \label{eq:BCmec}
\end{array} \end{equation*}
The above mechanical conditions are equivalent to
\begin{equation}\begin{array}{llll}
  \psi = 0, &
  \displaystyle \zeta = \PD{^2\psi}{z^2} = 0. &
  (x=0,2d, & z=0,d)
\end{array}\end{equation}
%All variables are periodically continuous in the horizontal direction
% with a period of the width of the computational domain, $l=2d$,
%\begin{equation} \begin{split}
%  \zeta(x,z,t) &= \zeta(x+l,z,t), \\
%  T    (x,z,t) &= T    (x+l,z,t), \\
%  \psi (x,z,t) &= \psi (x+l,z,t).
%\end{split} \end{equation}

The initial conditions are given as
\begin{equation} \begin{array}{llll}
  \zeta = 0, & T = T_0, & \psi = 0. & \textrm{($t=0$)}
\end{array} \end{equation}
Since the convectional current cannot occur without any horizontal
 temperature gradient as we have seen in (\ref{eq:zeta}),
 a small perturbation in temperature is given
 at the center of the domain,
\begin{equation}
  T(x={d}, z=\frac{d}{2}) = T_0 + 0.1\Delta T.
\end{equation}

In the calculations, the time step $\Delta t$ must satisfy the following
 condition to obtain stable solution,
\begin{equation}
  \Delta t < \textrm{min}
       \left(\frac{1}{   \nu(k_\textrm{max}^2+l_\textrm{max}^2)},
             \frac{1}{\kappa(k_\textrm{max}^2+l_\textrm{max}^2)}\right),
\end{equation}
 where $k_\textrm{max}=\pi/\Delta x, l_\textrm{max}=\pi/\Delta z$,
 and where $\Delta x$ and $\Delta z$ are the spatial grid sizes of
 the horizotal and vertical directions, respectively.

%----------
\subsection{Standard experiment}
For instance, the water in a pot heated from bottom 
 has such parameters as
\begin{equation*}\begin{array}{llllll}
%  \rho_0 = 1000 \textrm{ kg/m$^3$}, &
%  \eta   = 10^{-3} \textrm{ Pa$\cdot$s}, &
  \nu    = 10^{-6} \textrm{ m$^2$/s}, &
  g      = 9.8 \textrm{ m/s$^2$}, &
  \kappa = 1.4 \times 10^{-7} \textrm{ m$^2$/s}, \\ 
  \alpha = 2.1 \times 10^{-4} \textrm{ /K}, &
  \Delta T = 200 \textrm{ K}, &
  d = 0.005 \textrm{ m},
\end{array}\end{equation*}
for which the non-dimensional numbers are Ra$=1.5 \times 10^9$ and
 Pr$=7.1$.
Since Ra clearly exceeds Ra$_c$, the thermal convection occurs in the pot. 

\begin{figure}[t]
  \begin{center}
    \includegraphics[width=10cm,keepaspectratio,clip]{fig/CA5_01.epsf}
    \caption{Result of Example 1.}
    \label{fig:CA5_01}
  \end{center}
\end{figure}

%----------
\section{Appendix}
In the program of the thermal convection, the main contour plot shows
 the stream function $\psi$.
According to the definition (\ref{eq:defpsi}), the convection
 currents occur along the contour lines clockwise around the domain of 
 higher $\psi$, and counterclockwise around that of lower $\psi$.

The plot in the right shows the temperature averaged in horizontal 
 ($x$) direction.
When the thermal transport is performed mainly by the thermal
 conduction, the $x$-averaged temperature linearly increases with depth.
Once the convection takes place, the most part of the fluid layer has
 a constant $x$-averaged temperature, although there exist the vertical
 movements of hot and cold plumes.
Meanwhile, in the top and the bottom of the fluid layer, 
 the $x$-averaged temperature drastically varies to form
 the thermal boundary layers, of which thicknesses depend on
 the Prandtl number, Pr.

Change the parameters, Ra and Pr, and observe the development of 
 the velocity and the temperature fields.

%-----------------------------------------------------------------------
%\end{document}