Simplex space-time meshes in two-phase flow simulations

TL;DR

This paper introduces a space-time finite element method with local temporal refinement for two-phase flow simulations, enabling high-resolution modeling of interfaces with complex local effects.

Contribution

It extends existing methods to include adaptive temporal refinement based on interface position, creating tetrahedral and pentatope meshes for improved accuracy and efficiency.

Findings

Effective local temporal refinement improves interface resolution.

Method handles complex effects like surface tension and property jumps.

Enhanced computational efficiency in two-phase flow simulations.

Abstract

In this paper, we present the numerical solution of two-phase flow problems of engineering significance with a space-time finite element method that allows for local temporal refinement. Our basis is the method presented in [3], which allows for arbitrary temporal refinement in preselected regions of the mesh. It has been extended to adaptive temporal refinement that is governed by a quantity that is part of the solution process, namely, the interface position in two-phase flow. Due to local effects such as surface tension, jumps in material properties, etc., the interface can, in general, be considered a region that requires high flexibility and high resolution, both in space and in time. The new method, which leads to tetrahedral (for 2D problems) and pentatope (for 3D problems) meshes, offers an efficient yet accurate approach to the underlying two-phase flow problems.

Tables

Table 1. Table 1: Typical performance behavior of the prism- and simplex-based calculations.

	Time	Nodes	Elements	System	System
	Steps	per step	per step	formation ($s$)	Solution ($s$)
Prism, $\Delta t=0.001s$	1600	730	624	331.65	5317.04
Prism, $\Delta t=0.005s$	320	730	624	77.77	1260.55
Simplex, $\Delta t=0.001-0.005s$	320	$\sim 919$	$\sim 3003$	142.21	1874.36

Equations

$${\rho }_i\left(\frac{\partial u}{\partial t}+u\cdot \nabla u-f\right)-\nabla \cdot {\mathbf{\sigma }}_i=0\mathrm{in}{\Omega }_i(t),$$

$$\nabla \cdot u=0\mathrm{in}{\Omega }_i(t),$$

$${\mathbf{\sigma }}_i(u,p)=-p\mathbb{I}+2{\mu }_i\mathbf{\epsilon }(u)\mathrm{in}{\Omega }_i(t),i=1,2,$$

$$\mathbf{\epsilon }(u)=\frac{1}{2}(\nabla u+(\nabla u{)}^T),$$

$$n\cdot [\mathbf{\sigma }{]}_{{\Gamma }_{int}(t)}=\gamma \kappa n\mathrm{on}{\Gamma }_{int}(t),\forall t\in [0,T],$$

$$[u{]}_{{\Gamma }_{int}(t)}=0\mathrm{on}{\Gamma }_{int}(t),\forall t\in [0,T].$$

$$\frac{\partial \varphi }{\partial t}+u\cdot \nabla \varphi =0\mathrm{in}{\Omega }_i(t),\forall t\in [0,T],$$

$$({\mathcal{S}}_u^h{)}_n=\{u^h|u^h\in [H^{1h}(Q_n){]}^{n_{sd}},u^h\dot{=}{g}^h\mathrm{on}(P_n{)}_g\},$$

$$({\mathcal{V}}_u^h{)}_n=\{w^h|w^h\in [H^{1h}(Q_n){]}^{n_{sd}},w^h\dot{=}0\mathrm{on}(P_n{)}_g\},$$

$$({\mathcal{S}}_p^h{)}_n=({\mathcal{V}}_p^h{)}_n=\{p^h|p^h\in H^{1h}(Q_n)\},$$

$$({\mathcal{S}}_{\varphi }^h{)}_n=\{{\varphi }^h|{\varphi }^h\in (H^{1h}{)}^{n_{sd}},{\varphi }^h\dot{=}{\hat{\varphi }}^h\mathrm{on}(P_n{)}_{\varphi }\},$$

$$({\mathcal{V}}_{\varphi }^h{)}_n=\{v^h|v^h\in (H^{1h}{)}^{n_{sd}},v^h\dot{=}0\mathrm{on}(P_n{)}_{\varphi }\}.$$

