1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

16 
\section{Level Set Method} 
17 

18 
The Level Set Method is used for tracking interfaces between two different types of fluids, which may have different physical parameter values for density or viscosity. The interface is represented by a signed distance function, $\phi(x)$, where the isocontour at $\phi(x)=0$ is used to defined the interface. A point in the domain can then be determined on which side of the interface it resides, based on the local sign of $\phi(x)$; for example positive $\phi(x)$ on one side of the interface and negative $\phi(x)$ on the other. Parameters values such as density and viscosity can then be defined for the two different mediums. The Level Set Method consists of two procedures, the advection and reinitialization of the signed distance function, $\phi$. The LevelSet class can be used in conjunction with the StokesProblemCartesian class for solving computational fluid dynamics problems involving the tracking of the interface. The advantage of the Level Set Method is that it can be used to track surfaces that break apart or intersect. Also, the Level Set Method avoids the need for remeshing, which is required by the LagrangianEulerian (ALE) method. 
19 

20 
% An example of using the Level Set Method is described in the tutorial Chapter, Section %\ref{LEVELSET CHAP}. 
21 

22 
\subsection{Solution Method} 
23 

24 
The displacement of the interface at the zero isocontour of $\phi(x)$ is calculated each timestep by using the velocity field. This is achieved my solving the advection equation: 
25 
% 
26 
\begin{equation} 
27 
\frac{\partial \phi}{\partial t} + \vec{v} \cdot \nabla \phi = 0, 
28 
\label{ADVECTION MODELS} 
29 
\end{equation} 
30 
% 
31 
where $\vec{v}$ is the velocity field. The advection equation is solved using a TaylorGalerkin scheme with the presence of diffusion; by expanding $\phi$ into a Taylor series: 
32 
% 
33 
\begin{equation} 
34 
\phi^{+} \simeq \phi^{} + dt\frac{\partial \phi^{}}{\partial t} + \frac{dt^2}{2}\frac{\partial^{2}\phi^{}}{\partial t^{2}}, 
35 
\label{TAYLOR EXPANSION MODELS} 
36 
\end{equation} 
37 
% 
38 
then by inserting 
39 
% 
40 
\begin{equation} 
41 
\frac{\partial \phi^{}}{\partial t} =  \vec{v} \cdot \nabla \phi^{}, 
42 
\label{INSERT ADVECTION MODELS} 
43 
\end{equation} 
44 
% 
45 
and 
46 
% 
47 
\begin{equation} 
48 
\frac{\partial^{2} \phi^{}}{\partial t^{2}} = \frac{\partial}{\partial t}(\vec{v} \cdot \nabla \phi^{}) = \vec{v}\cdot \nabla (\vec{v}\cdot \nabla \phi^{}), 
49 
\label{SECOND ORDER MODELS} 
50 
\end{equation} 
51 
% 
52 
into Equation (\ref{TAYLOR EXPANSION MODELS}), the calculation of the level set function is given by: 
53 
% 
54 
\begin{equation} 
55 
\phi^{+} = \phi^{}  dt\vec{v}\cdot \nabla \phi^{} + \frac{dt^2}{2}\vec{v}\cdot \nabla (\vec{v}\cdot \nabla \phi^{}). 
56 
\label{TAYLOR GALERKIN MODELS} 
57 
\end{equation} 
58 

59 
If $\nabla \cdot \vec{v}=0$ is assumed, then the calculation of the second order derivatives in Equation (\ref{TAYLOR GALERKIN MODELS}) can be avoided. 
60 

