ProjectReport.tex

\documentclass[12pt,a4paper,english]{article}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc} %for å bruke æøå
\usepackage{graphicx}%for å inkludere grafikk
\usepackage{verbatim} %for å inkludere filer med tegn LaTeX ikke liker
\usepackage{hyperref}%To include web address 
\usepackage{babel,url}%Denne kan ha samme namespace som den over
\usepackage{amsmath}
\usepackage{listings}

%\include{pythonlisting}%matlablisting

\bibliographystyle{plain}


\title{INF5620 Project\\ Extended Navier-Stokes solver for Platelet Aggregation}
\author{Jarle Sogn\\ Institutt for matematikk\\
Universitetet i Oslo\\ \url{jarlesog@math.no}}

\begin{document}
\maketitle

\begin{abstract}
%hensikt
%Gjennomføring
%viktigste konklusjonene i oppgave
The aim of the project is to solve an extended set of Navier-Stokes equations by the finite element method and \textit{fenics}. The model I'm using is a simplification of model found in Aaron L. Fogelson article: \href{http://www.jstor.org/stable/2102193}{\textbf{Continuum models of platelet aggregation: Formulation and mechanical properties}}. I will focus mainly on equation 1 - 3 in this article. My solver is an extantion of \href{http://fenicsproject.org/documentation/dolfin/1.0.0/python/demo/pde/navier-stokes/python/documentation.html}{\textit{fenics}} Navier-Stoke demo.
%This project might be relevant for my master thesis. 

\end{abstract}

\section*{The Model}
%\subsection*{a)}

I focus on equation 1 - 3 from the article, witch are:

\begin{equation}
\rho \left(\textbf{u}_t + \textbf{u} \cdot \nabla \textbf{u} \right) = -\nabla p + \mu \Delta \textbf{u} + \textbf{f}^g 
\label{eq:MainNS}
\end{equation}
 
\begin{equation}
\nabla \cdot \textbf{u} = 0
\label{eq:MainIncompressible}
\end{equation}

\begin{equation}
\frac{\partial \phi}{\partial t} + \textbf{u} \cdot \nabla \phi = C\Delta \phi - R\left( c \right)\phi 
\label{eq:MainPhi}
\end{equation}

The first two equation are Navier-Stokes equations for incompressible fluid. $\textbf{u}\left(t, \textbf{x} \right)$ is the fluid velocity field, $p\left(t, \textbf{x} \right)$ is the pressure, $\rho$ is the fluid's mass density(I'll sett this equal to 1, for simplicity.) and $\mu$ is the fluids viscosity(assumed to be constant). $\epsilon^{-3} \phi$ is the concentration of non-activated platelet($\epsilon^{-3}$ is a scaling constant). $C$ is a diffusion constant\footnote{Fogelson article uses $D_n$ and $\phi_n$ for $C$ and $\phi$.}. The term $R\left( c \right)\phi$ is the rate in which the non-activated platelets are converted to active platelets. This depends on the ADP concentration $c$, which vary. Since I'm not using the equation for $c$, I will have to manufacture a suitable function $c\left( \textbf{x}, t\right)$\footnote{To understand the equations and symbols better, I advise you to take a look at the first few first pages of Fogelson article}.     	


\subsection*{Solving Navier-Stokes equation}
Now I'm in detail going to talk about one way of solving the Navier-Stokes equations(\ref{eq:MainNS} and \ref{eq:MainIncompressible}) with \textit{fenics} and the finite element method. I will use Chorin's projection method. The idea is first to compute a tentative velocity($\textbf{u}^\star$) by ignoring the pressure in equation \ref{eq:MainNS} and then project the velocity onto the space of divergence free vector fields. Equation \ref{eq:MainNS} then becomes:

$$
\frac{\partial}{\partial t} \textbf{u}^\star + \textbf{u} \cdot \nabla \textbf{u} - \mu \Delta \textbf{u}^\star = \textbf{f}^g
$$

Let $V = H^{1}_{0}\left( \Omega \right) $ be a Sobolev space and $\Omega$ the domain. I want to write the equation above as a weak formulation. I multiply with the function $\textbf{v}\left(\textbf{x} \right)$ and integrate the space.