$${\int }_{Q_n}$$

$$+{\int }_{Q_n}{q}^h\nabla \cdot u^{h}dQ+{\int }_{{\Omega }_n}(w^h{)}_n^{+}\cdot {\rho }_i((u^h{)}_n^{+}-(u^h{)}_n^{-})d\Omega$$

$$+\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\tau }_{MOM}\frac{1}{{\rho }_i}\left[{\rho }_i\left(\frac{\partial w^h}{\partial t}+u^h\cdot \nabla w^h\right)-\nabla \cdot {\mathbf{\sigma }}_i(w^h,q^h)\right]$$

$$\cdot \left[{\rho }_i\left(\frac{\partial u^h}{\partial t}+u^h\cdot \nabla u^h-f\right)-\nabla \cdot {\mathbf{\sigma }}_i(u^h,p^h)\right]dQ$$

$$+\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\tau }_{CONT}\nabla \cdot w^h{\rho }_i\nabla \cdot u^{h}dQ={\int }_{(P_{int}{)}_n}{w}^h\cdot \gamma \kappa ndP.$$

$${\int }_{Q_n}$$

$$+\sum _{e=1}^{n_{el}}{\int }_{Q_n^e}(\frac{\partial v^h}{\partial t}+u^h\cdot \nabla v^h){\tau }_{LEV}(\frac{\partial {\varphi }^h}{\partial t}+u^h\cdot \nabla {\varphi }^h)dQ=0.$$

$${\int }_{(P_{int}{)}_n}{w}^h\cdot \gamma \kappa ndP$$

$$=-{\int }_{(P_{int}{)}_n}\gamma \underline{\nabla }{id}_{(P_{int}{)}_n}:\underline{\nabla }{w}^{h}dP,$$

$$Re=\frac{{\rho }_2\sqrt{gd}d}{{\mu }_2}=35,$$

$$Eo=\frac{g{\rho }_2{d}^2}{\gamma }=10.$$

$${\mathbb{X}}_c=(x_c,y_c)=\frac{{\int }_{{\Omega }_2}xdx}{{\int }_{{\Omega }_2}1dx},$$

$${\mathbb{X}}_c=(x_c,y_c,z_c)=\frac{{\int }_{{\Omega }_2}xdx}{{\int }_{{\Omega }_2}1dx}.$$

Full text

Keywords discontinuous-in-time Galerkin  $\cdot$ space-time finite elements  $\cdot$ space-time  $\cdot$ simplex  $\cdot$ evolving front  $\cdot$ two-phase flow  $\cdot$ level-set

1 Introduction

The spatial discretization of a time-dependent problem is often performed by means of the Galerkin or Petrov-Galerkin finite element (FE) method. However, the time discretization is typically being based on an explicit or implicit finite difference (FD) temporal discretization, such as $\theta$-family schemes or Runge-Kutta-family of methods. In recent years, the space-time finite element method has been steadily applied to more and more problems, e.g., advective-diffusive systems [13, 12], elastodynamics [11], Navier-Stokes equations [20, 10] and Navier-Stokes equations with deforming domains [22, 21, 9].

The space-time approach utilizes subsets of the temporal domain called space-time slabs, which are more or less similar to time steps of the standard semi-discrete approach. In most space-time implementations so far, the meshes for the space-time slabs are simply extruded in the temporal direction from a spatial mesh, resulting in semi-unstructured domains, which can be either unstructured or structured in space but structured in time. That means that space-time method has inherent flexibility to admit completely unstructured meshes with varying levels of refinement only in spatial dimensions, but does not allow different temporal refinement in different parts of the domain.

