The numerical solution scheme of the level set equation is very dependant of the application. First the different time stepping schemes are presented, followed by a description of further level set related algorithms.

Translational Motion Function

For the translational motion of solid interfaces a semi-Lagrange solver [[guenther2014]] can be used to determine the level-set function. In this approach, the error of the geometric representation of the surface is reduced to the constant interpolation error between two reference locations.

Rotational Motion Function

Arbitrary Motion Function

Initialisation of Level Set Field

The signed distance field data is initialized as initial condition from a reference STL file.

Band Cell Approach

As pointed out, the interface is commonly represented by the zero level set $\varphi_0$ bounding a region $\Omega^{+} \subset R^n$ and separating $\Omega^{-} \subset R^n$ from $\Omega^{+}$. It is embedded in the $n$-dimensional scalar level-set function $\varphi=\varphi(\mathbf{x}, t)$

Let the subset Γ of cells which are adjacent to the zero level set be defined by

\begin{align} \Gamma=\left\{C_{i, j, k}:\left(\Pi_{i^{\prime}, j, k}^{i, j, k} \varphi \leq 0\right) \vee\left(\Pi_{i, j^{\prime}, k}^{i, j, k} \varphi \leq 0\right) \vee\left(\Pi_{i, j, k^{\prime}}^{i, j k} \varphi \leq 0\right)\right\} \end{align}

That is, all cells in Γ are located within a distance $\Delta X$ from the zero level set.

The computational costs of solving the level-set equation can be reduced by orders of magnitude when the level-set method is localized, i.e., a solution is sought only in a small region around the zero level set, while all other areas are assigned a constant value indicating the location in $\Omega^{+}$ or $\Omega^{-}$. We consider a localized computational domain $\Omega_\varphi \subset \Omega$ moving along with the zero level set. All cells outside $\Omega_\varphi$ are discarded and the level-set algorithms are localized as proposed in , such that the computational costs of the overall level-set method scale with $\mathcal{O}(\mathcal{N})$. Let $\mathcal{B}$ designate the subset of cells $C_{i, j, k}$ which are used in the localized solution the level-set equation forming a narrow band around $\varphi_0$ bounded by the boundary cells $\hat{C}_{i, j, k} \in \partial \mathcal{B}$. Let us furthermore define $\partial \mathcal{B} \cap \Omega_\varphi=\emptyset$ such that boundary cells are outside of and adjacent to $\Omega_\varphi$. $\mathcal{B}$ is created using an efficient marching algorithm, which is based on neighbor relations. The level-set solver is based on a hierarchical quadtree/octree data structure, such that all neighbor information can be directly accessed. The narrow band $\mathcal{B}$ is regenerated before each reinitialization step, while the subset $\Gamma$ is updated after each time step.

Time and spatial discretization

The level-set equation is integrated in time with a 3-step third-order accurate TVD Runge-Kutta scheme [Shu] denoted by $\mathrm{RK}_3$,

