\documentclass[12pt,twocolumn]{article}
\usepackage{latexsym,amsmath,amssymb,amsfonts}
\usepackage{graphicx}


\begin{document}

\noindent \textbf{Physics 567}\\
\noindent \textbf{Homework 5}\\
\noindent \textit{Advection of Atmospheric Aerosols}\\
\noindent Due Friday, April 26, 2024 \\

Convective atmospheric eddies (thermally driven vortices with the axis parallel
to the  ground) of 1\,km size often form just above the Earth's surface. These
eddies pick up dust  particles and water droplets (aerosols) and carry
(advect) them to overlying mesoscale (10\,km)  eddies that reach up to the
lower boundary of the stratosphere (10-20\,km altitude). In this  assignment,
you will construct an elementary 2D simulation of the trajectory of an aerosol 
particle (trapped?) in an idealized system of steady, 1\,km convective
eddies.\footnote{The real atmosphere has an evolving, complex chain of eddy
systems, including smaller scale (100\,m) turbulent eddies that  serve to
enhance particle diffusion. Aerosols also experience diffusion due to Brownian 
motion. We will ignore these effects.}

The aerosol particle experiences two forces: gravity and aerodynamic
drag. The equation of  motion is
\[
\mathbf{\ddot{x}}=-g\mathbf{\hat{y}} + 
   \frac{1}{2m} C_D A \rho \, \mathbf{\left| w-\dot{x} \right| 
   \, \left( w-\dot{x} \right)}
\]
where $\mathbf{x} = (x, y)$ is the particle position. The $y$-axis points up. The dots represent time derivatives, $A$ is the cross-sectional area of the (spherical) particle, $g$ is the acceleration of gravity, $C_D$ is the drag coefficient, $\rho$ is the density  of the atmosphere, and $m$ is the mass of the particle. Assume that the particle has the density of water. The coefficient of drag for a sphere is highly dependent on the Reynolds  Number (${Re}$):
\[
C_D \approx \frac{24}{Re} \left( 1+\frac{Re^{2/3}}{6} \right),
\qquad Re \lesssim 1000,
\]
\[
\mbox{where}\qquad Re \equiv \frac{\rho v d}{\eta}.
\] 
The speed $v$ is the magnitude of the velocity of the fluid relative to the particle ($v \equiv \left|
\mathbf{w-\dot{x}} \right|$); $d$ is the particle diameter, and $\eta$ is
the viscosity of air.

We need a plausible expression for the convective velocity field
$\mathbf{w}(x,y)$. In fluid flow, conservation of mass is embodied in the
{\it continuity equation},
\[
  \frac{\partial \rho}{\partial t} = -\nabla \cdot (\rho \mathbf{w}).
\]
In our simplified world, fluid flow will be steady  ($\frac{\partial}{\partial
t}=0$) and incompressible ($\rho=$ constant). The velocity field therefore has
zero divergence. This is easily accomplished in 2D by defining $\mathbf{w}$ in
terms of a stream function $\psi(x,y)$ :
\[
w_x = - \frac{\partial \psi}{\partial y},\qquad
w_y = \frac{\partial \psi}{\partial x}.
\]
We choose the stream function 
\[
\psi (x,y)=  
a w_0 \cos \left( \dfrac{\pi x}{2a} \right)\,
   \sin \left( \dfrac{\pi y}{2a}\right).
\]
The contours of  the stream function are the streamlines (field lines
of the velocity field):
\begin{center}
\includegraphics[width=\columnwidth]{streamfxn.png}
\end{center}
 Note that the vortex directions alternate in a checkerboard pattern. This is
(crudely) what convection cells look like in cross-section. Note that the period of circulation for a given air molecule depends on which streamline it is on. With that caveat, we somewhat arbitrarily define the  convective turnover time as
\[ 
\tau \equiv \frac{2 \pi a}{w_0}.
\]


\vspace{16pt} \noindent {\bf Your assignment:}

\begin{enumerate}
\item Run simulations for particle diameters $d = 1, 10, 100,
1000\,\mu$m (microns).  For each aerosol particle, produce a 
plot of the  trajectory $x,y$ for 10
convective turnover times. Use the following initial conditions:
\[
\begin{matrix}
x_0 = 0, & y_0 = 1.5\,\mbox{km},\\
\dot{x}_0=0, & \dot{y}_0=0.
\end{matrix}
\]
\item Find the largest size of particle that can be trapped in the vortex with the above initial conditions. We will deem the particle trapped if it does not reach ground ($y=0$)
within 10 convective turnover times. 
\end{enumerate}


\vspace{16pt}
\noindent Use the following parameter values in your simulation:
\begin{center}
\begin{tabular}{lcc}
\hline
particle density & $\ldots$ & 1000\,kg/m$^3$ \\
air density & $\rho$ & 1\,kg/m$^3$ \\
vortex speed & $w_0$ & 10\,m/s \\
vortex radius & $a$ & 1000\,m \\
air viscosity & $\eta$ & $2 \times 10^{-5}$\,N\,s/m$^2$ \\
gravitational acceleration & $g$ & 10\,m/s$^2$\\
\hline
\end{tabular} 
\end{center}
\vspace{16pt}

Turn in the following by email to the grader:
%\texttt{<kankel@solar.physics.montana.edu>}:
\begin{enumerate}
\item Well-commented source code.
\item Trajectory plots for each of the particle sizes specified, for 10 turnover times.
\item Answer to the question: what is the largest size particle that will remain trapped? Your answer should be within 0.1\,$\mu$m or better. Include program output to substantiate your answer.
% \item A back-of-the-envelope estimate for (3).
\end{enumerate}
The last part will be more tractable if you automate your search. Please be advised that whether a particle is trapped may not be a monotonic function of particle size.
% \textit{Hint:} If your ODE integrator has difficulty getting started, try slightly 
% adjusting the initial position of the particle. 

\end{document}