Nowadays, the extraordinary flexibility of the space-time FE is being exploited, when dealing with varying resolution of complicated domains not only in space dimensions but also in the time dimension. That leads to the use of fully unstructured meshes in both space and time. The construction of simplicial meshes suitable for space-time discontinuous Galerkin finite-element methods was introduced in [7] and relied on the “Tent Pitcher” algorithm of [23]. In [3], the generation of simplex space-time meshes was demonstrated and tested in the context of the advection-diffusion equation. Wang and Persson [24] showed a fully unstructured space-time mesh, which can cope with any type of domain deformations, even with topological changes. They also used local mesh operations in order to avoid remeshing. An algorithm for arbitrary finite element discretizations of the space-time cylinder was presented by Neumüller and Steinbach [17] as well. This method does not depend on the time-slabs, resulting in adaptive meshes, movable in time. In [17], the decomposition of a pentatope into smaller ones is also proposed. This decomposition relies on the Freudenthal algorithm [8]. In [16], a new procedure for the subdivision of four-dimensional prisms intersected by a moving front into simplices was illustrated. Such a subdivision method is important for a better resolution of the space-time interface.

In the present paper, we show a fully-unstructured space-time discretization of an interface-capturing finite element method, designed for two-phase incompressible flows including surface tension effects. We use $P1P1$ finite elements with least-squares stabilization. This approach is based on the discontinuous-Galerkin method in time (space-time elements), details of which can be found in [13, 12]. The variational formulation of the problem is written over the associated space-time domain. The interface is approximated by the level-set method. Level-set method describes implicitly the interface, meaning that the formulation is able to cope with extreme topological changes of the evolving front between the two phases. The benchmark cases reveal that the simulation results obtained with the fully-unstructured space-time discretization are equivalent to those obtained with the standard discretization, but offer a potential reduction in the number of degrees of freedom.

The structure of this paper is the following: In Section $2\text{\,}$ and $3\text{\,}$, the governing equations and their discretization are described, respectively. Section $4\text{\,}$ deals with the generation of simplex-type space-time meshes. In Section $5\text{\,}$, the numerical results are presented. The concluding remarks are drawn in Section $6\text{\,}$.

2 Governing equations

We focus on the incompressible two-phase flow which is governed by the transient, isothermal and incompressible Navier-Stokes equations. A computational domain $\Omega$ is considered, which is a subset of $\mathbb{R}^{n_{sd}}$; $n_{sd}$ is the number of space dimensions. This domain encloses two immiscible Newtonian phases $\Omega_{1}(t)$ and $\Omega_{2}(t)$, where $\Omega_{1}(t)\cup\Omega_{2}(t)=\Omega$. The boundary of the domain is denoted by $\Gamma=\partial{\Omega}$, whereas $\Gamma_{int}(t)$ stands for the interface between the two fluids $\partial{\Omega_{1}(t)}\cap\partial{\Omega_{2}(t)}$. Note that the spatial subdomains of the two phases and their interface are time-dependent.

At each instant $t\in[0,T]$, the velocity, $\mathbb{u}(\mathbb{x},t)$, and the pressure, $p(\mathbb{x},t)$, in each phase are governed by the following equations:

[TABLE]

[TABLE]

for $i=1,2$ number of phases. Here, $\rho_{i}$ is the density of the corresponding fluid. The stress tensor $\sigma_{i}$ and the rate of strain tensor are defined as

[TABLE]

[TABLE]

where $\mu_{i}$ is the dynamic viscosity of each fluid.

At the interface between the two phases, we impose a boundary condition which is based on the Laplace-Young equation:

[TABLE]

where $\mathbb{n}$ is the outward unit normal vector on $\Gamma_{int}(t)$, $\gamma$ the surface tension coefficient and $\kappa$ the local curvature of $\Gamma_{int}(t)$. Furthermore, the velocities are assumed to be continuous across the interface:

[TABLE]

In order to elucidate the evolution of the interface between the two phases, the level-set transport equation [18]:

[TABLE]

is solved, where $\phi$ is a signed-distance function and $\mathbb{u}$ is the velocity field obtained from the Navier-Stokes equations (1) – (2).

3 Solution technique

