High-resolution transport of regional level sets for evolving complex interface networks

TL;DR

This paper introduces a high-resolution regional level-set method that improves the accuracy of simulating complex interface network evolutions while maintaining comparable computational costs to existing methods.

Contribution

It presents a novel high-resolution transport formulation with local level set construction and reconstruction, enhancing the prediction of complex interface networks.

Findings

Accurate simulation of complex interface network evolutions.

Comparable computational cost to existing semi-Lagrangian methods.

Significantly improved prediction accuracy.

Abstract

In this paper we describe a high-resolution transport formulation of the regional level-set approach for an improved prediction of the evolution of complex interface networks. The novelty of this method is twofold: (i) construction of local level sets and reconstruction of a global regional level sets, (ii) locally transporting the interface network by employing high-order spatial discretization schemes for improved representation of complex topologies. Various numerical test cases of multi-region flow problems, including triple-point advection, single vortex flow, mean curvature flow, normal driven flow and dry foam dynamics, show that the method is accurate and suitable for a wide range of complex interface-network evolutions. Its overall computational cost is comparable to the Semi-Lagrangian regional level-set method while the prediction accuracy is significantly improved. The approach thus offers a \textbf{viable} alternative to previous interface-network level-set method.

Tables

Table 1. Table 1: Area conservation error ( % percent \% ) and CPU times (in seconds) of the Semi-Lagrangian regional level-set method (SL-RLS) and our high-resolution regional level-set method (HR-RLS) for simulating the single vortex flow. The area loss Δ ​ A ​ ( t 0 ) = ( Δ ​ A 1 ​ ( t 0 ) , Δ ​ A 2 ​ ( t 0 ) ) Δ 𝐴 subscript 𝑡 0 Δ superscript 𝐴 1 subscript 𝑡 0 Δ superscript 𝐴 2 subscript 𝑡 0 \Delta A(t_{0})=(\Delta A^{1}(t_{0}),\Delta A^{2}(t_{0})) and Δ ​ A ​ ( t 1 ) = ( Δ ​ A 1 ​ ( t 1 ) , Δ ​ A 2 ​ ( t 1 ) ) Δ 𝐴 subscript 𝑡 1 Δ superscript 𝐴 1 subscript 𝑡 1 Δ superscript 𝐴 2 subscript 𝑡 1 \Delta A(t_{1})=(\Delta A^{1}(t_{1}),\Delta A^{2}(t_{1})) are measured at physical time t 0 = 3.0 subscript 𝑡 0 3.0 t_{0}=3.0 and t 1 = 6.0 subscript 𝑡 1 6.0 t_{1}=6.0 , where the superscripts 1 1 1 and 2 2 2 indicate the outer and inner circles, respectively. The CPU time is measured at t = 3.0 𝑡 3.0 t=3.0 .

Method	$h=\frac{1}{64}$				$h=\frac{1}{128}$				$h=\frac{1}{256}$
	$\Delta A\left(t_0\right)$	$\Delta A\left(t_1\right)$	CPU time		$\Delta A\left(t_0\right)$	$\Delta A\left(t_1\right)$	CPU time		$\Delta A\left(t_0\right)$	$\Delta A\left(t_1\right)$	CPU time
SL-RLS	(80.2, 100)	(100, 100)	5.32		(24.9, 100)	(51.9, 100)	24.68		(7.6, 98.8)	(15.6, 100)	134.88
HR-RLS	(72.4, 100)	(72.8, 100)	6.44		(24.5, 10.5)	(41.0, 59.8)	44.52		(2.3, 15.3)	(5.2, 17.6)	147.30

Equations

$$\max {\mathcal{N}}_{\omega }=\min (2^d,\mathcal{N}).$$

$$V_s=\{C_{k,l}|i-1<k<i+1,j-1<l<j+1\},$$

$$V_c=\{C_{k,l}|i-t<k<i+t,l=j\}\cup \{C_{k,l}|k=i,j-t<l<j+t\},$$

