\documentclass[10pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage{natbib}
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage{float}
\begin{document}
\title{Incompressible Flows About an Object In a Driven Vertically Repeating
Rectangular Domain}
\author[1]{Forrest Bullard}%
\affil[1]{California State University, Chico}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\sloppy
\section*{Introduction}
{\label{668455}}
In this report we will discuss how the Navier-Stokes momentum equation
for incompressible flows can be adapted to a numerical solution so that
analysis of fluid flow around a fixed object can be made. We will
discuss a simplification of the Navier-Stokes equation by a change of
variable from velocity to vorticity, the curl of velocity, and the
stream function who's spatial derivatives will return the components of
our velocity for each section of flow. As well we will examine what
boundary conditions are necessary in order to create a flow through our
domain and other conditions which will define the boundaries of our
object. Finally we will discuss some of the limitations of this
examination.
\section*{Math}
{\label{724577}}
In our derivation we began with the incompressible Navier-Stokes
equation of the form,
\par\null
\begin{equation} \label{eq:1}
\frac{d{\bf u}}{dt} = -\frac{1}{\rho}\nabla{p} + {\bf g} + \nu\nabla^2{\bf u}
\end{equation}
we will then take the curl of this formula to change it into the
vorticity form.~
\begin{equation} \label{eq:2}
\nabla \times [\frac{d{\bf u}}{dt} = -\frac{1}{\rho}\nabla{p} + {\bf g} + \nu\nabla^2{\bf u}]
\end{equation}
The curl of the first two terms on the right side of this equation will
both be zero as the curl of the gradient of any scalar function is zero.
As long as~\textbf{g} is conservative it can be expressed as the
gradient of a scalar. Here we must also note that the gradient of the
curl of any vector field is zero. Expanding the left side of this
equation gives,
\begin{align*}
\nabla \times [\frac{\partial {\bf u}}{\partial t} + ({\bf u}\cdot \nabla){\bf u}] = \frac{\partial {\bf \omega}}{\partial t} + \nabla \times [({\bf u} \cdot \nabla){\bf u}] = \frac{\partial {\bf \omega}}{\partial t} + \nabla \times [\nabla(\frac{\bf u \cdot u}{2}) + {\bf \omega} \times {\bf u}]
\end{align}
\begin{align*}
= \frac{\partial {\bf \omega}}{\partial t} + \nabla \times ({\bf \omega} \times {\bf u})
\end{align}
So our equation reduces to,
\begin{align*}
\frac{\partial {\bf \omega}}{\partial t} + \nabla \times ({\bf \omega} \times {\bf u}) = \nu\nabla^2{\bf \omega}
\end{align}
Using a few vector identities gives us the form of our equation before
the use of the stream function. Namely,\(\)
\begin{equation}
\frac{d{\omega}}{dt} = (\omega \cdot \nabla){\bf u} + \nu\nabla^2\omega
\end{equation}
We now have the field equations for the vorticity in a fluid with
constant~\(\rho\).~ We can see that~~\(\nu\nabla^2u\)~ will
represent the rate of change of the~\(\omega\) caused by
diffusion of vorticity and~~\(\left(\omega\cdot\nabla\right)u\)~ representing the rate of
vorticity caused by the stretching and tilting of vortex lines. We can
also see that the central body forces for gravity and pressure are not
found in this equation as they cause no torques. To relate this equation
to the equation we will use to numerically solve our system we will
relate the vorticity and velocity functions by the stream function such
that,
\par\null
\begin{equation}
\frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2\psi}{\partial y^2} = -\omega
\end{equation}
and
\par\null
\begin{equation}
u = \frac{\partial \psi}{\partial y};\ \ \ v = -\frac{\partial \psi}{\partial x}
\end{equation}
We may now write the advection-diffusion equation for the vorticity.
\par\null
\begin{equation}
\frac{\partial \omega}{\partial t} = -\frac{\partial \psi}{\partial y}\frac{\partial \omega}{\partial x} + \frac{\partial \psi}{\partial x}\frac{\partial \omega}{\partial y} + \nu(\frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2})
\end{equation}
This coupled with equation 4, the Poission equation that relates the
stream function to vorticity, gives us the self-contained system we will
need to solve for the flow given the appropriate boundary and initial
conditions. The boundary values for the stream function can be found
along the wall by use of equation 5. This as well shows that for a
stationary object inside our rectangular domain we will need a region
where both the vorticity and stream function is zero. We will also
create repeated boundaries at the top and bottom by solving for the
stream function at the bottom as if the next increment down is the first
zone from the top and then by resolving the stream function at the top
layer with the solutions from the bottom.
\section*{Numerical Approximations}
{\label{769605}}
The process for our numerical approximation will be as stated in Fluid
Mechanics section 6.3
\textbf{1.~}Given~\(\omega_{i,j}\) at all interior points, solve for
\(\psi_{i,j}\).
\textbf{2.} Find the boundary vorticity,~\(\omega_{wall}\).
\textbf{3.} Calculate the vorticity at the new time,~
\(\omega_{i,j}^{n+1}\)~ for all interior points.
\textbf{4.} Set~ \(t\ =t\ +\ \Delta t\)~ and go back to first step.
For step one we will use the equation,\(\)
\begin{equation}
\frac{\psi_{i+1,j}^n+\psi_{i-1,j}^n+\psi_{i,j+1}^n+\psi_{i,j-1}^n-4\psi_{i,j}^n}{h^2}=-\omega_{i,j}^n
\end{equation}
along with a successive over relaxation method so we will converge to a
solution faster. The boundary conditions for the stream function will
not change during the progression of the program but we will have to
calculate for the vorticity function along the wall separately before we
may calculate the next time step for the vorticity function inside. The
formula for solving the vorticity function along the wall will involve a
Taylor expansion at the boundaries and will be of the form,
\(\)
\begin{equation}
\omega_{wall}\ =\ \left(\psi_{i,1}-\psi_{i,2}\right)\frac{2}{h^2}+U_{wall}\cdot\frac{2}{h}+O\left(h\right)
\end{equation}
Once we know the vorticity function on the wall we can use a forward
difference method to approximate the next time step of the vorticity
function inside using the following equation,
\par\null
\begin{align*}
\frac{\omega_{i,j}^{n+1}-\omega_{i,j}^n}{\Delta t}=-\left(\frac{\psi_{i,j+1}^n-\psi_{i,j-1}^n}{2h}\right)\left(\frac{\omega_{i+1,j}^n-\omega_{i-1,j}^n}{2h}\right)+\left(\frac{\psi_{i+1,j}^n-\psi_{i-1,j}^n}{2h}\right)\left(\frac{\omega_{i,j+1}^n-\omega_{i,j-1}^n}{2h}\right)
\end{align}
\begin{equation}
+\nu(\frac{\psi_{i+1,j}^n+\psi_{i-1,j}^n+\psi_{i,j+1}^n+\psi_{i,j-1}^n-4\psi_{i,j}^n}{h^2})
\end{equation}
Note that these equations have been simplified under the assumption that
the spacing between the grid points will be the same in the x and y
direction. It was also necessary to reapply the region of our object
having zero vorticity and zero stream after every cycle.
\section*{Code}
{\label{133107}}
The following was the code, written in python, used to examine the flow
of fluid with different viscosities about a circular object. It was
necessary to find appropriate time steps as related to the size of the
box and to the given viscosity so that convergence could be maintained
as low viscosity does not allow for enough dissipation of the vorticity
function.
\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_
``''" Created on Fri Jan 4 13:19:28 2019
Based on code from Fluid Dynamics sixth edition section 6.3 (code 3)
Solution for Unsteady Two-Dimensional Navier-Stokes equations in
vorticity-stream function form, using an explicit forward in time,
centered in space scheme
Wind tunnel simulation expanded by repeated boundary conditions on
horizontal boundary and net zero flux on vertical boundaries
Succesive Over Relaxation (SOR) is used for generation of stream
function, assuming uniform resolution
@author: fbullard ``''"
from scipy import zeros, linspace, meshgrid, sqrt from numpy import sum
as sumpy import matplotlib.pyplot as plt
\begin{verbatim}
#parameters for SOR
\end{verbatim}
scale=1 ; Ny=60 ; Nx=(Ny*scale)-scale+1 ; MaxStep=1000 ; Visc=.5 ;
dt=0.00005 ; t = 0.0; MaxIt=600 ; Beta=1.5 ; MaxErr=0.001 ;
\begin{verbatim}
#functions and intitial conditions
\end{verbatim}
sf=zeros({[}Ny, Nx{]}) ; vt=zeros({[}Ny, Nx{]}) ; vto=zeros({[}Ny,Nx{]})
; Ui = .1; height=2.0
\begin{verbatim}
#for contour mapping. h is also used for boundary conditions on stream function
\end{verbatim}
x=zeros({[}Ny,Nx{]}) ; y=zeros({[}Ny,Nx{]}) ; h=height/(Ny-1);
\begin{verbatim}
# for creation of shapes
\end{verbatim}
shape = zeros({[} Ny , Nx {]}) xi =
linspace(-(height\emph{scale)/2,(height}scale)/2, Nx , endpoint = True)
yj = linspace(-height/2, height/2, Ny, endpoint = True) xx, yy =
meshgrid(xi, yj, indexing = `xy')
for a in range(Nx): for b in range(Ny): if sqrt(xx{[}b,a{]}\textbf{2 +
yy{[}b,a{]}}2) \textgreater{} 0.3: shape{[}b,a{]} = 1
\begin{verbatim}
#could be replaced by meshgrid
\end{verbatim}
for i in range(Nx): for j in range(Ny): x{[}j,i{]} = h\emph{i y{[}j,i{]}
= h}j
\begin{verbatim}
#generate boundaries for stream function
\end{verbatim}
for j in range(Ny): sf{[}j,0{]} = Ui\emph{j}h sf{[}j,Nx-1{]} =
Ui\emph{j}h
for tstep in range(MaxStep): for SOR in range(MaxIt): vto = sf.copy()
for i in range(1,Nx-1): for j in range(1,Ny-1): sf{[}j,i{]} =
0.25\emph{Beta}(sf{[}j,i+1{]} + sf{[}j,i-1{]} + sf{[}j+1,i{]} +
sf{[}j-1,i{]} + h\emph{h}vt{[}j,i{]}) + (1.0-Beta)*sf{[}j,i{]}
\begin{verbatim}
#for repeated boundary condition solve bottom with solution of third from top then set bottom two rows to top two rows
for i in range(1,Nx-1):
sf[Ny-1, i] = 0.25*Beta*(sf[Ny-1,i+1] + sf[Ny-1,i-1] + sf[2,i] + sf[Ny-2,i] + h*h*vt[j,i]) + (1.0-Beta)*sf[j,i]
sf[:2,1:Nx-1] = sf[Ny-2:Ny,1:Nx-1]
sf *= shape
#check for convergence of the stream function given some maximum error and steps
error = sumpy(abs(sf - vto))
if error <= MaxErr:
break
#generate vorticity function for boundaries
vt[1:Ny-1,0] = (2.0*(sf[1:Ny-1,0] - sf[1:Ny-1,1]))/(h*h)
vt[1:Ny-1,Nx-1] = (2.0*(sf[1:Ny-1,Nx-1] - sf[1:Ny-1,Nx-2]))/(h*h)
vt[0,1:Nx-1] = (2.0*(sf[0,1:Nx-1] - sf[1,1:Nx-1]))/(h*h)
vt[Ny-1,1:Nx-1] = (2.0*(sf[0,1:Nx-1] - sf[1,1:Nx-1]))/(h*h)
vto = vt.copy()
#print(tstep)
for i in range(1,Nx-1):
for j in range(1,Ny-1):
vt[j,i] = vt[j,i]+dt*(-0.25*((sf[j+1,i]-sf[j-1,i])*(vto[j,i+1]-vto[j,i-1]) - (sf[j,i+1]-sf[j,i-1])*(vto[j+1,i]-vto[j-1,i]))/(h*h) \
+Visc*(vto[j,i+1]+vto[j,i-1]+vto[j+1,i]+vto[j-1,i]-4.0*vto[j,i])/(h*h) )
vt *= shape
t = t + dt
\end{verbatim}
plt.subplot(121) plt.contour(x,y,vt,40) plt.subplot(122)
plt.contour(x,y,sf) plt.show()
\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_
\section*{Results}
{\label{327120}}
This code is meant to eventually reach a steady state but the
requirements for this to happen will be held in the Reynolds number.
That is, the flow into the box must not be too high as compared to the
viscosity of the fluid inside the box. In the case that this is not true
we will not get enough dissipation in the vorticity function and our
solution will eventually diverge. To solve this one must only turn up
the viscosity. However we must also note that the time scale at which we
advance to the next solution is relative to the size of the box. If we
are going to increase the scale at which we are looking to determine a
solution we must also decrease the time increment, failure to do so will
also result in divergence. This code will determine the stream function
and the vorticity function but the often more useful result of the
velocity vectors can be easily generated from equations 5. We can see
from the next figure an appropriate choice of time step, size,
viscosity, and inflow that allows the program to come close to its
steady state. These figures show flow lines of constant velocity.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Screenshot-from-Vorticity-Animation-mp4/Screenshot-from-Vorticity-Animation-mp4}
\caption{{Viscosity = .1,~ Inflow=.1, grid: (20x20), time step = .0005
{\label{672880}}%
}}
\end{center}
\end{figure}
The following shows what happens just before instability given a poor
choice of parameters.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Screenshot-from-Vorticity-Animation2-mp4/Screenshot-from-Vorticity-Animation2-mp4}
\caption{{Viscosity=.1, Inflow=.1, grid: (60x60), time step =.00005%
}}
\end{center}
\end{figure}
One issue with this specific numerical method is the false creation of
curl at the corners do to the inability to generate our stream functions
in these locations. Part of the reason for the instability as seen in
the figure above. This also would affect our ability to determine the
true nature of the flow around the object in question but this can be
resolved some amount by pushing the boundaries far away as compared to
the size of the object in question.
\selectlanguage{english}
\FloatBarrier
\end{document}