For the discretization of equations (1), (2) and (7), we use a first-order interpolation of all degrees of freedom in combination with a Galerkin/Least-Squares (GLS) stabilization in space and DG in time [6]. For each space-time slab $Q_{n}$, the following finite element interpolation and weighting function spaces are defined for the velocity, the pressure and the level-set function:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $P_{n}$ denotes the space-time extruded boundary, whereas $\mathbb{g}$ and $\hat{\phi}$ the prescribed values for the velocity and the level-set function, respectively, on the appropriate subset of $P_{n}$.

The stabilized space-time formulation for the Navier-Stokes equations (1) and (2) can then be written as follows: Given $(\mathbb{u}^{h})_{n}^{-}$, find $\mathbb{u}^{h}\in(\mathscr{S}_{\mathbb{u}}^{h})_{n}$ and $p^{h}\in(\mathscr{S}_{p}^{h})_{n}$ such that $\forall\mathbb{w}^{h}\in(\mathscr{V}_{\mathbb{u}}^{h})_{n},\forall q^{h}\in(\mathscr{V}_{p}^{h})_{n}$:

[TABLE]

The stabilized space-time formulation for the level-set transport equation (7) can then be written as follows: Given $(\phi^{h})_{n}^{-}$, find $\phi^{h}\in(\mathscr{S}_{\phi}^{h})_{n}$ such that $\forall v^{h}\in(\mathscr{V}_{\phi}^{h})_{n}$:

[TABLE]

The notation $(\bullet^{h})^{+}_{n}$ and $(\bullet^{h})^{-}_{n}$ in the weak formulations (13) and (14) denote the upper and lower values, respectively, of the discontinuous variable at the lower surface $\Omega$ of the space-time slab $Q_{n}$. The problem is solved sequentially for each space-time slab, starting with $(\bullet^{h})^{+}_{n}=\bullet_{0}$. Details on the parameters $\tau_{MOM}$ and $\tau_{CONT}$ in (13) can be found in [2], whereas details on the parameter $\tau_{LEV}$ in (14) are given by Sauerland and Fries [19].

The Laplace-Beltrami technique, as proposed in [14], is employed to reformulate the surface tension term in the weak formulation (13) and results in:

[TABLE]

where $\underline{\Delta}$ is the Laplace-Beltrami operator, $\underline{\nabla}$ is tangential gradient and $\mathbb{id}$ is the identity mapping on the space-time evolving interface $(P_{int})_{n}$.

4 Simplex space-time meshes

The key to space-time simulations with varying degrees of temporal refinement is the generation of simplex-based space-time meshes. A straightforward and robust algorithm that produces this type of meshes was already described by Behr [3]. Here, we use an updated version of this algorithm that allows arbitrary temporal refinement in selected portions of space-time slabs based on the level-set field.

To start with, we need a spatial mesh in $n_{sd}$ dimensions, which can be generated by any of the freely or commercially available mesh generators. We confine ourselves to $n_{sd}$-simplex-based meshes, such as triangular meshes in 2D and tetrahedral meshes in 3D. In the traditional space-time implementation, the spatial mesh is at first extruded in the time dimension to fill the space-time slab contained between time levels $t_{n}$ and $t_{n+1}$. The extruded mesh is composed of prisms (6-noded space-time elements for 2D problems and 8-noded space-time elements for 3D problems). These elements, referred to as 3d6n and 4d8n in [3] and shown in Figures 1(a) and 1(b), respectively, are the basis of the traditional space-time approach and are considered here as a reference.

Our goal is to subdivide these prism-type elements into simplex-type elements, referred to as 3d4n (familiar tetrahedrons) and 4d5n (pentatopes) in [3] and illustrated in 1(c) and 1(d), respectively. The initial space-time mesh contains only two nodes for each of the nodes in the spatial mesh, one located at the bottom of the slab and the other one located at the top of the slab. The temporal refinement is accomplished by adding, in parts of the domain where temporal accuracy is to be increased, one or more nodes along the line connecting the original nodes in the temporal direction. The space-time faces of the prism-type elements will be later divided into ($n_{sd}-1$)-node simplices according to the Delaunay criteria, independently for each prism. However, we need to ensure the uniqueness of the Delaunay process and the compatibility of the ($n_{sd}-1$)-simplices between the neighboring space-time prisms. A solution proposed in [3] for the aforementioned problem is to perturb the time coordinates of some or all nodes randomly. Additionally, a posteriori sliver elimination is unavoidable but easily executed by computing 3D or 4D element volumes and by rejecting elements whose volume falls below a specified threshold.

