Vortexy_doc.tex

\documentclass[12pt]{article}
\usepackage[dvips]{graphicx}
\usepackage{color}
\usepackage[finnish, english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{wrapfig}
\usepackage{caption}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{fancyhdr}
\usepackage{titling}
\usepackage[top=52pt, bottom=2cm, left=2cm, right=2cm]{geometry}
\usepackage{float}
\usepackage{hyperref}
\usepackage{authblk}
\usepackage{bm}
\usepackage{comment}
\usepackage{mathtools}
\usepackage[linesnumbered,ruled]{algorithm2e}

\usepackage{rakmacros/rakmacros}

\pagestyle{fancy}

\title {
  Vortexy fluid simulator
}

\date{\today}

\def \firstauth{
  Roni Koitermaa
}

\def \email{roninkoi@iki.fi}

\author[1] {
  \firstauth\thanks{\href{mailto: \email}{\email}}
}

\fancyhf{}
\setlength{\headheight}{15pt}
\lhead{\thetitle}
\rhead{\firstauth}
\cfoot{\thepage}

\renewcommand\maketitlehooka{\vspace{0.2\textheight}}
% \renewcommand\maketitlehookd{\vfill}

% fluids
\newcommand{\Ren}{\x{Re}}
\newcommand{\fnb}{F_{\x{nb}}}
\newcommand{\sumnb}{\sum_{f \in \fnb}}

\begin{document}

\setlength{\belowcaptionskip}{10pt}

\selectlanguage{english}

\normalsize

\begin{titlingpage}
  \maketitle

  % \begin{abstract}
  % \end{abstract}
  \begin{center}
    Software documentation
  \end{center}
\end{titlingpage}

\newpage

\tableofcontents

\newpage

\section{Introduction}

{\bf Vortexy} is a computational fluid dynamics (CFD) simulation code. It is written in C and uses the finite volume method with the SIMPLE algorithm to calculate flow of incompressible fluids, namely liquids.

The simulator is based on irregular tetrahedral meshes. These meshes can be computed from surfaces using the program Tetgen. The simulator takes a configuration file as input that contains paths to the simulation mesh and boundary conditions in addition to other settings. The state of the system is periodically written to an output file specified in the config. Included is also a renderer that uses OpenGL to visualize results.

\section{Background}

\subsection{Navier-Stokes equations}

The \textit{Navier-Stokes equations} form the basis for all of fluid dynamics. The momentum equation is typically written as \cite[p.~59]{tri}
\begin{equation}
  \frac{\partial \vb u}{\partial t} + (\vb u \cdot \nabla) \vb u = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \vb u + g,
\end{equation}
where $\vb u = (u, v, w)$ is velocity in $\un{m} \ \up{s}{-1}$, t time in $\un{s}$, $\rho$ density in $\un{kg} \ \up{m}{-3}$, $p$ pressure in $\un{Pa}$, $\nu=\frac{\mu}{\rho}$ kinematic viscosity in $\up{m}{2} \ \up{s}{-1}$, $g$ gravity in $\un{m} \ \up{s}{-2}$.

The continuity equation must be satisfied for incompressible fluids that have no sinks of sources
\begin{equation}
  \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vb u) = 0 \rightarrow \nabla \cdot \vb u = 0,
\end{equation}
where $\vb u = (u, v, w)$ is velocity in $\un{m} \ \up{s}{-1}$.

\subsection{Turbulence}

A simple way of predicting onset of turbulence is the \textit{Reynolds number} \cite{tri}:
$$
\text{Re} = \frac{\mu u L}{\rho} = \frac{u L}{\nu}
$$
Turbulence models in simulations include RANS (Reynolds Averaged) e.g. $k$-$\varepsilon$, LES (Large Eddy) and DNS (Direct).

\section{Methods}

\trifig{plots/streamex.pdf}{Stream plot example.}{plots/pipe_vm.pdf}{Poiseuille, velocity \mbox{$x$-component}, cross section.}{plots/pipe_vdp.png}{Pipe flow velocity magnitude, 3D view.}

\subsection{Finite volume method}

The \textit{finite volume method} (FVM) is based on a simulation mesh with volume elements. This enables evaluation of partial differential equations (PDEs) prevalent in physics. The divergence theorem allows us to convert volume integrals to surface integrals
$$
\int_V \nabla \cdot \vb F \ \D V = \oint_S \vb F \cdot \D \vb S,
$$
so volume terms can be computed from fluxes at element faces.

Gradients need to be calculated for both fields. This can be achieved using the Green-Gauss gradient, which can be calculated as
$$
\nabla \phi = \frac{1}{V} \int_S \phi \ \D \vb S = \frac{1}{V} \sum_{f \in \fnb} \phi_f \ \vb S_f.
$$
Face values can be computed from the volume element values by linear interpolation. Second-order gradients are computed by iterative correction.