$$
\int_\Omega \frac{\partial \textbf{u}^\star}{\partial t} \cdot \textbf{v} + \left( \textbf{u} \cdot \nabla \textbf{u}\right) \cdot \textbf{v} -\mu \Delta \textbf{u}^\star \cdot \textbf{v} d\textbf{x} = \int_\Omega \textbf{f}^g \cdot \textbf{v} d\textbf{x} \quad \forall \textbf{v} \in V
$$
Now if I integrate by parts the lest term in the left hand side(or use Green's formula), I obtain:
$$
\int_\Omega \frac{\partial \textbf{u}^\star}{\partial t} \cdot \textbf{v} + \left( \textbf{u} \cdot \nabla \textbf{u}\right) \cdot \textbf{v} +\mu \nabla \textbf{u}^\star \cdot \nabla \textbf{v} d\textbf{x} = \int_\Omega \textbf{f}^g \cdot \textbf{v} d\textbf{x} \quad \forall \textbf{v} \in V
$$
Finely I discretesize the time($n$) and space($h$), $V_h\subset V$. Then $\textbf{u}$ becomes $\textbf{u}^{n}_{h}$. I define $\left\langle \cdot , \cdot \right\rangle$ as the $L^2\left( \Omega \right)$ norm. The above equation can be written as:

\begin{equation}
\left\langle D^{n}_{t} \textbf{u}^{\star}_{h} , \textbf{v} \right\rangle + \left\langle \textbf{u}^{n-1}_{h} \cdot \nabla \textbf{u}^{n-1}_{h} , \textbf{v} \right\rangle +  \left\langle \mu \nabla \textbf{u}^{\star}_{h} , \nabla \textbf{v} \right\rangle = \left\langle \textbf{f}^{g,n} , \textbf{v} \right\rangle \quad \forall \textbf{v} \in V_h
\label{eq:discreteNS}
\end{equation}
  
The projection give us these two equation:


\begin{equation}
\frac{\textbf{u}^{n}_{h} - \textbf{u}^{\star}_{h}}{k_n} + \nabla p^{n}_{h} = 0
\label{eq:projection1}
\end{equation}

\begin{equation}
\nabla \cdot \textbf{u}^{n}_{h} = 0
\label{eq:projection2}
\end{equation}

where $k_h$ is the size of the local time step. I now multiply \ref{eq:projection1} with $\nabla q$ where $q \in Q_h \subset L^2\left(\Omega \right)$ and integrate. 

$$
\frac{1}{k_n} \int_{\Omega} \textbf{u}^{n}_{h} \cdot \nabla q - \textbf{u}^{\star}_{h} \cdot \nabla q + k_n \nabla p^{n}_{h} \cdot \nabla q dx = 0 \quad \forall q \in Q_h
$$

I integrate by parts the first two terms. The first term becomes zero from equation \ref{eq:projection2} and $\textbf{u}\mid_{d\Omega}$ = 0.

$$
\int_{\Omega} \left( \nabla \cdot \textbf{u}^{\star}_{h} \right)q/k_n + \left( \nabla p^{n}_{h} \cdot \nabla q\right)  dx = 0 \quad \forall q \in Q_h
$$ 
This can be written as:
\begin{equation}
\left\langle \nabla p^{n}_{h} , \nabla q\right\rangle  = - \left\langle \nabla \cdot \textbf{u}^{\star}_{h} , q \right\rangle /k_n \quad \forall q \in Q_h
\label{eq:discreteProjection1}
\end{equation}

Finally I multiply equation \ref{eq:projection1} with $\textbf{v} \in V_h$ and integrate. I obtain:

\begin{equation}
\left\langle  \textbf{u}^{n}_{h} , v \right\rangle = \left\langle  \textbf{u}^{\star}_{h} , v \right\rangle - k_n\left\langle \nabla p^{n}_{h} , v\right\rangle \quad \forall v \in V_h
\label{eq:discreteProjection2}
\end{equation}

The three equation \ref{eq:discreteNS}, \ref{eq:discreteProjection1} and \ref{eq:discreteProjection2} is the Chorin's projection method scheme for solving Navier-Stokes equation. First solve for $\textbf{u}^{\star}_h$ in equation \ref{eq:discreteNS}, the solve the pressure $p^{n}_{h}$ in equation \ref{eq:discreteProjection1} and finally find the velocity $\textbf{u}^{n}_{h}$ in equation \ref{eq:discreteProjection2}. Repeat this for each time-step.
%Forklar rekke følgen og hva du løser for hver av ligningene


\subsection*{Stability}
When running the demo program I note that the speed is less then $1 mm/s$. My model is suppose to model blood flow throw arteries. The velocities throw arteries is about $1000 mm/s$, hence I need to turn up the pressure to acquire this velocity. \\
After increasing the velocity in the demo program, it became unstable. The velocity started fluxing and became non-physically high. To avoid this problem I can use a semi-implicit scheme instead of a explicit scheme. I simply change equation \ref{eq:discreteNS} to:

$$\left\langle D^{n}_{t} \textbf{u}^{\star}_{h} , \textbf{v} \right\rangle + \left\langle \textbf{u}^{\star}_{h} \cdot \nabla \textbf{u}^{n-1}_{h} , \textbf{v} \right\rangle +  \left\langle \mu \nabla \textbf{u}^{\star}_{h} , \nabla \textbf{v} \right\rangle = \left\langle \textbf{f}^{g,n} , \textbf{v} \right\rangle \quad \forall \textbf{v} \in V_h$$


Another problem I ran into, is the advection-diffusion. This happens when $\textbf{u} \cdot \nabla \phi \gg C \Delta \phi$. One way to fix this is to use a finer mesh. However this is not an option for me, since I'm reading the mesh from a file. This problem can be fixed with streamline-diffusion(Petrov-Galerkin method). I have chosen the easy way and increased $C$ to achieve stability. 

\section*{The extension}
As mention earlier I have extended an already existing program. Making it also solve equation \ref{eq:MainPhi}. In order to do this I have to rewrite \ref{eq:MainPhi} to variation form. I multiply the test function $\psi \in \Psi \subset H^{1}_{0}\left( \Omega \right)$ and integrate.


$$
\int_{\Omega} \left( D_t \phi \right)\psi + \textbf{u}\cdot\left( \nabla \phi \right)\psi - C\left( \Delta\phi\right) \psi dx = - \int_{\Omega} R(\textbf{x}, t) \phi\psi dx \quad \forall \psi \in \Psi
$$


after integrating by part and using backward euler as $D_t$, I obtain:
$$
\int_{\Omega} \left( \frac{\phi^n - \phi^{n-1}}{k_n}\right) \psi + \textbf{u}\cdot\left( \nabla \phi^n \right)\psi + C\left( \nabla\phi^n \cdot \nabla \psi \right) + R\left( \textbf{x}, t^n \right) \phi^n \psi dx = 0 \quad \forall \psi \in \Psi
$$
I end up with the scheme:
\begin{equation}
\left\langle \phi^n -\phi^{n-1}, \psi \right\rangle \frac{1}{k_n} + \left\langle \psi \nabla \phi^n, \textbf{u}^n\right\rangle + C\left\langle \nabla \phi^n, \nabla \psi \right\rangle + \left\langle R\phi^n , \psi \right\rangle = 0 \quad \forall \psi \in \Psi_h
\label{eq:discretephi}
\end{equation}

%Spesifiser problemt og hvilke konstanter, funcsjoner, BC, og IS jeg bruker. 
\subsection*{The specified problem}

I use no-slip boundary for the velocity $\textbf{u}$ and $\phi$. This means I sett there value equal to zero at boundary except at the inflow and outflow. $\phi_{in} = 1$ and $\phi_{out} = 0$. The pressure at the inflow is a sinus function that is zero if sinus have negative value(Pressure can't be negative). This initial value of $p$, $\textbf{u}$ and $\phi$ is zero. I use $R = 0$ and $\textbf{f}^g$, but I did one test run with $R\left( \textbf{x},t \right)  = 3x\left( 1-x \right)y\left(1-y \right)\left(t+1\right)$ to check if it was implemented correctly. $\nu = 0.01$ and $C = 10$, this is not physical realistic. \\
After running the program \textit{NSextended.py} and obtaining the results, I used paraview to create a animation of $\phi$, se \textit{moviePHI.avi}. From this animation we can see how the $\phi$ follow the flow $\textbf{u}$, with in turn is powered by the pressure. 



\section*{Verification}

To verify the Navier-Stoke solve I have to simplify the problem. I copyed the code from \textit{NSextended.py} to a new program \textit{verificationEasy.py}, stripped away the extension $\phi$ and created a unitsquare mesh. I force the pressure to be zero and want to make the solution to be $\textbf{u} = \left( Ae^{at}sin\left(\pi y\right) , 0\right)$. Then I have to chose $\textbf{f}$ such that I obtain this.
$$
\textbf{f} = \textbf{u}_t +\textbf{u} \cdot \nabla \textbf{u} - \Delta \textbf{u}
$$ 
I find $\textbf{f}$ by inserting in for $\textbf{u}$.
$$
\textbf{f} = \left( \left(a + \pi^2\right) Ae^{at}sin\left(\pi y\right), 0\right) 
$$

The goal of the verification is to see if the error converges with the correct rate as I make a finer mesh(or smaller time step). Note that the flow only moves in x direction. So If I want to increase accuracy by a finer mesh, I only have to make y axis finer. The Error should be:
$$E = O\left(\Delta t\right) + O\left(\Delta y^2 \right)$$
I now decrease $\Delta t$ such that the error from the time step can be neglected and than see if $E \approx O\left(\Delta y^2 \right)$. The result can be seen en the table under. 
\begin{center}
\begin{tabular}{|r|l|}
  \hline
  $N_y$ & $E$ \\
  \hline
  	4 & 0.0486353 \\
  \hline
    8 & 0.0123844 \\
  \hline
    16 & 0.003104 \\
  \hline
    32 & 0.000775 \\
  \hline
\end{tabular}
\end{center}

When I divide $E_4$ by $E_8$(and $E_8$ by $E_{16}$ etc) I get roughly the number 4(witch is expected).

\subsection*{Other ways of verification}
My previous verification of Navier-Stake was friarly simple. It only tests in one directions and completely ignores the pressure. I can use $\textbf{u} = \nabla \times sin(xy)sin(t)$ and chose a pressure $p$ and then find a function $\textbf{f}$ that fits. 

\end{document}