\input{aefm-macro}


\setcounter{chapter}{6}
\def\bu{{\mathbf u}}
\def\bq{{\mathbf q}}
\setcounter{section}{5}
\section{ Rayleigh-Darcy (or Horton-Rogers-Lapwood) instability in a porous layer}
\setcounter{equation}{0}
6-6-Lapwood.tex\hfil\today

Nield \& Bejan, Chapter 6 {\bf Convection in Porous Media}\\
Related: Rayleigh-Bernard Problem (Chandrsekhar, Chapter II, {\bf Hydrodynamic and Hydromagnetic Stability}) 

\vskip0.3cm
If  a layer of viscous fluid is heated from below, instability can occur and leads to convection cells important in meteorology. (Rayleigh-Benard Problem).\\
If a saturated porous layer is heated from below, similar instabilily and convection can occur. This is   of basic interest to geothermal convection and is relevant to the complex problem of heat transport due to the burial of nuclear waste in mountains or in a seabed. 

To give some visual ideas of what can happen in porous media, we shall borrow  some photographic evidence for the mathematically similar Rayleigh-Benard  problem of a pure fluid layer heated form below.  

 
\begin{figure} [h]\begin{center}
\includegraphics[scale=1]{geotherm.eps}\end{center}
\caption{A saturated porus medium in geothermal gradient }
 \label{fig:geotherm}
\end{figure} 


 
\underline{Static state}
\[\bu = 0, \quad \f{\p }{\p t }=0 \]
Hence \be T_S = T_o +\Delta T \lp 1-\f{z}{h} \rp\ee
Static equilibrium:
\[ \f{d p_S}{dz}=-\rho_og \lb T_o-{\beta \Delta T}\lp 1-\f{z}{h}\rp \rb \]
hence
\be p_S = p_o -\rho_o g   \lcb T_o z + \f{\beta  \Delta T}{2} \lp \f{z^2}{h} - 2z \rp
\rcb\ee
 
Consider the  \underline{ perturbed state of small disturbances}: 
\be \bu
= 0+ 
\bu',
\quad T=T_S + T',
\quad P= p_S + p'
\ee
then 
 \be \nabla \cdot \bu' = 0  \ee