Considering two-phase flow problems, the region of the domain that demands higher temporal accuracy is the propagating interface between the two fluids. The level-set function gives us the opportunity to refine the mesh at a narrow band around the front, i.e., the zero level-set field, as presented in Figure 2. Inside this narrow band, the temporal refinement can be either uniform or non-uniform. That means that the number of nodes, added along the edges connecting the original nodes in the temporal direction of the extruded elements, can be different among these edges. Furthermore, there are elements outside the narrow band, but still next to it, for which, some of their edges contain more than two nodes in the temporal direction. These elements comprise a transition zone between the refined narrow band and the unrefined region of the domain. The sample connectivity of these elements for $2$D problems is illustrated in Figure 3.

5 Numerical results

5.1 Rising bubble in 2D

As a first test case, a two-dimensional bubble, rising in an initially motionless liquid column due to the buoyancy effects, is considered and used for validating the unstructured space-time mesh solver for problems in two-space dimensions. The computational domain is rectangular and has the size of 1.0\text{\,}\mathrm{m}$\times$2.0\text{\,}\mathrm{m}. An initially circular bubble with diameter $d=$0.5\text{,}\mathrm{m}$$ is placed inside the domain, as illustrated in Figure 4.

The fluids have the following properties: $\rho_{1}=$1000\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\rho_{2}=$100\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\mu_{1}=$10\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s} and $\mu_{2}=$1\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s}. The applied gravitational acceleration is $f_{y}=-g=$-0.98\text{,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$ and the surface tension coefficient $\gamma$ is equal to $24.5\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{s}^{2}$. The characteristic dimensionless Reynolds and Eötvös numbers are the following:

[TABLE]

[TABLE]

The spatial resolution consists of triangular elements. The aforementioned triangulation is based on an equidistant grid of 80\text{\,}$\times$160\text{\,} quadrilateral elements. The time-slab size is $\Delta t=$0.01\text{,}\mathrm{s}$$. No-slip boundary condition is imposed on the top and bottom of the computational domain and slip boundary condition on the vertical walls. Zero pressure is defined at the upper wall, and the initial velocity field is set to zero. With the given configuration, the bubble shape should become ellipsoidal, as stated in [15]. That means, the surface tension effects are dominant enough to hold the bubble together and no break up should be expected.

The bubble rose first for $t=$3.0\text{,}\mathrm{s}$$ using the usual discontinuous-in-time Galerkin time stepping (prismatic space-time elements, as shown in Figure 5). These standard results were then compared with the results obtained with the tetrahedral-based space-time mesh discretization of a slab of the same thickness without applying any temporal refinement (cf. Figure 5). Figure 6 shows the rising bubble over time and compares the bubble position obtained using the aforementioned types of discretization.

In Figure 7 the bubble shape at $t=$3.0\text{,}\mathrm{s}$$ is compared with results published in [15]. As we can see, the results obtained with both types of the different space-time elements show a very good agreement with the reference data. The centroid of the two-dimensional bubble, which is defined by:

[TABLE]

is also tracked.

Figure 8 depicts the position of the center of mass of the bubble. This quantity is also in good agreement with the TP2D simulation results reported in [15]. We have compared our results only with TP2D simulation results, published in [15], in order that the visibility of the Figures 7 and 8 could be maintained. All the different approaches in [15] lead to much the same results for this test case as well.

5.2 Rising droplet in 3D

The benchmark case of the rising bubble in 2D is now repeated in 3D. This example serves as the initial validation of the unstructured space-time mesh solver for problems in three-space dimensions. The bubble rose first for $t=$3.0\text{,}\mathrm{s} using the usual discontinuous-in-time Galerkin time stepping (prismatic space-time elements, linear-in-time interpolation). The time step size was at first chosen as $\Delta t=$0.01\text{\,}\mathrm{s}. These standard results were then compared with the rising results obtained with the pentatope-based space-time mesh discretization of a slab $\Delta t=$0.01\text{,}\mathrm{s}$$ thick without any temporal refinement.

The computational domain occupies a cuboid tank with the dimensions 1.0\text{\,}\mathrm{m}$\times$2.0\text{\,}\mathrm{m}$\times$1.0\text{\,}\mathrm{m}. The initial position of the droplet is at the center of the tank at a distance of $0.5\text{\,}\mathrm{m}$ from the bottom wall and its initial shape is assumed to be spherical (cf. Figure 9). The properties of the fluids are given by: $\rho_{1}=$1000\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\rho_{2}=$100\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\mu_{1}=$10\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s}, $\mu_{2}=$1\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s} and $f_{y}=-g=$-0.98\text{,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}. The surface tension coefficient is $\gamma=$24.5\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{s}^{2}. No-slip conditions are applied on all outer boundaries. Zero pressure is applied at the top wall and the initial velocity field is set to zero. Due to buoyancy effects, the droplet will start rising, change its shape and become ellipsoidal-shaped.

