\documentclass[12pt]{article}

\renewcommand\thesubsection{\alph{subsection}}

\usepackage{float}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage{fancyhdr}
\usepackage{enumerate}
\usepackage{amsmath}
\usepackage{bm}
\usepackage{cancel}
\usepackage[justification=centering]{caption}
\usepackage{parskip}
\usepackage{array}
\usepackage{tabu}
\usepackage{subcaption}
\usepackage[style=ieee]{biblatex}
\usepackage{setspace}
\usepackage[section]{placeins}
\pagestyle{fancy}
\fancyhf{}
\rhead{CEE 262B Week 5}
\lhead{Section Notes}
\rfoot{\thepage}

\begin{document}
\section{TA Section Agenda}
\begin{itemize}
\item Derive stress distribution, law of the wall, log law, eddy viscosity for turbulent open channel flow, and Rouse profile
\item Lecture notes summary
\end{itemize}

\section{Derivations}


\subsection{Stress distribution}
Because we're still talking about turbulence, we'll start with the RANS equations:

\begin{equation*}
\frac{\partial \overline{u_i}}{\partial t} + \frac{\partial }{\partial x_j} (\bar{u_i} \bar{u_j}) = -\frac{1}{\rho_0} \frac{\partial \overline{p}}{\partial x_i} -\frac{g}{\rho_0} \overline{\rho} \delta_{i3} + \nu \frac{\partial^2 \overline{u_i}}{\partial x_j \partial x_j} + \frac{1}{\rho_0}\frac{\partial \tau_{ij}}{\partial x_j}
\end{equation*}

Here, we have defined the Reynolds stress tensor as:
\begin{equation*}
\tau_{ij} = -\rho_0 \overline{u_i' u_j'}
\end{equation*}


For 1D, open-channel flow, we can assume that $\overline{u_2} = \overline{u_3} = 0$, and that $\overline{u_1} = \overline{u_1} (x_3)$. So the RANS equations reduce to:
\begin{equation*}
\frac{\partial \overline{u_1}}{\partial t} + \frac{\partial }{\partial x_j} (\bar{u_1} \bar{u_j}) = -\frac{1}{\rho_0} \frac{\partial \overline{p}}{\partial x_1} + \nu \frac{\partial^2 \overline{u_1}}{\partial x_j \partial x_j} + \frac{1}{\rho_0} \frac{\partial \tau_{1j}}{\partial x_j}
\end{equation*}