\subsection{SIMPLE algorithm}

The \textit{Semi-Implicit Method for Pressure-Linked Equations algorithm} solves the Navier-Stokes equations using an iterative procedure \cite{pat}. First the momentum equation is solved, producing a momentum-conserving field $u^*$. However, this resultant field is not necessarily divergence free, meaning it does not satisfy the continuity equation. The SIMPLE algorithm solves this by calculating a correction $u'$ to the intermediate field by solving the \textit{pressure equation}. The pressure equation is derived from the continuity equation.

The SIMPLE algorithm can be summarized as follows \cite{mou, pat}:
\begin{enumerate}
\item Set boundary conditions, set $u$ and $p$
\item Compute gradients $\nabla u$ and $\nabla p$
\item Compute mass fluxes $j_m$, flow rate $\dot m = j_m \cdot A$
\item Solve \textit{momentum equation} using velocity guess $u^0$
\item Solve \textit{pressure correction equation} to get $p'$
\item Correct pressure $p = p^* + p'$ and velocity $v = v^* + v'$
\item Increment time $t^{(n+1)} = t^{(n)} + \Delta t$
\item Jump to step 2
\end{enumerate}
Rhie-Chow interpolation is performed as \cite{mou}
$$
\vb u_f = \overline{\vb u_f} - \underbrace{\overline{C_f} \lb \nabla p_f - \overline{\nabla p_f} \rb}_\text{Pressure correction} + \underbrace{(1 - \alpha_u) \lb \vb u_f^{(i)} - \overline{\vb u_f^{(i)}} \rb}_\text{Relaxation correction} + \underbrace{\frac{\overline{C_f}}{\overline{C_f^{(t)}}} \lb \vb u_f^{(t)} - \overline{\vb u_f^{(t)}} \rb}_\text{Transient correction}.
$$

\subsection{Discretization}

The momentum equation is written in a form conducive for discretization \cite{mou}:
$$
  \underbrace{\frac{\partial \vb u}{\partial t}}_\text{Transient} + \underbrace{(\vb u \cdot \nabla) \vb u}_\text{Convection} = \underbrace{\nu \nabla \cdot (\nabla \vb u)}_\text{Diffusion} + \underbrace{\nu \nabla \cdot (\nabla \vb u)^\T + \frac{1}{\rho} \nabla p + \vb g}_\text{Source}.
$$
The continuity equation for liquids is
$$
\nabla \cdot \vb u = 0.
$$
This saddle-point problem can be represented in matrix form as \cite{mou}
$$
\x A \vb u = \begin{pmatrix}
  \x F & \x B^\T \\
  \x B & 0
\end{pmatrix}
\begin{pmatrix}
  \vb u \\
  p
\end{pmatrix} =
\begin{pmatrix}
  \vb f_b \\
  0
\end{pmatrix}.
$$
Discretization for the one-dimensional momentum equation \cite{mou}:
\begin{equation}
  a_V u_V + \sum_{f \in \fnb} a_f u_f = b_V,
\end{equation}
where $a_V$ refers to a volume element and $a_f$ to a face element coefficient. This can be represented in matrix form as
$$
\begin{pmatrix}
  a_{11} & \dots & 0 & \dots & 0 \\
  \vdots & \ddots & &  & \vdots \\
  a_f & a_f & a_V & a_f & a_f \\
  \vdots & &  & \ddots & \vdots \\
  0 & \dots & 0 & \dots & a_{nn}
\end{pmatrix}
\begin{pmatrix}
  u_1 \\
  \vdots \\
  u_n
\end{pmatrix} =
\begin{pmatrix}
  b_{1} \\
  \vdots \\
  b_n
\end{pmatrix},
$$
where $u_1 \dots u_n$ represents volume element velocities. In the pressure calculation these are referred to as $u^*$, i.e. momentum conserving. The coefficients $a$ and $b$ are calculated for each volume element and assembled into a matrix. Components $x$, $y$ and $z$ are calculated one after another.
Face fluxes:
\begin{align}
  \phi_f = \max(\dot m_f, 0) + \mu \frac{E_f}{d_{vf}}, \\
  \Phi_f = -\max(-\dot m_f, 0) - \mu \frac{E_f}{d_{vf}}, \\
  \psi_f = - \mu (\nabla \vb u_f) \cdot \vb T_f + \dot m_f (\vb u_f^{\text{hr}} - \vb u_f^{\text{uw}}),
\end{align}
where $\vb S_f = A_f \hat n_f = \vb E_f + \vb T_f$ and $\dot m_f = \rho \vb u_f \cdot \vb S_f$.
Volume fluxes:
\begin{align}
  \phi_V = \frac{\rho V}{\Delta t}, \\
  \psi_V = \frac{\rho V}{\Delta t} \vb u_V - \vb f_b V.