61 
As the computation of the distance function progresses, it becomes distorted, and so it needs to be updated in order to stay regular \cite{SUSSMAN1994}. This process is known as the reinitialization procedure. The aim is to iteratively find a solution to the reinitialization equation: 
62 
% 
63 
\begin{equation} 
64 
\frac{\partial \psi}{\partial \tau} + sign(\phi)(1  \nabla \psi) = 0. 
65 
\label{REINITIALISATION MODELS} 
66 
\end{equation} 
67 
% 
68 
where $\psi$ shares the same level set with $\phi$, $\tau$ is pseudo time, and $sign(\phi)$ is the smoothed sign function. This equation is solved to meet the definition of the level set function, $\lvert \nabla \psi \rvert = 1$; the normalization condition. Equation (\ref{REINITIALISATION MODELS}) can be rewritten in a similar form to the advection equation: 
69 
% 
70 
\begin{equation} 
71 
\frac{\partial \psi}{\partial \tau} + \vec{w} \cdot \nabla \psi = sign(\phi), 
72 
\label{REINITIALISATION2 MODELS} 
73 
\end{equation} 
74 
% 
75 
where 
76 
% 
77 
\begin{equation} 
78 
\vec{w} = sign(\phi)\frac{\nabla \psi}{\nabla \psi}. 
79 
\label{REINITIALISATION3 MODELS} 
80 
\end{equation} 
81 
% 
82 
$\vec{w}$ is the characteristic velocity pointing outward from the free surface. Equation (\ref{REINITIALISATION2 MODELS}) can be solved by a similar technique to what was used in the advection step, using the TaylorGalerkin procedure. 
83 
When the distance function, $\phi$, is calculated, the physical parameters, density and viscosity, are updated using the sign of $\phi$. The region along the interface is assumed to be of finite thickness of $\alpha h$, where $h$ is the size of the elements in the computational mesh and $\alpha$ is a smoothing parameter. The parameters are updated by the following expression: 
84 
% 
85 
\begin{equation} 
86 
P = 
87 
\left \{ \begin{array}{l} 
88 
P_{1} \hspace{5cm} where \ \ \psi <  \alpha h \\ 
89 
P_{2} \hspace{5cm} where \ \ \psi > \alpha h \\ 
90 
(P_{2}  P_{1}) \psi/2\alpha h + (P_{1} + P_{2})/2 \ \ \ \ \ \ where \ \ \psi < \alpha h. 
91 
\end{array} 
92 
\right. 
93 
\label{UPDATE PARAMETERS MODELS} 
94 
\end{equation} 
95 
% 
96 
where the subscripts $1$ and $2$ denote the different fluids. 
97 

98 

99 
\subsection{Functions} 
100 

101 
\begin{classdesc}{LevelSet}{domain, func, reinit\_max, reinit\_each, tolerance, smooth} 
102 
opens the LevelSet \index{Level Set} on the \Domain domain. \var{func} defines the initial Level Set function representing the interface between two fluids. \var{reinit\_max} sets the maximum number of iterations to satisfy the normal condition, $\nabla \phi=1$, during the reinitialization of the Level Set function. \var{reinit\_each} sets the frequency of reinitialization for a number of timesteps. \var{tolerance} sets the convergence tolerance to satisfy the normal condition during the reinitialization of the Level Set function. \var{smooth} sets the bandwidth of size 2$\alpha h$ along the interface to smooth the physical parameters of density and viscosity; $h$ is the size of the elements in the mesh and $\alpha$ is the smoothing parameter, usually set to 1. 
103 
\end{classdesc} 
104 

105 
\begin{methoddesc}[LevelSet]{update\_parameter}{par1, par2} 
106 
updates the physical parameters using the sign of $\phi$. \var{par1} and \var{par2} are the physical parameter values for fluid1 and fluid2 respectively. Usually this method is called twice during each timestep to update the density and viscosity of the two fluids. 
107 
\end{methoddesc} 
108 

109 
\begin{methoddesc}[LevelSet]{update\_phi}{vel, dt, t\_step} 
110 
updates the Level Set function. It performs the advection and reinitialization procedures. \var{vel} is the velocity field of the fluid domain, \var{dt} is the timestep size, and \var{t\_step} is the current timestep to determine when to reinitialize. 
111 
\end{methoddesc} 