Expanding, and crossing out terms (let's justify all of them):

\begin{align*}
\cancel{\frac{\partial \overline{u_1}}{\partial t}} &+ \cancel{\frac{\partial }{\partial x_1} (\bar{u_1} \bar{u_1})} + \cancel{\frac{\partial }{\partial x_2} (\bar{u_1} \bar{u_2})}  + \cancel{\frac{\partial }{\partial x_3} (\bar{u_1} \bar{u_3})} = -\frac{1}{\rho_0} \frac{\partial \overline{p}}{\partial x_1} + \nu \frac{\partial^2 \overline{u_1}}{\partial x_3^2} +\\ 
&\cancel{ \frac{1}{\rho_0}\frac{\partial \tau_{11}}{\partial x_1}}  + \cancel{ \frac{1}{\rho_0}\frac{\partial \tau_{12}}{\partial x_2}} +  \frac{1}{\rho_0}\frac{\partial \tau_{13}}{\partial x_3}
\end{align*}

So we're left with:
\begin{equation*}
\frac{1}{\rho_0} \frac{\partial \overline{p}}{\partial x_1} = \nu \frac{\partial^2 \overline{u_1}}{\partial x_3^3} +  \frac{1}{\rho_0}\frac{\partial \tau_{13}}{\partial x_3}
\end{equation*}

Breaking the RHS into a viscous, and turbulent stress component, i.e. 
\begin{equation*}
\tau = \mu \frac{\partial \overline
{u_1}}{\partial x_3} + \tau_{13}
\end{equation*}

We obtain:
\begin{equation*}
\frac{\partial \overline{p}}{\partial x_1} = \frac{\partial \tau}{\partial x_3}
\end{equation*}

The goal now is to use this to find an expression for the stress, $\tau$. To do this, we first assume a form for the pressure. Let's try hydrostatic, i.e.
\begin{gather*}
\frac{\partial \overline{p}}{\partial x_3} = -\rho_0 g \\
\overline{p} = \rho_0 g (h-x_3)\\
\frac{\partial \overline{p}}{\partial x_1} = \rho_0 g \frac{\partial h}{\partial x_1}
\end{gather*}

So then the equation becomes:
\begin{equation*}
\frac{\partial \tau}{\partial x_3} = \rho_0 g \frac{\partial h}{\partial x_1}
\end{equation*}

Integrating both sides over the depth, from 0 to H:
\begin{gather*}
\cancel{\tau(H)} - \tau(0) = \rho_0 g H \frac{\partial h}{\partial x_1}\\
\end{gather*}

Rearranging:
\begin{gather*}
\frac{\tau (0)}{\rho_0} = u_*^2 = -g H \frac{\partial h}{\partial x_1}\\
\frac{\partial h}{\partial x_1} = - \frac{u_*^2}{gH}
\end{gather*}
So plugging this back in to solve for $\tau$:
\begin{gather*}
\tau = \rho_0 g \frac{\partial h}{\partial x_1} x_3 + C\\
\tau= -\rho_0 u_*^2 \frac{x_3}{H} + C
\tau(0) = C = \rho_0 u_*^2\\
\boxed{\tau = \rho_0 u_*^2 \left(1-\frac{x_3}{H} \right)}
\end{gather*}

\subsection{Velocity in near-wall region}
Close to the wall we expect turbulent stresses to vanish, i.e. $\overline{u_1' u_3'} \rightarrow 0$ as $x_3 \rightarrow 0$. This leaves us with only the viscous stress term: 
\begin{gather*}
\nu \frac{\partial \overline{u_1}}{\partial x_3} = u_*^2\\
\overline{u_1} = \frac{u_*^2 x_3}{\nu}\\
\frac{\overline{u_1}}{u_*} = \frac{u_* x_3}{\nu}\\
\end{gather*}

Now if we define a viscous sublayer thickness $\delta = \nu/u_*$, and a new ``wall unit" height $z^+ = z/\delta$, then we have our velocity distribution: 

\begin{equation*}
 \boxed{\frac{\overline{u_1}}{u_*} = z^+}
\end{equation*}

\subsection{Log law}
Further from the wall, we expect viscous stresses to vanish, and the stress to be given as:

\begin{equation*}
\frac{\tau_{13}}{\rho_0} = -\overline{u_1' u_3'} 
\end{equation*}

Now we invoke Prandtl's Mixing Length Hypothesis, which is: 
\begin{equation*}
|u_1'| \sim l_m \left|\frac{\partial u_1}{\partial x_3}\right| \qquad |u_3'| \sim l_m \left|\frac{\partial u_1}{\partial x_3}\right|
\end{equation*}

Or:

\begin{equation*}
-\overline{u_1' u_3'} = l_m^2 \left(\frac{\partial \overline{u_1}}{\partial x_3}\right)^2
\end{equation*}

Next, we assume that the mixing length is given by:
\begin{equation*}
l_m = \kappa x_3
\end{equation*}
where $\kappa = 0.41$ is the von Karman constant. Putting this all together yields: 

\begin{gather*}
\frac{\tau_{13}}{\rho_0} = \kappa^2 x_3^2 \left(\frac{\partial \overline{u_1}}{\partial x_3}\right)^2 = u_*^2\\
\kappa x_3 \left(\frac{\partial \overline{u_1}}{\partial x_3}\right) = u_*\\
\frac{\overline{u_1}}{u_*} = \frac{1}{\kappa} \ln x_3 + c
\end{gather*}
If we enforce that $\overline{u_1} = 0$ at some ``rougness height" $z_0$, then we end up with the standard form of the log law: 
\begin{equation*}
\boxed{\frac{\overline{u_1}}{u_*} = \frac{1}{\kappa} \ln\left(\frac{x_3}{z_0}\right)}
\end{equation*}

\subsection{Eddy Viscosity}
Finally, we'll derive the equation for eddy viscosity. This requires the ``eddy viscosity model":
\begin{equation*}
\frac{\tau_{13}}{\rho_0} = -\overline{u_1' u_3'} \approx \nu_T \left(\frac{\partial \overline{u_1}}{\partial x_3}\right) 
\end{equation*}
Where $\nu_T$ is the eddy viscosity. 

So, going back to the split version of our stress for channel flow:
\begin{gather*}
\frac{\tau}{\rho_0} = \nu \frac{\partial \overline{u_1}}{\partial x_3} + \frac{\tau_{13}}{\rho_0} = u_*^2\left(1 - \frac{x_3}{H}\right) \\
(\nu + \nu_T) \frac{\partial \overline{u_1}}{\partial x_3} = u_*^2\left(1 - \frac{x_3}{H}\right)
\end{gather*}

Now we assume a log-law for the velocity, and ignore the kinematic viscosity because viscous effects are assumed small when deriving the log law:

\begin{gather*}
\nu_T \left(\frac{u_*}{\kappa x_3}\right) = u_*^2\left(1 - \frac{x_3}{H}\right)\\
\boxed{\nu_T = \kappa u_* x_3 \left(1 - \frac{x_3}{H}\right)}
\end{gather*}
 
 If we average this over the depth, we end up with an equation you will soon be very familiar with: 
\begin{equation*}
\overline{\nu_T} = 0.067 u_* H
\end{equation*}

\subsection{Application: Rouse Profile} 
We now have (semi)-analytic equations for the velocity and eddy viscosity profiles, which should tell us almost everything we need to know about how things are mixed around in the water column. For an application of this, consider a water column full of sediment. We can write the full, 3d suspended sediment concentration (SSC) transport equation as:

\begin{equation*}
\frac{\partial C}{\partial t} + \frac{\partial}{\partial x}(uC) + \frac{\partial}{\partial y}(vC) + \frac{\partial}{\partial z}((w-w_s)C) = 0
\end{equation*}

Here, $w_s$ is the (assumed constant) sediment fall velocity. Because we're in a turbulent flow, we can separate the variables into mean and fluctuating components:

\begin{equation*}
\frac{\partial (\overline{C} + C')}{\partial t} + \frac{\partial}{\partial x}((\overline{u}+u')(\overline{C}+C')) + \frac{\partial}{\partial y}((\overline{v}+v')(\overline{C}+C'))  + \frac{\partial}{\partial z}(((\overline{w}+w')-w_s)(\overline{C}+C')) = 0
\end{equation*}

Expanding, and averaging each terms (i.e. standard RANS approach): 
\begin{gather*}
\frac{\partial \overline{ (\overline{C} + C')}}{\partial t} + \frac{\partial}{\partial x}(\overline{\overline{u}\overline{C}+u'\overline{C} + C'\overline{u} + u'C'}) + ...\\
+ \frac{\partial}{\partial y}(\overline{\overline{v}\overline{C}+v'\overline{C} + C'\overline{v} + v'C'}) ++ \frac{\partial}{\partial z}(\overline{\overline{w}\overline{C}+w'\overline{C} + C'\overline{w} + u'C' -w_s\overline{C} - w_sC'})  = 0
\end{gather*}

The average of a fluctuation is identically 0, and the average of an average is simply an average, so we can reduce all of this to:

\begin{equation*}
\frac{\partial \overline{C}}{\partial t} + \frac{\partial}{\partial x}(\overline{u}\overline{C} + \overline{u'C'}) + \frac{\partial}{\partial y}(\overline{v}\overline{C} + \overline{v'C'}) + \frac{\partial}{\partial z}(\overline{w}\overline{C} + \overline{w'C'} - w_s\overline{C}) = 0 
\end{equation*}

Rearranging so that the Reynolds flux terms are on the right hand side:

\begin{equation*}
\frac{\partial \overline{C}}{\partial t} + \frac{\partial}{\partial x}(\overline{u}\overline{C}) + \frac{\partial}{\partial y}(\overline{v}\overline{C}) + \frac{\partial}{\partial z}(\overline{C}(\overline{w} - w_s)) =  -\frac{\partial}{\partial x}(\overline{u'C'}) - \frac{\partial}{\partial y}(\overline{v'C'}) -\frac{\partial}{\partial z}(\overline{w'C'})
\end{equation*}

Directly applying this to steady, 1D, open channel flow, we can safely say that $\partial/\partial x = \partial/\partial y = \partial/\partial t = 0$. This reduces our equation to: 

\begin{equation*}
\frac{\partial}{\partial z}(\overline{C}(\overline{w} - w_s)) = -\frac{\partial}{\partial z}(\overline{w'C'})
\end{equation*}

Integrating both sides yields a balance of vertical settling and a vertical sediment flux:
\begin{equation*}
\overline{C}(\overline{w} - w_s) = -\overline{w'C'}
\end{equation*}

Here, we invoke the eddy diffusivity model, which is a direct analogue of the eddy viscosity model, but relating to mass fluxes rather than momentum fluxes. We can write the turbulent Reynolds sediment flux in terms of a turbulent eddy diffusivity $\epsilon_t$ as:

\begin{equation*}
\overline{C}(\overline{w} - w_s) = \epsilon_t \frac{\partial \overline{C}}{\partial z}
\end{equation*}

Letting $(\overline{w} - w_s) = -W_s$, we can drop all of the overbars and write this equation as:

\begin{equation*}
\epsilon_t \frac{\partial C}{\partial z} + W_s C = 0
\end{equation*}

Here is where we apply what we derived in the previous section. The turbulent eddy diffusivity can be written in terms of the turbulent eddy viscosity and the turbulent Schmidt number, $\sigma_t =\nu_t/\epsilon_t$. This is a ratio of the diffusivity of sediment to the diffusivity of momentum. The Schmidt number is to mass as the Prandtl number is to heat. Using this definition:

\begin{equation*}
\frac{\nu_t}{\sigma_t} \frac{\partial C}{\partial z} + W_s C = 0
\end{equation*}

And rearranging:

\begin{equation*}
\frac{\partial C}{\partial z} + \frac{W_s \sigma_t}{\nu_t} C = 0
\end{equation*}

Now we plug in our definition of the parabolic eddy viscosity for turbulent open channel flow:
\begin{equation*}
\frac{\partial C}{\partial z} + \frac{W_s \sigma_t}{\kappa u_* z(1-\frac{z}{h})} C = 0
\end{equation*}

Now we can define the Rouse Number, Ro = $\frac{W_s \sigma_t}{\kappa u_*}$, and rewrite this as:
\begin{equation*}
\frac{\partial C}{\partial z} + \frac{\textrm{Ro}}{z(1-\frac{z}{h})} C = 0
\end{equation*}

This is a linear first-order ODE (the best kind!) with solution:

\begin{equation*}
C(z) = C_a\left(\frac{H-z}{z}\frac{a}{H-a}\right)^{\textrm{Ro}}
\end{equation*}

Here, $C_a$ is a reference concentration at height $a$. Note that this is usually valid for sediment concentrations below 100 g/L. 




\newpage
\section{Turbulence models and physics}
This set of notes is a good reference when solving problems related to turbulent flow. There are quite a few equations, and many of them are semi-empirical. This can be frustrating, but they are still based on physical assumptions, so focus on those and things will make sense! And note our word choice here: these are just turbulence ``models". We will be solving idealized problems, with methods that work surprisingly well considering how complex turbulence is.

\subsection{Eddy Viscosity and the Law of the Wall}

Here we will denote a new parameter, $\nu_t$ [$m^2/s$], called the eddy viscosity. Note that the units are the same as a diffusion coefficient, or kinematic viscosity. The analogy here is that the eddy viscosity describes how momentum is diffused through the water column via turbulent eddies. With this parameter, we can now state the eddy viscosity turbulence model:

\begin{equation*}
-\overline{u'w'} = \nu_t \frac{\partial U}{
\partial z}
\end{equation*}

This says that the $\overline{u'w'}$ Reynolds Stress is equal to the vertical velocity gradient times the eddy viscosity. For open channel flow, eddy viscosity takes the form:

\begin{equation*}
\nu_t = \kappa u_* z \left(1 - \frac{z}{h}\right)
\end{equation*}

Where $u_*$ is the friction velocity [$m/s$], defined by: 
\begin{equation*}
u_* = \sqrt{\frac{\tau_b}{\rho}}, \qquad \tau \equiv \textrm{Bottom Shear Stress}
\end{equation*}

$\kappa$ is the Von Karman constant, usually taken to be $\approx 0.4$, and $h$ is the water column depth. All told, this gives us the classic parabolic eddy viscosity profile, which reaches a maximum in the center of the water column, and goes to 0 at the top and bottom. 

However, we still need the vertical velocity profile. The notes contain a derivation, but I'll simply write out the result here; it's the Law of the Wall:
\begin{equation*}
U = \frac{u_*}{\kappa}\textrm{ln}(z) + C
\end{equation*}

Where C is a constant determined by wall roughness.

Note that differentiating our velocity profile gives:
\begin{equation*}
\frac{\partial U}{\partial z} = \frac{u_*}{\kappa z}
\end{equation*}

And plugging this back into our eddy viscosity model, we obtain: 
\begin{equation*}
-\overline{u'w'} = u_*^2 \left(1 - \frac{z}{h}\right)
\end{equation*}

We can define things a little differently for a \textbf{smooth wall}. In the near-wall region, our flow is entirely dominated by viscous forces, so we have:

\begin{equation*}
\nu \frac{\partial U}{\partial z} = u_*^2  \qquad U = \frac{u_*^2 z}{\nu}
\end{equation*}

And thus, 
\begin{equation*}
u^+ = z^+
\end{equation*}

Where $u^+ = U/u_*$, and $z^+ = zu_*/\nu$. We usually say that this estimate holds until $z^+ \approx 12$. 

However, in general, we don't have smooth walls. For walls of some roughness, it is useful to define a ``roughness length", denoted by $z_0$, which is the length at which the \textit{mean velocity} is 0. This is often determined empirically. Using this roughness length, we can write the Law of the Wall as: 

\begin{equation*}
U = \frac{u_*}{\kappa}\textrm{ln}\left(\frac{z}{z_0}\right)
\end{equation*}

Note that this the same as what we had before, we simply account for the integration constant C with the empirical roughness length. 

Another way that we define roughness is via a bottom drag coefficient, $C_D$:

\begin{equation*}
C_D = \frac{u_*^2}{U_{ref}^2}
\end{equation*}

So rearranging the law of the wall at a reference height $z_{ref}$, we obtain: 

\begin{equation*}
C_D = \frac{\kappa^2}{\textrm{ln}^2\left(\frac{z_{ref}}{z_0}\right)}
\end{equation*}

One more thing to mention here. Recall from the last lecture that we often assume turbulent production and dissipation are in balance, i.e. 
\begin{equation*}
P = -\overline{u'w'} \frac{\partial U}{\partial z} \approx \epsilon
\end{equation*}
Near a wall, we have  $-\overline{u'w'} \approx u_*^2$, and $\frac{\partial U}{\partial z} = u_*/\kappa z$, which gives:

\begin{equation*}
P \approx \epsilon \approx \frac{u_*^3}{\kappa z} 
\end{equation*}

\subsection{Complications}

\begin{itemize}
\item Waves: Surface waves enhance bottom stress. We can see these effects by looking at at the bottom drag coefficient $C_D$ as a function of mean velocity/orbital velocity (see plot in lecture slides)
\item Stratification: Stratification damps/suppresses turbulence, i.e. reduces the effective eddy viscosity. We can see these effects through the Richardson Number 
\end{itemize}

Recall that the Richardson Number is defined by:

\begin{equation*}
\textrm{Ri} = \frac{N^2}{S^2}
\end{equation*}

Where we define the square of the Buoyancy Frequency as:
\begin{equation*}
N^2 = -\frac{1}{\overline{\rho}}\frac{\partial \overline{\rho}}{\partial z} g
\end{equation*}

And the mean shear strength as: 
\begin{equation*}
S^2 = \left(\frac{\partial U}{\partial z}\right)^2
\end{equation*}

So, a higher Ri means more density stratification, and thus, more suppression of turbulence. We'll examine this in much greater detail when we talk about mixing and circulation in estuaries. 
\newpage

\section{SFD with Turbulence}

Here we will combine the last two topics, shear flow dispersion and turbulence, and look at the case where they are occuring together. This is a pretty straightforward extension, though there are some modifications to the velocity profile (and as always, notation) that we need to take into account. 

\subsection{Notation and Governing Equations}

First, notation. When talking about turbulent SFD, we use $\epsilon$, [$m^2/s$] to denote the turbulent eddy diffusivity. In addition, variables with overbars (e.g. $\overline{u}, \overline{C}$) denote cross-sectional averages of Reynolds averaged variables, instead of the Reynolds averaged variables themselves. Similarly, $u'$ and $C'$ refer to deviations from those cross-sectional averages. Now, moving onto the governing equations: 

\begin{equation*}
u' \frac{\partial \overline{C}}{\partial \xi} = \frac{\partial}{\partial y} \epsilon (y) \frac{\partial C'}{\partial y}
\end{equation*}

Here, $\xi$ is our spatial variable in a translated reference frame, defined by: 
\begin{equation*}
\xi = x - \overline{u}t
\end{equation*}

And our equation for the turbulent dispersion coefficient:

\begin{equation*}
K = -\frac{1}{h}\int_0^h u'(y) \int_0^y \frac{1}{\epsilon (y)} \int_0^y u'(y)\, dy\, dy\, dy
\end{equation*}
Note that this is very similar to the laminar SFD equation, except our turbulent eddy diffusivity is now a function of our traverse directional variable, $y$, so we can't pull it out of the integral. 

\subsection{Timescales}
Recall from the Laminar SFD section we were able to calculate a time after which the above governing equations hold. This time scale will be dependent on the geometry of our problem, and the turbulent diffusivity, $\epsilon$. In general we take the time scale at which SFD applies as:

\begin{equation*}
T = 0.4 \frac{d^2}{\overline{\epsilon}}
\end{equation*}

Where $d$ is some relevant geometric length.

For open channel flow, we use the relation:
\begin{equation*}
\overline{\epsilon_z} = 0.07 u_* H 
\end{equation*}

Where $u_*$ is the friction velocity, and H is the channel depth. 

For pipe flow, we use the relation:

\begin{equation*}
\overline{\epsilon} = 0.054 u_* a
\end{equation*}

Where $a$ is the pipe radius. 

\subsection{Idealized Cases}

Luckily, we don't have to compute that integral every time we want to estimate a dispersion coefficient, because other people have done it for us. Here are a few idealized cases that you may find useful

\subsubsection{Turbulent pipe flow}
\begin{equation*}
K = 10.1 a u_*
\end{equation*}

Where $a$ is the pipe radius.

\subsubsection{Laminar pipe flow}
\begin{equation*}
K = \frac{1}{192}a^2 \frac{u_0^2}{D}
\end{equation*}

Where $u_0$ is the maximum velocity along the centerline, and $D$ is the diffusion coefficient. 

\subsubsection{Turbulent flow, depth d, down inclined plane}
This is often called ``Elder's Result", and is given by:
\begin{equation*}
K = 5.93 d u_*
\end{equation*}

This is derived assuming a logarithmic velocity profile.

\subsubsection{General Longitudinal Dispersion, Unidirectional Shear Flow}

We can generalize Elder's result with the following formula:

\begin{equation*}
K = \frac{h^2 \overline{u'^2}}{E} I
\end{equation*}

Here, $h$ is the dimension in the direction of velocity variation, e.g. the transverse y-direction of we have a velocity profile $u(y)$. We define $E$ as the cross sectionally averaged $\epsilon$:
\begin{equation*}
E = \frac{1}{h} \int_0^h \epsilon (y) dy
\end{equation*}

And $I$ is the dimensionless integral:
\begin{equation*}
I = - \int_0^1 u'' \int_0^{y'} \frac{1}{\epsilon'}\int_0^{y'} u'' \,dy\,dy\,dy
\end{equation*}

Just a couple more definitions now...
\begin{equation*}
u'' = u'/\sqrt{\overline{u'^2}}
\end{equation*}

Think of $u''$ as the velocity perturbation normalized by the intensity of the velocity perturbations. 

And finally, 

\begin{equation*}
\epsilon'=\epsilon/E
\end{equation*}

Often times we just take \textbf{$I=0.1$}, but there is also a table in your notes and the textbook that gives slightly more accurate values depending on the circumstances. As you can see in the example in the notes, this method may give similar results to Elder, but doesn't assume the logarithmic velocity profile \textit{a priori}.

\end{document}