\end{align}
Coefficients:
\begin{align}
  a_V = \phi_V + \sum_{f \in \fnb} \phi_f = \frac{\rho V}{\Delta t} + \sum_{f \in \fnb} \lb \max(\dot m_f, 0) + \mu \frac{E_f}{d_{vf}} \rb, \\
  a_f = \Phi_f = -\max(-\dot m_f, 0) - \mu \frac{E_f}{d_{vf}},
\end{align}
\begin{equation}
  \vb b_V = -\psi_v - \sum_{f \in \fnb} \psi_f + \sum_{f \in \fnb} \lsb \mu (\nabla \vb u_f)^\T \cdot \vb S_f \rsb - V \nabla p_v,
\end{equation}
$$
= - \frac{\rho V}{\Delta t} \vb u_V - \vb f_b V - \sum_{f \in \fnb} \lsb - \mu (\nabla \vb u_f) \cdot \vb T_f + \dot m_f \lb \vb u_f - \vb u_f^{\text{uw}}\rb \rsb + \sum_{f \in \fnb} \lb \mu (\nabla \vb u_f)^\T \cdot \vb S_f \rb - V \nabla p_v.
$$
For interpolation we define:
$$
\vb C_V = \frac{V}{a_V}.
$$
Similarly, the pressure correction equation can be discretized as \cite{mou}
\begin{equation}
  a_V p_v' + \sum_{f \in \fnb} a_f p_f' = b_V
\end{equation}
with coefficients
\begin{align}
  a_f = -\rho \mathcal{D}_f, \\
  a_V = \sum_{f \in \fnb} \rho \mathcal{D}_f,\\
  b_V = - \sum_{f \in \fnb} \dot m_f^* = - \sum_{f \in \fnb} \rho \vb u_f^* \cdot \vb S_f,
\end{align}
where $\mathcal{D}_f$ is calculated e.g. using the minimum correction approach:
$$
\mathcal{D}_f = \frac{d_{vf}^x S_f^x \overline{C_f^x} + d_{vf}^y S_f^y \overline{C_f^y} + d_{vf}^z S_f^z \overline{C_f^z}}{(d_{vf}^x)^2 + (d_{vf}^y)^2 + (d_{vf}^z)^2},
$$
where the overlined values have been interpolated from volume elements to faces.

All of these are $n \times n$ diagonally dominant matrices which can be solved using a numerical matrix solver, such as Gauss-Seidel. The off-diagonal coefficients represent neighbouring elements, so they will be mostly zero.

\subsection{Gauss-Seidel method}

The Gauss-Seidel method is an iterative procedure to solve a linear system of equations. It belongs to a class of successive over-relaxation methods (SOR). The linear system is represented as a diagonally dominant square matrix $\x A$ and given a vector $\vb b$ the solution vector $\vb x$ to the equation
$$
\x A \vb x = \vb b
$$
is found by an iterative algorithm. GS has better efficiency for large systems than other methods such as matrix inversion.

The iteration for the GS method can be defined as \cite[p.~510]{golub}

\begin{equation}
  x_i^{(k + 1)} = \frac{b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k + 1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)}}{a_{ii}}.
\end{equation}

\subsection{Boundary conditions}

\begin{table}[H]
  \centering
  \begin{tabular}{| l | l |}
    \hline
    Type & Specified quantity \\
    \hline
    No-slip wall & $u$ \\
    \hline
    Slip wall & $u$ \\
    \hline
    Inlet & $u$ or $p$ \\
    \hline
    Outlet & $\dot m$ or $p$\\
    \hline
  \end{tabular}
  \caption{Four main types of boundary conditions}
\end{table}
The boundary conditions at each face define the matrix coefficients $a_f$. For example, a slip wall would have a coefficient $a_f = 0$ at the corresponding column in the matrix $\x A$.

\subsection{Simulation mesh}

The simulation mesh is an irregular, non-orthogonal Delaunay tetrahedralization that is generated from a closed boundary surface. A triangular surface can be used to generate tetrahedra that fill up the bounded volume with constraints such as angle, volume and edge length. Each volume element has four faces and up to four neighbouring volume elements.

\section{Software}

Vortexy consists of two parts: the simulator and the renderer. The renderer can be used to view output of the simulator in real time and from files. It can also be used to set boundary conditions, but this functionality is limited. The simulator is implemented in the folder \verb|src/phys|. This contains code for finite volume mesh computation (\verb|volume.c|), gradient computation using the Green-Gauss theorem (\verb|gradient.c|) and interpolation (\verb|interpolate.c|). Finally, there is the solver itself that assembles the matrix equations (\verb|solver.c|). The SIMPLE algorithm is implemented in file \verb|sys.c| function \verb|p_sysTick()|. This is the main loop that solves the matrix equations at every time step. The equations are solved iteratively until the desired residual is reached.