\begin{align} \left\{\begin{array}{l} \varphi^{(0)}=\varphi^w \\ \varphi^{(k)}=\alpha_k \varphi^{(0)}+\beta_k \varphi^{(k-1)}-\gamma_k \Delta t L\left(\varphi^{(k-1)}\right) \\ \varphi^{w+1}=\varphi^{(N)} \end{array}\right. \end{align}

where $N=3$ and the coefficients $\boldsymbol{\alpha}=\left(0, \frac{3}{4}, \frac{1}{3}\right), \boldsymbol{\beta}=\left(1, \frac{1}{4}, \frac{2}{3}\right)$, and $\boldsymbol{\gamma}=\left(1, \frac{1}{4}, \frac{2}{3}\right)$ are used. The superscript $k$ denotes the Runge-Kutta step, while the superscript $w$ counts the time steps $\Delta t$. The operator $L(\varphi)$ denotes the numerical approximation of the term $v_{\Gamma} \cdot \boldsymbol{\nabla} \varphi$ in the level-set equation, which is specified in the following. To spatially discretize the level-set equation, unlimited third- and fifth-order upstream central schemes denoted by $\mathrm{UC}_3$ and $\mathrm{UC}_5$ are used, respectively. The upwind-biased spatial derivatives read as follows:

\begin{align} {UC}_3 : \left\{\begin{array}{l} D_{i, j, k}^{x, L}=\frac{1}{6 \Delta x}\left(\varphi_{i-2, j, k}-6 \varphi_{i-1, j, k}+3 \varphi_{i, j, k}+2 \varphi_{i+1, j, k}\right) \\ D_{i, j, k}^{x, R}=\frac{1}{6 \Delta x}\left(-\varphi_{i+2, j, k}+6 \varphi_{i+1, j, k}-3 \varphi_{i, j, k}-2 \varphi_{i-1, j, k}\right) \end{array}\right. \end{align}

\begin{align} {UC}_5 : \left\{\begin{array}{c} D_{i, j, k}^{x, L}=\frac{1}{60 \Delta x}\left(-2 \varphi_{i-3, j, k}+15 \varphi_{i-2, j, k}-60 \varphi_{i-1, j, k}+20 \varphi_{i, j, k}+30 \varphi_{i+1, j, k}-3 \varphi_{i+2, j, k}\right) \\ D_{i, j, k}^{x, R}=\frac{1}{60 \Delta x}\left(2 \varphi_{i+3, j, k}-15 \varphi_{i+2, j, k}+60 \varphi_{i+1, j, k}-20 \varphi_{i, j, k}-30 \varphi_{i-1, j, k}+3 \varphi_{i-2, j, k}\right) \end{array}\right. \end{align}

Likewise derivatives with respect to the $y$ and $z$ direction are obtained by exchanging the respective subscripts. Let us define $D_{i, j, k}^{x, \pm} \equiv D_{i, j, k}^{x, L} \pm D_{i, j, k}^{x, R}$ and introduce the vector operator $\mathfrak{D}_{i, j, k}^{\pm}=\left(D_{i, j, k}^{x, \pm}, D_{i, j, k}^{y, \pm}, D_{i, j, k}^{z, \pm}\right)^T$. Then, $L(\varphi)$ can be computed by

\begin{align} L\left(\varphi_{i, j, k}\right)=\frac{1}{2}\left\{\left(\sum_\zeta\left|f_\zeta\right| \mathbf{e}_\zeta\right) \cdot \mathfrak{D}_{i, j, k}^{+}+v_{\Gamma} \cdot \mathfrak{D}_{i, j, k}^{-}\right\}, \zeta=\{x, y, z\} \end{align}

Specification of the extension velocity

The motion of the zero level set $\varphi_0$ is governed by the extension velocity $v_{\Gamma}=v_{\Gamma}(\mathbf{x}, t)$,

\begin{align} v_{\Gamma}=\mathbf{v}+s \mathbf{n}, \end{align}

with the components $v_{\Gamma}=\left(f_x, f_y, f_z\right)^T$. It comprises the advection by an underlying flow velocity field $\mathbf{v}=\mathbf{v}(\mathbf{x}, t)$ and the propagation of the front relative to the flow field in the normal direction to $\varphi_0$ by $s$. The normal direction is defined by the outward normal vector

\begin{align} \mathbf{n}=-\frac{\boldsymbol{\nabla} \varphi}{|\boldsymbol{\nabla} \varphi|} \end{align}

pointing into $\Omega^{-}$, where $\boldsymbol{\nabla}=\left(\partial_x, \partial_y, \partial_z\right)^T$ denotes the vector operator of spatial derivatives. The local speed of propagation $s$ may be induced by several effects such as curvature and, in the case of premixed, level-set based combustion, the flame propagation into the unburnt gas. For geometric motion, the movement function of the boundary surface would also appear as a part of the translational term as a function of time.

The extension velocity vector $v_{\Gamma}$ is often only defined at the interface $\varphi_0$. Furthermore, one would like to move the entire level-set field along with the zero level set rather than move each level set individually according to the local velocity. It is therefore useful to compute $v_{\Gamma}$ or alternatively $f_n$ only on the cells in $\Gamma$ and then propagate this information along the normal direction to the interface. The first-order partial differential equation

\begin{align} \partial_\tau f+S(\varphi) \frac{\nabla \varphi}{|\nabla \varphi|} \cdot \nabla f=0 \end{align}

is formulated in artificial time $\tau$ for a generic quantity $f$, which is to be propagated from the interface. The given is a hyperbolic equation of Hamilton-Jacobi type. That is, the characteristics flow away from the zero level set along the normal vectors $\mathbf{n}$ and $\mathbf{-n}$.

Multiple Level Set Approach

Let $\Gamma(t)$ be the boundary of the fluid domain $\Omega(t)$, which is composed of a set of $N$ surface elements $\Gamma_i(t)$. Each element represents various parts of the geometry, e.g., individual STL files. Each surface element can be expressed as the zero contour of the signed-distance function $\varphi_{i+1}(\mathbf{x},t)$, which may evolve in time by individual motion functions. The overall surface $\Gamma(t)$ is constructed from all surfaces $\Gamma_i(t)$ using a combined level set $\varphi_0(\mathbf{x},t)$ as described in [[guenther2014]]. Where $\varphi_0(\mathbf{x},t)$ is constructed as:

\begin{align} \end{align}

In this approach, the distance from each surface can be easily evaluated by the absolute value of the level set value $\varphi_i(\mathbf{x},t)$. Thus, the the gap width between interfaces approaching

This approach can be used to determine

each other, i.e., when the engine valves are close to the liner valve seat. When the gap width is smaller than a specified threshold, the combined $\varphi_0(\mathbf{x},t)$ level-set field in the vicinity of the gap is modified to an artificially closed gap [[guenther2014]]. Since only the combined level-set is changed, no surface location information is lost through this modification.