Figure 10 shows the rising droplet over time and compares the simulation results of the droplet position obtained with the prismatic and the simplex space-time mesh discretization. We also track the position of the centroid of the three-dimensional droplet, in order to acquire a more accurate comparison. The center of the mass is defined by:

[TABLE]

The position of the center of mass of the droplet is presented in Figure 11. This quantity is also in good agreement with the DROPS simulation results reported in [1]. We have compared our results only with DROPS simulation results, published in [1], again for maintaining the visibility of the Figure 11. NaSt3DGPF, reported in [1], leads to very similar results for this test case as well. However, the position of the droplet obtained with OpenFOAM is shifted in the vertical direction. One explanation for this difference might be that OpenFOAM employs the VOF method as an interface capturing method, as stated by Adelsberger et al. [1].

5.3 Step cavity in 2D with time refinement

As the last example for verifying our numerical approach, we simulate the filling stage of a step cavity, which has a lower and an upper step, considered by Dhatt et al. [5] and Cruchaga et al. [4]. The molten material enters the mold with uniform velocity and displaces the air, which is initially quiescent. This inflow velocity is deliberately chosen to be low. As a consequence, the spreading impact of gravity on the evolving interface is apparent.

The computational domain is illustrated in Figure 12. The spatial discretization of the domain consists of $365\text{\,}$ triangular elements. The time-slab size is $\Delta t=$0.005\text{,}\mathrm{s}$$. Slip boundary conditions are assumed on the horizontal and vertical walls, except for those of the inflow (leftmost vertical boundary) and outflow (rightmost vertical boundary). A uniform velocity is imposed at the inflow boundary, whereas traction-free boundary conditions are used at the outflow boundary. We consider isothermal condition, so natural convection and phase-change effects are disregarded.

The material properties correspond to those of [5] and are as follows: $\rho_{1}=$0.1\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\rho_{2}=$100\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}, $\mu_{1}=$0.02\text{,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s}, $\mu_{2}=$0.2\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{/}\mathrm{s}. The gravitational acceleration is equal to $f_{y}=-g=$-9.80\text{,}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$. The surface tension effects and the wall friction are neglected.

The molten material filled the mold first for $t=$1.6\text{,}\mathrm{s}$$ using the usual discontinuous-in-time Galerkin time stepping (prismatic space-time elements). These standard results were then compared with the filling results obtained with a hybrid tetrahedral-based space-time mesh discretization of the slab. This hybrid mesh was generated using the technique presented in [3], with each of the 320 time slabs being 0.005 thick and discretized differently with one to five elements in the time direction (cf. Figure 13). The part of the domain, where temporal accuracy is increased, covers the area close to the propagating interface, as shown in Figure 13.

Figure 14 illustrates the front position of the molten material at various time instances. The results obtained with the usual prismatic space-time discretization and with the hybrid mesh are also compared with those reported in Reference [4]. As we can see from Figure 14, the results show a good agreement with the reference data [4].