\subsection{Rendering}

OpenGL 3.1, GLSL 150.

\subsection{Compilation}

The software is compiled using CMake. The \verb|CMakeLists.txt| file is found in the root directory. The simulator can be configured either in graphical or non-graphical mode. Rendering can be disabled by passing the \verb|-DRENDER=OFF| flag to CMake. The non-graphical version is linked entirely statically, so static C development libraries are needed. The graphical version has its \verb|RPATH| set to \verb|lib/|, so required dynamically linked libraries can be placed there.

The root directory contains shell scripts that can be used for compilation. To build the simulator without the renderer use the script \verb|build.sh 0|, which sets \verb|RENDER_ENABLED| to 0 in \verb|sim.h|.\\
Building using default options (dependencies: libc, libm, OpenGL, GLEW, GLFW):
\begin{verbatim}
  cmake .
  make
\end{verbatim}

\subsection{Configuration}

\subsubsection{Simulation config}

\begin{verbatim}
# simulator options
render 0 # use renderer
simulating 1 # start simulation right away
autoquit 1 # quit when not simulating
divhalt 0 # halt simulation upon divergence

fluid data/lid16/lid100.cfg # fluid config file
endt 0.15 # end time / sim duration
maxit 10000 # max solver iterations
epsilon 1.0e-6 # solver target
dtmaxit 15 # max time-step iterations
residual 5.0e0 # solution residual
transient 0 # include Chie-Chow transient term
convsch 0 # convection scheme, 0 = upwind, 1 = second-order upwind
gradit 2 # gradient iterations, if 0 -> first-order
relaxm 1.0 # relax mass flow

# IO options
file data/lid16/out100.dat # input/output file
printitn 1 # print time-step iteration number
outputting 1 # writing to output file?
outputt 0 # outputting time data
outputf 1 # frequency of output (1/cycle)
inputf 20 # frequency of input (1/cycle)
inputram 0 # bytes to load into memory, 0 for direct file access

# rendering options
# f  v  t  l  w  m
# 64 32 16 8  4  2
rmode 10 # rendering mode (see above)
rz -4.0 # camera z
rs 0.3 # render scale
bgcol 1.0 1.0 1.0 # background color

vs 1.0 # velocity scale (m/s)
ps 100000.0 # pressure scale (Pa)

end # optional
\end{verbatim}

\subsubsection{Fluid config}

\begin{verbatim}
mu 2.0e-1 # viscosity
rho 1.0e-1 # density

dt 0.001 # time step

vr 0.7 # velocity urf
pr 0.3 # pressure under-relaxation factor

mesh data/lid16/plane_16x16_s2.1.mesh

f 0.0 0.0 0.0

bp 1.0e3 # base pressure, p = bp + ip

ebc 1 # boundary condition for all edges

# bc <f id> <bc> = boundary condition
# 0 = open, 1 = no-slip wall, 2 = slip wall
# 3 = velocity inlet, 4 = pressure inlet
# 5 = mass flow outlet, 6 = pressure outlet
# iv <f id> <v> = initial velocity, ip <f id> <p> = initial pressure
# cv <f id> <v> = constant velocity, cp <f id> <p> = constant pressure

bc 244 3 # velocity inlet for face id 244
cv 244 10.0e2
iv 244 10.0e2

end # optional
\end{verbatim}

\subsection{Output data}

\begin{verbatim}
s <sim tick>
o <object id> t <time in s> f <face id> x <face centroid x in m> <y> <z> 
v <face velocity x in m/s> <y> <z> p <face pressure in Pa>
e
\end{verbatim}

\addcontentsline{toc}{section}{References}
\begin{thebibliography}{9}
\bibitem{tri} D. J. Tritton. (1977). Physical Fluid Dynamics. Oxford University Press.  
\bibitem{mou} F. Moukalled, L. Mangani \& M. Darwish. (2016). The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM and Matlab. Fluid Mechanics and Its Applications Volume 113. Springer International Publishing, Cham. \url{https://doi.org/10.1007/978-3-319-16874-6}
\bibitem{golub} G. H. Golub \& C. F. Van Loan. (1996). Matrix Computations. The John Hopkins University Press.
\bibitem{pat} S. V. Patankar. (1980). Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation.
\end{thebibliography}

% See also:
% https://www.cfd-online.com/Wiki/SIMPLE_algorithm
% https://www.cfd-online.com/Wiki/Gradient_computation

\end{document}

% Local Variables:
% coding: utf-8-unix
% TeX-engine: luatex
% End: