\documentclass[11pt,b5paper,DIV=calc,BCOR1.3cm,headings=small,%
               footinclude=false,headsepline]{scrbook}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsfonts,amsthm,amssymb,mathtools}
\usepackage{textcomp}
\usepackage[slantedGreek,sc]{mathpazo}
\usepackage{avant}
\usepackage{url}
\usepackage{fix-cm}
\usepackage{siunitx}
\usepackage{scrpage2}
\usepackage{pgfplots}
\pgfplotsset{compat=1.6}
\usepackage{tikz}
\usepackage[grey]{quotchap}
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage[bf,small,margin=10pt]{caption}
\usepackage[margin=1em]{subcaption}
\usepackage{bm}
\usepackage{csquotes}
\usepackage[backend=biber,
  maxnames=2,
  maxbibnames=9,
  firstinits,
  natbib=true]{biblatex}
\usepackage{booktabs}
\usepackage{microtype}
\usepackage{pdfpages}
\usepackage[pdfborder={0 0 0},hyperfigures]{hyperref}
\usepackage{cleveref}
\usepackage{cancel}
\usepackage{centernot}
\usepackage[perpage,symbol*]{footmisc}

\addbibresource{thesis.bib}
\addbibresource{.reference.bib}
\DeclareNameAlias{default}{last-first}

% Set PDF file information
\hypersetup{
  pdfinfo={
    Title={Calculation of interface curvatures with the level-set method for
      two-phase flow simulations and a second-order diffuse-domain method for
      elliptic problems in complex geometries},
    Author={Karl Yngve Lervåg},
    Subject={PhD Thesis},
  }
}

% TikZ settings
\usetikzlibrary{calc}
\usetikzlibrary{arrows}
\usetikzlibrary{positioning}
\usetikzlibrary{intersections}
\usetikzlibrary{decorations.pathreplacing}
\usetikzlibrary{decorations.pathmorphing}
\usetikzlibrary{shapes.geometric}

% Figure paths
\graphicspath{{figures/}}

% Clever reference package
\crefname{equation}{Equation}{Equations}
\creflabelformat{equation}{#2(#1)#3}
\crefname{pluraleq}{equations}{equations}
\creflabelformat{pluraleq}{#2(#1)#3}
\crefname{section}{Section}{Sections}
\crefname{chapter}{Chapter}{Chapters}
\crefname{figure}{Figure}{Figures}
\crefname{table}{Table}{Tables}
\crefname{subfigure}{Figure}{Figures}

% Don't allow small figures to appear alone on a page
\renewcommand{\textfraction}{0.20}
\renewcommand{\topfraction}{0.80}
\renewcommand{\bottomfraction}{0.80}
\renewcommand{\floatpagefraction}{0.80}

% Fancy differential from Claudio Beccari, TUGboat
% * No need for manual tweak of spaces.
% * Copied from Svend Tollak Munkejord.
\makeatletter
\newcommand*{\dif}{\@ifnextchar^{\DIfF}{\DIfF^{}}}
\def\DIfF^#1{\mathop{\mathrm{\mathstrut d}}\nolimits^{#1}\gobblesp@ce}
\def\gobblesp@ce{\futurelet\diffarg\opsp@ce}
\def\opsp@ce{%
  \let\DiffSpace\!%
  \ifx\diffarg(%
    \let\DiffSpace\relax
  \else
    \ifx\diffarg[%
      \let\DiffSpace\relax
    \else
      \ifx\diffarg\{%
        \let\DiffSpace\relax
      \fi\fi\fi\DiffSpace}
\makeatother

% Define the jmp macro
\newbox\bokstav
\newdimen\hoyde
\newcommand{\bgl}{{\hbox{$\left\lbrack\vbox to 8.5pt{}\right.\nOspace$}}}
\newcommand{\Bgl}{{\hbox{$\left\lbrack\vbox to 11.5pt{}\right.\nOspace$}}}
\newcommand{\bggl}{{\hbox{$\left\lbrack\vbox to 14.5pt{}\right.\nOspace$}}}
\newcommand{\Bggl}{{\hbox{$\left\lbrack\vbox to 17.5pt{}\right.\nOspace$}}}
\newcommand{\bgr}{{\hbox{$\left\rbrack\vbox to 8.5pt{}\right.\nOspace$}}}
\newcommand{\Bgr}{{\hbox{$\left\rbrack\vbox to 11.5pt{}\right.\nOspace$}}}
\newcommand{\bggr}{{\hbox{$\left\rbrack\vbox to 14.5pt{}\right.\nOspace$}}}
\newcommand{\Bggr}{{\hbox{$\left\rbrack\vbox to 17.5pt{}\right.\nOspace$}}}
\newcommand{\nOspace}{\nulldelimiterspace=0pt \mOth}
\newcommand{\mOth}{\mathsurround=0pt}
% ordinary
\newcommand{\Ljmp}{\mathopen{\lbrack\!\lbrack}}
\newcommand{\Rjmp}{\mathclose{\rbrack\!\rbrack}}
% big
\newcommand{\bgLjmp}{\mathopen{\Bgl\!\!\Bgl}}
\newcommand{\bgRjmp}{\mathclose{\Bgr\!\!\Bgr}}
% Big
\newcommand{\BgLjmp}{\mathopen{\Bgl\!\!\Bgl}}
\newcommand{\BgRjmp}{\mathclose{\Bgr\!\!\Bgr}}
% bigg
\newcommand{\bggLjmp}{\mathopen{\bggl\!\!\bggl}}
\newcommand{\bggRjmp}{\mathclose{\bggr\!\!\bggr}}
% Bigg
\newcommand{\BggLjmp}{\mathopen{\Bggl\!\!\Bggl}}
\newcommand{\BggRjmp}{\mathclose{\Bggr\!\!\Bggr}}
\newcommand{\jmp}[1]{%
\setbox\bokstav=\hbox{$ \left. #1\right. $}
\hoyde=\ht\bokstav
\advance\hoyde by \dp\bokstav%
% \showthe\hoyde
\hbox{$
  \ifinner
    \ifdim\hoyde<10pt
      \Ljmp #1 \Rjmp%
    \else
      \ifdim\hoyde <11pt
        \Ljmp #1 \Rjmp%
      \else
        \ifdim\hoyde <14pt
          \bgLjmp #1 \bgRjmp%
        \else
          \ifdim\hoyde <20pt
            \BgLjmp #1 \BgRjmp%
          \else
            \bggLjmp #1 \bggRjmp%
          \fi
        \fi
      \fi
    \fi
  \else
    \ifdim\hoyde<8.5pt
      \Ljmp #1 \Rjmp%
    \else
      \ifdim\hoyde <11.5pt
        \bgLjmp #1 \bgRjmp%
      \else
        \ifdim\hoyde <14.5pt \Ch
          \BgLjmp #1 \BgRjmp%
        \else
          \ifdim\hoyde <17.5pt
            \bggLjmp #1 \bggRjmp%
          \else
            \BggLjmp #1 \BggRjmp%
          \fi
        \fi
      \fi
    \fi
  \fi
$}}

% New commands
\newcommand*{\ifrac}[2]{\ensuremath{#1/#2}}
\newcommand*{\cd}[1]{\ensuremath{\pdt{#1} + \vct u\cdot\grad #1}}
\newcommand*{\dt}[1]{\ensuremath{\frac{\dif #1}{\dif t}}}
\newcommand*{\dti}[1]{\ifrac{\dif #1}{\dif t}}
\newcommand*{\pdt}[1]{\ensuremath{\frac{\partial #1}{\partial t}}}
\newcommand*{\pdti}[1]{\ifrac{\partial #1}{\partial t}}
\newcommand*{\td}[1]{\frac{\mathrm{D} #1}{\mathrm{D} t}}
\newcommand*{\tdi}[1]{\mathrm D #1 /\mathrm Dt}
\newcommand*{\od}[2]{\ensuremath{\frac{\dif#1}{\dif{#2}}}}
\newcommand*{\odi}[2]{\ensuremath{{\dif#1}/{\dif{#2}}}}
\newcommand*{\pd}[2]{\ensuremath{\frac{\partial #1}{\partial{#2}}}}
\newcommand*{\pdi}[2]{\ensuremath{{\partial #1}/{\partial{#2}}}}
\newcommand*{\pdd}[2]{\ensuremath{\frac{\partial^2 #1}{\partial{#2}^2}}}
\newcommand*{\vct}[1]{\ensuremath{\boldsymbol{#1}}}
\newcommand*{\normal}{\ensuremath{\vct n}}
\newcommand*{\del}{\boldsymbol\nabla}
\renewcommand*{\div}{\del\cdot}
\newcommand*{\divs}{\del_{\text s}\cdot}
\newcommand*{\grad}{\del}
\newcommand*{\grads}{\del_{\text s}}
\newcommand*{\lapl}{\boldsymbol\Delta}
\newcommand*{\lapls}{\boldsymbol\Delta_{\text s}}
\newcommand*{\sint}[2]{\ensuremath{\int_{#1}#2\dif #1}}
\newcommand*{\surint}[2]{\ensuremath{\int_{#1}#2\dif S}}
\newcommand*{\chr}{\chi}
\newcommand*{\bigo}[1]{\ensuremath{\mathcal O\left(#1\right)}}
\newcommand*{\smallo}[1]{\ensuremath{o\left(#1\right)}}
\newcommand*{\eint}[1]{\ensuremath{\int_{-\infty}^\infty #1\dif z}}
\newcommand*{\zint}[1]{\ensuremath{\int #1\dif z}}
\newcommand*{\ejmp}[1]{\ensuremath{\left[#1\right]_{-\infty}^\infty}}
\newcommand*{\zlimp}{\ensuremath{\lim_{z\to\infty}}}
\newcommand*{\zlimm}{\ensuremath{\lim_{z\to-\infty}}}
\newcommand*{\zlimpm}{\ensuremath{\lim_{z\to\pm\infty}}}
\newcommand*{\ndot}{\ensuremath{\vct n\cdot}}
\newcommand*{\vext}{\ensuremath{V_{\text{ext}}}}
\newcommand*{\txin}{\ensuremath{\textup{in}\ }}
\newcommand*{\txon}{\ensuremath{\textup{on}\ }}
\newcommand*{\tenstr}{\vct T}
\newcommand*{\tenI}{\vct I}
\newcommand*{\tenD}{\vct D}
\newcommand*{\funsfd}{\vct f_{\text{sfd}}}
\newcommand*{\funssf}{\vct f_{s}}
\newcommand*{\intoepdx}[1]{\int_{\Omega_\epsilon}#1\dif\vct x}
\newcommand*{\intdoepds}[1]{\int_{\partial\Omega_\epsilon}#1\dif s}
\newcommand{\set}[2]{\ensuremath{\left\{#1\,:\,#2\right\}}}
\newcommand*{\e}[1]{\ensuremath{\times 10^{-#1}}}

% Declare math operators
\DeclareMathOperator\trace{tr}
\DeclareMathOperator\sech{sech}

% Define header and footer
\setkomafont{pageheadfoot}{\normalfont\sffamily\bfseries}
\setkomafont{pagenumber}{\normalfont\bfseries}
%\setheadsepline{0.4pt}
\ohead[]{\pagemark}
\ihead[]{\headmark}
\ifoot[]{}
\ofoot[]{}
\cfoot[\pagemark]{}
\pagestyle{scrheadings}

% Remove overfull hbox warnings
\hbadness=10000
\hfuzz=50pt

\begin{document}

\frontmatter

% Fakechapter: Titlepage
\begin{titlepage}
  \setlength{\parindent}{0cm}
  \addtolength{\parskip}{\baselineskip}
  {\large\sffamily Karl Yngve Lervåg}

  \vspace{2cm}%{\vspace{0.5cm}}

  {\sffamily\large\bfseries Calculation of interface curvatures with the
    level-set method for two-phase flow simulations and a second-order
    diffuse-domain method for elliptic problems in complex geometries}

  \vfill

  {\large Doctoral thesis\\
    for the degree of Philosophiae Doctor

    Trondheim, June 2013}

  \vspace{2cm}

  \begin{tabular}{m{0.18\textwidth}m{0.80\textwidth}}
    \includegraphics[width=2.0cm]{ntnublaa}
    & \parbox{\linewidth}{%
      {\Large\textbf{NTNU}}\\
      \textbf{Norwegian University of Science and Technology}\\
      Faculty of Engineering Science and Technology\\
      Department of Energy and Process Engineering}
  \end{tabular}

  \newpage
  \thispagestyle{empty}
  {~}
  \vfill
  {\scriptsize
    \textbf{NTNU}\\
    Norwegian University of Science and Technology

    Thesis for the degree of Philosophiae Doctor

    Faculty of Engineering Science \& Technology\\
    Department of Energy and Process Engineering

    \copyright\ 2013 Karl Yngve Lervåg

    ISBN 978-82-471-4544-9 (printed version)\\
    ISBN 978-82-471-4545-6 (electronic version)\\
    ISSN 1503-8181

    Doctoral theses at NTNU, 2013:214

    Printed by Skipnes kommunikasjon
  }
\end{titlepage}

% Fakechapter: Dedication
\thispagestyle{empty}
\vspace*{7cm}
\begin{center}
  Dedicated to Jon Vegard Lervåg (1979--2011)
\end{center}

\chapter*{Abstract}
\addcontentsline{toc}{chapter}{Abstract}

This thesis considers in the first part the mathematical modelling of
incompressible two-phase flow, in particular the calculation of interface
curvatures and normal vectors with the level-set method.  The main contribution
is the development of two new numerical methods that enable a more robust
calculation of the curvature and normal vectors in areas where the gradient
of the level-set method is discontinuous.

Incompressible two-phase flow is in this thesis modelled by the Navier-Stokes
equations with a singular source term at the interface between the phases.  The
singular source term leads to a set of interface jump conditions.  These jump
conditions are used in the ghost-fluid method to solve two-phase flow in
a sharp manner.  The interface position is captured and evolved in time with
the level-set method.  The Navier-Stokes equations for two-phase flow are
solved with projection methods and discretized by finite differences in space
and Runge-Kutta methods in time.  The advective terms in the governing
equations are discretized by a weighted essentially non-oscillatory scheme.

In the second part, the thesis considers the more general problem of solving
partial-differential equations (PDEs) in complex geometries.  An extension of
a diffuse-domain method is presented, where the accuracy is improved by adding
a correction term.  The extension is derived for elliptic problems with Neumann
and Robin boundary conditions.  One of the advantages of the diffuse-domain
methods is that they allow the use of standard tools and methods because they
are based on solving PDEs reformulated in larger and regular domains.

\chapter*{Preface}
\addcontentsline{toc}{chapter}{Preface}

This thesis is submitted to the Norwegian University of Science and Technology
(NTNU) for partial fulfilment of the requirements for the degree of
philosophiae doctor.  The doctoral work has been performed at the Department of
Energy and Process Engineering, NTNU, Trondheim, with Professor Bernhard Müller
as main supervisor and with Svend Tollak Munkejord, chief scientist at SINTEF
Energy Research, as co-supervisor.  The work was carried out in the period from
September 2010 to June 2013.

The project was financed through the research project ``Enabling low emission
LNG systems'', performed under the Petromaks program and coordinated by SINTEF
Energy Research.  I gratefully acknowledge the support from the project
partners: Statoil, GDF SUEZ, and the Research Council of Norway (contract
number 193062/S60).

I am very thankful to both of my supervisors, Bernhard Müller and Svend Tollak
Munkejord.  In the regular meetings throughout the PhD project  they have given
me very helpful and encouraging comments and feedback on my work.  They have
allowed me freedom to pursue my own ideas, but at the same time they have
ensured that I was on track so that I finished my PhD work on time.

I am also indebted to Professor John Lowengrub for inviting me to stay at the
University of Irvine, California.  My stay at UC Irvine was a very enlightening
and enjoyable experience, both on a personal and a professional level.  I also
want to thank Esteban Meca at Lowengrub's lab for many helpful and inspiring
discussions.

I am very grateful to the Fulbright Foundation, both for the financial support
for the stay in Irvine and for the invaluable aid in the practical matters of
living abroad.  I am particularly grateful to Ann Kerr for organising several
very interesting seminars and events in the Los Angeles area that allowed both
me and my wife to meet a lot of wonderful people.

I would like to thank all my colleagues and friends from the Department of
Energy and Process Engineering, NTNU.  In particular, I would like to thank
Claudio Walker for helpful discussions about the modelling of two-phase flow
with the level-set method.  Further, I would like to thank Åsmund Ervik for
many fruitful meetings and discussions, both with regard to our paper about the
LOLEX method, and about the intricacies and complications with our numerical
code.  Last, but not least, I thank my office mate Halvor Lund.  It has been
a great pleasure to share the office with him these last three years.

I also want to thank Halvor Lund, Frode Bjørdal, and Lars Eivind Lervåg for
proofreading my manuscript.

Finally, I extend my deepest gratitude to my wife for her unconditional love
and support.

\vspace{2em}
\begin{flushright}
  Trondheim, June 2013 \\
  Karl Yngve Lervåg
\end{flushright}

% Fakechapter: Table of contents
\tableofcontents

\mainmatter

\begin{savequote}[8.2cm]
  ``The scientific man does not aim at an immediate result.  He does not expect
  that his advanced ideas will be readily taken up.  His work is like that of
  the planter -- for the future.  His duty is to lay the foundation for those
  who are to come, and point the way.''
  \qauthor{--- Nikola Tesla (1856--1943)}
\end{savequote}
\chapter{Introduction}
This thesis considers the mathematical modelling and numerical computation of
two-phase flows.  It focuses on developing more robust numerical methods to
calculate the curvature and the normal vector of the interface between the two
phases.  In addition it considers a diffuse-domain method for solving partial
differential equations in complex geometries, and it derives an asymptotically
second-order method for elliptic problems.

\section{Background and motivation}
Two-phase flows are particular examples of multiphase flows of gas and liquid
with an interface that separates the two phases.  In the pedantic sense,
two-phase flow is flow of a single fluid that occurs as two different phases,
for example steam and water.  However, it is common to be more general, and in
this thesis we use the term two-phase flow also for immiscible mixtures of
different fluids, such as water and oil.

Two-phase flows are crucial to a large amount of processes, both in nature and
in industry.  Examples range from weather phenomena, such as rain drops falling
through air, to industrial processes, for instance the separation of water from
oil.  In general, multiphase flow phenomena influence any process where liquids
and gases are involved.  In the oil and gas industry, most processes are
two-phase or multiphase flow processes.  Needless to say, the understanding of
these phenomena is fundamental in the development of new or improved processes.

Consider the international trade of liquefied natural gas (LNG), which is
a particular branch of the oil and gas industry that has undergone an
exceptional growth in the last decades.  There is a strong focus both in Norway
and internationally on producing LNG on large floaters
(FLNG)\footnote{FLNG facilities do not yet exist, although a facility is under
  development by Royal Dutch Shell~\cite{Prelude1}.  The construction of this
  facility was started in 2012, and the first drilling is stated to commence in
  2013~\cite{Prelude2}.}.  There are a number of both environmental and
economic advantages of such FLNG facilities.  In particular, FLNG facilities
would remove the need of long pipelines from the gas fields to the shore, there
would be no requirement for compression units to pump the gas to the shore, and
one would not need to construct onshore production facilities.  This would
significantly reduce the environmental footprint, and would help preserve
marine and coastal environments.  Since an FLNG facility can be moved to a new
location when a field has been depleted, it would make it economically viable
to open up new business opportunities to develop offshore gas fields that would
otherwise remain stranded.

Moving the LNG production to an offshore facility presents a demanding set of
challenges.  A particular challenge is that the elements of a conventional LNG
facility need to fit into an area roughly one quarter the original size.  Heat
exchangers are among the main challenges in the design and operation of LNG
plants \cite{Hesselgreaves01}.  Compact and efficient heat exchangers are
needed to obtain an energy efficient plant with low emissions.  More optimized
designs require more accurate tools for design and operation.  Such tools can
only come as a result of an improved physical understanding of the complex
two-phase flows occurring in the heat exchangers.  This can be achieved by more
detailed mathematical modelling, together with dedicated laboratory
observations, cf.\ \cite{Olsen07}.

\section{Overview of methods for two-phase flow simulations}
In the following, a brief overview is given of different methods for modelling
two-phase flows where the interface location is known.  In particular, we focus
on methods that handle the interface evolution.  For more in-depth reviews, see
for instance \cite{Cristini04,Scardovelli99,Sethian03}.

Interface propagation methods comprise a range of methods that are often
categorised as either front tracking or front capturing.  Front-tracking
methods use Lagrangian particles to track the interface
explicitly~\cite{Tryggvason01}, while front-capturing methods use an Eulerian
approach to capture the interface implicitly.  Examples of the latter are
volume-of-fluid methods~\cite{Scardovelli99,Tryggvason11}, phase-field
methods~\cite{Anderson98,Jacqmin99,Lowengrub98}, and level-set
methods~\cite{Osher03}.

The main advantages of front-tracking methods are their inherent accuracy and
that topological changes do not occur without explicit action.  Hence,
unphysical numerical reconnection does not occur.  This means that if front
tracking is used for the simulations of two drops that collide, these drops
will not coalesce due to numerical reconnection or merging.  However, the
handling of topological changes are challenging, in particular in three
dimensions~\cite{Shin02}.  Also, there are some issues of numerical
instabilities, as discussed by \citet{Sethian85} and \citet{Osher88}.

The volume-of-fluid method utilizes a volume-fraction function whose values
represent the characteristic function of one of the fluid domains.  Its values
are zero or one, except in those cells cut by the interface.  A considerable
advantage of the volume-of-fluid method is that it conserves the mass of both
fluids well.  However, the reconstruction of the interface from volume
fractions is not simple, and computation of geometric quantities such as the
interface curvature is not straightforward.  Also, spurious bubbles and drops
may be created, cf.\ \cite{Lafaurie94}.

The phase-field methods treat the interface in a diffuse manner, where the
fluid properties, such as density and viscosity, change rapidly but smoothly
across the interface.  These methods typically solve the coupled
Cahn-Hilliard/Navier-Stokes equations, where the Cahn-Hilliard equation is
based on the free energy of an interface~\cite{Cahn57}.  Through this energy
formulation, one can model more advanced interface physics such as van der
Waals interactions, electrostatic forces, and fluids with varying
miscibilities.  However, the phase-field methods require that one resolves very
small length scales at the interfaces.  This poses a severe restriction on the
applicability of phase-field methods for two-phase flow, where the length
scales of the flow are generally much larger than those of the interface.

In this thesis, we have used the level-set method~\cite{Osher88}, which
implicitly captures the interface as an isocontour of a function defined in the
entire domain.  The main motivation of this choice is both that the level-set
method handles topological changes of evolving interfaces automatically,
cf.~\citet{Sethian03}, and that it is relatively straightforward to implement.

It should be noted that the automatic handling of topological changes is not
based on physical principles.  For instance, when two interfaces approach each
other and their distance becomes less than the spacing of the grid, the
level-set method can no longer resolve both interfaces and so they are merged.
An important consequence is that the level-set method does not model the
physics involved in the coalescence process, and in particular in the smaller
scales of the film-drainage process.  This process involves a wide range of
length scales, varying from the nanometer scale where van der Waals force are
active to the length scales of the external flow.  Some effort has been made to
include the effects of the smaller scales in front-capturing methods,
cf.~\cite{Nobari96,Tryggvason10}.  Recently, \citet{Kwakkel13} presented
a level-set/volume-of-fluid method that is coupled with a film-drainage model
that predicts if and when two colliding droplets will coalesce.  To prevent the
numerical coalescence, each droplet has its own locally defined level-set
function.

The level-set method has been used to model several diverse phenomena,
such as tumour growth~\cite{Macklin06,Macklin05,Macklin08}, propagation of
wildland fire~\cite{Mallet09}, and computer RAM production~\cite{Melicher08}.
For a good and thorough introduction to the level-set method, see
\cite{Osher03}.

A weakness of the level-set method is that it does not conserve the mass of the
two fluids, in particular in areas of low resolution and/or high curvature.
Different approaches have been developed to overcome this disadvantage, for
example the conservative level-set method~\cite{Olsson05,Olsson07}, the
particle level-set method~\cite{Enright02}, or the coupled
level-set/volume-of-fluid method~\cite{Sussman00}.

When we use the level-set method to capture the interface for immiscible and
incompressible two-phase flow simulations, there will be a sharp change in the
density and viscosity across the interface.  In this thesis we use the
ghost-fluid method~\cite{Fedkiw99}, which is a sharp-interface method where the
jumps are included in the spatial discretizations in a sharp manner.  An
alternative method is the continuous surface-force method introduced by
\citet{Brackbill92}, where the density and the viscosity are smeared out across
the interface through a smoothed Heaviside function.

\clearpage
\section{Goal and contribution of the present thesis}
The main goal of the PhD project has been to develop fundamental knowledge of
two-phase flow phenomenon that are relevant for compact heat exchangers.

In order to do detailed theoretical studies of phenomena that are relevant for
compact heat exchangers, we need to consider two-phase flows with mass and heat
transfer in confined and complex geometries.  One such relevant phenomenon is
the drop-film collision process \cite{Zhao09}.  Even when restricted to
isothermal and immiscible flows, this simple phenomenon remains a challenge.
When the level-set method is used to capture the interface, one must be
particularly careful about how one calculates the interface curvature and
normal vector.  For instance, when two drops are in near collision there is
a kink region in the level-set function between the drops, where the gradient
is discontinuous, see \cref{fig:curvature-problem}.  This discontinuity may
lead to large errors in the curvature and the normal vector if it is not taken
into account in the discretization stencils.
\begin{figure}[tbp]
 \centering
 \begin{subfigure}[t]{0.47\textwidth}
   \centering
   \begin{tikzpicture}
     [ scale=0.6,
       drop/.style={shading=ball, ball color=black!10}, ]

     % Drops and a slice
     \shadedraw[drop] (-1.9,0) circle (1.5);
     \shadedraw[drop] ( 1.9,0) circle (1.5);
     \draw[thick] (-3.8,0.0) node[left] {$\phi(x)$} -- (3.8,0.0);
     \draw[thick, dotted] (0,1.6) -- (0,-1.6);
   \end{tikzpicture}
   \caption{Drops in near contact}
 \end{subfigure}
 \begin{subfigure}[t]{0.47\textwidth}
   \centering
   \begin{tikzpicture}
     [ scale=0.6,
       >=stealth ]

     % Axes and the level-set function
     \draw[<-] (-3.8,2.5) node[above] {$\varphi(x)$} -- (-3.8,-0.5);
     \draw[->] (-3.8,-0.5) -- ( 3.8,-0.5) node[right] {$x$};
     \draw[thin,black] (-3.8,1.5) node[left] {$0$} -- (3.8,1.5);
     \coordinate (a) at (-3.8,1.9);
     \coordinate (b) at (-1.9,0);
     \coordinate (c) at ( 0.0,1.9);
     \coordinate (d) at ( 1.9,0);
     \coordinate (e) at ( 3.8,1.9);
     \draw[thick] (a) -- (b) -- (c) -- (d) -- (e);

     % Marking the kinks
     \fill ( 0.0,1.85) circle (2.5pt);
     \fill (-1.9,0.05) circle (2.5pt);
     \fill ( 1.9,0.05) circle (2.5pt);
   \end{tikzpicture}
   \caption{Slice of the level-set function}
 \end{subfigure}
 \caption{(a) Two drops in near contact.  The dotted line marks a region
   where the derivative of the level-set function, $\varphi(x)$, is not
   defined.  (b) A one-dimensional slice of the level-set function.  The dots
   mark points where the derivative of $\varphi(x)$ is discontinuous.}
 \label{fig:curvature-problem}
\end{figure}

As a step towards computing two-phase flow simulations in confined geometries,
the thesis has considered the diffuse-domain method~\cite{Li09}.  This is
a method where partial-differential equations in complex domains are extended
into larger, regular domains with the use of diffuse approximations of the
physical boundaries.  The approximations converge asymptotically to the
original problem when the width of the diffuse boundary is reduced.  With this
method one can use standard numerical methods to solve equations that
incorporate complex boundaries.

The main contributions of the present thesis are two new methods to calculate
the curvature and normal vector in a robust manner with the level-set method.
These methods are shown to yield more accurate calculations of drop-film and
drop-drop collisions.  The methods are compared with standard methods and with
other methods from the literature.

In addition, the thesis presents an extension of a diffuse-domain method by
a high-order correction term for the solution of elliptic problems in complex
geometries.  New analysis provided in the thesis improves the understanding of
the diffuse-domain method, and the derived method is shown to be more accurate
than the existing diffuse-domain method.

\section{Outline of the thesis}
The thesis is organised as follows: \Cref{chap:two-phase-flow} gives a brief
overview of the derivation of the governing equations for two-phase flow.  It
includes a consideration of the fluid-fluid interface conditions and the
derivation of a simplified jump tensor for the viscous term at the interface.

\Cref{chap:numerical-methods} gives a detailed description of the numerical
methods that are used to solve the two-phase flow equations.  In particular, it
describes the level-set method, which is used to capture the interfaces,
projection methods that are used to solve the Navier-Stokes equations, and the
spatial and temporal discretization schemes.  At the end of the chapter, an
overview of the novel discretization methods for the curvature and the normal
vector is given.

\Cref{chap:diffuse-domain} gives a short introduction to the diffuse-domain
method for an elliptic problem with Neumann boundary conditions.  It introduces
the high-order correction term, and shows that the new method converges
asymptotically with second-order to the original problem.  The chapter presents
new analysis that shows that the correction term is not necessary for
second-order convergence.

In \cref{chap:contributions} the main results of the contributed papers A--E
are summarized, and the author's contributions are highlighted.  Finally,
\cref{chap:conclusions} gives concluding remarks and provides an outlook for
future work.

Full-text versions of the research papers A--E are provided in the Appendices
at the end of the thesis.

% Fakesection: Quote
\begin{savequote}[8.4cm]
  ``But it is just this characteristic of simplicity in the laws of nature
  hitherto discovered which it would be fallacious to generalize, for it is
  obvious that simplicity has been a part cause of their discovery, and can,
  therefore, give no ground for the supposition that other undiscovered laws
  are equally simple.''
  \qauthor{--- Bertrand Russel (1872--1970)}
\end{savequote}
\chapter{Governing equations for two-phase flow}
\label{chap:two-phase-flow}
In this chapter we will give a brief overview of the derivation of the
governing equations for two-phase flow.  We begin with a consideration of the
Navier-Stokes equations for single-phase flow.  We then introduce a singular
surface-force term and derive the Navier-Stokes equations for two-phase flow.
Finally, we use the Navier-Stokes equations for two-phase flow to derive jump
conditions at the interface between the phases.

\section{The Navier-Stokes equations for single-phase flow}
The following is a brief derivation of the Navier-Stokes equations for
single-phase flows.  For a more thorough derivation of these equations, see for
instance \citet[§4 and §5]{Aris89} or \citet[Chapter~4]{White03}.

We consider a single-phase, viscous flow in some domain $\Omega$ with boundary
$\partial\Omega$.  When temperature effects are neglected the flow is described
by the Cauchy equation
\begin{equation}
  \label{eq:cauchy_equation}
  \rho\left(\cd{\vct u}\right) = \div\tenstr+\rho\vct f,
\end{equation}
and the mass conservation equation,
\begin{equation}
  \label{eq:cons_of_mass}
  \pdt\rho + \div(\rho\vct u) = 0.
\end{equation}
Here $\rho$ is the fluid density, $\vct u$ is the flow velocity, $t$ is time,
$\tenstr$ is the stress tensor, $\vct f$ denotes body forces, and $\del$ is
used to denote the gradient and divergence operators.  The stress tensor for
Newtonian fluids with zero bulk viscosity is
\begin{equation}
  \tenstr = -p\tenI + 2\mu\tenD - \frac 2 3\mu (\trace\tenD) \tenI,
  \label{eq:viscous_stress_tensor}
\end{equation}
where $p$ is the pressure, $\tenI$ is the identity tensor, $\mu$ is the dynamic
viscosity, and $\trace\tenD$ denotes the trace of the strain-rate tensor
$\tenD$,
\begin{equation}
  \tenD = \frac{1}{2}\left( \grad\vct u + (\grad\vct u)^T \right).
  \label{eq:deformation_tensor}
\end{equation}
For incompressible flow, the mass-conservation equation reduces to
\begin{equation}
  \label{eq:navier-stokes-1}
  \div\vct u = 0,
\end{equation}
that is, the velocity field must be divergence free.  If we further assume that
the viscosity is constant, then it follows that the divergence of the stress
tensor reduces to
\begin{equation}
  \label{eq:stress_tensor}
  \div\tenstr = -\grad p + \mu\lapl\vct u,
\end{equation}
and so the Cauchy equation \eqref{eq:cauchy_equation} becomes
\begin{equation}
  \rho\left(\cd{\vct u}\right) = -\grad p + \mu\lapl\vct u + \rho\vct f.
  \label{eq:navier-stokes-2}
\end{equation}
Equations \eqref{eq:navier-stokes-1} and \eqref{eq:navier-stokes-2} are the
incompressible Navier-Stokes equations for single-phase flow with constant
viscosity.

\section{The Navier-Stokes equations for two-phase flow}
We now consider an immiscible two-phase flow of two Newtonian fluids, each with
its own viscosity and density.  We let $\Omega_1$ and $\Omega_2$ denote the
domains occupied by fluid 1 and fluid 2, respectively, and let the interface
between the fluids be denoted by $\Gamma$.  Then $\Omega
= \Omega_1\cup\Omega_2$ and $\partial\Omega
= (\partial\Omega_1\cup\partial\Omega_2)\setminus\Gamma$ are the fluid domain
and its boundary, respectively.  See \cref{fig:domain} for an illustration of
a two-phase flow domain.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [
    scale=0.8,
    >=stealth,
    wall/.style={
      decoration={border,angle=45,segment length=4},
      postaction={decorate,draw}},
    ]
    \draw[wall] (0,0) rectangle(10,5);
    \draw[fill=black!1] (3,1)
      .. controls (2,1) and (2,4) .. (4,3.5)
      .. controls (6,3) and (5,5) .. (7,4)
      .. controls (9,3) and (9,1.5) .. (8,1.5) node[below] {$\Gamma$}
      .. controls (4,1.5) and (4,1) .. (3,1);
    %\node (g2) [below=0.4cm of g1] {$\Gamma$};
    %\draw[->] (g2.west) to [out=150, in=240] (g1.west);
    \node at (1,1) {$\Omega_2$};
    \node at (4,2) {$\Omega_1$};
  \end{tikzpicture}
  \caption{An illustration of a two-phase flow domain.  The interface $\Gamma$
    separates the two phases $\Omega_1$ and $\Omega_2$.}
  \label{fig:domain}
\end{figure}

The extension of the single-phase model to account for two fluids can be made
by adding a singular surface-force term to represent the effects of surface
tension between the fluids.  We assume that the surface tension is constant, in
which case the singular surface force can be defined as
\begin{equation}
  \funssf(\vct x,t)
  = \int_{\Gamma}\sigma\kappa\vct n\delta(\vct x - \vct x_I(\vct s))\dif\vct s,
  \label{eq:ssf}
\end{equation}
where $\sigma$ is the surface tension, $\kappa$ is the local curvature, $\vct
n$ is the normal vector, $\delta$ is the Dirac delta function, and $\vct
x_I(\vct s)$ is a parametrisation of the interface.  Thus the Navier-Stokes
equations for immiscible and incompressible two-phase flow are
\begin{align}
  \label{eq:ns1}
  \div\vct u &= 0, \\
  \label{eq:ns2}
  \rho\left(\cd{\vct u}\right) &= -\grad p + \mu\lapl\vct u
    + \rho\vct f + \funssf.
\end{align}

\section{Fluid-fluid interface conditions}
The surface-tension force and the jump in viscosity and density across the
interface $\Gamma$ lead to a set of interface conditions that must be satisfied
along the interface $\Gamma$.  The following is a brief derivation of these
conditions.

First, we only consider flows where there is no mass transfer, which implies
that
\begin{equation}
  \jmp{\vct u}\cdot\vct n = 0,
\end{equation}
where $\jmp{\cdot}$ denotes the jump at the interface, for instance $\jmp{\mu}
= \mu_2 - \mu_1$.  Further, for viscous flows there is no slip at the
interface, and thus the tangential velocity component of the two fluids must be
equal at the interface,
\begin{equation}
  \jmp{\vct u}\cdot\vct t = 0.
\end{equation}
It follows that
\begin{align}
  \label{eq:jmpu}
  \jmp{\vct u} &= 0, \\
  \label{eq:jmpJut}
  \jmp{\grad \vct u}\cdot\vct t &= 0.
\end{align}
The latter is a direct consequence of the former.  Both fluids are
incompressible, which gives the trivial identity $\jmp{\div\vct u} = 0$.  If we
use the identity
\begin{equation}
  \div\vct u
    = \vct n\cdot\grad\vct u\cdot\vct n + \vct t\cdot\grad\vct u\cdot\vct t
  \label{eq:divuident}
\end{equation}
together with \eqref{eq:jmpJut} we get
\begin{equation}
  \vct n\cdot\jmp{\grad\vct u}\cdot\vct n = 0,
  \label{eq:jmpngradun}
\end{equation}
which means that the normal component of the normal derivative of the velocity
field is continuous across the interface.  Note that in \eqref{eq:divuident}
and in similar expressions in the following, the nabla operator is only applied
to $\vct u$, not the normal and tangential vectors that follow.

\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}

    \clip (-4.4,-3.7) rectangle (4.6,3.4);

    % Draw the domains
    \def\domain{(-3,2) .. controls (-5, 0) and ( 0,-5) .. ( 3,-2)
                       .. controls ( 6, 1) and (-1, 4) .. (-3, 2)}
    \fill[gray!5] [scale=1.2] \domain;
    \fill[white]  [scale=0.8] \domain;
    \draw [thick]          \domain;
    \draw [dotted,scale=1.2] \domain;
    \draw [dotted,scale=0.8] \domain;

    % Draw annotations
    \node at ( 1.0,-1.0) {$\Omega_2$};
    \node at (-3.0,-2.5) {$\Omega_1$};
    \draw[very thin,->,>=stealth]
        (-1,-0.5) node[right] {$\Gamma$} -- (-2.1,-1.5);
    \draw[very thin,<->,>=stealth]
        (1.1,-3) -- (1.1,-3.5) node[midway,right] {$\epsilon$};
    \clip (-3.5,1.5) rectangle (-2.5,2.5);
    \draw [densely dashed,very thick,gray!5] \domain;
    \node[fill=gray!5,fill opacity=0.5,text opacity=1.0]
        at (-3,2) {$\Omega_\epsilon$};

  \end{tikzpicture}
  \caption{A sketch of the control volume $\Omega_\epsilon$ that covers the
    interface, $\Gamma$.}
  \label{fig:domain_omega_epsilon}
\end{figure}%
Next, we consider the conservation of momentum.  We define a control volume
\begin{equation}
  \Omega_{\epsilon} = \set{\vct x \in \Omega}{\min_{\vct x_I \in \Gamma}|\vct
    x - \vct x_I| \leq \epsilon},
\end{equation}
where $\epsilon > 0$, see \cref{fig:domain_omega_epsilon}.  We then integrate
the Cauchy equation \eqref{eq:cauchy_equation} with the added singular
surface-force term over the domain~$\Omega_\epsilon$,
\begin{multline}
  \intoepdx{\rho\td{\vct u}}
    = \intoepdx{\div\tenstr} + \intoepdx{\rho\vct f} \\
      + \intoepdx{ \int_{\Gamma}\sigma\kappa\vct n\delta(\vct x - \vct
        x_I(\vct s))\dif\vct s}.
\end{multline}
Here $\tdi{\vct u} = \partial\vct u/\partial t + \vct u\cdot\grad\vct u$
denotes the convective derivative.  Now we apply the Gauss theorem and change
the order of integration in the last term to obtain
\begin{equation}
  \intoepdx{\rho\td{\vct u}} = \intdoepds{\tenstr\cdot\vct n}
    + \intoepdx{\rho\vct f} + \int_\Gamma\sigma\kappa\vct n\dif\vct s.
\end{equation}
If we let $\epsilon$ go to zero, the left-hand side and the second term on the
right-hand side vanish, and we get
\begin{equation}
  0 = \int_\Gamma\left( \jmp{\tenstr}\cdot\vct n + \sigma\kappa\vct
    n \right)\dif s.
  \label{eq:stressjumpint}
\end{equation}
Since the derivation above is also valid for any subset of $\Omega_\epsilon$
containing a part of $\Gamma$, \eqref{eq:stressjumpint} must hold for any part
of $\Gamma$.  Therefore
\begin{equation}
  0 = \jmp{\vct T}\cdot\vct n + \sigma\kappa\vct n.
\end{equation}
With $\vct T = -p\vct I + 2\mu\vct D$ we finally get the interface condition
for the stresses,
\begin{equation}
  \label{eq:stress_jump}
  \jmp{p}\vct n - \jmp{2\mu\tenD}\cdot\vct n = \sigma\kappa\vct n.
\end{equation}
The surface tension force is seen to introduce a discontinuity in the normal
stresses across the interface.  The tangential stresses are continuous.

\section{Interface conditions for the pressure and the viscous term}
The previous section gave a brief derivation of the interface conditions for
immiscible and incompressible two-phase flow without mass transfer.  In order
to use these conditions for the discretization of the Navier-Stokes equations
\eqref{eq:ns1} and \eqref{eq:ns2}, we need to rewrite them into a more suitable
form.  In particular, we want to find explicit expressions for the jump in the
pressure $\jmp p$ and the viscous term $\jmp{\mu\grad\vct u}$.  The following
derivation is based on \cite{Hansen05} and \cite{Kang00}.

First, a jump condition for the pressure is obtained by taking the inner
product of \eqref{eq:stress_jump} with the normal vector $\vct n$,
\begin{equation}
  \jmp p = \vct n\cdot\jmp{2\mu\tenD}\cdot\vct n + \sigma\kappa
      = 2\jmp{\mu}\vct n\cdot \grad\vct u \cdot\vct n
        + \sigma\kappa,
  \label{eq:jump_pressure}
\end{equation}
where \eqref{eq:jmpngradun} was used for the second equality.

To find the jump in the viscous term, it is first decomposed into an interface
normal coordinate system as
\begin{equation}
  \begin{split}
    \jmp{\mu\grad\vct u} &=
  \left(\vct n\cdot\jmp{\mu\grad\vct u}\cdot\vct n\right)\vct n\otimes\vct n +
  \left(\vct t\cdot\jmp{\mu\grad\vct u}\cdot\vct t\right)\vct t\otimes\vct t
        \\ &\quad+
  \left(\vct n\cdot\jmp{\mu\grad\vct u}\cdot\vct t\right)\vct n\otimes\vct t +
  \left(\vct t\cdot\jmp{\mu\grad\vct u}\cdot\vct n\right)\vct t\otimes\vct n,
  \end{split}
\end{equation}
where $\otimes$ denotes the dyadic product.  We already showed that
$\jmp{\grad\vct u}\cdot\vct t = 0$ and $\vct n\cdot\jmp{\grad\vct u}\cdot\vct
n = 0$, cf.\ \eqref{eq:jmpJut} and \eqref{eq:jmpngradun}, which gives
\begin{align}
  \vct n\cdot\jmp{\mu\grad\vct u}\cdot\vct n &=
    \jmp{\mu}\,\vct n\cdot\grad\vct u\cdot\vct n, \\
  \vct t\cdot\jmp{\mu\grad\vct u}\cdot\vct t &=
    \jmp{\mu}\,\vct t\cdot\grad\vct u\cdot\vct t, \\
  \vct n\cdot\jmp{\mu\grad\vct u}\cdot\vct t &=
    \jmp{\mu}\,\vct n\cdot\grad\vct u\cdot\vct t.
\end{align}
We then take the inner product of \eqref{eq:stress_jump} with $\vct t$, which
gives
\begin{equation}
  \vct t\cdot\jmp{\mu\left(\grad\vct u + (\grad\vct u)^T \right)}\cdot\vct n
  = \vct t\cdot\jmp{\mu\grad\vct u} \cdot\vct n
    + \jmp{\mu}\vct n\cdot\grad\vct u \cdot\vct t
  = 0,
\end{equation}
or
\begin{equation}
  \vct t\cdot\jmp{\mu\grad\vct u}\cdot\vct n =
    - \jmp{\mu}\vct n\cdot\grad\vct u \cdot\vct t.
\end{equation}
The jump in the viscous term becomes
\begin{multline}
  \jmp{\mu\grad\vct u} = \jmp{\mu} \Big(
       (\vct n\cdot\grad\vct u\cdot\vct n)\vct n\otimes\vct n
     + (\vct n\cdot\grad\vct u\cdot\vct t)\vct n\otimes\vct t \\
     - (\vct n\cdot\grad\vct u\cdot\vct t)\vct t\otimes\vct n
     + (\vct t\cdot\grad\vct u\cdot\vct t)\vct t\otimes\vct t \Big).
\end{multline}
This can be simplified further by noting that
\begin{equation}
  (\grad\vct u\cdot t)\otimes\vct t
  = (\vct n\cdot\grad\vct u\cdot\vct t)\vct n\otimes\vct t
    + (\vct t\cdot\grad\vct u\cdot\vct t)\vct t\otimes\vct t.
\end{equation}
Thus we obtain the expression for the jump in the viscous term that has been
used in the present work,
\begin{multline}
  \jmp{\mu\grad\vct u} = \jmp{\mu} \Big(
       (\vct n\cdot\grad\vct u\cdot\vct n)\vct n\otimes\vct n
     - (\vct n\cdot\grad\vct u\cdot\vct t)\vct t\otimes\vct n \\
     + (\grad\vct u\cdot\vct t)\otimes\vct t \Big).
  \label{eq:jumptens}
\end{multline}

\section{Summary}
This chapter has given a brief derivation of the immiscible and incompressible
Navier-Stokes equations for two-phase flow \eqref{eq:ns1} and \eqref{eq:ns2},
where the viscosity is assumed constant in each phase.  The effect of surface
tension is included as a singular source term.

Fluid-fluid interface conditions have been discussed, and it has been shown how
the singular source term and a jump in viscosity lead to a set of interface
jump conditions for the pressure \eqref{eq:jump_pressure} and the viscous term
\eqref{eq:jumptens}.

% Fakesection: Quote
\begin{savequote}[7cm]
  ``As long as you recognize your sinful ways, correct your behaviour and adopt
  the right method and the right step size, your misdemeanours will be forgiven
  and your solution will prosper.''
  \qauthor{--- Arieh Iserles (1947)}
\end{savequote}
\chapter{Numerical methods}
\label{chap:numerical-methods}
This chapter describes the numerical methods that have been used in this thesis
to solve the Navier-Stokes equations for two-phase flow.  First, a brief
introduction to the level-set method is given, which is followed by an outline
of the two projection methods that have been used to solve the Navier-Stokes
equations.  The spatial and temporal discretization methods are then
summarized, and at the end of the chapter the new methods to calculate the
curvature and the normal vectors are presented.

\section{The level-set method}
\label{sec:level-set}
In order to solve the Navier-Stokes equations for two-phase flow, we need to
know the location of the interface.  The level-set method proposed by
\citet{Osher88} allows us to capture the interface location as the zero level
set of the level-set function $\varphi(\vct x,t)$.  The level-set function is
typically defined as a signed-distance function,
\begin{equation}
  \varphi(\vct x,t) = \begin{cases}
     d(\vct x,t) & \text{if } \vct x\in\Omega_2, \\
    -d(\vct x,t) & \text{if } \vct x\in\Omega_1, \\
  \end{cases}
\end{equation}
where $d(\vct x,t)$ is the shortest distance to the interface $\Gamma$,
\begin{equation}
  d(\vct x,t) = \min_{\vct x_I \in \Gamma} |\vct x - \vct x_I|.
\end{equation}
Thus the interface can be defined implicitly as
\begin{equation}
  \Gamma(t) = \set{\vct x\in\Omega}{\varphi(\vct x,t)=0},
              \quad t\in \mathbb R^+.
\end{equation}

The position of the interface is updated by solving an advection equation for
$\varphi$,
\begin{equation}
  \label{eq:ls_adeq}
  \pdt{\varphi} + \vct{\hat u}\cdot\grad\varphi = 0,
\end{equation}
where $\vct{\hat u}$ is the velocity at the interface extended to the entire
domain.  We extend the interface velocity through solving
a velocity-extrapolation equation,
\begin{equation}
  \label{eq:ls_velextr}
  \pd{\vct{\hat u}}{\tau} + S(\varphi)\vct n\cdot\grad\vct{\hat u} = 0,
  \quad \vct{\hat u}_{\tau=0} = \vct u,
\end{equation}
to steady state, cf.\ \cite{Adalsteinsson99,Zhao96}.  Here $\tau$ is
a pseudo-time and $S$ is a smeared sign function which is equal to zero at the
interface,
\begin{equation}
  S(\varphi) = \frac{\varphi}{\sqrt{\varphi^2+2\Delta x^2}}.
\end{equation}

When we solve the level-set equation \eqref{eq:ls_adeq}, the non-uniform
advection will distort the signed-distance property of the level-set function,
and numerical dissipation error adds to this distortion.  The level-set
function is therefore reinitialized regularly by solving
\begin{equation}
  \label{eq:ls_reinit}
  \begin{split}
    \pd{\varphi}{\tau} + S(\varphi_0)(|\grad\varphi|-1) &= 0, \\
    \varphi(\vct x,0) &= \varphi_0(\vct x),
  \end{split}
\end{equation}
to steady state as proposed by \citet{Sussman94}.  The level-set function just
before initialization is used as the initial condition $\varphi_0$.

The level-set equations \eqref{eq:ls_adeq}, \eqref{eq:ls_velextr}, and
\eqref{eq:ls_reinit} are discretized in time and space as described in the
following sections.  The method presented by \citet{Adalsteinsson95} is used to
improve the computational speed.  The method is often called the narrow-band
method, since the level-set function is only updated in a narrow band across
the interface at each time step.

One of the advantages of the level-set method is that normal vectors and
curvatures can be readily calculated from the level-set function, that is,
\begin{align}
  \label{eq:norm}
  \vct n &= \frac{\grad\varphi}{|\grad\varphi|}, \\
  \label{eq:curv}
  \kappa &= \div\left(\frac{\grad\varphi}{|\grad\varphi|}\right).
\end{align}

\section{Projection methods}
\label{sec:projection-method}
We have employed two different projection methods in this thesis to solve the
incompressible Navier-Stokes equations for two-phase flow \eqref{eq:ns1} and
\eqref{eq:ns2}.  Papers A, B, and C used the direct projection
scheme~1~(DP1)~\cite{Hansen05}, and Paper D used the more standard Chorin
projection method \cite{Chorin68}.

The projection methods are a family of methods for the solution of the
incompressible Navier-Stokes equations that are based on the Helmholtz-Hodge
theorem.  This theorem states that an arbitrary vector field can be decomposed
into a divergence-free part and a rotation-free part.  That is, any vector
field $\vct a$ can be written as
\begin{equation}
  \vct a = \vct a' + \grad\psi,
\end{equation}
where $\vct a'$ is a vector with $\div\vct a' = 0$ and $\psi$ is a scalar
potential.  The proof of this theorem can be found in for instance
\cite[§3.44]{Aris89} or \cite{Chorin00}.

\subsection{The direct projection scheme 1}
The DP1 was developed by \citet{Hansen05} and is based on a direct application
of the Helmholtz-Hodge theorem.  We assume that the velocity $\vct u$ is
sufficiently smooth, and then rewrite \eqref{eq:navier-stokes-2} to get
\begin{equation}
  \label{eq:divFreeVel}
  \pdt{\vct u} = \vct w - \frac{\grad p}{\rho},
\end{equation}
where
\begin{equation}
  \vct w = - (\vct u\cdot\grad)\vct u + \nu\grad^2\vct u + \vct f.
  \label{eq:continuousRightHandVector}
\end{equation}
The divergence of \eqref{eq:divFreeVel} yields a Poisson equation for the
pressure,
\begin{equation}
  \div\left(\frac{\grad p}{\rho}\right) = \div\vct w.
  \label{eq:PoissonDP1}
\end{equation}
The DP1 scheme follows from a direct numerical discretization of the above
equations.  First, $\vct w$ is calculated with
\eqref{eq:continuousRightHandVector}.  Then the pressure is found by solving
the Poisson equation \eqref{eq:PoissonDP1}.  Finally, an Euler step is used to
solve \eqref{eq:divFreeVel}, that is,
\begin{equation}
  \vct u^{n+1} = \vct u^n + \Delta t\left(\vct w - \frac{\grad p}{\rho}\right).
  \label{eq:dp1-update}
\end{equation}

\subsection{The Chorin projection method}
The Chorin projection method was presented by \citet{Chorin68}, and is today
one of the standard methods for solving the Navier-Stokes equations, e.g.\
\cite{Walker12}.  The method is briefly presented in the following.

Let $\Delta t$ be the time step, and consider the discretization of the
momentum equation \eqref{eq:navier-stokes-2} with the forward Euler method,
\begin{equation}
  \frac{\vct u^{n+1} - \vct u^n}{\Delta t}
  = -\vct u^n\cdot\grad\vct u^n
    - \frac{\grad p^{n+1}}{\rho} + \frac\mu\rho\lapl\vct u^n + \vct f,
  \label{eq:discmomeq}
\end{equation}
where $\vct u^n\equiv\vct u(\vct x,n\Delta t)$ and $p^n\equiv p(\vct x,n\Delta
t)$ are assumed known at time level $n$.  Note that the pressure gradient is
evaluated at time level $n+1$.  Next, to solve \eqref{eq:discmomeq} in two
steps, we introduce the intermediate velocity field $\vct u^\star$,
\begin{equation}
  \frac{\vct u^{n+1} - \vct u^\star + \vct u^\star - \vct u^n}{\Delta t}
  = -\vct u^n\cdot\grad\vct u^n
    - \frac{\grad p^{n+1}}{\rho} + \frac\mu\rho\lapl\vct u^n + \vct f,
\end{equation}
which is chosen such that
\begin{align}
  \label{eq:ustar}
  \frac{\vct u^\star - \vct u^n}{\Delta t}
    &= -\vct u^n\cdot\grad\vct u^n + \frac\mu\rho\lapl\vct u^n + \vct f, \\
  \label{eq:unext}
  \frac{\vct u^{n+1} - \vct u^\star}{\Delta t} &= - \frac{\grad p^{n+1}}{\rho}.
\end{align}
From \eqref{eq:ustar} we get an explicit expression for $\vct u^\star$,
\begin{equation}
  \vct u^{\star} = \vct u^n
    + \Delta t\left(-\vct u^n\cdot\grad\vct u^n
        + \frac\mu\rho\,\lapl\vct u^n + \vct f\right).
\end{equation}
Next, the divergence of \eqref{eq:unext} and $\div\vct u^{n+1} = 0$ yields
a Poisson equation for the pressure,
\begin{equation}
  \div\left(\frac{\grad p^{n+1}}{\rho}\right)
    = \frac{\div\vct u^\star}{\Delta t}.
  \label{eq:PoissonChorin}
\end{equation}
Finally, we obtain for \eqref{eq:unext}
\begin{equation}
  \vct u^{n+1} = \vct u^\star - \Delta t\frac{\grad p^{n+1}}{\rho}.
  \label{eq:velcorrection}
\end{equation}

\section{Spatial discretization}
\label{sec:x-discr}
The governing equations \eqref{eq:ns1}, \eqref{eq:ns2}, \eqref{eq:ls_adeq},
\eqref{eq:ls_velextr}, and \eqref{eq:ls_reinit} are discretized on a staggered
grid~\cite{Harlow65}, where the scalar values are located at the cell centres
and the vector values are located at the cell edges, see
\cref{fig:staggered-grid}.  The domain boundary coincides with cell edges, and
the fixed grid spacing is $\Delta x$ in the $x$ direction and $\Delta y$ in the
$y$ direction.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [scale=1.5]

    \draw[very thin] (-2,-2) grid (1,1);
    \node[below] at (-1.5,-2.1) {$i-1$};
    \node[below] at (-0.5,-2.1) {$i$};
    \node[below] at ( 0.5,-2.1) {$i+1$};
    \node[left]  at (-2.1,-1.5) {$j-1$};
    \node[left]  at (-2.1,-0.5) {$j$};
    \node[left]  at (-2.1, 0.5) {$j+1$};
    \foreach \x in {-1.5,-0.5,0.5}
      \foreach \y in {-1.5,-0.5,0.5}
        \fill (\x,\y) circle(1pt);
    \foreach \x in {-2,-1,0,1}
      \foreach \y in {-1.5,-0.5,0.5} {
        \fill (\x,\y)+(-1pt,-1pt) rectangle +(1pt,1pt);
        \draw[very thin] (\y,\x)+(-1pt,-1pt) rectangle +(1pt,1pt); };

    \draw[decorate,decoration=brace]
      ( 1.1, 0.0) -- node[right=0.2em] {$\Delta y$} (1.1,-1.0);
    \draw[decorate,decoration=brace]
      (-1.0, 1.1) -- node[above=0.2em] {$\Delta x$} (0.0, 1.1);
  \end{tikzpicture}
  \caption{An illustration of a small part of a uniform staggered grid, where
    the fixed grid spacings $\Delta x$ and $\Delta y$ are indicated.  The
    scalar values are stored at the cell centres (filled circles), the
    $x$-component of the vector values is stored at the cell edges $(i+\frac
    1 2,j)$ (filled squares), and the $y$-component of the vector values is
    stored at the cell edges $(i,j+\frac 1 2)$ (open squares).}
  \label{fig:staggered-grid}
\end{figure}

The $x$- and $y$-derivatives, divergence, and Laplacian operators are
discretized by the second-order central-difference scheme,
\begin{align}
  \left. G_x p \right|_{i+\frac 1 2,j}
    &= \frac{p_{i+1,j} - p_{i,j}}{\Delta x}, \\
  \left. G_y p \right|_{i,j+\frac 1 2}
    &= \frac{p_{i,j+1} - p_{i,j}}{\Delta y}, \\
  \left. \vct D \cdot \vct u \right|_{i,j}
    &=   \frac{u_{i+1/2,j} - u_{i-1/2,j}}{\Delta x}
       + \frac{v_{i,j+1/2} - v_{i,j-1/2}}{\Delta y}, \\
  \left. Lp \right|_{i,j}
    &=   \frac{p_{i+1,j} - 2p_{i,j} + p_{i-1,j}}{(\Delta x)^2}
       + \frac{p_{i,j+1} - 2p_{i,j} + p_{i,j-1}}{(\Delta y)^2},
\end{align}
respectively.  Here $\vct u = (u,v)$ is a vector and $p$ is a scalar.  Note
that the gradient of a scalar is a vector and has values located at the cell
edges.  These definitions are consistent with the staggered grid, since the
gradient returns a vector with components defined on the cell faces, while the
divergence and the Laplacian returns a scalar defined at the cell centres.

\subsection{Advective terms}
The advective terms, $\vct u\cdot\grad\vct u$ in the momentum equation
\eqref{eq:ns2} and $\vct{\hat u}\cdot\grad\varphi$ in the level-set equation
\eqref{eq:ls_adeq}, and the normal derivative $\vct n\cdot\grad\vct{\hat u}$ in
the velocity-extrapolation equation \eqref{eq:ls_velextr} are discretized with
the weighted essentially non-oscillatory (WENO) scheme, cf.~\cite{Fedkiw99}
and~\cite{Jiang96}.  The WENO scheme is a high-order upwind scheme that is
fifth-order accurate in smooth regions.  In nonsmooth regions the accuracy is
reduced to a minimum of third order.  The following is a brief outline of the
WENO scheme.

First, when we calculate the advective term for the velocity, the vector
components are required on all cell edges.  That is, we must interpolate the
$x$-component of the velocity to the location of the $y$-component of the
velocity and vice versa.  To do this we use a linear interpolation,
\begin{align}
  u_{i,j+\frac 1 2} &= \frac 1 4\left(
      u_{i-\frac 1 2,j  } + u_{i-\frac 1 2,j+1}
    + u_{i+\frac 1 2,j+1} + u_{i+\frac 1 2,j  } \right), \\
  v_{i+\frac 1 2,j} &= \frac 1 4\left(
      v_{i,  j-\frac 1 2} + v_{i+1,j-\frac 1 2}
    + v_{i+1,j+\frac 1 2} + v_{i  ,j+\frac 1 2} \right).
\end{align}

Similarly, when we calculate the advective term for the level-set function
$\varphi$, we first interpolate the velocity to the scalar grid,
\begin{align}
  u_{i,j} &= \frac 1 2\left( u_{i+\frac 1 2,j} + u_{i-\frac 1 2,j} \right), \\
  v_{i,j} &= \frac 1 2\left( v_{i,j+\frac 1 2} + v_{i,j-\frac 1 2} \right).
\end{align}

We now consider the WENO scheme for the advective operator in 1D at the point
$x_{i+\frac 1 2}$.  The scheme extends naturally to higher dimensions.  First,
if $u_{i+\frac 1 2} = 0$, then
\begin{equation}
  \left. u\pd u x \right|_{i+\frac 1 2} = 0.
\end{equation}
Otherwise we need to calculate a set of five differences that depend on the
upwind direction.  The differences are denoted $\Delta v_k$ for $k=1,\dots,5$.
If $u_{i+\frac 1 2}>0$, then we calculate
\begin{equation}
  \Delta v_k = \frac{u_{i+\frac{2k-5} 2} - u_{i+\frac{2k-7} 2}}{\Delta x}.
\end{equation}
Otherwise if $u_{i+\frac 1 2}<0$, we calculate
\begin{equation}
  \Delta v_k = \frac{u_{i-\frac{2k-9} 2} - u_{i-\frac{2k-7} 2}}{\Delta x}.
\end{equation}
Next we calculate expressions for the smoothness of three substencils,
\begin{align}
  S_1 &= \frac{13}{12}(\Delta v_1 - 2\Delta v_2 + \Delta v_3)^2
    + \frac 1 4 (\Delta v_1 - 4\Delta v_2 + 3\Delta v_3)^2, \\
  S_2 &= \frac{13}{12}(\Delta v_2 - 2\Delta v_3 + \Delta v_4)^2
    + \frac 1 4 (\Delta v_2 - \Delta v_4)^2, \\
  S_3 &= \frac{13}{12}(\Delta v_3 - 2\Delta v_4 + \Delta v_5)^2
    + \frac 1 4 (3\Delta v_3 - 4\Delta v_4 + \Delta v_5)^2,
\end{align}
where a small $S$ indicates a smooth substencil.  These smoothness factors are
then used to compute weights for the substencils,
\begin{equation}
  w_k = \frac{b_k}{b_1+b_2+b_3},
\end{equation}
for $k = 1,2,3$, where
\begin{equation}
  b_1 = \frac 1 {10}\frac 1 {(\epsilon+S_1)^2},\quad
  b_2 = \frac 6 {10}\frac 1 {(\epsilon+S_2)^2},\quad
  b_3 = \frac 3 {10}\frac 1 {(\epsilon+S_3)^2}.
\end{equation}
Here $\epsilon$ is a regularization parameter that is used to avoid division by
zero.  We have used $\epsilon = 10^{-6}$ in this work.  Finally, the WENO
scheme for the gradient becomes
\begin{multline}
  \left. \pd u x \right|_{i+\frac 1 2} \simeq
      w_1\left(\frac 1 3\Delta v_1 - \frac 7 6\Delta v_2
        + \frac{11} 6 \Delta v_3\right) \\
    + w_2\left(-\frac 1 6\Delta v_2 + \frac 5 6\Delta v_3
        + \frac 1 3 \Delta v_4\right) \\
    + w_3\left(\frac1 3\Delta v_3 + \frac 5 6\Delta v_4
        - \frac 1 6 \Delta v_5\right).
\end{multline}

\subsection{The viscous term}
The viscous term $\mu\lapl\vct u$ in the Navier-Stokes equations is discretized
by standard second-order central differences using the ghost-fluid method
(GFM)~\cite{Fedkiw99,Kang00,Liu00}.  This method includes the jump in the
viscous term at the interface \eqref{eq:jumptens} in the discretization stencil
in a sharp manner.

In the following we present the GFM scheme in the one-dimensional case.  The
extension to higher dimensions is straightforward.  In order to simplify the
notation, we omit the half indices and let $k\equiv i+\frac 1 2$.  Further, we
consider a general case where
\begin{equation}
  \jmp u = a_\Gamma, \qquad \jmp{\mu\pd u x} = b_\Gamma.
\end{equation}
For the viscous term, $a_\Gamma = 0$, cf.\ \eqref{eq:jmpu}, and $b_\Gamma$ is
given by the jump tensor~\eqref{eq:jumptens}.  The jump tensor is calculated at
the cell centres near the interface with the second-order central-difference
scheme.  It is then interpolated linearly to the cell edges, from where it is
again interpolated linearly to the interface when needed.

If the interface does not cross the stencil, then the GFM scheme reduces to the
standard second-order central-difference stencil.  Else there are four
stencils, depending on the location of the interface.  One such interface
configuration is sketched in \cref{fig:gfm-scheme-example}.  In all the
stencils, $\theta$ is defined as the relative distance from the interface to
the node on the left, for instance
\begin{equation}
  \theta = \frac{|\varphi_{k}|}{|\varphi_{k}| + |\varphi_{k+1}|},
\end{equation}
where $\varphi_k$ and $\varphi_{k+1}$ are the level-set function values
linearly interpolated to the vector grid.  The four stencils are given below,
where $h\equiv\Delta x$.
\begin{enumerate}
  \item Phase 1 is to the left and interface lies between $k$ and $k+1$:
    \begin{align}
      \begin{split}
        \left.\pd{} x\left(\mu\pd u x\right)\right|_{x_k}
          &= \frac{\hat\mu\left(u_{k+1} - u_{k}\right)
                 - \mu_1\left(u_{k} - u_{k-1}\right)}{h^2} \\
          &\qquad
             - \frac{\hat\mu a_\Gamma}{h^2}
             - \frac{\hat\mu b_\Gamma (1-\theta)}{\mu_2 h},
      \end{split} \\
       \hat\mu &= \frac{\mu_1\mu_2}{\theta\mu_2 + (1-\theta)\mu_1}.
    \end{align}
  \item Phase 1 is to the left and interface lies between $k-1$ and $k$:
    \begin{align}
      \begin{split}
        \left.\pd{} x\left(\mu\pd u x\right)\right|_{x_k}
          &= \frac{\mu_2\left(u_{k+1} - u_{k}\right)
                - \hat\mu\left(u_{k} - u_{k-1}\right)}{h^2} \\
          &\qquad
            + \frac{\hat\mu a_\Gamma}{h^2}
            - \frac{\hat\mu b_\Gamma \theta}{\mu_1 h},
      \end{split} \\
      \hat\mu &= \frac{\mu_1\mu_2}{\theta\mu_2 + (1-\theta)\mu_1}.
    \end{align}
  \item Phase 2 is to the left and interface lies between $k$ and $k+1$:
    \begin{align}
      \begin{split}
        \left.\pd{} x\left(\mu\pd u x\right)\right|_{x_k}
          &= \frac{\hat\mu\left(u_{k+1} - u_{k}\right)
                - \mu_2\left(u_{k} - u_{k-1}\right)}{h^2} \\
          &\qquad
            + \frac{\hat\mu a_\Gamma}{h^2}
            + \frac{\hat\mu b_\Gamma (1-\theta)}{\mu_1 h},
      \end{split} \\
      \hat\mu &= \frac{\mu_1\mu_2}{\theta\mu_1 + (1-\theta)\mu_2}.
    \end{align}
  \item Phase 2 is to the left and interface lies between $k-1$ and $k$:
    \begin{align}
      \begin{split}
        \left.\pd{} x\left(\mu\pd u x\right)\right|_{x_k}
          &= \frac{\mu_1\left(u_{k+1} - u_{k}\right)
                - \hat\mu\left(u_{k} - u_{k-1}\right)}{h^2} \\
          &\qquad
            - \frac{\hat\mu a_\Gamma}{h^2}
            + \frac{\hat\mu b_\Gamma \theta}{\mu_2 h},
      \end{split} \\
      \hat\mu &= \frac{\mu_1\mu_2}{\theta\mu_1 + (1-\theta)\mu_2}.
    \end{align}
\end{enumerate}
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}[scale=2.0]
    \fill[black!2] (-0.9,0) rectangle (0.45,1);
    \draw[dotted] (-0.9,0) grid (2.5,1);
    \draw[thick] (0.45,0) -- (0.45,1) node[above] {$\Gamma$};
    \fill (0.0,0.5) circle(1pt) node[above left] {$k-1$};
    \fill (1.0,0.5) circle(1pt) node[above left] {$k$};
    \fill (2.0,0.5) circle(1pt) node[above left] {$k+1$};
  \end{tikzpicture}
  \caption{One-dimensional case where the interface $\Gamma$ separates
    fluid 1 on the left from fluid 2 on the right.}
  \label{fig:gfm-scheme-example}
\end{figure}

\subsection{Pressure Poisson equation}
\label{sec:pressure-poisson}
The Laplace operator in the Poisson equations \eqref{eq:PoissonDP1} and
\eqref{eq:PoissonChorin} is also discretized by standard second-order central
differences using the GFM~\cite{Kang00,Liu00}.  The GFM was discussed in the
previous section, where stencils were given for the viscous term in the
one-dimensional case.  For the pressure, the stencils are the same except we
use $p$ instead of $u$, and $1/\rho$ and $1/\hat\rho$ instead of $\mu$ and
$\hat\mu$.  To find the jump in the pressure at the interface, $\jmp
p = a_\Gamma$, the jump condition \eqref{eq:jump_pressure} is calculated at the
cell centres and then interpolated to the interface.  In addition, we use that
$b_\Gamma = 0$, which is justified by \citet[Section~3.7]{Kang00}.

The resulting linear system for the pressure is solved with a solver from the
Portable, Extensible Toolkit for Scientific Computation (PETSc)~\cite{PETSc}.
PETSc makes available a large selection of solvers.  In most cases, we used the
direct solver based on LU factorisation or the conjugate gradient method with
incomplete Cholesky factorisation.  In some cases, we used the GMRES solver
with a preconditioner, either incomplete LU factorisation or an algebraic
multigrid method.  See \cite{PETSc} for more details and a list of available
methods, and see for example~\cite{Iserles09} or~\cite{Saad03} for an
introduction to linear solvers.

Finally, the gradient of the pressure, used in \eqref{eq:divFreeVel} and
\eqref{eq:velcorrection}, is also calculated with the GFM.  In one dimension,
the stencil is
\begin{align}
  \left. \frac 1 \rho \pd p x\right|_{i+\frac 1 2}
           &= \frac 1 {\hat\rho} \frac{(p_{i+1} - a_\Gamma) - p_i} h, \\
  \hat\rho &= \theta\rho_1 + (1-\theta)\rho_2,
\end{align}
if $\phi_i\leq 0$ and $\phi_{i+1}>0$, or else
\begin{align}
  \left. \frac 1 \rho \pd p x\right|_{i+\frac 1 2}
           &= \frac 1 {\hat\rho} \frac{(p_{i+1} + a_\Gamma) - p_i} h, \\
  \hat\rho &= \theta\rho_2 + (1-\theta)\rho_1,
\end{align}
where
\begin{equation}
  \theta = \frac{|\varphi_{i}|}{|\varphi_{i}| + |\varphi_{i+1}|}.
\end{equation}

\section{Temporal discretization}
\label{sec:t-discr}
The temporal discretization is done with the explicit strong
stability-preserving Runge-Kutta (SSP-RK) schemes, see
\cite{Ketcheson05,Shu88}.  The idea behind the SSP-RK schemes is to preserve
the stability of a low-order method when it is extended to higher order.  It is
argued in \cite{Gottlieb01} that if one extends the forward Euler method, which
is total variation diminishing (TVD) for a suitable discretization of a scalar
conservation law, to a non-SSP higher-order method, then overshoots may occur
at discontinuities.

An SSP-RK method can be written as a convex sum of explicit Euler steps with
time step $\Delta t$,
\begin{equation}
  \mathcal E\left(x^n\right) = x^n + \Delta t \mathcal F(x^n,t^n).
\end{equation}
where $\mathcal F$ is the residual of the PDE to be solved.  In this thesis we
have used two methods:  The third-order three-stage SSP-RK method~\cite{Shu88},
\begin{equation}
  \begin{split}
    x^{(1)} &= \mathcal E\left(x^n\right), \\
    x^{(2)} &= \frac 3 4 x^n + \frac 1 4\mathcal E\left(x^{(1)}\right), \\
    x^{n+1} &= \frac 1 3 x^n + \frac 2 3\mathcal E\left(x^{(2)}\right),
  \end{split}
  \label{eq:ssprk33}
\end{equation}
and the third-order four-stage SSP-RK method~\cite{Kraaijevanger91},
\begin{equation}
  \begin{split}
    x^{(1)} &= \frac 1 2 x^n + \frac 1 2\mathcal E\left(x^n\right), \\
    x^{(2)} &= \frac 1 2 x^{(1)} + \frac 1 2\mathcal E\left(x^{(1)}\right), \\
    x^{(3)} &= \frac 2 3 x^n + \frac 1 6 x^{(2)}
               + \frac 1 6\mathcal E\left(x^{(2)}\right), \\
    x^{n+1} &= \frac 1 2 x^{(3)} + \frac 1 2\mathcal E\left(x^{(3)}\right).
  \end{split}
  \label{eq:ssprk34}
\end{equation}

The semi-discretized velocity-extrapolation equation \eqref{eq:ls_velextr} and
level-set reinitialization equation \eqref{eq:ls_reinit} can be written as
systems of ordinary differential equations (ODEs) of the form
\begin{equation}
  \od{\vct \psi} \tau = \vct F(\vct\psi,\tau),
\end{equation}
where $\vct\psi$ is a vector containing the discrete variables and $\vct F$
contains the spatially discretized terms.  The explicit Euler step becomes
\begin{equation}
  \mathcal E(\vct\psi^n) = \vct\psi^n + \Delta\tau \vct F(\vct\psi^n,\tau^n).
\end{equation}
In this thesis, these equations are solved with the third-order four-stage
SSP-RK method \eqref{eq:ssprk34}.  The four-stage method is used because it is
more accurate than the three-stage method.  We wanted to make sure that the
error is mainly dominated by the spatial discretization.

The semi-discretized Navier-Stokes equations \eqref{eq:navier-stokes-1} and
\eqref{eq:navier-stokes-2} are solved together with the level-set advection
equation \eqref{eq:ls_adeq} with the third-order three-stage SSP-RK method
\eqref{eq:ssprk33}.  Here an explicit Euler step $\mathcal E(\vct
u^n,\varphi^n)$ consists of the following steps:
\begin{enumerate}
  \item Calculate the curvature and the normal vectors from $\varphi^n$ with
    one of the methods described in \cref{sec:curvature}.
  \item Calculate the intermediate velocity field, either $\vct w$
    \eqref{eq:continuousRightHandVector} or $\vct u^\star$ \eqref{eq:ustar}.
  \item Solve the Poisson equation \eqref{eq:PoissonDP1} or
    \eqref{eq:PoissonChorin} for the pressure as described in
    \cref{sec:pressure-poisson}.
  \item Correct the velocity field with either \eqref{eq:dp1-update} or
    \eqref{eq:velcorrection}.
  \item Advect the level-set function $\varphi^n$ with \eqref{eq:ls_adeq}.
\end{enumerate}

\section{Time step restriction}
\label{sec:deltaT}
We employ the Courant-Friedrich-Lewy~(CFL) condition to allow adaptive time
stepping and to enforce stability.  The CFL condition that is used in this
thesis is
\begin{equation}
  \Delta t = \frac{C}{\frac{C_c + C_v} 2
             + \sqrt{(C_c + C_v)^2 + 4C_g^2 + 4C_s^2}},
\end{equation}
where the CFL restriction $C<1$, and $C_c$, $C_v$, $C_s$, and $C_g$ represent
the contributions from the convective term, the viscous stresses, the surface
tension, and the gravity, respectively,
\begin{align}
  C_c &=  \frac{\max_{i,j}|u_{i,j}|}{\Delta x}
        + \frac{\max_{i,j}|v_{i,j}|}{\Delta y}, \\
  C_v &= 2\max\left(\frac{\mu_1}{\rho_1},\frac{\mu_2}{\rho_2}\right)
         \left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right), \\
  C_s &= \sqrt{\frac{\sigma\max_{i,j}|\kappa_{i,j}|}
                    {\max(\rho_1,\rho_2)\min(\Delta x^2,\Delta y^2)}}, \\
  C_g &= \sqrt{\frac{|g_x|}{\Delta x} + \frac{|g_y|}{\Delta y}}.
\end{align}
This CFL condition is discussed in more detail by~\citet{Lervag08}, and it is
based on the condition by \citet{Kang00}.

\section{Axisymmetry}
For axisymmetric flow, the governing equations \eqref{eq:ns1} and
\eqref{eq:ns2} become
\begin{equation}
  \frac 1 r \pd{} r\left(r u\right) + \pd v z = 0,
\end{equation}
\begin{multline}
  \rho\left(\pdt u + u\pd u r + v\pd u z\right) = -\pd p r \\
    + \mu\left(\frac 1 r \pd{} r\left(r \pd u r\right)
          + \pdd u z - \frac u {r^2}\right) + \rho f_r + {f_s}_r,
\end{multline}
and
\begin{multline}
  \rho\left(\pdt v + u\pd v r + v\pd v z\right) = -\pd p z \\
    + \mu\left(\frac 1 r \pd{} r\left(r \pd v r\right)
      + \pdd v z\right) + \rho f_z + {f_s}_z.
\end{multline}
Here $u$ and $v$ are the radial and axial velocity components, $f_r$ and $f_z$
are the radial and axial body-force components, and ${f_s}_r$ and ${f_s}_z$ are
the radial and axial components of the singular surface force, respectively.
The equations are solved as explained in the previous sections.

\section{Discretization of the curvature and the normal vector}
\label{sec:curvature}
As stated in \cref{sec:level-set}, the normal vector \eqref{eq:norm} and the
curvature \eqref{eq:curv} can be calculated from the level-set function as
\begin{align*}
  \vct n &= \frac{\grad\varphi}{|\grad\varphi|}, \\
  \kappa &= \div\left(\frac{\grad\varphi}{|\grad\varphi|}\right).
\end{align*}
They are typically discretized with the standard second-order
central-difference scheme, cf.\ \cite{Kang00,Sethian03,Xu06}.  The normal
vector is calculated at the cell edges, and the curvature is calculated at the
grid nodes.  The curvature is then interpolated to the interface where needed
with linear interpolation, for instance with
\begin{equation}
  \kappa_\Gamma = \frac{|\varphi_{i,j}|\kappa_{i+1,j}
  + |\varphi_{i+1,k}|\kappa_{i,j}}{|\varphi_{i,j}| + |\varphi_{i+1,j}|}.
  \label{eq:curvature-gamma}
\end{equation}

If the level-set method is used to capture non-trivial geometries, then it will
contain kink regions, that is, areas where the gradient of the level-set
function is discontinuous.  \Cref{fig:level-set-function} shows a simple
example of such a kink region for a level-set function in a one-dimensional
domain that captures two interfaces, one on each side of $x_i$.  The kink at
$x_i$ may lead to large errors both for the curvature and the normal vector if
one is not careful.  Errors in the curvature lead to errors in the surface
tension force and in the pressure, which in turn lead to errors in the
interface evolution and in the two-phase flow.  Errors in the normal vector
affect both the calculation of the viscous jump condition and the advection of
the interface.  If the level-set method is used to study for example
coalescence and breakup of drops, these errors may severely affect the
simulations.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [ % Define styles
    axes/.style={thick,gray!150,>=stealth},
    phi/.style={black},
    scale=0.75,
    inner sep=0mm,
    filledcircle/.style={minimum size=3pt,fill=black,circle},
    ]

    % Create grid
    \node[anchor=south] (posphi) at ( 0.00, 2.0 ) {$\varphi>0$};
    \node[anchor=north] (negphi) at ( 0.00,-2.0 ) {$\varphi<0$};
    \node[anchor=east]  (zerphi) at (-0.25, 0.0 ) {$\varphi=0$};
    \node[anchor=north] (xi)     at ( 5.00,-0.1 ) {$x_i$};
    \node[anchor=west]  (x)      at (10.00, 0.0 ) {$x$};
    \draw[<-,axes] (posphi) -- (negphi);
    \draw[->,axes]  (zerphi) -- (x);
    \foreach \x in {0,...,9}
      \draw (\x,2pt) -- (\x,-2pt);

    % Create phi function
    \draw[phi] (1.10,-1.45) -- (5.00, 0.40) node[filledcircle] {}
                            -- (8.90,-1.45) ;
  \end{tikzpicture}
  \caption{A level-set function with a gradient that is discontinuous at
    $x_i$.}
  \label{fig:level-set-function}
\end{figure}

This problem was to our knowledge first described by \citet{Smereka03}, who
increase the numerical smoothing in the curvature discretization to lessen the
effect.  Several non-smearing approaches have subsequently been developed.
\citet{Macklin05} used the level-set method to study tumor growth, and they
present a one-sided direction-difference scheme for the discretization of the
normal vector and the curvature.  Later, \citet{Macklin06} presented an
improved geometry-aware curvature discretization, where the curvature is
calculated based on a local least-squares parametrisation of the interface.

A different approach to avoid the kinks was presented by \citet{Salac08}.  They
used a level-set extraction technique, where an extraction algorithm was used
to reconstruct separate level-set functions for each distinct \emph{body}.  The
term \emph{body} is used here to denote a subset of a given phase or fluid.
For example, in the case of two drops of water colliding in air, the water
drops would make two distinct bodies.  One can also avoid the extraction
algorithm altogether by use of multiple marker functions for different bodies,
see for instance \cite{Coyajee09,Kwakkel12}.  Note, however, that the latter
approach means that the different bodies will not coalesce unless explicit
action is made, cf.~\cite{Kwakkel13}.  Also, both of the approaches mentioned
here fails to handle the problem of kinks from a single body.  That is, there
may still be kink areas due to deformed bodies, for instance bodies with thin
filaments or tails, or bodies shaped like horse shoes.

In the following, we first present the direction-difference
scheme~\cite{Macklin05}.  We then describe the curve-fitting discretization
method and the local level-set extraction method.

\subsection{Direction-difference scheme}
\label{sec:dds}
The direction-difference scheme (DDS) was introduced by \citet{Macklin05}.  It
uses a quality function to ensure that the difference stencils never cross any
kink regions.  The DDS is used in Papers A--C to calculate the normal vectors.
However, in Paper C it is shown that the DDS does not always yield an accurate
approximation of the normal vector.

The basic strategy is to use a combination of central differences and one-sided
differences based on the values of a quality function,
\begin{equation}
  Q(\vct x) = \left| 1 - |\grad\varphi(\vct x)|\right|.
  \label{eq:quality}
\end{equation}
The quality function is approximated with central differences, and is used to
detect the areas where the level-set function differs from the signed-distance
function.  In the following, let $Q_{i,j} \equiv Q(\vct x_{i,j})$ and define
a parameter $\eta>0$.  This threshold parameter is tuned such that the quality
function will detect all the kinks.

The quality function is used to define a direction function,
\begin{equation}
  \vct D(\vct x_{i,j}) = (D_x(\vct x_{i,j}), D_y(\vct x_{i,j})),
\end{equation}
where
\begin{equation}
  D_x(\vct x_{i,j}) = \begin{cases}
    -1 & \text{if $Q_{i-1,j}<\eta$ and $Q_{i+1,j}\geq\eta$,} \\
     1 & \text{if $Q_{i-1,j}\geq\eta$ and $Q_{i+1,j}<\eta$,} \\
     0 & \text{if $Q_{i-1,j}<\eta$ and $Q_{i,j}<\eta$ and $Q_{i+1,j}<\eta$,} \\
     0 & \text{if $Q_{i-1,j}\geq\eta$ and $Q_{i,j}\geq\eta$ and
         $Q_{i+1,j}\geq\eta$,} \\
     4 & \text{otherwise.}
  \end{cases}
  \label{eq:dirfunc}
\end{equation}
$D_y(\vct x_{i,j})$ is defined in a similar manner.  If $D_x(\vct x_{i,j})
+ D_y(\vct x_{i,j}) > 2$, then $\vct D(\vct x_{i,j})$ is chosen as the vector
normal to $\grad\varphi(\vct x_{i,j})$.  It is normalized, and the sign is
chosen such that it points in the direction of the best quality.  See
\cite{Macklin05} for more details.

The DDS is then defined as
\begin{equation}
  \partial_x f_{i,j} = \begin{cases}
    \frac{f_{i,j} - f_{i-1,j}}{\Delta x}    & \text{if $D_x(x_i,y_j) = -1$,} \\
    \frac{f_{i+1,j} - f_{i,j}}{\Delta x}    & \text{if $D_x(x_i,y_j) =  1$,} \\
    \frac{f_{i+1,j} - f_{i-1,j}}{2\Delta x} & \text{if $D_x(x_i,y_j) =  0$,}
  \end{cases}
  \label{eq:Ddifference}
\end{equation}
and similarly for $\partial_y f_{i,j}$, where $f_{i,j}$ is a piecewise smooth
function.  The DDS is equivalent to using central differences in smooth areas
and one-sided differences in areas close to the kinks.

\subsection{Curve-fitting discretization method}
\label{sec:curvature-curve-fitting}
The curve-fitting discretization method (CFDM) was first presented
in~\cite{Lervag11} and is based on the method by \citet{Macklin06}.  The main
idea is to identify kink regions with the quality function \eqref{eq:quality},
and to use a curve parametrisation of the closest interface to calculate the
curvature in regions where the quality function is larger than the threshold
parameter, $\eta$.

The CFDM applied to the curvature or the normal vector at the grid point $\vct
x_{i,j}$ can be summarized as follows.  See also \cref{fig:cfdm}, which shows
an example of the CFDM used at $\vct x_{i,j}$.
\begin{enumerate}
  \item If the quality of the level-set function in the neighbourhood of $\vct
    x_{i,j}$ is good, that is
    \[
      Q(\vct{x}_{n,m}) \le \eta \; \forall (n,m) \in [i-1,i+1]\times[j-1,j+1],
    \]
    then we use a standard discretization.  Otherwise continue to the next
    step.
  \item Locate the closest interface, $\Gamma$.
  \item Find a set of points on the located interface, $\vct x_1, \dots, \vct
    x_n \in \Gamma$.
  \item Create a parametrisation $\vct \gamma(s)$ of the points $\vct x_1,
    \dots, \vct x_n$.
  \item Use the parametrisation $\vct \gamma(s)$ to calculate a local level-set
    function.
  \item Use a standard discretization of the local level-set function to
    calculate the curvature or the normal vector.
\end{enumerate}
\begin{figure}[tbp]
  \centering
  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \draw[thick]   (-2.25, 2.35) arc  ( 80: 44:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \fill[black!2] (-2.25, 2.35) arc  ( 80: 44:9.0cm) -- (-1,-1) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \foreach \i in {1,2,3}
        \foreach \j in {1,2,3}
          \fill (\i,\j) circle(1pt);

      \fill (2.83, 2.83) circle(2pt) node[right] {$\vct x_3$};
    \end{tikzpicture}
    \caption{First locate the closest interface, here represented with $\vct
      x_3$.}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \draw[thick]   (-2.25, 2.35) arc  ( 80: 44:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \fill[black!2] (-2.25, 2.35) arc  ( 80: 44:9.0cm) -- (-1,-1) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \foreach \i in {1,2,3}
        \foreach \j in {1,2,3}
          \fill[gray] (\i,\j) circle(1pt);

      \fill (1.52, 5.00) circle(2pt) node[right] {$\vct x_1$};
      \fill (2.00, 4.04) circle(2pt) node[right] {$\vct x_2$};
      \fill (2.83, 2.83) circle(2pt) node[right] {$\vct x_3$};
      \fill (3.63, 2.00) circle(2pt) node[right] {$\vct x_4$};
      \fill (4.92, 1.00) circle(2pt) node[right] {$\vct x_5$};
    \end{tikzpicture}
    \caption{Then find a set of points along the closest interface.}
  \end{subfigure}

  \vspace{1em}

  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \draw[thick]   (-2.25, 2.35) arc  ( 80: 44:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \fill[black!2] (-2.25, 2.35) arc  ( 80: 44:9.0cm) -- (-1,-1) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \foreach \i in {1,2,3}
        \foreach \j in {1,2,3}
          \fill[gray] (\i,\j) circle(1pt);

      \fill (1.52, 5.00) circle(2pt);
      \fill (2.00, 4.04) circle(2pt);
      \fill (2.83, 2.83) circle(2pt);
      \fill (3.63, 2.00) circle(2pt);
      \fill (4.92, 1.00) circle(2pt);
      \draw[very thick, dashed]    (1.52, 5.00) -- (2.00, 4.04)
                                -- (2.83, 2.83) -- (3.63, 2.00)
                                -- (4.92, 1.00);
    \end{tikzpicture}
    \caption{Construct a curve parametrisation from the points $\vct
      x_1,\dots,\vct x_5$.}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \draw[thick]   (-2.25, 2.35) arc  ( 80: 44:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \fill[black!2] (-2.25, 2.35) arc  ( 80: 44:9.0cm) -- (-1,-1) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \foreach \i in {1,2,3}
        \foreach \j in {1,2,3}
          \fill (\i,\j) circle(2pt);

      \draw[very thick, dashed]    (1.52, 5.00) -- (2.00, 4.04)
                                -- (2.83, 2.83) -- (3.63, 2.00)
                                -- (4.92, 1.00);
    \end{tikzpicture}
    \caption{Calculate a new local level-set function at the grid points around
      and including $\vct x_{i,j}$.}
  \end{subfigure}
  \caption{Example of the CFDM at a grid point $\vct x_{i,j}$.  First five
    points are found on the interface closest to $\vct x_{i,j}$.  Then a curve
    parametrisation (dashed line) is calculated, and the parametrisation is
    used to calculate a new local level-set function at the grid points
    surrounding and including $\vct x_{i,j}$.}
  \label{fig:cfdm}
\end{figure}

\subsection{Local level-set extraction method}
\label{sec:lolex}
The local level-set extraction (LOLEX) method is based on the method presented
by \citet{Salac08}, here called the SLM.  It was found that the latter method
was insufficient, because it did not treat all the kink problems.  The LOLEX
method is therefore a further development of the SLM, in that it handles the
kink regions in a more general manner.

The LOLEX method applied to the curvature or the normal vector at a grid point
$\vct x_{i,j}$ is summarized by the following algorithm.  The algorithm is
presented with 2D notation for clarity, and extends easily to 3D.  See also
\cref{fig:lolex}, which gives a simple example of the procedure.
\begin{enumerate}
  \item If the quality in the neighbourhood of $\vct x_{i,j}$ is good, that is
    \[
      Q(\vct{x}_{n,m}) \le \eta \; \forall (n,m) \in [i-1,i+1]\times[j-1,j+1],
    \]
    then we use a standard discretization.  Otherwise continue to the next
    step.
  \item Copy a small, local square centred around $\vct x_{i,j}$ from the
    level-set function $\varphi$ into a local array $\varphi_\text{loc}$.
  \item Identify and enumerate all the \emph{bodies} in the local array
    $\varphi_\text{loc}$.  A \emph{body} is here defined as a set of
    neighbouring points where $\varphi_\text{loc}<0$, see \cref{fig:lolex}.
    There will be $n\geq 0$ bodies in any given $\varphi_\text{loc}$ array.
  \item If no body is identified, that is, if $n=0$, then use a standard
    discretization with the global level-set function $\varphi$.  Otherwise
    continue to the next step.
  \item For each body $n$, extract the relevant parts of $\varphi_\text{loc}$
    into an array $\varphi_\text{loc}^n$.  If necessary, extrapolate values to
    ghost cells.
  \item For each body, $n$, reinitialize $\varphi_\text{loc}^n$.
  \item At this step, all the bodies in the local grid have their own
    local level-set functions that have been reinitialized to proper
    signed-distance functions.  Due to the separation of the bodies, there
    are no longer any kinks.
  \item Use the standard discretization of the curvature and the normal
    vector at the local level-set function that represents the body that is
    closest to $\vct x_{i,j}$.
\end{enumerate}
\begin{figure}[tbp]
  \centering
  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \draw[thick]   (-2.25, 2.35) arc  ( 80: 44:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \fill[black!2] (-2.25, 2.35) arc  ( 80: 44:9.0cm) -- (-1,-1) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \fill (2,2) circle(1pt);

      \node at (0.5,0.5) {1};
      \node at (4.5,4.5) {2};
    \end{tikzpicture}
    \caption{First enumerate the bodies.  In this case there are $n=2$ bodies.
      Then the bodies are extracted into separate level-set functions.}
  \end{subfigure}
  \begin{subfigure}[t]{0.49\textwidth}
    \centering
    \begin{tikzpicture}[scale=0.9]
      \clip (-0.25,-0.25) rectangle (5.5,5.25);
      \draw[thick]   ( 0.80, 8.25) arc  (182:257:9.0cm);
      \fill[black!2] ( 0.80, 8.25) arc  (182:257:9.0cm) -- ( 6, 6) -- cycle;
      \draw[dotted]  (-0.25,-0.25) grid (5.25,5.25);

      \node[below right] at (2,2) {$\vct x_{i,j}$};
      \foreach \i in {1,2,3}
        \foreach \j in {1,2,3}
          \fill (\i,\j) circle(1pt);
    \end{tikzpicture}
    \caption{Reinitialize the separated level-set functions, then use the
      function that represents the closest interface.}
  \end{subfigure}
  \caption{An example of the LOLEX method at a grid point $\vct x_{i,j}$.
    First the bodies are identified and enumerated.  Then they are extracted
    into separate level-set functions, which are reinitialized.  Finally, the
    curvature or normal vector is calculated based on the level-set function
    for the closest body.}
  \label{fig:lolex}
\end{figure}

\section{Summary}
In this chapter we have described the numerical methods that have been used to
solve the Navier-Stokes equations for two-phase flow \eqref{eq:ns1}
and~\eqref{eq:ns2}.

We first gave a brief introduction to the level-set method, which is used to
capture the interface.  We then presented the spatial and temporal
discretization methods, including a brief overview of the projection methods
that were used to decouple the pressure from the Navier-Stokes equations.  In
the final section, we presented the new methods for calculating the curvature
and normal vector.  These are the curve-fitting discretization scheme (CFDM)
and the local level-set extraction (LOLEX) method.

% Fakesection: Quote
\begin{savequote}[6cm]
  ``If we want to solve a problem that we have never solved before, we must
  leave the door to the unknown ajar.''
  \qauthor{--- Richard P. Feynman (1918--1988)}
\end{savequote}
\chapter{The diffuse-domain approach}
\label{chap:diffuse-domain}
In the previous chapters, we have considered the modelling of two-phase flows,
and in particular methods for calculating the curvature and normal vector with
the level-set method in a reliable manner.  In this chapter, we consider
a different problem of a more general nature:  How to solve partial
differential equations (PDEs) in complex domains.  In particular, we consider
an extension of a diffuse-domain method (DDM) by a high-order correction term
that gives increased accuracy with respect to interface-width refinements.

We begin with a short introduction to the diffuse-domain approach.  We then
outline how it can be used to derive a DDM for the steady reaction-diffusion
equation with Neumann boundary conditions.  Next, we continue with a brief
introduction to the method of matched asymptotic expansions.  Finally, we
introduce the high-order correction term derived in Paper~E and show that the
resulting DDM converge with second order in the diffuse-interface width to the
original problem.  The analysis also shows that the correction term is not
necessary for second-order convergence.

\section{Introduction}
There exist several methods for solving PDEs in complex domains.  Most of them
have in common that they require tools or methods that are not frequently
available in standard finite-element or finite-difference software packages.
Examples of such methods include the immersed-interface method
\cite{LeVeque94}, the matched interface and boundary method \cite{Zhou06}, the
extended and composite finite-element method \cite{Dolbow09}, embedded boundary
methods \cite{Johansen98}, cut-cell methods \cite{Ji06}, and ghost-fluid
methods \cite{Fedkiw99}.  A different approach, known as the fictitious domain
method~\cite{Glowinski94,Glowinski96} or the domain imbedding
method~\cite{Buzbee71}, either augments the original system with equations for
Lagrange multipliers to enforce the boundary conditions, or use the penalty
method to enforce the boundary conditions weakly.  For a more complete list of
references, see Paper~E.

The DDM is an alternative method for solving PDEs in complex domains.  The main
idea is to use an implicit representation of the boundary, where the sharp
boundary is replaced by a diffuse layer.  The PDEs are then reformulated on
a larger, regular domain, and the boundary conditions are incorporated via
source terms in the diffuse layer.  When the thickness of the diffuse layer is
reduced, these source terms tend towards singular source terms.  The resulting
PDEs can then be solved with the use of standard tools and methods.

The diffuse-domain approach was first introduced by \citet{Kockelkoren03} to
study diffusion inside a cell with homogeneous Neumann boundary conditions at
the cell boundary.  It was later used by \citet{Li09} to develop a DDM for
solving PDEs in complex evolving domains with Dirichlet, Neumann and Robin
boundary conditions, which is hereafter called the DDM1.  The DDM1 has been
used by \citet{Teigen09-b}, who modelled bulk-surface coupling of material
quantities on a deformable interface.  It was also used by \citet{Aland10} to
simulate incompressible two-phase flows in complex domains in 2D and 3D, and by
\citet{Teigen11} to study two-phase flows with soluble surfactants.

An analysis of the error behaviour of the diffuse-domain approach was done by
\citet{Franz12} for a diffuse-domain approximation of an elliptic problem with
Dirichlet boundary conditions.  They considered the infinity norm of the
difference of the approximated solution and the exact solution, and their
analysis shows that the approximation quality is of order one in the interface
width.

In Paper~E, we present the DDM2, which is an extension of the DDM1 by
a high-order correction term.  The DDM2 is derived for elliptic problems with
Neumann and Robin boundary conditions, and it is shown to be asymptotically
second-order accurate in the interface width.  However, the analysis in Paper~E
is somewhat lacking in that it assumes that the DDM1 is only first-order
accurate.  In the following sections, we extend the analysis of Paper~E and
show that the DDM1 is also second-order accurate.  The analysis is shown for
the steady reaction-diffusion equation with Neumann boundary conditions, but
the same technique also applies for the corresponding Robin problem.

\section{The DDM for a Neumann problem}
Consider the steady reaction-diffusion equation with Neumann boundary
conditions,
\begin{equation}
  \begin{alignedat}{2}
    \lapl u - u  &= f        & \txin D, \\
    \ndot\grad u &= g \qquad & \txon \partial D,
  \end{alignedat}
  \label{eq:ddan}
\end{equation}
where $f$ and $g$ are given.  Let $\chi_D$ be the characteristic function of
$D$,
\begin{equation}
  \chi_D = \begin{cases}
    1 & \text{if $x\in D$,} \\
    0 & \text{if $x\notin D$.}
  \end{cases}
\end{equation}
The main idea with the diffuse-domain approach is to extend the original
equation \eqref{eq:ddan} into a larger and regular domain $\Omega\supset D$, as
depicted in \cref{fig:ddadomain}.  The extension can be written as
\begin{equation}
  \div(\chi_D\grad u) - \chi_D u + \text{BC} = \chi_D f,
  \label{eq:ddm}
\end{equation}
where BC is a singular source term that represents the physical boundary
condition on $\partial\Omega$.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [
    scale=0.8,
    wall/.style={
      decoration={border,angle=45,segment length=4},
      postaction={decorate,draw}},
    ]
    \draw[wall] (0,0) rectangle(9,5);
    \draw (2,1) .. controls (1,1) and (1,4) .. (3,3.5)
                .. controls (5,3) and (4,5) .. (6,4)
                   node[above right] {$\partial D$}
                .. controls (8,3) and (8,1.5) .. (7,1.5)
                .. controls (3,1.5) and (3,1) .. (2,1);
    \node at (2.3,2.3) {$D$};
    \node at (7.6,0.8) {$\Omega$};
    \node at (4.8,2.4) {$\chi_D=1$};
    \node at (5.4,0.5) {$\chi_D=0$};
  \end{tikzpicture}
  \caption{A regular domain $\Omega$ that contains a complex domain $D$.}
  \label{fig:ddadomain}
\end{figure}

The characteristic function is typically approximated by the phase-field
function,
\begin{equation}
  \chi_D \simeq \phi(\vct x,t) = \frac{1}{2} \left( 1
  - \tanh \left( \frac{3r(\vct x,t)}{\epsilon} \right) \right),
  \label{eq:characteristic}
\end{equation}
where $\epsilon$ is the interface width and $r(\vct x,t)$ is the
signed-distance function with respect to the boundary $\partial D$, which is
taken to be negative inside $D$.

The main difficulty with the diffuse-domain approach is the derivation of
approximations for the boundary condition term BC.  \citet{Li09} give four
approximations that are shown to converge asymptotically with first order in
$\epsilon$ to the original equation when $\epsilon$ is decreased.  In the
following we consider the approximation
\begin{equation}
  \text{BC} \simeq |\grad\phi| g.
  \label{eq:bc1}
\end{equation}
If we combine the above approximations \eqref{eq:characteristic} and
\eqref{eq:bc1}, we get a DDM1 equation for \eqref{eq:ddan},
\begin{equation}
  \div\left(\phi\grad u\right) - \phi u +  |\grad\phi| g = \phi f.
  \label{eq:ddm1}
\end{equation}

\section{The method of matched asymptotic expansions}
The following is a brief introduction to the method of matched asymptotic
expansions, which is used to show that a given diffuse-domain approximation
converges to the original problem when the interface width is decreased.  More
details can be found in Paper~E and in \cite{Pego88}.

Let $u$ be some diffuse-domain variable.  The asymptotic convergence of a given
diffuse-domain approximation can be shown through expansions of the
diffuse-domain variables in powers of the interface thickness $\epsilon$ in
regions close to and far from the interface.  For example, the expansions of
$u$ are
\begin{align}
  u(\vct x)        &= \sum_{k=0}^\infty \epsilon^k      u^{(k)}(\vct x), \\
  \hat u(z,\vct s) &= \sum_{k=0}^\infty \epsilon^k \hat u^{(k)}(z,\vct s),
\end{align}
where $u(\vct x)$ and $\hat u(\vct s,z)$ denote the outer and inner expansions,
respectively.  Here $z$ is a stretched variable,
\begin{equation}
  z = \frac{r(\vct x)}{\epsilon},
\end{equation}
where $r$ is the signed distance from the point $\vct x$ to $\partial D$ and is
taken to be negative inside $D$.  Further, $z$ and $\vct s$ form a local
coordinate system such that
\begin{equation}
  \vct x(\vct s,z) = \vct X(\vct s) + \epsilon z\vct n(\vct s),
\end{equation}
where $\vct X(\vct s)$ is a parametrisation of the interface, $\vct n(\vct s)$
is the interface normal vector, and $z$ is a stretched variable.

When the inner and outer expansions are found, they are matched in a region
where both solutions are valid and where $\epsilon z = \bigo 1$, see
\cref{fig:regions}.  The outer solution is then evaluated in the inner
coordinates, which leads to a set of matching conditions that must hold when
$\epsilon\to 0$.  If we consider $\epsilon$ to be fixed and let $z \to
\pm\infty$, we get the following asymptotic matching conditions:
\begin{equation}
  \label{eq:match1}
  \zlimpm \hat u^{(0)}(z,\vct s) = u^{(0)}(\vct s),
\end{equation}
and as $z \to \pm\infty$,
\begin{align}
  \label{eq:match2}
  \hat u^{(1)}(z,\vct s) &= u^{(1)}(\vct s)
    + z\ndot\grad u^{(0)}(\vct s) + \smallo 1, \\
  \begin{split}
    \hat u^{(2)}(z,\vct s) &= u^{(2)}(\vct s)
      + z\ndot\grad u^{(1)}(\vct s) \\
      &\quad + \frac{z^2}{2} (\ndot\grad)\grad u^{(0)}(\vct s)\cdot\vct n
      + \smallo 1.
    \label{eq:match3}
  \end{split}
\end{align}
\begin{figure}[b!p]
  \centering
  \begin{tikzpicture}
    [
    yscale=0.8,
    interface/.style={thick},
    inner/.style={fill=gray,dotted,fill opacity=0.2},
    outer/.style={fill=gray,dashed,fill opacity=0.3},
    labels/.style={above right, font=\small},
    wall/.style={
      decoration={border,angle=45,segment length=4},
      postaction={decorate,draw}},
    ]

    % Domains
    \begin{scope}[scale=0.8]
      \draw[wall] (0,0) rectangle(9,5);
      \draw (2,1) .. controls (1,1) and (1,4) .. (3,3.5)
                  .. controls (5,3) and (4,5) .. (6,4)
                  .. controls (8,3) and (8,1.5) .. (7,1.5) node (g1) {}
                  .. controls (3,1.5) and (3,1) .. (2,1);
      \node at (7.4,0.4) {$\Omega$};
      \node at (2.0,2.3) {$D$};
      \coordinate (a) at (2.4,1.6);
      \coordinate (b) at (3.4,1.6);
      \coordinate (c) at (2.4,0.6);
      \coordinate (d) at (3.4,0.6);
      \draw[very thin] (c) rectangle (b);
    \end{scope}

    % Closeup
    \begin{scope}[xshift=2.5cm,yshift=4.4cm]
      % Zoom lines
      \draw[very thin] (a) -- (0, 3);
      \draw[very thin] (b) -- (8, 3);
      \draw[very thin] (c) -- (0,-3);
      \draw[very thin] (d) -- (8,-3);
      \fill[white] (0,-3) rectangle (8,3);
      \draw[very thin, densely dotted] (a) -- (0, 3);
      \draw[very thin, densely dotted] (b) -- (8, 3);

      % Interface
      \draw[interface] (0,0)
        .. controls (2, 0.5) and (3, 0.5) .. (4,0)
        .. controls (5,-0.5) and (6,-1.0) .. (8,0);

      % Inner region
      \draw[inner] (0,2.0)
        .. controls (2,2.5) and (3,2.5) .. (4,2.0)
        .. controls (5,1.5) and (6,1.0) .. (8,2.0) -- (8,-2.0)
        .. controls (6,-3.0) and (5,-2.5) .. (4,-2.0)
        .. controls (3,-1.5) and (2,-1.5) .. (0,-2.0) -- cycle;

      % Outer region
      \draw[outer] (0,3.0) -- (8,3.0) -- (8,1.0)
        .. controls (6,0.0) and (5,0.5) .. (4,1.0)
        .. controls (3,1.5) and (2,1.5) .. (0,1.0) -- cycle;
      \draw[outer] (0,-3.0) -- (8,-3.0) -- (8,-1.0)
        .. controls (6,-2.0) and (5,-1.5) .. (4,-1.0)
        .. controls (3,-0.5) and (2,-0.5) .. (0,-1.0) -- cycle;

      % Labels
      \node[labels] at (0.1,2.25) {Outer region};
      \node[labels] at (0.1,1.25) {Overlapping region};
      \node[labels] at (0.1,0.25) {Inner region};
      \draw[decorate,decoration=brace] (8.1, 3.0) --
        node[right=0.5em] {$D$} (8.1, 0.1);
      \draw[decorate,decoration=brace] (8.1,-0.1) --
        node[right=0.5em] {$\Omega$} (8.1,-3.0);
      \node[right=0.5em] at (8.1,0) {$\partial D$};
    \end{scope}
  \end{tikzpicture}
  \caption{A sketch of the regions used for the matched asymptotic expansions.
    The inner region is marked with a light gray color and the outer region
    with a slightly darker gray color.  The overlapping region is marked with
    the darkest gray color.}
  \label{fig:regions}
\end{figure}

We remark that the inner expansion is used to obtain the boundary condition on
$\partial D$, and that the outer solution is used to obtain the sharp-interface
equation inside the physical domain $D$.

To show that a given DDM approximation converges with second order, one must
show that the order-one term of the outer solution of the DDM equation is zero
in $D$.  As an example, we consider the outer solution of the DDM1 equation
\eqref{eq:ddm1}, which is
\begin{equation}
  \begin{split}
    \lapl u^{(0)} - u^{(0)} &= f, \\
    \lapl u^{(1)} - u^{(1)} &= 0, \\
    \lapl u^{(k)} - u^{(k)} &= 0,\qquad k = 2,3,\dots.
  \end{split}
  \label{eq:ddm_outer}
\end{equation}
For the solution to be asymptotically second order, that is $u = u^{(0)}
+ \bigo{\epsilon^2}$, we must have that $u^{(0)}$ satisfies the original
problem \eqref{eq:ddan} and that $u^{(1)}=0$.  Thus the inner expansion must
yield a boundary condition for $u^{(1)}$ to enforce $u^{(1)}=0$.

\section{Asymptotic analysis of the DDM1 and the DDM2}
\label{sec:DDM2}
In Paper~E we present the DDM2, which extends the DDM1 \eqref{eq:ddm1} with
a high-order correction term,
\begin{equation}
  \div\left(\phi\grad u\right) - \phi u +  |\grad\phi| g
    + r|\grad\phi|\left(f - \kappa  g - \lapls u + u\right) = \phi f.
  \label{eq:ddm2}
\end{equation}
Here $r|\grad\phi|\left(f - \kappa  g - \lapls u + u\right)$ is the
correction term, $\kappa$ is the curvature of the boundary $\partial D$, and
$\lapls u$ is the surface Laplacian of $u$, which can be defined as
\begin{equation}
  \lapls u \equiv
  \left( I - \vct n\otimes\vct n \right)
        \div \left( I - \vct n\otimes\vct n \right)\grad u,
\end{equation}
where $I$ is the identity matrix and $\vct n$ is the normal vector.  The
curvature can be calculated from the phase-field function
\eqref{eq:characteristic} as
\begin{equation}
  \kappa = -\div\frac{\grad\phi}{|\grad\phi|}.
\end{equation}

In the following, we use the method of matched asymptotic expansions to show
that the DDM1 \eqref{eq:ddm1} and the DDM2 \eqref{eq:ddm2} are both
second-order approximations of \eqref{eq:ddan} in $\epsilon$.  First, it is
easy to see that the outer expansions of both approximations are given by
\eqref{eq:ddm_outer}.

Next, we consider the inner expansion of the DDM2 \eqref{eq:ddm2}, which is
\begin{multline}
  \frac{1}{\epsilon^2}\left(\phi\hat u_z\right)_z
  + \frac \kappa \epsilon \phi\hat u_z
  + \phi\lapls \hat u
  - \phi\hat u \\
  - \frac 1 \epsilon \phi_z g
  - z\phi_z\left(\hat u + \hat f - \kappa g - \lapls\hat u\right)
  = \phi \hat f.
  \label{eq:ddm_inner}
\end{multline}
We expand $\hat u(z,\vct s)$ in powers of $\epsilon$ and collect the lowest
order terms,
\begin{equation}
  \left(\phi\hat u_z^{(0)}\right)_z = 0.
\end{equation}
If we integrate over all $z$, we get that $\hat u_z^{(0)} = 0$.  The next order
terms of~\eqref{eq:ddm_inner} then give
\begin{equation}
  \left(\phi\hat u_z^{(1)}\right)_z = \phi_z g,
\end{equation}
and again we integrate, which gives that
\begin{equation}
  \phi\hat u_z^{(1)} = \phi g + C
  \label{eq:uz1}
\end{equation}
where the constant $C$ must be zero, since $\lim_{z\to\infty}\phi(z) = 0$.  Now
consider the limit $z\to-\infty$ and use the matching condition
\eqref{eq:match2} to get
\begin{equation}
  \ndot\grad u^{(0)} = g.
\end{equation}
Thus $u^{(0)}$ satisfies the original problem at least to first order in
$\epsilon$.  This shows that both DDM1 and DDM2 are first-order approximations
of the sharp-interface problem.

To obtain the result for the next order, we need to apply the derivative of
the matching condition \eqref{eq:match3},
\begin{equation}
  \hat u_z^{(2)} = \ndot\grad u^{(1)}
      + z (\ndot\grad)\grad u^{(0)}\cdot\vct n.
  \label{eq:match3_z}
\end{equation}
Further, we use that $u^{(0)}$ satisfies
\begin{equation}
  \lapl u^{(0)} - u^{(0)} = f^{(0)},
\end{equation}
and that the Laplacian may be decomposed as
\begin{equation}
  \lapl u = (\ndot\grad)\grad u\cdot\vct n + \kappa\ndot\grad u + \lapls u,
\end{equation}
to get
\begin{equation}
  (\ndot\grad)\grad u^{(0)}\cdot\vct n
    = u^{(0)} + f^{(0)} - \kappa g - \lapls u^{(0)}.
  \label{eq:ngradun}
\end{equation}
Now insert \eqref{eq:ngradun} into \eqref{eq:match3_z} and use the matching
condition \eqref{eq:match1} to obtain a modified matching condition,
\begin{equation}
  \hat u_z^{(2)} - z\left(\hat u^{(0)} + \hat f^{(0)}
        - \kappa g - \lapls \hat u^{(0)}\right)
  = \ndot\grad u^{(1)}.
  \label{eq:match3_z2}
\end{equation}

We are now ready to consider the zeroth order terms,
\begin{multline}
  \label{eq:o1}
  \left(\phi\hat u^{(2)}_z\right)_z
  + \phi\kappa\hat u_z^{(1)}
  + \phi\lapls \hat u^{(0)}
  - \phi\hat u^{(0)} \\
  - z\phi_z\left(\hat u^{(0)}
      + \hat f^{(0)} - \kappa g - \lapls\hat u^{(0)}\right)
  = \phi \hat f^{(0)}.
\end{multline}
The modified matching condition \eqref{eq:match3_z2} motivates that we subtract
and add the term
\[
  \bigg(z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z
\]
to \eqref{eq:o1}, which gives
\begin{multline}
  \bigg(\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z \\
  + \bigg(z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z \\
  + \phi\kappa\hat u_z^{(1)}
  + \phi\lapls \hat u^{(0)}
  - \phi\hat u^{(0)} \\
  - z\phi_z\left(\hat u^{(0)}
      + \hat f^{(0)} - \kappa g - \lapls\hat u^{(0)}\right)
  = \phi \hat f^{(0)}.
\end{multline}
We expand the terms and use \eqref{eq:uz1},
\begin{multline}
  \bigg(\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z \\
  + \cancel{\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)} \\
  + z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)_z \\
  + \cancel{\phi\kappa g
  + \phi\lapls \hat u^{(0)}
  - \phi\hat u^{(0)}}
  = \cancel{\phi \hat f^{(0)}},
  \label{eq:o2}
\end{multline}
or
\begin{multline}
  \bigg(\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z \\
  + z\phi\left(\cancel{\hat u_z^{(0)}} + \hat f_z^{(0)}
      - \cancel{(\kappa g)_z}
      - \cancel{\lapls \hat u_z^{(0)}}\right) = 0.
\end{multline}
If we assume that $\hat f_z^{(0)}$, we get
\begin{equation}
  \bigg(\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z = 0.
  \label{eq:o3}
\end{equation}
Finally, we integrate the left-hand side and take the limit,
\begin{align}
  \nonumber
  \eint{\bigg(\phi\hat u_z^{(2)}
    &- z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
     - \lapls \hat u^{(0)}\right)\bigg)_z} \\
  \nonumber
  &= \ejmp{\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)} \\
  \nonumber
  &= -\zlimm{\left(\hat u_z^{(2)}
    - z\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\right)} \\
  &= -\ndot\grad u^{(1)},
\end{align}
thus
\begin{equation}
  \ndot\grad u^{(1)} = 0.
  \label{eq:ngradu1}
\end{equation}
Combined with \eqref{eq:ddm_outer}, this shows that $u^{(1)}=0$, and so DDM2
converges asymptotically with second order to the original problem.

The analysis above also holds for DDM1, except instead of \eqref{eq:o3} we get
\begin{multline}
  \bigg(\phi\hat u_z^{(2)}
    - z\phi\left(\hat u^{(0)} + \hat f^{(0)} - \kappa g
    - \lapls \hat u^{(0)}\right)\bigg)_z \\
  = - z\phi_z \left(\hat u^{(0)} + \hat
        f^{(0)} - \kappa g - \lapls \hat u^{(0)}\right)
  = - z\phi_z D,
  \label{eq:o4}
\end{multline}
where $D$ is independent of $z$.  Now we use that
\begin{equation}
  \phi_z = -(3 \sech^2 3z)/2,
\end{equation}
which follows from the definition of the phase-field function
\eqref{eq:characteristic}.  We integrate the right-hand side, which gives
\begin{equation}
  D\eint{z\phi_z} = -D\frac 3 2\eint{z\sech^2 3z} = 0.
\end{equation}
Thus DDM1 is also second order in $\epsilon$.

The difference between the DDM1 and the DDM2 is therefore that the correction
term with the DDM2 directly cancels the term on the right-hand side in
\eqref{eq:o4}.  This should give an increase of accuracy, but the convergence
order remains the same.

The analysis for the corresponding Robin problem is essentially the same as the
above.  In Paper~E, the DDM1 and DDM2 are compared for several elliptic
problems with both Neumann and Robin boundary conditions.  The results of
Paper~E show that the correction term in the DDM2 leads to an increase of
accuracy and that both DDM1 and DDM2 converge with second-order accuracy.

\section{Summary}
In this chapter, we have given a brief introduction to the diffuse-domain
method (DDM).  We considered a steady reaction-diffusion equation with Neumann
boundary conditions \eqref{eq:ddan} and two DDM approximations:  DDM1
\eqref{eq:ddm1} and DDM2 \eqref{eq:ddm2}.  The DDM2 is an extension of DDM1 by
a high-order correction term, and was first derived in Paper~E.

Next, we gave an outline of the method of matched asymptotic expansions, and we
used it to show that both the DDM1 and the DDM2 converged asymptotically with
second order in the diffuse-interface width to the original equation
\eqref{eq:ddan}.  The analysis shows that the correction term in the DDM2 leads
to a cancellation in the asymptotic expansions.  By doing the integration, we
see that this cancellation is not necessary for obtaining the second order
convergence.

% Fakesection: Quote
\begin{savequote}[7cm]
  ``Count what is countable, measure what is measurable, and what is not
  measurable, make measurable.''
  \qauthor{--- Galileo Galilei (1564--1642)}
\end{savequote}
\chapter{Summary of contributions}
\label{chap:contributions}
This chapter presents summaries of the papers that constitute parts of this
thesis.  Each summary gives a brief discussion of the results of each paper,
and the contribution of the author is highlighted for each paper.

\section[Paper A]{Paper A:  Calculation of interface curvature with the
  level-set method}
Karl Yngve Lervåg.  Published in \emph{MekIT'11 - 6th National Conference on
  Computational Mechanics, Trondheim}, 2011.  ISBN: 978-82-519-2798-7.

In this paper I address a problem with the calculation of the interface
curvature with the level-set method, cf.\ \cref{sec:curvature}.  The curvature
can be calculated from the level-set function, $\phi$, as
\begin{equation}
  \kappa = \div\vct n = \div\frac{\grad\phi}{|\grad\phi|}.
\end{equation}
It is typically discretized by standard methods such as the second-order
central-difference scheme (CD-2), and interpolated to the interface where
needed \cite{Osher88,Kang00,Osher03}.  However, the level-set function as
a signed-distance function will tend to have kinks where its gradient is
discontinuous.  The standard methods may lead to large errors in the curvature
close to these regions, which in turn may lead to errors in the surface tension
force.

The main contribution of this paper is a new curve-fitting discretization
method (CFDM) for the curvature (see \cref{sec:curvature-curve-fitting}).  The
method is based on the approach developed by \citet{Macklin06}.  It differs in
that it uses a cubic Hermite spline parametrisation of the interface, and that
the curvature values are calculated on the grid and then interpolated to the
interface as opposed to using a localised grid centred at the interface as in
\cite{Macklin06}.

\begin{figure}[tbp]
  \centering
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \includegraphics[width=\textwidth]{simple_old}
    \caption{CD-2}
    \label{fig:A1a}
  \end{subfigure}
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \includegraphics[width=\textwidth]{simple_new}
    \caption{CFDM}
    \label{fig:A1b}
  \end{subfigure}
  \caption{A comparison of curvature calculations between standard
    discretization (CD-2) and the improved method (CFDM).  The standard
    discretization leads to large errors in the curvatures in areas that are
    close to two interfaces.}
  \label{fig:A1}
\end{figure}

The CFDM is tested and compared with the CD-2 for two test cases, and it is
shown to yield better results in both cases.  \Cref{fig:A1} shows one of these
results, where the calculated curvature values are compared.  In the example,
a cylindrical drop impacts on a liquid film.  The figure shows that the CD-2
leads to large errors in the curvature calculations in the kink regions, that
is, the red and dark blue regions in \cref{fig:A1a} near the liquid film.
These errors are not present with the new method, cf.\ \cref{fig:A1b}.

\paragraph{My contribution:}  I developed the method and implemented it into
our in-house finite-difference code for two-phase flow based on the methods of
\cref{sec:deltaT,sec:t-discr,sec:x-discr,sec:projection-method,sec:curvature}.
I ran the numerical simulations.  I wrote the paper and presented the work at
the conference.

\section[Paper B]{Paper B:  Curvature calculations for the level-set method}
Karl Yngve Lervåg and Åsmund Ervik.  Published in \emph{ENUMATH 2011}
proceedings volume, Springer, 2013.  ISBN: 978-3642331336.

This paper is a continuation of Paper A.  The main contribution in this paper
is a comparison of different methods for calculating the curvature in a robust
manner with the level-set method in the kink regions.  In particular, the
CFDM\footnote{The method is called LM in the paper.  Here CFDM is used, in
  order to be consistent with the rest of the thesis.} that was presented in
Paper A is compared with Macklin and Lowengrub's method (MLM)~\cite{Macklin06}.
In addition, the method is compared with the second-order central-difference
scheme (CD-2) and the more recent method presented by \citet{Salac08}, here
called Salac and Lu's method (SLM).

The main result in the paper is shown in \cref{fig:B1}, which shows
a comparison of the methods for a case where two drops collide in a 2D shear
flow.  In particular, it shows snapshots of the evolution of the interfaces and
the curvature at times $t=\SI{2.30}{s}$, $t=\SI{2.75}{s}$, and
$t=\SI{3.10}{s}$.  The results show that all of the improved methods, that is
MLM, CFDM, and SLM, handle the kink region in a more reliable manner than CD-2.
The reason that the result with MLM differs from those with CFDM and SLM might
be that it uses a localised grid centred at the interface to calculate the
curvature, which means that it does not need to use interpolation of the
curvature from the grid to the interface.  Note that the difference is mainly
that the MLM results in slightly earlier coalescence in the given case.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [
    time/.style={fill=white,text width=0.24\textwidth,inner sep=1pt},
    caption/.style={font=\small},
    ]

    \node (standard1) at (0,0)
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_nolc_230}};
    \node (standard2) [below=-0.15cm of standard1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_nolc_275}};
    \node (standard3) [below=-0.15cm of standard2]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_nolc_310}};
    \node (macklin1) [right=-0.15cm of standard1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_macklin_230}};
    \node (macklin2) [below=-0.15cm of macklin1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_macklin_275}};
    \node (macklin3) [below=-0.15cm of macklin2]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_macklin_310}};
    \node (lervag1) [right=-0.15cm of macklin1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_lc_230}};
    \node (lervag2) [below=-0.15cm of lervag1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_lc_275}};
    \node (lervag3) [below=-0.15cm of lervag2]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_lc_310}};
    \node (salac1) [right=-0.15cm of lervag1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_salac_230}};
    \node (salac2) [below=-0.15cm of salac1]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_salac_275}};
    \node (salac3) [below=-0.15cm of salac2]
      {\includegraphics[width=0.24\textwidth]{grey_ycf84_salac_310}};
    \node[time,above=-10pt of standard1] {$t=\SI{2.30}{s}$};
    \node[time,above=-10pt of standard2] {$t=\SI{2.75}{s}$};
    \node[time,above=-10pt of standard3] {$t=\SI{3.10}{s}$};

    % Captions
    \node[caption,below=3pt of standard3] {{\bf(a)} CD-2};
    \node[caption,below=3pt of macklin3]  {{\bf(b)} MLM};
    \node[caption,below=3pt of lervag3]   {{\bf(c)} CFDM};
    \node[caption,below=3pt of salac3]    {{\bf(d)} SLM};

    % Legend
    \node (legend) [right=0.0cm of salac1.south east]
      {\includegraphics[width=1.8cm]{grey_ycf84_legend}};
    \node [above=-0.3cm of legend] {$\kappa\ [\si{\per\metre}]$};

  \end{tikzpicture}
  \caption{A comparison between the different discretization schemes of the
    interface evolution and the curvature $\kappa$ of drop collision in shear
    flow.}
  \label{fig:B1}
\end{figure}

\paragraph{My contribution:}  I wrote the manuscript, implemented CFDM and MLM,
and produced the results with CD-2, CFDM, and MLM.  Åsmund Ervik implemented
the SLM and ran the simulations that used the SLM.  He also gave feedback on
the manuscript.  I presented the work at the conference.

\section[Paper C]{Paper C:  Calculation of the interface curvature and normal
  vector with the level-set method}
Karl Yngve Lervåg, Bernhard Müller, and Svend Tollak Munkejord.  Published in
Computers and Fluids, volume 84 (2013), 218--230.

Paper A presented the curve-fitting discretization method (CFDM) for the
calculation of the curvature with the level-set method.  The method was
designed to be robust in the calculation of the curvature in kink regions, that
is regions where the gradient of the level-set function is not smooth.  This
paper presents the details of the CFDM and applies it to the calculation of
both the curvature and the normal vector.

In the paper we compare the CFDM with the second-order central-difference
scheme (CD-2) for several test cases.  In the first case, we consider the
curvature calculations for a nontrivial geometry that includes some kink
regions.  This is a static test case with no flow, and the results show that
the CD-2 leads to large errors for the curvature calculations in areas close to
kink regions and that these errors are not present with the CFDM.

In the following two cases, we consider the collision of two drops in a 2D
shear flow and in an axisymmetric flow.  These cases show that the errors in
the curvature calculations in the kink regions with the CD-2 lead to errors in
the pressure that prevents coalescence.  These errors are prevented with the
CFDM.  The curvature and the evolution of the interfaces for the axisymmetric
case are shown in \cref{fig:C1,fig:C2}.  As in the earlier results of Papers
A and B, the figures show that the errors in the curvature calculation with
CD-2 prevent coalescence, in this case leading to a slower coalescence process.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [
    plot/.style={inner sep=0em},
    axes/.style={->,>=stealth',thick},
    time/.style={above left=1em, fill=white},
    kapp/.style={below left},
    ]

    \node[plot] (p1) at (0,0)
      {\includegraphics[width=0.26\textwidth]{ref_0}};
    \coordinate (dy) at ($0.25*(p1.north) - 0.25*(p1.south)$);
    \draw let \p1=(dy) in (\y1,0) coordinate (dx);
    \draw[axes] (p1.south west) -- (p1.south east) node[right=0pt] {$r$};
    \draw[axes] (p1.south west) -- (p1.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/-0.5, 2/0.0, 3/0.5}
      \draw ($(p1.south west) + \n*(dy) + (0.1,0)$) -- +(-0.2,0)
        node[left] {\label};

    \node[plot] (p2) [right=1.5em of p1]
      {\includegraphics[width=0.26\textwidth]{ref_1}};
    \draw[axes] (p2.south west) -- (p2.south east) node[right=0pt] {$r$};
    \draw[axes] (p2.south west) -- (p2.north west) node[above=0pt] {$z$};

    \node[plot] (p3) [right=1.5em of p2]
      {\includegraphics[width=0.26\textwidth]{ref_2}};
    \draw[axes] (p3.south west) -- (p3.south east) node[right=0pt] {$r$};
    \draw[axes] (p3.south west) -- (p3.north west) node[above=0pt] {$z$};

    \node[plot] (p4) [below=2em of p1]
      {\includegraphics[width=0.26\textwidth]{ref_3}};
    \draw[axes] (p4.south west) -- (p4.south east) node[right=0pt] {$r$};
    \draw[axes] (p4.south west) -- (p4.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/-0.5, 2/0.0, 3/0.5}
      \draw ($(p4.south west) + \n*(dy) + (0.1,0)$) -- +(-0.2,0)
        node[left] {\label};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p4.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};

    \node[plot] (p5) [right=1.5em of p4]
      {\includegraphics[width=0.26\textwidth]{ref_4}};
    \draw[axes] (p5.south west) -- (p5.south east) node[right=0pt] {$r$};
    \draw[axes] (p5.south west) -- (p5.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p5.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};

    \node[plot] (p6) [right=1.5em of p5]
      {\includegraphics[width=0.26\textwidth]{ref_5}};
    \draw[axes] (p6.south west) -- (p6.south east) node[right=0pt] {$r$};
    \draw[axes] (p6.south west) -- (p6.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p6.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};
    \node[kapp] at (p6.north east) {$\kappa$};

    \node[time] at (p1.south east) {$t=\SI{0.30}{s}$};
    \node[time] at (p2.south east) {$t=\SI{0.42}{s}$};
    \node[time] at (p3.south east) {$t=\SI{0.43}{s}$};
    \node[time] at (p4.south east) {$t=\SI{0.48}{s}$};
    \node[time] at (p5.south east) {$t=\SI{0.49}{s}$};
    \node[time] at (p6.south east) {$t=\SI{0.60}{s}$};

  \end{tikzpicture}
  \caption{Drop collision in axisymmetric flow calculated with the CD-2.  The
    legend for the colour contours of the curvature $\kappa$ is shown in the
    last image.  The velocity vectors are displayed to show the evolution of
    the flow during the collision.}
  \label{fig:C1}
\end{figure}
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    [
    plot/.style={inner sep=0em},
    axes/.style={->,>=stealth',thick},
    time/.style={above left=1em, fill=white},
    kapp/.style={below left},
    ]

    \node[plot] (p1) at (0,0)
      {\includegraphics[width=0.26\textwidth]{locurv_0}};
    \coordinate (dy) at ($0.25*(p1.north) - 0.25*(p1.south)$);
    \draw let \p1=(dy) in (\y1,0) coordinate (dx);
    \draw[axes] (p1.south west) -- (p1.south east) node[right=0pt] {$r$};
    \draw[axes] (p1.south west) -- (p1.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/-0.5, 2/0.0, 3/0.5}
      \draw ($(p1.south west) + \n*(dy) + (0.1,0)$) -- +(-0.2,0)
        node[left] {\label};

    \node[plot] (p2) [right=1.5em of p1]
      {\includegraphics[width=0.26\textwidth]{locurv_1}};
    \draw[axes] (p2.south west) -- (p2.south east) node[right=0pt] {$r$};
    \draw[axes] (p2.south west) -- (p2.north west) node[above=0pt] {$z$};

    \node[plot] (p3) [right=1.5em of p2]
      {\includegraphics[width=0.26\textwidth]{locurv_2}};
    \draw[axes] (p3.south west) -- (p3.south east) node[right=0pt] {$r$};
    \draw[axes] (p3.south west) -- (p3.north west) node[above=0pt] {$z$};

    \node[plot] (p4) [below=2em of p1]
      {\includegraphics[width=0.26\textwidth]{locurv_3}};
    \draw[axes] (p4.south west) -- (p4.south east) node[right=0pt] {$r$};
    \draw[axes] (p4.south west) -- (p4.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/-0.5, 2/0.0, 3/0.5}
      \draw ($(p4.south west) + \n*(dy) + (0.1,0)$) -- +(-0.2,0)
        node[left] {\label};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p4.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};

    \node[plot] (p5) [right=1.5em of p4]
      {\includegraphics[width=0.26\textwidth]{locurv_4}};
    \draw[axes] (p5.south west) -- (p5.south east) node[right=0pt] {$r$};
    \draw[axes] (p5.south west) -- (p5.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p5.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};

    \node[plot] (p6) [right=1.5em of p5]
      {\includegraphics[width=0.26\textwidth]{locurv_5}};
    \draw[axes] (p6.south west) -- (p6.south east) node[right=0pt] {$r$};
    \draw[axes] (p6.south west) -- (p6.north west) node[above=0pt] {$z$};
    \foreach \n/\label in {1/0.5, 2/1.0}
      \draw ($(p6.south west) + \n*(dx) + (0,0.1)$) -- +(0,-0.2)
        node[below] {\label};
    \node[kapp] at (p6.north east) {$\kappa$};

    \node[time] at (p1.south east) {$t=\SI{0.30}{s}$};
    \node[time] at (p2.south east) {$t=\SI{0.42}{s}$};
    \node[time] at (p3.south east) {$t=\SI{0.43}{s}$};
    \node[time] at (p4.south east) {$t=\SI{0.48}{s}$};
    \node[time] at (p5.south east) {$t=\SI{0.49}{s}$};
    \node[time] at (p6.south east) {$t=\SI{0.60}{s}$};

  \end{tikzpicture}
  \caption{Drop collision in axisymmetric flow calculated with the CFDM.  The
    legend for the colour contours of the curvature $\kappa$ is shown in the
    last image.  The velocity vectors are displayed to show the evolution of
    the flow during the collision.}
  \label{fig:C2}
\end{figure}

In the final test case, we consider the calculation of the normal vector, and
we compare the CD-2, the CFDM, and the direction-difference scheme (DDS)
presented by \citet{Macklin05}, cf.\ \cref{sec:dds}.  The results show that
both the CFDM and the DDS generally lead to good results, see \cref{fig:C3}.
Here the red and green vectors depict the DDS and CFDM results, respectively.
The red vectors are plotted below the green vectors, and since the results
agree well at most points, the red vector is often covered by its corresponding
green vector.  However, at the point in the middle between the drops, the DDS
completely fails to calculate the normal vector.  Here the CFDM still gives
a reasonable result.
\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.70\textwidth]{two-disc-normal-vectors.pdf}
  \caption{A comparison of the DDS and the CFDM for calculating normal vectors.
    The thick black lines depict the interfaces, the green vectors are the
    results with the CFDM, and the red vectors are the results with the DDS.
    The red vectors are covered by the green vectors at most points, because
    the results agree well at those points.}
  \label{fig:C3}
\end{figure}

\paragraph{My contribution:}  I designed the new method, implemented it into
our in-house finite-difference code, ran the simulations, and wrote the paper
manuscript.  The co-authors contributed with feedback on the manuscript and
discussions of the results.

\clearpage
\section[Paper D]{Paper D:  A robust method for calculating interface curvature
  and normal vectors using an extracted local level set}
Åsmund Ervik, Karl Yngve Lervåg, and Svend Tollak Munkejord.  Submitted to
Journal of Computational Physics, 2013.

In this paper we present an alternative method for the calculation of the
curvature and the normal vector of an interface with the level-set method in
kink regions, hereafter called the local level-set extraction (LOLEX) method.
The method is based on a method presented by \citet{Salac08} (SLM), who handle
the kink region by extracting different bodies of a domain into separate
level-set functions.  This procedure removes most of the kink regions, but it
does not handle kink regions that are due to complex interfaces of single
bodies.  Our method extends the SLM by making it local.  That is, we only
consider the local area around the point for which we are calculating the
curvature or the normal vector.  This leads to a method that is more generally
applicable, as shown in \cref{fig:D1}.  The figure shows a comparison between
the LOLEX method, the SLM, and the standard central differences (CD-2).  CD-2
leads to curvature spikes at the kink regions, as explained in
\cref{sec:curvature}.  The SLM gives a better result for the kink regions
around the rightmost disc.  However, since the other two discs are connected to
each other and to the film, they are considered to be the same body and are
extracted into the same level-set function.  Several kink regions are therefore
not removed.  Since the LOLEX method only considers the local area, as
explained in \cref{sec:lolex}, it is able to handle all the kink regions in
a robust manner.
\begin{figure}[tbp]
  \centering
  \begin{subfigure}[t]{\textwidth}
    \centering
    \includegraphics[width=0.6\textwidth]{init-test-lolex-2.png}
    \caption{LOLEX method}
  \end{subfigure} \\
  \begin{subfigure}[t]{\textwidth}
    \centering
    \includegraphics[width=0.6\textwidth]{init-test-salac-2.png}
    \caption{SLM}
  \end{subfigure} \\
  \begin{subfigure}[t]{\textwidth}
    \centering
    \includegraphics[width=0.6\textwidth]{init-test-standard-2.png}
    \caption{CD-2}
  \end{subfigure}
  \caption{Comparison of curvature calculation methods for three discs and
    a film with an angle at the right-hand side.  The film is connected to the
    leftmost disc, which is connected to the middle disc.  The rightmost disc
    is disjoint.  The color indicates the curvature; white is zero, blue is
    negative and red is positive.}
  \label{fig:D1}
\end{figure}

The LOLEX method has proven to be a good alternative to the CFDM presented in
Paper A.  Its main advantages are that it does not rely on complex algorithms
as used in the CFDM or by \citet{Macklin06}, and that the method easily extends
to 3D as demonstrated in the paper in Section~4.4.

In the previous papers A--C, we used the DP1 projection method by
\citet{Hansen05}.  In this paper we instead used the more standard Chorin
projection method.  These methods differ in that the DP1 assumes
\begin{equation}
  \div\left(\pdt{\vct u}\right) = 0.
\end{equation}
When compared with the Chorin method, this assumption becomes equivalent to
assuming that $\div\vct u^{n} = 0$ in \eqref{eq:PoissonChorin}.  That is, the
DP1 assumes that the initial velocity field is divergence free.  We have found
that the DP1 works well in most cases, but that it is less robust than the
Chorin method.  In particular, the Chorin method is not equally affected by
errors in the curvature calculations in kink regions.  In other words, the
difference between using a standard discretization and an improved
discretization of the curvature is smaller with the Chorin method than with the
DP1.

The LOLEX method is used for several test cases and compared with CD-2.  The
results indicate that even though we use the Chorin projection method, the
LOLEX method outperforms CD-2 in all cases.  Further, the results agree well
with experiments, except for time instants, as shown in \cref{fig:D2}.  The
exact reason why the time instants do not match is not known, but one reason
may be that the initial condition of the numerical simulation did not match the
corresponding state of the experiment.
\begin{figure}[tbp]
  \centering
  \begin{subfigure}[t]{\textwidth}
    \centering
    \includegraphics[width=\linewidth]{comp-water-partial_exp}
    \caption{Experimental result}
  \end{subfigure}
  \begin{subfigure}[t]{\textwidth}
    \centering
    \includegraphics[width=\linewidth]{comp-water-partial_sim}
    \caption{Simulation result}
  \end{subfigure}
  \caption{Experimental results (top) and simulation results (bottom) for
    a $\SI{0.18}{mm}$ water drop falling through air and impacting a deep
    pool of water at $\SI{0.29}{m/s}$.  Figure (a) is reprinted
    from~\cite{Zhao11}, Copyright (2011), with permission from Elsevier.}
  \label{fig:D2}
\end{figure}

When we study the drop-film collision processes, an important consequence of
the error in the curvature calculation is a loss of kinetic energy.  This can
be seen in \cref{fig:D3}, which compares the LOLEX method with CD-2 at two
different stages of the collision process.  The figure also compares two
different frequencies of reinitialization of the level-set function:  Every
7 time steps ((a) and (c)) and every single time step ((b) and (d)).  The
figure shows that the error in the curvature calculation with CD-2 leads to
a shorter neck, as seen in \cref{fig:D3}, (c) and (d).  The error is larger for
the higher frequency of reinitialization.  The results indicate that CD-2 leads
to a loss of kinetic energy during the collision process when compared with the
LOLEX method.  The LOLEX method is not significantly affected by the amount of
reinitialization, and the kink region does not affect the curvature
calculation, cf.\ \cref{sec:lolex}.  Thus the pressure field is more sensible,
as seen in \cref{fig:D3} (a) and (b).  Finally, we remark that some authors
have noted \cite{Blanchette06} that the height of the neck and the dynamics of
the capillary waves are important factors for the partial coalescence
mechanism, which implies that the correct calculation of the curvature is
important to capture the correct physical behavior.
\begin{figure}[tbp]
  \tikzstyle{label}=[text width=2.5cm, font=\scriptsize]
  \tikzstyle{label1}=[label, above right=5pt]
  \tikzstyle{label2}=[label, above=2pt]
  \centering
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \begin{tikzpicture}
      \node (img) at (0,0)
        {\includegraphics[width=\textwidth]{pres_coll_i7}};
      \node[font=\footnotesize,below right=5pt]
        at (img.north west) {$p$ [\si{Pa}]};
      \node[label1] (lbl1) at (img.south west) {$t=\SI{1.327e-4}{s}$};
      \node[label2] at (lbl1)                  {LOLEX};
      \node[label1] (lbl2) at (img.south)      {$t=\SI{1.386e-4}{s}$};
      \node[label2] at (lbl2)                  {CD-2};
    \end{tikzpicture}
    \caption{Reinitialization every 7 time steps}
  \end{subfigure}
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \begin{tikzpicture}
      \node (img) at (0,0)
        {\includegraphics[width=\textwidth]{pres_coll_i1}};
      \node[label1] (lbl1) at (img.south west) {$t=\SI{1.323e-4}{s}$};
      \node[label2] at (lbl1)                  {LOLEX};
      \node[label1] (lbl2) at (img.south)      {$t=\SI{1.342e-4}{s}$};
      \node[label2] at (lbl2)                  {CD-2};
    \end{tikzpicture}
    \caption{Reinitialization every time step}
  \end{subfigure}
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \begin{tikzpicture}
      \node (img) at (0,0)
        {\includegraphics[width=\textwidth]{pres_neck_i7}};
      \node[font=\footnotesize,below right=4pt]
        at (img.north west) {$p$ [\si{Pa}]};
      \node[label1] (lbl1) at (img.south west) {$t=\SI{2.408e-4}{s}$};
      \node[label2] at (lbl1)                  {LOLEX};
      \node[label1] (lbl2) at (img.south)      {$t=\SI{2.459e-4}{s}$};
      \node[label2] at (lbl2)                  {CD-2};
    \end{tikzpicture}
    \caption{Reinitialization every 7 time steps}
  \end{subfigure}
  \begin{subfigure}[t]{0.47\textwidth}
    \centering
    \begin{tikzpicture}
      \node (img) at (0,0)
        {\includegraphics[width=\textwidth]{pres_neck_i1}};
      \node[label1] (lbl1) at (img.south west) {$t=\SI{2.412e-4}{s}$};
      \node[label2] at (lbl1)                  {LOLEX};
      \node[label1] (lbl2) at (img.south)      {$t=\SI{2.544e-4}{s}$};
      \node[label2] at (lbl2)                  {CD-2};
    \end{tikzpicture}
    \caption{Reinitialization every time step}
  \end{subfigure}
  \caption{Water drop falling onto a pool, a comparison between the LOLEX
    method and CD-2.  The interfaces are shown as solid black lines and the
    pressure field is shown as colored contours.  (a) and (b): just before the
    interfaces merge.  (c) and (d):  when the neck reaches its highest
    position.}
  \label{fig:D3}
\end{figure}

\paragraph{My contribution:}  The manuscript was written by Åsmund Ervik.  The
new method was developed and implemented into our in-house finite-difference
code by Åsmund Ervik, and most of the numerical results are due to Åsmund
Ervik.  I contributed with discussions during the development of the new
method, designed the test case in Section~4.1, ran simulations for Section~5.2,
created some of the result figures, and gave feedback on the manuscript.
I also assisted in some of the programming efforts for initializing the test
cases in Chapter~5.  Svend Tollak Munkejord contributed with discussions of the
manuscript and some code testing.

\clearpage
\section[Paper E]{Paper E:  Towards a second-order diffuse-domain approach for
  solving PDEs in complex geometries}
Karl Yngve Lervåg and John Lowengrub.  Submitted to Communications in Math.
Sciences, 2013.

\citet{Li09} developed a diffuse-domain method for solving partial-differential
equations (PDEs) inside complex, dynamic geometries with Dirichlet, Neumann,
and Robin boundary conditions.  This method is in the following referred to as
DDM1.  They use the diffuse-domain approach~\cite{Kockelkoren03}, where the
geometry is represented implicitly and the sharp boundary is replaced by
a diffuse layer with a fixed interface width.  The original governing equations
are then reformulated on a larger, regular domain and the boundary conditions
are incorporated via singular source terms.  The method of matched asymptotic
expansions is used to show that the reformulated problem converges
asymptotically to the original problem.

In the present paper, we use the method of matched asymptotic expansions to
extend the DDM1 with include a high-order correction term in the diffuse
formulation, cf.\ \cref{sec:DDM2}.  The extension is derived for elliptic
problems with Neumann and Robin boundary conditions, where the correction term
is shown to yield an asymptotically second-order accurate approximation of the
original problem.  The new method is referred to as the DDM2.

The DDM1 and DDM2 are compared for a selection of test problems.  The resulting
equations were discretized by standard second-order central-difference schemes
on uniform grids, and solved by a multigrid method.  A red-black Gauss-Seidel
type iterative method was used as a smoother, see \cite{Wise07}.

In addition to the comparison of DDM1 and DDM2, we compared two different
approximations of the boundary conditions.  These correspond to different
diffuse-interface surface delta functions, and for the Neumann boundary
conditions they are
\begin{equation}
  \text{BC1} = |\grad\phi| g,
\end{equation}
and
\begin{equation}
  \text{BC2} = \epsilon|\grad\phi|^2 g.
\end{equation}
The approximations are similar for the Robin boundary condition, although here
it was found that only BC1 resulted in a valid asymptotically second-order
accurate DDM2.

Our results show that the global accuracy and convergence of DDM2 is better
than DDM1, however both methods perform well and the global convergence rate is
around two for each.  The error was measured as the L$_2$ norm of the
difference between an analytic solution and the solution $u_\epsilon$ for
a given interface width, $\epsilon$.

\Cref{fig:E1,fig:E2} show two of the results, the first one with the Neumann
boundary conditions and the second with the Robin boundary condition.  The
results indicate that DDM2 performs slightly better than DDM1.  They also show
that the approximation BC1 gives more accurate results than BC2.  In
particular, we found that with BC2 we needed much finer grids to obtain
convergence for a given $\epsilon$.  For the smallest values of $\epsilon$ we
were not able to refine the grids enough to obtain valid results.
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    \begin{loglogaxis}[
      xlabel={Interface width, $\epsilon$},
      ylabel={$E_\epsilon$},
      x dir=reverse,
      width=0.8\textwidth,
      legend entries={DDM1 BC1,
                      DDM2 BC1,
                      DDM1 BC2,
                      DDM2 BC2},
      legend cell align=left,
      legend style={column sep=0.5em,draw=white},
      ]

    \addplot[solid,mark=*,black] plot coordinates {
      (8.000e-01, 2.458e-01)
      (4.000e-01, 7.303e-02)
      (2.000e-01, 1.935e-02)
      (1.000e-01, 5.161e-03)
      (5.000e-02, 1.580e-03)
    };

    \addplot[dashed,mark=square*,mark options=solid,black] plot coordinates {
      (8.000e-01, 1.960e-01)
      (4.000e-01, 2.582e-02)
      (2.000e-01, 5.205e-03)
      (1.000e-01, 1.201e-03)
      (5.000e-02, 3.700e-04)
    };

    \addplot[solid,mark=*,green] plot coordinates {
      (8.000e-01, 2.221e-01)
      (4.000e-01, 6.989e-02)
      (2.000e-01, 1.952e-02)
      (1.000e-01, 5.879e-03)
    };

    \addplot[dashed,mark=square*,mark options=solid,green] plot coordinates {
      (8.000e-01, 3.530e-01)
      (4.000e-01, 5.101e-02)
      (2.000e-01, 1.081e-02)
      (1.000e-01, 2.991e-03)
    };

    \end{loglogaxis}
  \end{tikzpicture}
  \caption{Errors for the Neumann problem with respect to $\epsilon$ for Case
    2, as labelled.}
  \label{fig:E1}
\end{figure}
\begin{figure}[tbp]
  \centering
  \begin{tikzpicture}
    \begin{loglogaxis}[
      xlabel={Interface width, $\epsilon$},
      ylabel={$E_\epsilon$},
      x dir=reverse,
      width=0.8\textwidth,
      legend entries={DDM1 BC1,
                      DDM2 BC1,
                      DDM1 BC2,
                      DDM2 BC2},
      legend cell align=left,
      legend style={column sep=0.5em,draw=white},
      ]

    \addplot[solid,mark=*,black] plot coordinates {
      (8.000e-01, 1.323e-01)
      (4.000e-01, 2.754e-02)
      (2.000e-01, 5.271e-03)
      (1.000e-01, 1.018e-03)
      (5.000e-02, 2.082e-04)
      (2.500e-02, 4.876e-05)
    };

    \addplot[dashed,mark=square*,mark options=solid,black] plot coordinates {
      (8.000e-01, 2.747e-02)
      (4.000e-01, 8.767e-03)
      (2.000e-01, 2.199e-03)
      (1.000e-01, 5.472e-04)
      (5.000e-02, 1.354e-04)
      (2.500e-02, 3.431e-05)
    };

    \addplot[solid,mark=*,green] plot coordinates {
      (8.000e-01, 1.137e-01)
      (4.000e-01, 2.536e-02)
      (2.000e-01, 8.536e-03)
      (1.000e-01, 3.860e-03)
      (5.000e-02, 2.005e-03)
    };

    \addplot[dashed,mark=square*,mark options=solid,green] plot coordinates {
      (8.000e-01, 9.798e-02)
      (4.000e-01, 3.034e-02)
      (2.000e-01, 8.098e-03)
      (1.000e-01, 3.053e-03)
      (5.000e-02, 1.689e-03)
    };
    \end{loglogaxis}
  \end{tikzpicture}
  \caption{Errors for the Robin problem with respect to $\epsilon$ for Case 2,
    as labelled.}
  \label{fig:E2}
\end{figure}

\paragraph{My contribution:}  I took a leading role in analysing the equations,
designing the numerical algorithms, selecting the test cases, and performing
the numerical simulations.  I wrote the manuscript.  John Lowengrub contributed
with important insights in the analysis, feedback on the manuscript, and
discussions of the results.

% Fakesection: Quote
\begin{savequote}[5cm]
  ``What is research but a blind date with knowledge?''
  \qauthor{--- Will Harvey (1967)}
\end{savequote}
\chapter{Conclusions and outlook}
\label{chap:conclusions}
This thesis has considered two different problems:  The discretization of
interface curvature and normal vectors with the level-set method and
a diffuse-domain approach for solving partial-differential equations (PDEs) in
complex domains.  In the following, some general concluding remarks and
recommendations for further work are given.

\section*{Conclusions}
The first part of the thesis considered the modelling of incompressible
two-phase flow.  The main motivation was to study two-phase flow phenomena that
are relevant for compact heat exchangers, such as drop-drop and drop-film
collisions.  In particular, the thesis addressed a challenge with the
calculation of interface curvature and normal vectors with the level-set
method.  Two methods were presented to handle the discretization in the kink
regions: A curve-fitting discretization method (CFDM) and a local level-set
extraction (LOLEX) method.  Both methods were shown to be robust in the kink
regions, where the standard central-difference scheme (CD-2) fails.  Of these
methods, the LOLEX method is the preferred method, because it relies on a less
complicated algorithm that easily extends to 3D.

CD-2 and the LOLEX method were used for simulations of drop-film collisions
that were compared with experiments.  The results showed that CD-2 leads to
errors in the curvature that cause unphysical pressure spikes during the
coalescence.  The errors are shown to lead to a dissipation of kinetic energy
during the collision and to a slower coalescence process.  The LOLEX method was
shown to prevent these unphysical pressure spikes, and to produce to more
accurate results.

The second part of the thesis considered the diffuse-domain approach for
solving PDEs in complex domains.  The main contribution, presented in Paper~E,
was the derivation of an asymptotically second-order diffuse-domain method
(DDM2) for solving elliptic problems in complex geometries with Neumann and
Robin boundary conditions.  The new method is an extension by a high-order
correction term of the method presented in \cite{Li09}, here called DDM1.  The
DDM1 and DDM2 were compared, and the results indicated that the DDM2 was
slightly better than the DDM1.

The thesis has expanded on the results of Paper~E with a new asymptotic
analysis that shows that DDM1 is in fact also asymptotically second order.
This analysis helps to explain why the performance of DDM2 presented in Paper~E
is only slightly better than that of DDM1.  As such, the new analysis leads to
a better understanding of the DDM1 and the DDM2.

\section*{Outlook}
The following gives an outline of some possibilities for future work.
\begin{itemize}
  \item A natural continuation of this work is to perform more in-depth
    comparisons of simulations with experiments for the drop-film collision
    phenomenon that was started in Paper~D.  In addition, it would be
    interesting to study other two-phase flow phenomena that are relevant for
    heat-exchanger processes, for instance drop-drop collisions or flow across
    tube bundles.  The latter requires the treatment of more complex boundaries,
    which can be handled either with the diffuse-domain method or other methods
    from the literature.
  \item Paper~C gave a short comparison of the standard discretization method
    (CD-2), the direction-difference scheme (DDS), and the curve-fitting
    discretization scheme (CFDM) applied to the calculation of the normal
    vectors.  The results showed that the CD-2 leads to inaccurate results in
    the kink regions.  The DDS gives robust and accurate results in most cases,
    but an example is given where only the CFDM yields an accurate result.
    However, the impact of inaccurate calculations of the normal vector should
    be investigated in more detail.  The normal vector is used both for the
    solution of the level-set equations \eqref{eq:ls_adeq},
    \eqref{eq:ls_velextr}, and \eqref{eq:ls_reinit}, and for the calculation of
    the interface jumps \eqref{eq:jump_pressure} and \eqref{eq:jumptens}, so
    one can expect that large errors in the calculation of the normal vector
    may lead to large errors in the numerical solution.  This warrants
    a further study.
  \item Mass and heat transfer are obviously an important part of the
    heat-exchanger processes.  The models that have been used in this thesis
    should therefore be expanded with additional models for heat transfer and
    mass transfer to enable the simulation of more relevant phenomena.
  \item The DDM2 was only derived for problems with Robin and Neumann boundary
    conditions.  If possible, it should be extended to also work with Dirichlet
    boundary conditions.
  \item As explained in Paper~E, we found that we were unable to solve the
    discrete system of equations when the surface Laplacian part of the
    correction term for the DDM2 was included.  A further investigation of this
    problem should be done, and stable numerical methods to solve the full DDM2
    equations should be developed.
  \item The diffuse-domain approach is a promising method for solving problems
    in complex and confined geometries with standard tools and methods.
    \citet{Aland10} provide a diffuse-domain formulation of the Navier-Stokes
    Cahn-Hilliard equations for incompressible two-phase flow and use it to
    compute two-phase flow in both complex and confined geometries.  However,
    they do not provide an asymptotic analysis of the reformulated equations to
    show that the system converges.  Such an analysis would be interesting, in
    particular to verify that the equations do converge to the original problem
    when the interface width is decreased.
  \item Finally, it would be interesting to develop a second-order
    diffuse-domain formulation for the incompressible Navier-Stokes equations.
    A starting point would be to use the results and techniques from this
    thesis and Paper~E.
\end{itemize}

% Fakechapter:  Bibliography
\printbibliography
\addcontentsline{toc}{chapter}{Bibliography}

% Fakesection: Appendix
\appendix

\chapter{Calculation of interface curvature with the level-set method}
K.\ Y.\ Lervåg\\
Published in \emph{MekIT'11 - 6th National Conference on Computational
  Mechanics, Trondheim}, 2011.  ISBN: 978-82-519-2798-7.

\cleardoublepage
\includepdf[pages=-]{papers/2011-04-27_lervag_mekit2011.pdf}

\chapter{Curvature calculations for the level-set method}
K.\ Y.\ Lervåg and Å.\ Ervik\\
Published in \emph{ENUMATH 2011} proceedings volume, Springer, 2013.  ISBN:
978-3642331336.

\cleardoublepage
\includepdf[pages=-]{papers/2012-03-22_lervag_enumath.pdf}

\chapter{Calculation of the interface curvature and normal vector with the
  level-set method}
K.\ Y.\ Lervåg, B.\ Müller, and S.\ T.\ Munkejord\\
Published in Computers and Fluids, volume 84 (2013), 218--230.

\cleardoublepage
\includepdf[pages=-]{papers/2013-04-18_lervag_CAF.pdf}

\chapter{A robust method for calculating interface curvature and normal vectors
  using an extracted local level set}
Å.\ Ervik, K.\ Y.\ Lervåg, and S.\ T.\ Munkejord\\
Submitted to Journal of Computational Physics, 2013

\cleardoublepage
\includepdf[pages=-]{papers/2013-02-27_ervik_LOLEX.pdf}

\chapter{Towards a second-order diffuse-domain approach for solving PDEs in
  complex geometries}
K.\ Y.\ Lervåg and J.\ Lowengrub\\
Submitted to Communications in Mathematical Sciences, 2013

\cleardoublepage
\includepdf[pages=-]{papers/2013-01-29_lervag_DDA.pdf}

\end{document}