\be 0=-\nabla p' -\f{\mu}{k}\bu' + \beta g\rho_o T' {\bf k}
\label{2}\ee
 \be  (\rho C)_m \f{\p T'}{\p t} + (\rho C)_f \bu' \cdot \nabla T_S= K_m
\nabla^2 T' \label{3}\ee

\subsection{Non-dimensionalization}.  Let $K_m, \kappa_m =K_m/(\rho C)_f$  be the conductivity and diffusivity of the mixture, and $k$ the permeability. Define
\begin{eqnarray}
&&(x,y,z) \to h(x^*,y^*,z^*), \, \,  t\to\f{\sigma h^2}{ \kappa_m}t^*
,\nonumber \\ &&\bu' \to \f{\kappa_m}{h} \bu^*,\, \, 
T'=
\Delta T \theta , \, \, p'\to  \f{\mu \kappa_m}{k} p^*
\end{eqnarray}
Then, after omitting ${}^*$ for
brevity, we get 
  \be \nabla \cdot \bu  = 0\label{1} \ee
\be 0=-\nabla p  - \bu  + Ra\, \theta {\bf k}  \label{2}\ee
 \be   \boxed{\f{\p \theta }{\p t} -   w
=  \nabla^2 \theta } \label{3}\ee
where  
\be Ra= \f{\rho_f g k \beta \Delta T h}{\mu \kappa_m} ~~~\mbox{Rayleigh  number in a porous medium} 
\ee is the Rayleigh number (ratio of buoyancy force to  diffusive resistence) of the porous medium. In Benard's problem, Rayleigh number is defined as \be Ra = \f{\rho g \Delta T h^4}{\mu \kappa},~~~\mbox{Rayleigh  number in a pure fluid} 
\ee
 Sometimes one calls 
\be D=\f{k}{h^2}, \mbox{Darcy number}\ee   Darcy number so that Rayleigh number of a porous medium is the product of the traditional Rayleigh number and Darcy number.


Taking the curl of (\ref{2}),
\be  \nabla\times \bu =  Ra ({\bf i}\theta_{y}- {\bf j}\theta_{x} )\label{vorticity}\ee
Note that the z compoent of the vorticity vector is zero.

Taking the curl again and using 
\[\nabla \times \nabla \times \bu = \nabla
(\nabla \cdot \bu ) - \nabla ^2 \bu \]
we get \[ \nabla^2 \bu = -Ra\lb {\bf i } \theta_{xz} + {\bf j}\theta_{yz} - {\bf
k} (\theta_{xx} + \theta_{yy} )\rb\]
Taking the $z$ component, we get
 \be \boxed{\nabla^2 w  = Ra \nabla^2_H\theta }\label{4}
\ee
where \be \nabla^2_H=\f{\p^2}{\p x^2}+\f{\p^2}{\p y^2}\ee is the horizontal Laplacian. 

Equations (\ref{3}) and (\ref{4})  couple the two unknowns $w$ and $\theta$. The boundary conditions are 
\be w = \theta = 0, \quad z = 0, 1\ee
After they  are solved the other velocity components and pressure can be found.


\subsection{Solution for sinusoidal disturbances}
 Let 
\be (w,\theta ) =(W(z),\Theta(z)) \exp (ilx + imy -i\omega t) \ee
and \[ D = \f{d}{ d z }\]
then from (\ref{3}),
\be -i \omega \Theta -   W = (D^2 -a^2) 
\Theta\label{11} \ee
and from (\ref{4})
 \be  (D^2-a^2)W = -  a^2 Ra \Theta \label{10}\ee
where \be a^2=\ell^2+m^2 \ee
The boundary conditions  are :
\be W=\Theta = 0, \quad z =
0,1\ee 
Because the equations and the boundary conditons are
homogeneous, the problem for $W$ and $\theta$ is an eigenvalue
problem. 

Note that $l, m,  a$  are related to dimensional wave numbers by
\be l=k_x h=\f{2\pi h}{L_x}, ~~m=k_y h=\f{2\pi h}{L_y}, a=kh={2\pi h\lp
\f{1}{L_x^2}+\f{1}{L_y^2} \rp^{1/2}} ~~\mbox{with}~~k=\sqrt{k_x^2+k_y^2}\ee
  

\subsection{Principle of exchange of stabilities}
We shall first show  that
$\omega$  must be purely imaginary. 
Multiplying (\ref{10}) by $W^*$ and integrating from $z=0$ to $z=1$,
we get, after partial  integration and using the boundary conditions,
\be -\int_0^1 (|DW|^2 + a^2 |W|^2) dz = -a^2Ra \int_0^1 W^*\Theta dz
\label{eq31}
\ee Similarly we multipy (\ref{11}) by $\Theta^*$ and integrating from
$z=0$ to
$z=1$, and  get 
\be
-\int_0^1 (|D\Theta|^2 + a^2 |\Theta|^2) dz = -i\omega\int |\Theta|^2 dz -
\int_0^1 W \Theta^* dz
\ee
Taking the complex conjugate of the second equation 
\be -\int_0^1 (|D\Theta|^2 + a^2 |\Theta|^2) dz = i\omega^*\int |\Theta|^2 dz -
\int_0^1 W^* \Theta dz\label{eq33}
\ee
Eqs (\ref{eq31}) and (\ref{eq33}) can be combined by eliminating the
cross product terms, ,  
 \be -\int_0^1 (|D\Theta|^2 + a^2 |\Theta|^2) dz = -i\omega^*\int |\Theta|^2 dz  
-\f{1}{ a^2Ra}\int_0^1 (|DW|^2 + a^2 |W|^2) dz
\ee
Since all integrals above are real, $-i\omega* = -\omega_i-i\omega_r$ must also be real. We conclude that
\be \omega_r = 0, ~~\mbox{hence}~~
 -i\omega= -\omega_i \ee

 Marginal stability (the threshold of instability) occurs at
$\omega_r=\omega_i = 0$.  If $\omega _i >0$,  the static state is unstable;
$\omega_i = 0$,  marginally stable;  if $\omega_i < 0$, stable.   A problem where the eigenfrequency is real so that marginal instabililty occurs when $\omega=0$ is said to obey the {\it principle of exchange of stabilities}. 

\subsection{Solution to eigenvalue problem}
Consider the situation at marginal stability : $\omega = 0$,
\be  (D^2-a^2)W = -  a^2 Ra\Theta \ee
\be  -   W = (D^2 -a^2) 
\Theta \ee
Eliminating $\Theta$, we get
\be (D^2-a^2)^2 W =    a^2 Ra W \label{gov-w}\ee
subject to 
\be W=0,~~ D^2 W= 0, \quad z = 0, 1\ee
Expanding (\ref{gov-w}) 
\be   D^4 W-2a^2D^2W+a^4 W=a^2RaW\label{gov-w-2} \ee
               Clearly $D^4W=0$ on $z=0,1$. Differentiating  (\ref{gov-w-2}) twice we see that $D^6=0$ on $z=0,1$. Repeating the process we find 
\be D^{(2m)}=0, ~~m = 1,2,3,\cdots,~~~\mbox{on} ~~z=0,1\ee
Therefore the eigensolution  must be
\be W\sim \sin j \pi z\ee
  To satisfy (\ref{gov-w})
it is necessary that 
\be   Ra = \f{[j^2\pi^2 + a^2 ]^2}{a^2}, ~~\mbox{for} ~ j = 1,2,3...,  
  \ee
 which is the eigenvalue condition. For any $j$, $Ra$ becomes unbounded for both $a^2\to 0$ and $a^2\to \infty$ and is curve concave upward in the plane of $a^2$(abscissa)  vs,  $Ra$ (ordinate).  


The lowest threshold 
occurs at
$j=1$, and
\[
\f{\p Ra}{\p a^2 }= 0\] i.e., 
\be a^2 = \pi^2 \ee
or
\be  Ra_c = 4\pi^2 = 39.48\ee


\subsection{Possible convection patterns} This is similar to Benard's problem which has been exhaustively studied theoretically and experimentally. There are 
many possibilities. Let us consider the lowest mode only with $j=1$.

\underline{2-Dimensinal Rolls} : ($\ell = \pi, m = 0$) 

Take 
\be  w = \cos \pi x \sin \pi z \ee 
then from mass conservation,
\be  u = -\sin \pi x \cos \pi z \ee 
The dimensionless wavelength is $L_x=L_y=2$. Along lines $x=0, \pm n, n=1,2,3,4, ...$,  $u=0$ but $w\neq 0$.  Along $x=0, \pm 2m$, $w>0$, hence fluid rises vertically. Along $x=\pm 2m-1$, $w<0$ hence fluids sinks vertically. Along $z=0 \, \mbox{and}\, 1,  w=0$. On the bottom ($z=0$),  $u>0$ while on the top ($z=1$),  if $0<x<  1, 2<x<3, 4<x<5, ...$ . The streamlines are shown in Figure \ref{fig:rolls}.

\begin{figure} [h]\begin{center}
\includegraphics[scale=0.75]{rolls.eps}\end{center}
\caption{Rolls in a period }
 \label{fig:rolls}
\end{figure} 

Experimetal images of rolls in a pure fluid (Benard problem) is shown in Figure \ref{fig:rolls-vd}.
\begin{figure} [h]
\begin{center}
\includegraphics[scale=0.5]{rolls-vd.eps}\end{center}
\caption{Rolls in a pure fluid layer heated from below.  From van Dyke: {\bf An Album of Fluid Motion}}
 \label{fig:rolls-vd}
\end{figure} 

\underline{Rectangular  cells}: ($\ell = m=\pi/\sqrt{2}$). 
\be
w=\cos\f{\pi x}{\sqrt{2}} \cos\f{\pi y}{\sqrt{2}} \sin
\pi z \ee
From the $z$ component of the vorticity equation (\ref{vorticity}), 
\be \f{\p v}{\p x} - \f{\p u}{\p y} = 0 \ee
and from continuity
\be \f{\p u}{\p x}+\f{\p v}{\p y} = -\f{\p w}{\p z}\ee 
By cross differentiation, we get 
\be \f{\p^2 u}{\p x^2}+\f{\p^2 u}{\p y^2} = -\f{\p^2 w}{\p z \p y} = \f{\pi^2}{2\sqrt{2}}\sin\f{\pi x}{\sqrt{2}}\cos \f{\pi y}{\sqrt{2}}\cos \pi z\ee
and 
\be \f{\p^2 v}{\p x^2}+\f{\p^2 v}{\p y^2} = -\f{\p^2 w}{\p z \p x} = \f{\pi^2}{2\sqrt{2}}\cos\f{\pi x}{\sqrt{2}}\sin \f{\pi y}{\sqrt{2}}\cos \pi z\ee
These are easily solved to give 
\be  u  = -\f{1}{\sqrt{2}}\sin\f{\pi x}{\sqrt{2}} \cos\f{\pi y}{\sqrt{2}} \cos
\pi z \ee
and 
\be  v  = -\f{1}{\sqrt{2}}\cos\f{\pi x}{\sqrt{2}} \sin\f{\pi y}{\sqrt{2}} \cos\pi z \ee
The streamlines in a  horizontal plane is shown in Figure \ref{fig:rect-cell}.

\begin{figure} [h]
\begin{center}
\includegraphics[scale=1]{rectangular-cell.eps}\end{center}
\caption{Rectangular and square cells in a pure fluid heated from below.  From Chandrasekhar}
 \label{fig:rect-cell}
\end{figure} 


 \underline{Hexagonal cells}: See Chandrasekhar. 
See Figure \ref{fig:benard-vd}. 
\begin{figure} [h]\begin{center}
\includegraphics[scale=0.75]{benard-vd.eps}\end{center}
\caption{Hexagonal cells in a pure fluid headted from below.  From van Dyke}
 \label{fig:benard-vd}
\end{figure} 

\end{document}