6 Efficiency Aspects

The example in Section 5.3 is used to provide reliable timing measurements, although it can be considered small and not too complex. Table 1 summarizes the typical performance behavior. The number of time steps, the number of nodes, the number of elements, as well as the total time required to form and to solve the equation systems are listed in Table 1, for

(a) a prism-based space-time finite element formulation with small time step ($\Delta t=$0.001\text{,}\mathrm{s}$$),

(b) a prism-based space-time finite element formulation with large time step ($\Delta t=$0.005\text{,}\mathrm{s}$$), as used in Figure 14 (top row), and

(c) a simplex-based space-time finite element formulation with variable time step ($\Delta t=$0.001\text{,}$-$0.005\text{,}\mathrm{s}$$), as used in Figure 14 (bottom row).

The results of the prism-based space-time finite element formulation with small time step ($\Delta t=$0.001\text{,}\mathrm{s}$$) are not presented in Subsection 5.3, but look similar to the reference data, shown in Figure 14 (middle row).

According to the above table, the total time required to obtain a solution with time step $\Delta t=$0.001\text{,}\mathrm{s} is four times that required to obtain a solution with time step $\Delta t=$0.005\text{\,}\mathrm{s}, assuming fixed linear solver parameters and a prism-based space-time finite element formulation. The average number of nodes in the hybrid mesh with variable temporal refinement is ca. $919$, whereas the maximum and the minimum number of nodes is $730$ and $998$, respectively. That means that the solution time is also higher in the case of a simplex-based space-time finite element formulation with variable temporal refinement than the one required for a prismatic mesh with the same maximum time step size when using iterative solvers with linear scaling properties. It should be mentioned that at every time step, Newton-Raphson iterations are performed to solve the nonlinear discretized system resulting from the nonlinear Navier-Stokes equations and strong coupling iterations are also executed because of the mutual dependence of the level set field on the fluid velocity field and vice versa. A GMRES solver is used to solve the resulting linear system of equations. However, the conditioning of the linear systems arising from both types of elements was not examined in detail. The convergence of the iterative GMRES solver was similar in all the aforementioned test cases, though.

Furthermore, as the average number of elements in a simplex-based mesh is significantly higher than that in the prism-based meshes (ca. $3003$ on average versus $730$), even though the elements are simpler, the system formation time takes almost twice as long. This time remains lower, however, than the time required to repeatedly form the system for a prism-based space-time finite element formulation with the small time step ($\Delta t=$0.001\text{,}\mathrm{s}$$).

To sum up, the use of simplex space-time meshes without any temporal refinement has some disadvantages. The equation system size remains the same in comparison with a prismatic mesh of the same time step size, but the number of elements is significantly increased leading to increased system formation times. However, the efficiency of the discretization is much improved when using local temporal refinement, while providing a better resolution of the space-time evolving interface.

7 Concluding remarks

We have presented an updated version of the straightforward method for generating simplex space-time meshes, which was already introduced by Behr [3]. This version is based on the level-set method and allows arbitrary temporal refinement of the space-time slabs in the vicinity of evolving fronts. We have tested the resulting unstructured space-time meshes in the context of two-phase flow problems. The benchmark cases of the static bubble and the rising bubble/droplet have served as the initial validation of the unstructured space-time mesh solver for the Navier-Stokes equations and the level-set equation in two and three space dimensions. The benchmark case of a two-dimensional step cavity filling was used for checking the reliability of the results obtained with a hybrid mesh. Future work includes the extension of the arbitrary temporal refinement to more complicated problems, such as the filling of three-dimensional complex molds. The numerical behavior of stabilized FE formulations on such highly unstructured space-time meshes still needs to be examined. Likewise, the efficiency aspects are yet to be analyzed for more complex geometries.

Acknowledgements

The authors gratefully acknowledge the support of the German Research Foundation (DFG) under program SFB 1120. The computations were conducted on computing clusters provided by the RWTH Aachen University IT Center and by the Jülich Aachen Research Alliance (JARA).