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. It turned out to be advection-diffusion that cause the unstable behaver: $\textbf{u} \cdot \nabla \textbf{u} \gg \mu \Delta \textbf{u}$. Similar problem will probably happen for the extension when: $\textbf{u} \cdot \nabla \phi \gg \mu \Delta \phi$.  
%nevn at du søtte på advection diffusion problem. explici/simi implicit/ implicit
%Nevn også at enten kan dt bli liten eller så kan petro -galerkin metoden.

\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. 


%\verbatiminput{cinnforing.py}
%\begin{equation}
%erf^{-1}\left(\frac{T\left(x_i,t\right)-T_0}{\Delta T}\right) = D^{-1/2}\cdot \gamma \hspace{0.5cm} \gamma = %\frac{x_i}{\sqrt{4t}}
%\label{eq:Dligning}
%\end{equation}

%\ref{eq:Dligning}

%\begin{center}
%174 333 371 902 042 752
%\end{center}
\end{document}