$${\varphi }_{k,l}^r={\mathbf{C}}_r\left({\phi }_{k,l}^{\chi }\right)=\{\begin{array}{l}{\phi }_{k,l}\mathrm{if}\text{leavevmode}{\chi }_{k,l}={\chi }_r\\ -{\phi }_{k,l}\mathrm{otherwise}\end{array},$$

$$\mathbf{n}=\frac{\nabla {\varphi }^r}{|\nabla {\varphi }^r|},$$

$$\kappa =\nabla \cdot \frac{\nabla {\varphi }^r}{|\nabla {\varphi }^r|}.$$

$${\varphi }_{k,l}^r={\mathbf{C}}_r^{\star }\left({\phi }_{k,l}^{\chi }\right)=\{\begin{array}{l}{\mathbf{C}}_r\left({\phi }_{k,l}^{\chi }\right)\mathrm{if}\text{leavevmode}{\chi }_{k,l}={\chi }_r\\ {\mathbf{C}}_r\left(\min \left(\Delta x,d_s(k,l)\right)\right)\mathrm{otherwise}\end{array},$$

$${\phi }^{\chi }(x,y)=\mathbf{R}\left({\varphi }^r(x,y),r\in {\mathrm{X}}_s\right)=\left(\left|\max {\varphi }^r(\mathbf{x})\right|,\arg \underset{{\chi }_r}{\max }{\varphi }^r(\mathbf{x})\right).$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{v}\cdot \nabla \varphi =0,$$

$${\varphi }_{i,j}^{1,(s+1)}={\alpha }_s{\varphi }_{i,j}^{1,(s)}+(1-{\alpha }_s)\left[{\varphi }_{i,j}^{1,(s)}-\Delta t{\mathbf{v}}_{i,j}^n\cdot {\left(\nabla {\varphi }^1\right)}_{i,j}^{(s)}\right],s=0,\dots ,m,$$

$${\varphi }_{i,j}^{1,(0)}={\varphi }_{i,j}^{1,n},{\varphi }_{i,j}^{1,(m+1)}={\varphi }_{i,j}^{1,n+1}$$

$${\phi }_{i,j}^{\chi ,(s+1)}=\{\begin{array}{l}({\varphi }_{i,j}^{1,(s+1)},{\chi }_1)\mathrm{if}\text{leavevmode}{\varphi }_{\mathrm{i},\mathrm{j}}^{1,(\mathrm{s}+1)}\ge 0\\ (-{\varphi }_{i,j}^{1,(s+1)},{\chi }_2)\mathrm{otherwise}\end{array}.$$

$${\varphi }_{i,j}^{m,(s+1)}={\alpha }_s{\varphi }_{i,j}^{m,(s)}+(1-{\alpha }_s)\left[{\varphi }_{i,j}^{m,(s)}-\Delta t{\mathbf{v}}_{i,j}^n\cdot {\left(\nabla {\varphi }^m\right)}_{i,j}^{(s)}\right],{\chi }_m\in V_s,$$

$${\epsilon }_1={\int }_0^T\Vert \phi {\Vert }_1\mathrm{d}x,{\epsilon }_{\infty }={\int }_0^T\Vert \phi {\Vert }_{\infty }\mathrm{d}x,{\epsilon }_d={\int }_0^T{d}_H({\Gamma }^n,{\Gamma }^e)\mathrm{d}x,$$

$$\Vert \phi {\Vert }_1$$

$$B$$

$$d_H({\Gamma }^n,{\Gamma }^e)=\max (\underset{\text{y}\in {\Gamma }^e}{\sup }\underset{\text{x}\in {\Gamma }^n}{\inf }\Vert \text{x}-\text{y}{\Vert }_2,\underset{\text{x}\in {\Gamma }^e}{\sup }\underset{\text{y}\in {\Gamma }^n}{\inf }\Vert \text{x}-\text{y}{\Vert }_2).$$

$$B=\{(i,j)|{\phi }_{i,j}^n<1.2\Delta x\cap \Vert {\text{x}}_{i,j}-{\text{x}}_s{\Vert }_2<0.05L\},$$

$$\{\begin{array}{rl}& u=-2{\sin }^2\left(\pi x\right)\sin \left(2\pi y\right)\\ & v=2{\sin }^2\left(\pi y\right)\sin \left(2\pi x\right)\end{array}.$$

$${\mathbf{u}}_{{\Gamma }_{ab}}={\mathbf{n}}_b,{\mathbf{u}}_{{\Gamma }_{bc}}={\mathbf{n}}_c,{\mathbf{u}}_{{\Gamma }_{ca}}={\mathbf{n}}_a,$$

$$\frac{dA}{dt}=2\pi \gamma (\frac{n}{6}-1),$$

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot {\mathbf{F}}_c=\nabla \cdot {\mathbf{F}}_v+\sigma \kappa \delta \mathbf{n},$$

Full text

and

1 Introduction

Multi-region problems can occur when the motion of more than two immiscible fluids is to be described. In this case the interface network, separating the different fluid regions, evolves in time due to interactions of the different fluids across interface segments. These interactions often can be described by local fluid properties. Important applications include shock-driven multiphase flows [1, 2, 3], astrophysical events [4, 5, 6], foam dynamics [7, 8, 9, 10], multi-cellular tissue dynamics [11, 12, 13, 14, 15], and grain coarsening in polycrystalline materials [16, 17, 18, 19].

A range of numerical models have been proposed to compute the evolution of interface networks for multi-region problems. Generally, they can be classified as Lagrangian or Eulerian methods according to the representation of the interface. With Lagrangian methods, such as front-tracking [20], immersed-boundary [21], or arbitrary Lagrangian-Eulerian (ALE) [22] methods, the interface is represented explicitly by conforming discretization elements. Although these methods have been extended to multi-region systems [23, 24, 25, 26, 27, 10], it is difficult to handle complex topological changes during interface-network evolution, especially in three dimensions. With Eulerian methods, such as volume-of-fluid (VOF) [28] and level-set methods [29], the interface is reconstructed from scalar fields, i.e., volume fraction or level-set field. The interface is represented implicitly and captured by solving the corresponding transport equations. These methods generally can handle complex interface evolution with topological changes and are straightforward to implement in three dimensions. However, they often exhibit low accuracy due to numerical dissipation introduced by the transport-equation discretization.

Two additional problems are encountered when Eulerian methods are applied to multi-region problems [30, 31, 32]. One is that the number of scalar fields increases with the number of regions and entails additional computational operations and memory cost [33]. The other is that the interface reconstruction can produce voids and overlaps where more than two regions meet [31, 33, 32]. There are two main level-set-based approaches for capturing the evolution of multiple regions. One is to define multiple level-set functions (referred to as multiple level-set method in this paper) and to assign these to corresponding regions [30, 31, 32], followed by solving separate level-set transport equations for each region. Different numerical procedures may be employed to prevent the generation of voids and overlaps during the interface reconstruction. For example, Starinshak et al. [32] propose $N(N-1)/2$ level-set functions to represent all interfaces of $N$ regions and an additional pairwise voting strategy to remove overlaps and voids. The algorithm copes with the interface reconstruction problem but requires that a larger number of different level-set fields are stored and evolved in time than with multiple level-set methods [30, 31]. The method of Vese and Chan [34] reduces the number of level-set fields from $N$ to $\log_{2}N$ [33] and it avoids the generation of voids and overlaps. Naturally, using a single level-set function to represent an arbitrary number of regions is the optimal strategy to address the memory overhead [33]. Lie et al. [35] and Chung and Vese [36] develop such methods for image segmentation. Also, the regional level-set method [37] addresses the problem of multiple level-set fields [38] by employing a combination of a single unsigned level-set field and an integer region indicator function. Interface reconstruction for more than two regions, is handled by employing a low-order Semi-Lagrangian scheme. We note that this low-order Semi-Lagrangian scheme is more dissipative than high-order finite-difference schemes. As further development of the regional level-set method, the Voronoi implicit interface method [39] uses a transported $\epsilon$-level-set and reconstructs the interface network by a reinitialization step based on Voronoi diagrams.

Our objective is to develop an improved regional level-set method that has computational cost comparable to Semi-Lagrangian regional level-set methods [37, 38], but significantly improves prediction accuracy. Specifically, the method inherits the advantages of the original level-set method and the Semi-Lagrangian regional level-set method but is modified in such a way that it is suitable for high-resolution finite-difference discretization of the level-set transport equations. We define local level-set fields to capture the evolution of the interface network by a simple construction operator, followed by a reconstruction operator to obtain the global regional level-set field. The multi-region system and the definition of the regional level-set field are revisited in Sec. 2. The proposed multi-region method is detailed in Sec. 3. Accuracy and robustness are assessed in Sec. 4, followed by a brief conclusion in Sec. 5.

2 Multi-region system and regional level set

First we introduce the representation of a multi-region system by implicit functions. Let the domain $\Omega$ be an open set in $\mathbb{R}^{d}$, let $\mathbf{x}\in\Omega$ be an interior point of $\Omega$, and let $\partial\Omega$ be the boundary of the domain, where $d$ is the spatial dimension. Assuming that there are $\mathcal{N}$ regions within this domain, $\Omega$ is the union of a family of disjoint subsets, $\Omega=\underset{\chi\in\mathrm{X}}{\bigcup}\Omega^{\chi}$, where $\Omega^{\chi}$ is the subdomain of the region $\chi$, and $\mathrm{X}=\{\chi\in\mathbb{N}|1\leq\chi\leq\mathcal{N}\}$ is the index set for all regions. The entire multi-region system consists of the following elements:

$\bullet$ The region domain $\Omega^{\chi}$ which contains all interior points of the region $\chi$.

$\bullet$ The region boundary $\partial\Omega^{\chi}$ which contains the set of boundary points of $\Omega^{\chi}$.

$\bullet$ The pairwise interfaces $\Gamma_{\alpha\beta}=\partial\Omega^{\alpha}\cap\partial\Omega^{\beta}$ ($\alpha,\beta\in\mathrm{X}$ and $\alpha\neq\beta$), each as codimension manifold in $\mathbb{R}^{d}$ that separates two connected regions.

$\bullet$ The interface network $\Gamma=\underset{\alpha,\beta\in\mathrm{X},\alpha\neq\beta}{\bigcup}\Gamma_{\alpha\beta}$ is the union of all pairwise interfaces.

$\bullet$ Multiple junctions (high order junctions) $\mathrm{J}=\underset{\alpha,\beta\in\mathrm{X},\alpha\neq\beta}{\bigcap}\Gamma_{\alpha\beta}$ are intersections of pairwise interfaces. In most practically relevant applications they are triple points (2D) or triple lines (3D).

If there are only two regions ($\mathcal{N}=2$) in the system, the above elements can be implicitly defined by a level-set field (or signed distance function) [29] $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}$. The two region domains, $\Omega^{1}$ (or $\Omega^{+}$) and $\Omega^{2}$ (or $\Omega^{-}$), are identified by the sign of the level-set function: $\Omega^{1}=\{\mathbf{x}\in\mathbb{R}^{d}|\phi(\mathbf{x})>0\}$ and $\Omega^{2}=\{\mathbf{x}\in\mathbb{R}^{d}|\phi(\mathbf{x})<0\}$. The region boundaries of $\Omega^{1}$ and $\Omega^{2}$ coincide and are identical to the pairwise interface and interface network, $\partial\Omega^{1}=\partial\Omega^{2}=\Gamma_{\alpha\beta}=\Gamma=\{\mathbf{x}\in\mathbb{R}^{d}|\phi(\mathbf{x})=0\}$. There are no multiple junctions in this system. When more than two regions ($\mathcal{N}>2$) exist, as shown in Fig.1(a), the regional level-set method [37, 39] can be used to represent the system implicitly. For the regional level-set method, a mapping $\varphi^{\chi}:\Omega\subset\mathbb{R}^{d}\rightarrow\mathbb{R}\times\mathbb{N}$ is defined as $\varphi^{\chi}(\mathbf{x})=(\varphi(\mathbf{x}),\chi(\mathbf{x}))$, where $\varphi(\mathbf{x})\geq 0$ is the unsigned distance function and $\chi(\mathbf{x})$ is a positive integer region indicator. A subdomain $\Omega^{\alpha}$ is identified by region indicators, $\Omega^{\alpha}=\{\mathbf{x}\in\Omega|\chi(\mathbf{x})=\alpha\}$ and the interface network is defined as the zero distance contour $\Gamma=\{\mathbf{x}\in\mathbb{R}^{d}|\varphi(\mathbf{x})=0\}$. However, the region boundary, the pairwise interface and multiple junctions can not be identified directly by regional level-set fields. An additional operation, such as the Voronoi method described in [39] is required for this purpose. In this paper we propose as an alternative a simple construction operator to define these surface contours, see Sec. 3.

3 Numerical Method

3.1 Global and local index sets of regions

In this section, we develop the numerical method for the representation of the multi-region system and the evolution of the interface network. For this purpose we consider the multi-region system defined by $\varphi^{\chi}_{i,j}=(\varphi_{i,j},{\chi}_{i,j})$ at the centers of finite-volume cells with indices $i$ in $x_{1}$ and $j$ in $x_{2}$ coordinate directions on a two-dimensional uniform Cartesian grid, as shown in Fig. 1(a). An extension of the two-dimensional definitions to three dimensions is straightforward. We recall that globally there are $\mathcal{N}$ regions in the considered domain. Although $\mathcal{N}$ may be very large, locally the number of regions which occur in a neighborhood of a cell is limited to a small integer number. Consequently, in this local small subdomain $\omega$, the complexity of the system is significantly reduced, with the maximum possible number of regions $\mathcal{N}_{\omega}$ being

[TABLE]

Interface networks typically exhibit only very few finite-volume cells which contain multiple junctions. Thus, multiple signed level-set functions defined only locally to capture the evolution of the interface network allow to reduce the computational effort significantly. The required global and local definitions on a two-dimensional uniform Cartesian grid are as follows:

$\bullet$ $\mathcal{N}$ global regions are labelled by the global index set $\mathrm{X}$. The center $\mathbf{x}_{i,j}$ of each finite-volume cell $C_{i,j}$ is assigned to one region by an indicator $\chi_{i,j}$ even though $C_{i,j}$ may contain more than one region. Here $i,j$ are the global indices for cells.

$\bullet$ A local subdomain is defined as $\omega\subset\Omega$. The $\mathcal{N}_{\omega}$ regions contained in this local subdomain are labelled by a local index set $\mathrm{X}_{\omega}=\{r\in\mathbb{N}|1\leq r\leq\mathcal{N}_{\omega}\}$. We map the local region indicators into the global index set $\mathrm{X}$, and identify the corresponding global region indicator as $\chi_{r}$. A local subdomain $\Omega^{\alpha}_{\omega}$ of the global region $\Omega^{\alpha}$ is defined by $\Omega^{\alpha}_{\omega}=\{\mathbf{x}\in\omega|\chi_{r}(\mathbf{x})=\alpha,r\in\mathrm{X}_{\omega}\}\subset\Omega^{\alpha}$. Here, $\omega$ can be a finite-volume cell $C_{i,j}$ or a neighborhood, i.e., a set of cells $C_{k,l}$) of $C_{i,j}$, which will be defined in Sec. 3.2. $(k,l)$ denotes the local index pair for cells in such a neighborhood of $C_{i,j}$.

3.2 Cell neighborhoods and cell types

For a finite-volume cell $C_{i,j}$ one can define two types of neighborhoods, i.e. cell sets: a square-shaped near neighborhood

[TABLE]

and a cross-shaped stencil neighborhood

[TABLE]

where $2t+1$ is the required width of the stencil for the level-set transport discretization schemes. The number of separate regions $\mathcal{N}_{C}$ contained in cell ${C_{i,j}}$ and the number of region indicators in the square-shaped near neighborhood $\mathcal{N}_{s}$ have the relation $1\leq\mathcal{N}_{C}\leq\mathcal{N}_{s}\ll\mathcal{N}$ when $\mathcal{N}$ is large. Accordingly, the local index sets $\mathrm{X}_{\omega}$ of Sec. 3.1 specifically for $V_{s}$ and $V_{c}$ are $\mathrm{X}_{s}=\{r\in\mathbb{N}|1\leq r\leq\mathcal{N}_{s}\}$ and $\mathrm{X}_{c}=\{r\in\mathbb{N}|1\leq r\leq\mathcal{N}_{c}\}$, respectively. Note that $V_{c}$, $V_{s}$, $\mathcal{N}_{C}$ and $\mathcal{N}_{s}$ are defined for each cell $C_{i,j}$ and thus depend on the index pair $(i,j)$. We identify $\chi_{1}=\chi_{i,j}$ as “primary indicator”, and other region indicators $\chi_{r}$, $r\in\mathrm{X}_{s},r>1$, if they exist, are called “secondary indicators”.

A cell is denoted as a “full cell” if it is not intersected by an interface, i.e. $\mathcal{N}_{C}=1$, as a “two-region cut cell” if it is intersected by an interface segment shared by two regions, i.e. $\mathcal{N}_{C}=2$, or as a “complex-region cut cell” if the interface segment is shared by more than two regions, i.e. $\mathcal{N}_{C}>2$, see Fig.1(a). For details on identifying cell types, please refer to Algorithm 1. Analogously, according to the value of $\mathcal{N}_{s}$, a cell $C_{i,j}$ is categorized as a cell with full type $V_{s}$ ($\mathcal{N}_{s}=1$), as a cell with two-region type $V_{s}$ ($\mathcal{N}_{s}=2$), or as a cell with complex-region type $V_{s}$ ($\mathcal{N}_{s}>2$). We clarify that the square-shaped near neighborhood $V_{s}$ serves as a search stencil to identify all local regions that may share the cell $C_{i,j}$ and to generate the local index set $\mathrm{X}_{s}$. On the cross-shaped stencil neighborhood $V_{c}$, $\mathcal{N}_{s}$ auxiliary local level-set fields (see Sec.3.3) corresponding to the local regions identified in $V_{s}$ are constructed in order to evolve the interface network by discretized evolution equations and represent the region boundaries from the zero level set.

3.3 Description of the local construction and reconstruction operators

For any cell $C_{k,l}$ in the neighborhood of the cell $C_{i,j}$, we can apply a construction operator for generating the local multiple signed level-set fields and a reconstruction operator for reconstructing the global regional level-set field from the local level-set fields, see Fig. 2.

3.3.1 The local construction operator

The construction operator is used to generate the $\mathcal{N}_{s}$ local level-set fields $\phi^{r}_{k,l}:=\phi_{k,l}(\chi_{r})$ from the regional level-set field $\varphi^{\chi}_{k,l}$ at the center of cell $C_{k,l}$ for each region indicator $\chi_{r}\in\mathrm{X}_{s}$. The construction operator $\mathbf{C}_{r}:\mathbb{R}\times\mathbb{N}\rightarrow\mathbb{R}$ is defined as

[TABLE]

and the local level-set function $\phi^{r}(\mathbf{x}):=\phi(\mathbf{x},\chi_{r})$ at any point can be obtained by interpolation.

Upon construction of the local multiple signed level-set fields, the normal direction and curvature are obtained at the finite-volume cells,

[TABLE]

[TABLE]

Unlike with the original level-set method, the normal direction always points away from the interface network. The local construction operator $\mathbf{C}_{r}$ may be inaccurate for curvature calculation when the signed level set is not strictly a distance function, which may happen near multiple junctions, as demonstrated in 4.1.4. This issue can be handled by a construction operator $\mathbf{C}_{r}^{\star}$ on $V_{s}$ which additionally invokes a re-initialization procedure for recovering the signed-distance property of the local multiple signed level-set fields

[TABLE]

where $d_{s}(k,l)$ is the distance of the cell center $\mathbf{x}_{k,l}$ from the local interface segment in $V_{s}$.

Note that discrete derivatives in the local level-set advection equation Eq. (3.4.1), see below, which operate on constructed local level-set functions in the stencil neighborhood $V_{c}$ now can be calculated by high-resolution spatial discretization schemes.

3.3.2 The local reconstruction operator

The regional level-set field at an arbitrary point can be reconstructed from the multiple auxiliary local level sets for all region indicators in $V_{s}$ by the reconstruction operator $\mathbf{R}:\mathbb{R}^{\mathcal{N}_{s}}\rightarrow\mathbb{R}\times\mathbb{N}$ given by

[TABLE]

Note that if the interface network is static, application of the construction operation followed by the reconstruction operation leaves the original regional level-set field invariant. Consider the example in Fig. 2. If the regional level set at the central cell (colored by red) is $\varphi^{\chi}=(0.2,4)$, upon local construction, the $3$ generated local level-set data $\{\phi^{r}|r=1,2,3\}$ at this point are $\{0.2,-0.2,-0.2\}$. Upon applying the reconstruction, we obtain the original regional level-set data $\varphi^{\chi}=(0.2,4)$.

This reconstructed field is unique. The reconstruction generates topologies without artificial overlaps or voids. Thus it is suitable for determining the regional level set when the interface network has been evolved in time through updating multiple local level-set fields.

3.4 Evolution of the interface network

The evolution of the interface network in a multi-region system is equivalent to that of a signed level-set field whose interface advection is determined by the advection equation

[TABLE]

where $\mathbf{v}$ is the advection velocity. The advection equation in our method is formulated locally for each region, identified by the region indicator, and recovers the original level-set method for cells that are sufficiently far away from a multiple junction. More detailed considerations are necessary to predict the evolution of the interface for cells in which more than two regions meet. Consider the situation in Fig.1. The cell contains three regions colored by red, blue, and yellow. Initially, the largest fraction of the cell is occupied by region $\Omega^{\chi_{1}}$, with $\chi_{1}$ being the primary indicator. After advection by one time step, the interface has three possible configurations, as illustrated in Fig.1(a). The first is that the cell center still resides within $\Omega^{\chi_{1}}$, so that the primary indicator of the cell does not change. In the other two cases, the primary indicator changes to the secondary indicator, either $\chi_{2}$ or $\chi_{3}$. After solving an advection equation for the local level-set field of the three regions separately, each region boundary may shrink (dashed line) or expand (dash-dotted line), see Fig.1(b). For example, when the boundary $\partial\Omega^{\chi_{1}}$ of $\Omega^{\chi_{1}}$ moves away from the cell center its local level-set $\phi^{\chi_{1}}$ value increases (positive). For the converse case, $\partial\Omega^{\chi_{1}}$ moves across the cell center, so that it does not belong to $\Omega^{\chi_{1}}$ any longer, resulting in a corresponding sign change to $\phi^{\chi_{1}}<0$. That the cell center is not in $\Omega^{\chi_{1}}$ does not necessarily imply that it is in $\Omega^{\chi_{2}}$ or $\Omega^{\chi_{3}}$. A direct combination of all three independently advected region boundaries may introduce an overlap or a void. To address this issue, we apply the reconstruction operator $\mathbf{R}$ to determine the regional level set $\varphi^{\chi}$ from the three candidate local level-set fields $\{\phi^{\chi_{1}},\phi^{\chi_{2}},\phi^{\chi_{3}}\}$. This is physically reasonable because $\mathbf{R}$ identifies the most likely indicator corresponding to the region domain in which the current cell center will be located after one advection time step.

The update of the regional level set of each finite-volume cell $\varphi^{\chi,n}_{i,j}=(\varphi_{i,j},\chi_{1})^{n}$ at time-step $n$, by a sub-step of an explicit time-integration scheme, such as a strongly stable Runge-Kutta scheme [40], consists of three sub-steps: (1) construction of the local multiple signed level-set fields in the stencil neighborhood $V_{c}$; (2) computation of the new intermediate multiple local level-set fields by updating the locally constructed advection equations, see below; (3) reconstruction of the regional level set from the multiple local level-set fields. The local construction of the advection equations depends on the type of $V_{s}$ (full, two-region or complex-region). Note that we assume that under the Courant-Friedrichs-Lewy (CFL) condition, i.e. $\mathrm{CFL}\leq 1.0$, the type of $V_{s}$ does not change and the region indicators within the square-shaped near neighborhood remain unchanged during one sub-step of the time-integration scheme, whereas the actual finite-volume cell may change its type. Unlike the “repairing” (or “modification”) procedure in Ref. [30], no a posteriori operations are required after advection.

3.4.1 Updating a cell with full type $V_{s}$

For a cell with full type $V_{s}$, the local level set $\phi^{1,(n)}_{i,j}=\phi^{n}_{i,j}(\chi_{1})$ and the intermediate data $\phi^{1,(s)}_{k,l}=\phi^{(s)}_{k,l}(\chi_{1})$ in the stencil neighborhood are constructed referring to the primary indicator $\chi_{1}$. The intermediate-step update $\phi^{1,(s+1)}_{i,j}$ is obtained from

[TABLE]

where $m$ is the number of sub-steps, $\alpha_{s}$ is the parameter of Runge-Kutta sub-step $s$, and $\left(\nabla\phi^{1}\right)^{(s)}_{i,j}$ is the finite difference approximation of the spatial derivative at the center of a finite-volume cell $C_{i,j}$. The updated regional level set at the sub-step $s+1$ is $\varphi^{\chi,(s+1)}_{i,j}=(\phi^{1,(s+1)}_{i,j},\chi_{1})$.

3.4.2 Updating a cell with two-region type $V_{s}$

For a cell with two-region type $V_{s}$, the construction of the multiple local level-set fields and the regional level-set field updates are obtained essentially by the same operations as for a full cell. The difference is that the intermediate regional level-set field update depends on the sign of $\phi^{1,(s+1)}_{i,j}$

[TABLE]

3.4.3 Updating a cell with complex-region type $V_{s}$

For a cell with complex-region type $V_{s}$, the local level-set field is constructed and updated referring to the primary indicator $\chi_{1}$. If $\phi^{1,(s+1)}_{i,j}$ does not change sign, the new intermediate local level set is $\varphi^{\chi,(s+1)}_{i,j}=(\phi^{1,(s+1)}_{i,j},\chi_{1})$. Otherwise, the local level-set field is constructed and updated referring to the secondary indicators

[TABLE]

and the regional level set is reconstructed by the operator $\mathbf{R}$.

The reason for using different advection strategies based on the cell types is the following. A direct transport of the global unsigned level-set function which exhibits a discontinuity across the interface requires numerical diffusion for stabilization and thus produces a smeared interface. Consequently, the unsigned scalar function does not maintain the distance function property. Advection with a high-order, low dissipation scheme of the local signed level-set functions, however, maintains the sharp interface and the distance property much more accurately. We advect the unsigned level-set function for full type cells and the local constructed signed level-set functions for other cells. This essentially is an adaptive algorithm which first determines the cell type as the indicator of the local smoothness of the unsigned level-set. Then the unsigned level-set field is advected wherever the unsigned level-set is smooth or the signed level-set field is advected wherever the unsigned level-set is singular. The re-initialization is not subject to regularization constraints, and our numerical examples show that it can be chosen as infrequent as with the original level-set method.

Computational efficiency can be further increased by employing the narrow band technique [41]. For simple test cases, such as in Sec. 4.1.2 and Sec. 4.1.3, this is not necessary. For more complex cases, such as in Sec. 4.3 and Sec. 4.5, it is applied. The operation count of our method per time step is $\mathcal{O}(n^{2})$ (same as with the Semi-Lagrangian regional level-set method and original level-set method) and can be reduced to $\mathcal{O}(kn)$ by employing the narrow band technique [41], where $k$ is the band width and $n$ is the number of cells in any direction. For the entire system the construction and reconstruction operators are applied $d(2t+1)N_{1}+2d(2t+1)N_{2}+\mathcal{N}_{s}d(2t+1)N_{3}$ and $N_{3}$ times respectively, where $N_{1}$, $N_{2}$ and $N_{3}$ are the number of cells with full type, two-region type and complex-region type $V_{s}$ in the narrow band of interface network, respectively. At every Runge-Kutta sub-step we solve the local advection equations $N_{1}+N_{2}+\mathcal{N}_{s}N_{3}$ times.

3.5 Re-initialization

It is important to re-initialize the regional level set when necessary for maintaining its distance-function property with respect to the interface network. We employ two re-initialization methods. One widely used method is that of Sussman et al. [42] where a re-initialization equation is solved iteratively until steady-state is reached. The implementation for the regional level-set method involves the following steps region by region. For each subdomain $\Omega^{\chi}$, $\chi\in\mathrm{X}$, first the local level set $\phi^{r}$ is constructed and the re-initialization equation according to Ref. [42] is iterated until a steady state is reached. Alternatively, the explicit one-step method developed by Fu et al.[43] can be used which is significantly faster than iterative methods.

3.6 Summary of the numerical method

Here we summarize our high-resolution regional level-set method and comment on the description above. The entire numerical method contains the following steps:

1. The initialization step. According to the description in Sec. 3.1, we define an initial regional level set $\varphi^{\chi}_{i,j}=(\varphi_{i,j},{\chi}_{i,j})$ at each finite-volume cell $C_{i,j}$ center $\mathbf{x}_{i,j}$, where $\varphi_{i,j}$ is the distance from the interface network and ${\chi}_{i,j}$ indicates which region domain the center of $C_{i,j}$ is located at.

2. The evolution step contains three sub-steps for $C_{i,j}$: (1) construction of the $\mathcal{N}_{s}$ local level-set fields $\phi^{r,n}_{k,l}$ for the current time-step $n$ at the center of each cell $C_{k,l}$ which belongs to the stencil $V_{c}$ of $C_{i,j}$; (2) computation of $\phi^{r,n+1}_{i,j}$ at the next time-step $n+1$ by solving the $\mathcal{N}_{s}$ local advection equations; (3) reconstruction of the new regional level set $\varphi^{\chi}_{i,j}$ at the center of $C_{i,j}$ from the $\mathcal{N}_{s}$ new local level-set fields $\phi^{r,n+1}_{i,j}$ by the reconstruction operator $\mathbf{R}$.

3. The re-initialization step. Enforce the distance function property of the regional level set $\varphi^{\chi}_{i,j}$ as the distance from the interface network.We emphasize that the interface-network transport re-initialization not always is necessary, see Sec. 4.

4. The postprocessing step. If necessary, extract the interface network by the triangulation method of Ref. [44].

It should be mentioned that except for a postprocessing step the method does not need to extract explicitly the interface network, so it is a fully implicit method for multi-region problems. Another important feature of our method is that unlike using neighbouring $\epsilon$-level-set contours to reconstruct the interface network [39], our method directly captures the evolution of the interface network. The method is not a hybrid of the multiple level-set method and the regional level-set method. It rather can be viewed as a regional level-set method employing locally signed level-set fields.

4 Numerical validation

In this section, we assess accuracy and efficiency, the present method by a range of numerical examples. Our intention is to show that we recover the high computational efficiency of the Semi-Lagrangian (SL) method, improves, however, significantly the prediction accuracy. First, we compare different high-order discretizations of the level-set transport equation with SL results for constant rotation motion of a three-region case. Afterwards, two simple mean curvature flows are considered to verify the construction operators and re-initialization methods. Subsequently, the single vortex flow is used to demonstrate the ability of the present method to resolve long thin filaments. The computations of normal driven flow, mean curvature flow and their combination serve to assess the accuracy of the present method. Finally, we couple the present method with Navier-Stokes equations applied to dry-foam dynamics which undergo sudden breakups, in order to demonstrate the capability for coping with complex configurations.

4.1 Simple test cases

Three simple cases are considered to test suitable high-resolution finite-difference schemes, construction operators, and re-initialization methods for the simulation of multi-region problems. We consider the 5th-order weighted essentially non-oscillatory scheme (WENO) [45] and the central-upwind weighted essentially non-oscillatory scheme (WENO$\_$CU6) [46]. Both are compared with the SL scheme of the regional level-set method [37].

4.1.1 Circle expansion

We start with a 2-region expansion case in Ref. [44]. A circle with radius of $0.2$ expands with a uniform speed to $t=0.2$. The computational domain is $\left[0,1\right]\times\left[0,1\right]$. The explicit Euler scheme is used for time marching with a CFL number of $0.5$.The 5th-order WENO scheme is employed for spatial discretization. Fig. 3(a) shows three error measures

[TABLE]

with

[TABLE]

The Hausdorff distance is

[TABLE]

The superscripts ‘e’ and ‘n’ stand for the exact and the numerical solutions. Our method achieves the expected 5th-order convergence rate for all three error measures. The SL method is 1st-order, as expected, as is the Voronoi implicit interface method, see Fig. 4 of Ref. [44]. Fig. 3(b) shows the interfaces of our method and SL method at $t=0.2$, compared the exact solution.

4.1.2 Triple point advection

We assess the prediction accuracy of our method by triple point advection cases which have analytical solutions. The first configuration contains only one triple point which initially is located at $(0.2,0.5)$. The second is a circle of radius $r_{0}=0.3$ divided into two parts. The computational domain is $\left[0,1\right]\times\left[0,1\right]$. The velocity field is given by $(u,v)=(1.0,0.0)$. The explicit Euler scheme is used for time marching with a CFL number of $0.6$. The 5th-order WENO scheme is used for spatial discretization. The simulations are performed until time $t=0.4$. Re-initialization is not employed. Error measures of Eq. (13) are computed within a narrow band

[TABLE]

where $\textbf{x}_{i,j}$ and $\textbf{x}_{s}$ are the locations of cell center and triple points, respectively. Due to the smoothness properties of the exact solution at most first order accuracy in $\varepsilon_{d}$ can be achieved. This is reproduced by our method, as shown in Figs. 4(a) and 4(b). The benefit of the high-resolution discretization becomes evident for the global measures $\varepsilon_{1}$ and $\varepsilon_{\infty}$. The interface recovers the exact solution as the resolution increases. A comparison with the SL regional level-set method, shown in Fig. 4(b), demonstrates that the error magnitude of the SL regional level-set method is significantly larger than that of our method.

4.1.3 Constant rotation of three regions

A two-dimensional circle of radius $r_{0}=0.3$ is divided into two equal parts which undergo constant rotation. The computational domain is $\left[0,1\right]\times\left[0,1\right]$. Symmetry conditions are employed at all boundaries. The grid spacing is $h=\frac{1}{64}$. The explicit Euler scheme is used for time marching with a CFL number of 0.6. Two high-resolution schemes (5th-order WENO and WENO$\_$CU6), see the last two rows of Fig.5, are compared with respect to their ability to capture the two triple points, in comparison with the results of the 1st-order SL scheme, see the first row of Fig.5. Three regional level-set contours are shown in Fig.5 at $t=\frac{1}{8}\pi$, $t=\frac{1}{4}\pi$, and $t=\frac{1}{2}\pi$. Apparently, in each simulation, the deflection angles of the separation interface at all times agree with the theoretical values. However, the level-set contours surrounding the interface alter slightly. The SL method results in more smeared contours near the triple points due to numerical dissipation. The two high-resolution schemes have less numerical dissipation and thus preserve sharp corners. The relative radius difference $|\triangle r|/r_{0}$ of the circle for the WENO and WENO$\_$CU6 results at $t=0.5\pi$ are $0.211\%$ and $0.285\%$, respectively, and are much smaller than for the SL result which is $3.42\%$, indicating better area conservation with WENO and WENO$\_$CU6. An error analysis for this case, shown in Fig. 6(a), indicates 1st-order convergence of the interface location. The extracted interface network converges to the exact solution asymptotically as the resolution increases, as shown in Fig. 6(b).

4.1.4 Two dimensional mean curvature flows

The interface network of three-region and five-region systems is evolved under mean curvature flow with $\mathbf{u}=\kappa\mathbf{n}$, where $\mathbf{n}$ and $\kappa$ are calculated by Eqs. (5) and (6). In our simulation we employ the 5th-order WENO scheme for advection and consider different construction operators and re-initialization methods. Since the mean curvature uses 2nd-order derivatives of $\phi^{r}$, it is more sensitive to the local constructed level-set fields. The computational domain extent is $\left[0,1\right]\times\left[0,1\right]$ in x and y directions. Symmetry boundary conditions are employed at all domain boundaries. The explicit Euler scheme is used for time marching, and the time step $\Delta t$ is the same as that in [44], $\Delta t=\frac{h^{2}}{4}$. As depicted in Fig.7, the T junction transforms to a Y junction under the effect of mean curvature. The different columns in this figure show results for different re-initialization methods and construction operators. When we use the construction operator $\mathbf{C}_{r}$, the displacement of the triple point is overestimated in both cases, see Fig.7(b). Significant improvement is observed by the operator $\mathbf{C}_{r}^{\star}$ where $\chi(\mathbf{x})\neq\chi_{r}$, irrespective of the employed re-initialization method, see Figs. 7(c) and 7(d). Note that noise in the contours far away from the interface visible for the explicit re-initialization is irrelevant to our method as only the smooth inner contours are used for extracting the interface.

We conclude that the construction operator $\mathbf{C}_{r}^{\star}$ is more suitable for mean curvature flows. Both re-initialization methods are suitable for capturing the interface network. We obtain convergened interface locations, as shown in Fig. 8(b), where the interface networks are extracted by the triangulation method of Ref. [44]. As shown in Fig. 8(b), first order convergence is achieved for the triple-point locations.

4.2 Single vortex flow

For the single vortex flow case, we use the setup introduced by Bell et al. [47] to test the ability and accuracy on resolving thin filaments under the deformation by the velocity field:

[TABLE]

Again, we employ the 5th-order WENO scheme for advection. The initial circle deforms into a filament wrapping around the center of the domain. This structure wraps back into the initial circle upon reversing the velocity at $t=3$. The centers of concentric circles are $(0.5,0.75)$ and their radii are $0.08$ and $0.22$, respectively. Symmetric boundary conditions and a second-order strongly stable Runge-Kutta scheme [48] with a CFL number of 0.6 are employed. The grid size is refined from $h=\frac{1}{256}$ to $h=\frac{1}{1024}$.

It can be seen from Fig.9 that two spirals of long filaments are successfully captured at $t=1.5$ and $3.0$. With increasing grid resolution the filaments become longer, especially the inner one. At $t=6.0$ with $h=\frac{1}{256}$ two additional unresolvable triple points are generated such that the two initial circle contours are connected. This phenomenon can be attributed to the fact that for extremely stretched filaments the two interfaces may reside in the same cell, and a small numerical perturbation leads to a topology change while reverse rotation does not disconnect the interface. We emphasize that although we use $256\times 256$ cells overall, the number of cells across the initial inner circle is only $40\times 40$. Increasing mesh resolution removes such artifacts, as shown in Figs.9(b) and (c). Both circles recover their initial shapes with good area conservation.

In Fig. 10, we compare the results of our method with that of the SL regional level-set method. We observe that the SL method shows large numerical dissipation and can not reproduce the long filament. The area conservation errors are listed in Table. 1. Note that the inner circle is very poorly resolved, i.e., for $h=\frac{1}{64}$, we have only $12\times 12$ cells across the initial inner circle. With such low resolution $100\%$ mass loss of the inner circle is inevitable for any method. At all three grid resolutions, our method exhibits significantly smaller area loss compared to the SL method. The CPU time measurement indicates that our method is nearly as fast as the SL regional level-set method even though we employ the high-order spatial schemes. We conclude that our approach achieves improved accuracy and comparable computational cost compared with the SL regional level-set method.

4.3 Normal driven flow

The constant normal driven flow has been studied in [44]. The pairwise $\Gamma_{ab}$ interface separating region domains $\Omega^{a}$ and $\Omega^{b}$ moves in its normal direction with a constant speed. As shown in Fig.11, two initially neighboring region domains $\Omega^{b}$ (colored by green) and $\Omega^{c}$ (colored by blue) are surrounded by a background region domain $\Omega^{a}$ (colored by white). We define the normal velocity of the interface of each region as

[TABLE]

where $\mathbf{u}_{\Gamma_{ab}}$ is the velocity at the pairwise interface $\Gamma_{ab}$ and $\mathbf{n}_{b}$ is the normal direction of the region boundary $\partial\Omega_{b}$ at $\Gamma_{ab}$. Thus the flow wraps $\Omega^{a}$ into $\Omega^{b}$ which in turn is wrapped into $\Omega^{c}$. The computation is carried out on a unit square with different resolutions employing the 5th-order WENO schemes for advection. Symmetry boundary conditions and the explicit Euler scheme are employed with a CFL number of $0.6$.

For a comparison between the results of our method and those in [44] we plot five snapshots at the same time instants as those of [44], from $t=0$ to $t=0.288$. During evolution, the number of spirals increases quickly, and we observe that our results actually exhibit more visible spirals compared to that in [44], both on the same $256\times 256$ grid, see third row of Fig.11. This observation can be attributed to the high-resolution scheme employed in our method. Another reason is that we directly advect the interface network. The number of spirals increases linearly with time and proportionally to grid resolution, see Figs 13 and 13, respectively. In order to demonstrate the validity of our method in three dimensions, a 3D example of normal driven flow on a $128\times 128\times 128$ grid, is shown in Fig. 12.

4.4 Von Neumann-Mullins’ law validation

We consider a case involving more regions and which serves to verify our method in a configuration with multiple triple points within a multi-region system. Initially $15$ regions are randomly placed in a domain $\left[0,1\right]\times\left[0,1\right]$ and evolve under a mean curvature generated velocity $\mathbf{u}=\kappa\mathbf{n}$. Periodic boundary conditions are imposed at the domain boundaries. The explicit Euler scheme is employed for temporal discretization. The time-step size is $\Delta t=\frac{h^{2}}{4}$ and the grid spacing is $h=\frac{1}{128}$. We want to verify that the von Neumann-Mullins’ law is reproduced, which states that the rate of area $A$ growth or decay is a function of the number of edges $n$ of the phase [44]

[TABLE]

where we set $\gamma=1.0$. The initial number of regions $15$ is successively reduced to $6$ under the effect of mean curvature, as shown in Fig.14. According to Fig.14, region “$h$” initially has four edges and shrinks under mean curvature, leading to region destruction, consistent with von Neumann-Mullins’ law. More specifically, the temporal growth and decay of the selected $8$ regions is shown in Fig. 15 which exhibits a piecewise linear profile for each region, in agreement with von Neumann-Mullins’ law which is indicated by a colored line in Fig. 15.

4.5 Dry foam dynamics

As demonstration for a complex application we couple the high-resolution regional level-set method with the Navier-Stokes equation Eq. (20) to simulate dry-foam dynamics.

[TABLE]

where $\mathbf{U}=(\rho$, $\rho u$, $\rho v$, $\rho w$, $\rho E)^{T}$, in which $\rho$, $u$, $v$, $\rho w$, and $\rho E$ are the density, the three velocity components, and the total energy, respectively. $\mathbf{F}_{c}$ and ${\mathbf{F}}_{v}$ are the convective and viscous flux tensor, respectively. The surface tension term $\sigma\kappa\delta\mathbf{n}$ on the right-hand side describes foam dynamics, with $\sigma$ being the surface tension and $\delta$ being a smoothed Dirac delta function. The surface tension force is calculated by Eq. (14) of Ref. [44]. The Ohnesorge number is $\mathrm{Oh}=\mu/\sqrt{\rho\sigma L}=0.00024$, i.e., the flow is surface-tension dominated. We note that the above surface-tension model may be imroved by Eq. (16) in Ref. [44], which is, however, not our main concern. Symmetry boundary conditions are applied. We use the 5th-order WENO scheme and a 4th-order central scheme for spatial discretization and a second-order strongly stable Runge-Kutta scheme for temporal discretization, both for solving N-S equation and level-set advection. The CFL number is $0.6$ and the grid size is $h=\frac{1}{256}$.

In a dry-foam cluster the gas bubbles are separated by thin liquid films which corresponds to the interface network of our method. As we are not concerned with a specific physical problem, inertia and gravity effects of liquid membranes, gas exchange across permeable membranes, and Marangoni forces at the liquid-gas interface are neglected here for simplicity. We consider the interconnected membranes as massless and infinitely thin. Initially, four bubbles evolve from an artificial initial configuration to an equilibrium state, as shown in the first row of Fig.16. Following the breakup of $\Omega^{a}$, a second equilibrium state is attained under the surface tension force, see the second row of Fig.16. Then $\Omega^{b}$ breaks up, and a similar process is observed for the remaining two bubbles. Note that the bubble breakup is triggered explicitly, which is unphysical, it serves, however, our purpose of demonstrating that our method has the capability to capture the interface evolution of a foam cluster subjected to stimulated bubble breakup.

5 Concluding remarks

The proposed method employs locally constructed level-set fields and the regional level-set method to overcome typical problems encountered with the numerical simulation of multi-region problems. The proposed numerical approach for multi-region problems and its algorithmic formulation have the following main properties: (1) As the proposed regional level-set method employs local signed level-set fields generated from the regional level-set field with a simple construction operator, it permits the implementation of high-resolution schemes for level-set transport in a straightforward way. (2) Instead of explicitly constructing the interface at every time step, we use a reconstruction operator to assemble the regional level set from multiple local level-set fields. This way we can ensure that the implicitly defined topology has no void or overlap artifacts. Moreover, a simple algorithm allows to distinguish different types of cells. It can be concluded from a range of test cases that the proposed method is more accurate than the Semi-Lagrangian regional level-set method. High order accuracy is demonstrated for some simple test cases where analytical results are known. Several increasingly complex configurations serve to demonstrate that improved accuracy and efficiency transfer to such test cases. A region-deconstruction example in foam dynamics demonstrates the feasibility of the method for complex applications.

Acknowledgment

This work is partially supported by China Scholarship Council (No. 201306290030), National Natural Science Foundation of China (No. 11628206) and Deutsche Forschungsgemeinschaft (HU 1527/6-1). The project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (No. 667483).