Simplex Space-Time Meshes in Compressible Flow Simulations

TL;DR

This paper introduces the use of simplex space-time meshes in compressible flow simulations, enabling flexible, efficient, and topology-changing simulations with local temporal refinement, demonstrated through complex 2D and 3D cases.

Contribution

It presents a finite element formulation for compressible flows on simplex space-time meshes, allowing topology changes and local temporal refinement in simulations.

Findings

Successful simulation of a valve with topology change

Parallel 3D simulation of piston ring blow-by

Computational time savings from local temporal refinement

Abstract

Employing simplex space-time meshes enlarges the scope of compressible flow simulations. The simultaneous discretization of space and time with simplex elements extends the flexibility of unstructured meshes from space to time. In this work, we adopt a finite element formulation for compressible flows to simplex space-time meshes. The method obtained allows, e.g., flow simulations on spatial domains that change topology with time. We demonstrate this with the two-dimensional simulation of compressible flow in a valve that fully closes and opens again. Furthermore, simplex space-time meshes facilitate local temporal refinement. A three-dimensional transient simulation of blow-by past piston rings is run in parallel on 120 cores. The timings point out savings of computation time gained from local temporal refinement in space-time meshes.

Figures

Tables

Table 1. Table 1: Mesh data and timings for piston ring blow-by simulations.

	$n_{\mathrm{ts}}$	$n_{\mathrm{el}}$ [Mio]	$n_{\mathrm{dof}}$[k]	$t_A$ [s]	$t_S$ [s]	$t_T$ [s]
Coarse FST	100	0.45	683	362	540	946
SST	100	2.91	894	812	1003	1898
Fine FST	400	0.45	683	1456	2175	3761

Equations

$$\tilde{\mathbf{U}}={\left[\rho ,\rho u_1,\rho u_2,\rho u_3,\rho e_{tot}\right]}^{\top },$$

$$\mathbf{Y}={\left[p,u_1,u_2,u_3,T\right]}^{\top },$$

$$p=\rho RT,\text{and}e=\frac{R}{\gamma -1}T.$$

$${\tilde{\mathbf{U}}}_{,t}+{\tilde{\mathbf{F}}}_{i,i}^{\text{adv}}-{\tilde{\mathbf{F}}}_{i,i}^{\text{diff}}-\tilde{\mathbf{S}}=0.$$

$$\mathbf{\tilde{F}}^{\text{adv}}_{i}=\left[\begin{array}[]{c}\rho u_{i}\\ \rho u_{i}u_{1}+p\delta_{1i}\\ \rho u_{i}u_{2}+p\delta_{2i}\\ \rho u_{i}u_{3}+p\delta_{3i}\\ \rho u_{i}e_{tot}+pu_{i}\end{array}\right],\qquad\mathbf{\tilde{F}}^{\text{diff}}_{i}=\left[\begin{array}[]{c}0\\ \tau_{1i}\\ \tau_{2i}\\ \tau_{3i}\\ \tau_{ij}u_{j}-q_{i}\end{array}\right].$$

$${\tau }_{ij}=\mu \left(u_{i,j}+u_{j,i}\right)+\lambda u_{k,k}{\delta }_{ij},i,j=1,\dots ,n_{\mathrm{sd}},$$

$$q_i=-\kappa T_{,i},$$

$${\mathbf{U}}_{,t}+{\mathbf{F}}_{i,i}^{\text{adv p}}+{\mathbf{F}}_{i,i}^{\text{p}}+{\mathbf{F}}^{\text{sp}}-{\mathbf{F}}_{i,i}^{\text{diff}}-\mathbf{S}=0.$$

$$\mathbf{U}={\left[\rho ,\rho u_1,\rho u_2,\rho u_3,\rho e\right]}^{\top },$$

$$\mathbf{F}^{\text{adv \textbackslash p}}_{i}=\left[\begin{array}[]{c}\rho u_{i}\\ \rho u_{i}u_{1}\\ \rho u_{i}u_{2}\\ \rho u_{i}u_{3}\\ \rho u_{i}e\end{array}\right],\qquad\mathbf{F}^{\text{p}}_{i}=\left[\begin{array}[]{c}0\\ p\delta_{1i}\\ p\delta_{2i}\\ p\delta_{3i}\\ 0\end{array}\right],\qquad\mathbf{F}^{\text{sp}}=\left[\begin{array}[]{c}0\\ 0\\ 0\\ 0\\ pu_{i,i}-\tau_{i,j}u_{j,i}\end{array}\right],\qquad\mathbf{F}^{\text{diff}}_{i}=\left[\begin{array}[]{c}0\\ \tau_{1i}\\ \tau_{2i}\\ \tau_{3i}\\ -q_{i}\end{array}\right].$$

$${\mathbf{A}}_0:=\frac{\partial \mathbf{U}}{\partial \mathbf{Y}},{\mathbf{A}}_i^{\text{adv p}}:=\frac{\partial {\mathbf{F}}_i^{\text{adv p}}}{\partial \mathbf{Y}},{\mathbf{A}}_i^{\text{p}}:=\frac{\partial {\mathbf{F}}_i^{\text{p}}}{\partial \mathbf{Y}}.$$

$${\mathbf{A}}_i^{\text{sp}}{\mathbf{Y}}_{,i}={\mathbf{F}}^{\text{sp}},{\mathbf{K}}_{ij}{\mathbf{Y}}_{,j}={\mathbf{F}}_i^{\text{diff}}.$$

$$\text{Res}(\mathbf{Y}):={\mathbf{A}}_0{\mathbf{Y}}_{,t}+\left({\mathbf{A}}_i^{\text{adv p}}+{\mathbf{A}}_i^{\text{p}}+{\mathbf{A}}_i^{\text{sp}}\right){\mathbf{Y}}_{,i}-({\mathbf{K}}_{ij}{\mathbf{Y}}_{,j}{)}_{,i}-\mathbf{S}=0.$$

$$\mathbf{x}=\left[\begin{array}[]{c}x_{1}\\ \vdots\\ x_{d}\\ \end{array}\right].$$

$$\mathbf{x}=\Phi (\mathbf{\xi })=\sum _{q=1}^{d+1}{\mathbf{x}}_q{\chi }_q(\mathbf{\xi }),$$

$$x(\xi ,\eta )$$

$$y(\xi ,\eta )$$

$${\chi }_1(\xi ,\eta )$$

$${\chi }_2(\xi ,\eta )$$

$${\chi }_3(\xi ,\eta )$$

$${\chi }_1(\mathbf{\xi })$$

$${\chi }_j(\mathbf{\xi })$$

$$\mathbf{J}=\left(\begin{array}[]{cc}\frac{\partial x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\ \frac{\partial y}{\partial\xi}&\frac{\partial y}{\partial\eta}\\ \end{array}\right).$$

$$\mathbf{J}_{k}=\frac{\partial(x,y)}{\partial(\xi,\eta)}=\left(\begin{array}[]{cc}x_{2}-x_{1}&x_{3}-x_{1}\\ y_{2}-y_{1}&y_{3}-y_{1}\\ \end{array}\right).$$

$$\mathbf{J}_{k}=\frac{\partial\mathbf{x}}{\partial{\boldsymbol{\xi}}}=\left(\begin{array}[]{cccc}\mathbf{x}_{2}-\mathbf{x}_{1},&\mathbf{x}_{3}-\mathbf{x}_{1},&\dots,&\mathbf{x}_{d+1}-\mathbf{x}_{1}\\ \end{array}\right).$$

$${\hat{G}}_{ij}=\sum _{k=1}^d\frac{\partial {\xi }_k}{\partial x_i}\frac{\partial {\xi }_k}{\partial x_j},i,j=1,\dots ,d.$$

$$\text{vol(△∗)}=\frac{1}{d!}\stackrel{!}{=}\frac{l^d}{d!}\sqrt{\frac{d+1}{2^d}}=\text{vol(△r)}$$

$$l=\sqrt{2(d+1{)}^{-\frac{1}{d}}}.$$

$$\mathbf{L}=\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_r}=\frac{\partial \mathbf{x}}{\partial \mathbf{\xi }}\frac{\partial \mathbf{\xi }}{\partial {\mathbf{x}}_r}.$$

$${\mathbf{G}}^{-1}=\mathbf{L}{\mathbf{L}}^{\top }=\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_r}{\left(\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_r}\right)}^{\top }=\frac{\partial \mathbf{x}}{\partial \mathbf{\xi }}\frac{\partial \mathbf{\xi }}{\partial {\mathbf{x}}_r}{\left(\frac{\partial \mathbf{\xi }}{\partial {\mathbf{x}}_r}\right)}^{\top }{\left(\frac{\partial \mathbf{x}}{\partial \mathbf{\xi }}\right)}^{\top },$$

$$\mathbf{G}={\left({\mathbf{G}}^{-1}\right)}^{-1}={\left(\frac{\partial \mathbf{x}}{\partial \mathbf{\xi }}\right)}^{-\top }{\left(\frac{\partial \mathbf{\xi }}{\partial {\mathbf{x}}_r}\right)}^{-\top }{\left(\frac{\partial \mathbf{\xi }}{\partial {\mathbf{x}}_r}\right)}^{-1}{\left(\frac{\partial \mathbf{x}}{\partial \mathbf{\xi }}\right)}^{-1}={\left(\frac{\partial \mathbf{\xi }}{\partial \mathbf{x}}\right)}^{\top }{\left(\frac{\partial {\mathbf{x}}_r}{\partial \mathbf{\xi }}\right)}^{\top }\frac{\partial {\mathbf{x}}_r}{\partial \mathbf{\xi }}\frac{\partial \mathbf{\xi }}{\partial \mathbf{x}}.$$

$$\mathbf{G}={\left(\frac{\partial \mathbf{\xi }}{\partial \mathbf{x}}\right)}^{\top }\mathbf{M}\frac{\partial \mathbf{\xi }}{\partial \mathbf{x}},\text{with}\mathbf{M}={\left(\frac{\partial {\mathbf{x}}_r}{\partial \mathbf{\xi }}\right)}^{\top }\frac{\partial {\mathbf{x}}_r}{\partial \mathbf{\xi }}={\mathbf{J}}_r^{\top }{\mathbf{J}}_r.$$

Full text

1 Introduction

In this work, we combine two constituents. A finite element formulation for compressible Navier–Stokes equations and simplex space-time meshes.

Regarding the first constituent, our formulation can be traced back to the extension of the streamline upwind/Petrov-Galerkin (SUPG) method [1] to compressible flow problems.[2] The formulation based on conservation variables as primary unknowns posed a breakthrough in handling the advective character of compressible flows on equal-order interpolation elements. Building on this, Shakib in his thesis [3] developed and analyzed space-time finite element algorithms for compressible flow simulations using entropy variables. Entropy variables as primary unknowns lead to a symmetric from of the compressible flow equations. Guidance in the question of which variable set to use as primary unknowns can be found in the comparative study of Hauke.[4] The performance of four varialbe sets, namely conservation, entropy, density-, and pressure-primitive variables, has been evaluated for various test cases. Pressure-primitive variables were found to be well behaved in the incompressible limit and the most convenient to prescribe boundary conditions. Recently, Xu et al. further developed the formulation based on pressure-primitive variables and extended it for moving domains, weakly enforced essential boundary conditions, and sliding interfaces, with a generalized $\alpha$-method for time integration.[5] In the following sections, we will adopt essential parts of this formulation for simplex space-time finite elements.

Simplex space-time meshes as the second constituent of our work provide a simultaneous discretrization of space and time. The meshes are characterized by an unstructured discretization of the spatial and temporal extent of the computational domain. In 1988, the potential of unstructured space-time meshes has been first postulated by Hughes and Hulbert.[6] They proposed solving a one-dimensional elastodynamics problem efficiently on an unstructured space-time mesh. Later, Maubach numerically solved the two-dimensional advection-diffusion equations on unstructured space-time meshes and provided a mathematical analysis of the method.[7] Unstructured space-time meshes were also applied in other fields, such as viscoelastic problems[8] for example. However, the spatial domains of the transient problems remained two-dimensional for roughly 20 years.

The key step in using unstructured space-time meshes for transient three-dimensional problems is the mesh generation. A robust and simple procedure to construct four-dimensional simplex-based space-time meshes has been presented by Behr.[9] The full space-time cylinder is first subdivided into space-time slabs which are in turn discretized with prism-type elements, followed by a Delaunay triangulation of each prism-type element. Neumüller and Steinbach proposed an alternative approach to four-dimensional simplex mesh generation, which includes the refinement of space-time meshes based on the Freudenthal algorithm.[10] A third four-dimensional meshing strategy based on combinatorics was described in the thesis of Wang.[11]

With mesh generation procedures available, recently three-dimensional transient problems were solved harnessing the advantages of unstructured space-time meshes. We describe here three examples in the field of fluid dynamics. A high-order discontinuous Galerkin method for compressible flows on domains with large deformations [12] was extended to three dimensions.[11] Also in conjunction with a discontinuous Galerkin method, Lehrenfeld employed unstructured space-time discretizations to solve two-phase mass transport problems in three dimensions.[13] As a third example, the incompressible two-phase flow simulations of Karyofylli et al. were boosted by adaptive temporal refinement along the moving fluid interface using simplex finite elements.[14] Also linear algebra aspects of unstructured space-time discretizations were recently explored. Steinbach and Yang studied the performance of algebraic multigrid preconditioned GMRES methods in solving linear systems arising from adaptive three- and four-dimensional space-time discretizations of the heat equation.[15]

Now, with both constituents introduced, we proceed with the combination as follows. In Section 2, we briefly review the Navier–Stokes equations of compressible flows and recast them as generalized advective-diffusive system. In Section 3, we discuss the numerical treatment and present a formulation suitable for simplex space-time finite elements. Section 4 aims at experimentally validating the formulation, as well as, demonstrating the particular potential of compressible flow simulations on simplex space-time meshes. Section 5 contains concluding remarks.

2 Governing equations of compressible flows

2.1 Preliminaries

As the Navier–Stokes equations of compressible flows state the pointwise conservation of mass, momentum, and energy, they are naturally expressed in terms of the conserved quantities, called conservation variables $\mathbf{\tilde{U}}$. For a three-dimensional problem the conservation varialbes are

[TABLE]

where $\rho$ is the density, $u_{i}$ are the velocity components in each space direction, $i=1,...,{n_{\mathrm{sd}}}$, collected in $\mathbf{u}$, and $e_{tot}=e+\frac{1}{2}\|\mathbf{u}\|^{2}$ is the total energy per unit mass of the fluid, being the sum of the internal energy $e$ and the kinetic energy $\frac{1}{2}\|\mathbf{u}\|^{2}$.

An alternative set of five independent variables are the pressure-primitive variables

[TABLE]

where $p$ is the pressure and $T$ is the temperature. Throughout this paper we will consider an ideal, calorically perfect gas with

[TABLE]

$R$ is the specific gas constant, and $\gamma$ the ratio of specific heats.

2.2 Strong form

Expressed in conservation variables, the strong form of the compressible Navier–Stokes equations is given as

[TABLE]

Therein, $\mathbf{\tilde{S}}$ is the vector of source terms. The vector of advective fluxes $\mathbf{\tilde{F}}^{\text{adv}}_{i}$ and the vector of diffusive fluxes $\mathbf{\tilde{F}}^{\text{diff}}_{i}$ are defined as

[TABLE]

Here, $\delta_{ij}$ is the Kronecker delta. The viscous stress tensor $\boldsymbol{\tau}$ models the stress contribution caused by the strain rate; its entries $\tau_{ij}$ are defined as

[TABLE]

where $\mu$ is the dynamic viscosity and $\lambda=-\frac{2}{3}\mu$ the bulk viscosity. The entries $q_{i}$ of the heat flux vector are

[TABLE]

where $\kappa$ is the coefficient of thermal conductivity. In the following, we denote partial time derivatives with $(\cdot)_{,t}$, partial derivatives in spatial directions with $(\cdot)_{,i}$. The Einstein summation convention applies to repeated indices.

2.3 Reduced form of the energy equation

By algebraic manipulation of the equation system (4), the energy equation of the compressible Navier–Stokes equations can be cast in a reduced form,[5] leading to the following system

[TABLE]

The conservation variables with reduced energy equation,

[TABLE]

differ from $\mathbf{\tilde{U}}$ in the fifth entry. $\rho e$ is solely the internal energy per unit mass instead of $\rho e_{tot}$ composed of internal and kinetic energy per unit mass.

To facilitate a separate treatment in the weak form, the advective fluxes are split into the contribution to the advective fluxes without the pressure contribution $\mathbf{F}^{\text{adv \textbackslash p}}_{i}$ and the pressure contribution $\mathbf{F}^{\text{p}}_{i}$. Further, the balance laws with reduce energy equation contain the contribution of stress power to the energy equation $\mathbf{F}^{\text{sp}}$ and the diffusive fluxes $\mathbf{F}^{\text{diff}}_{i}$. The four flux vectors read

[TABLE]

$\mathbf{S}$ can be used to model source terms such as body forces in the momentum equations.

2.4 Generalized advective-diffusive system

A formulation with $\rho$ as variable is not well-defined in the incompressible limit.[4] Also, the specification of boundary conditions is cumbersome when using the conservation variables $\mathbf{U}$ as primary unknowns. Hence, we consider the conservation variables $\mathbf{U}$ as function of the pressure-primitive variables $\mathbf{Y}$, leading to the change of variables $\mathbf{U}=\mathbf{U}(\mathbf{Y})$.[4] For convenience, the following mappings are introduced

[TABLE]

Further, the matrices $\mathbf{A}^{\text{sp}}_{i}$ and $\mathbf{K}_{ij}$ are constructed, such that

[TABLE]

Explicit forms of the matrices are given in the appendix of the publication by Xu et al. [5]

Using the change of variables $\mathbf{U}=\mathbf{U}(\mathbf{Y})$, and the mappings above, Equation (8) can be written as generalized advective-diffusive system

[TABLE]

3 Numerical treatment with simplex space-time finite elements

3.1 Simplex finite elements

A point $\mathbf{x}\in\mathbb{R}^{d}$ is characterized by its coordinates $x_{i},i=1,...,d$. As commonly done, we arrange the coordinates in a column vector

[TABLE]

In the case $d=2$, we denote $x_{2}$ by $y$. To discretize any finite-dimensional polytope region $R\subset\mathbb{R}^{d}$, without gaps, $d$-simplices can be used. Simplices for $d=1,2,3$, and $4$ are shown in Figure 1. The 1-simplex is a line, the 2-simplex a triangle, the 3-simplex a tetrahedron, and the 4-simplex a pentatope. Each $d$-simplex has $d+1$ vertices. We refer to the points where the simplex vertices of our $d$-dimensional tesselation meet as nodes.

We further define the isoparametric local-global mapping $\Phi$ from a reference element $\bigtriangleup_{*}$ with coordinates ${\boldsymbol{\xi}}$ to a physical element $\bigtriangleup_{k}$ with coordinates $\mathbf{x}$ as

[TABLE]

where the index $q$ runs over the nodes of the physical element $\bigtriangleup_{k}$.[16] See Figure 2 for a two-dimensional visualization. In two dimensions, the mapping reads

[TABLE]

with the $\mathbb{P}_{1}$ basis functions defined on the reference element $\bigtriangleup_{*}$

[TABLE]

Generalized to $d$-dimensions, the $\mathbb{P}_{1}$ basis functions are

[TABLE]

We follow the usual convention and define the Jacobian matrix of a vector-valued function by varying the differential operator, e.g., $\frac{\partial}{\partial\xi}$, between columns. In two dimensions, the Jacobian matrix reads

[TABLE]

Inserting the $\mathbb{P}_{1}$ basis function definitions (17) into (16) and that into (19), $\mathbf{J}_{k}$ can be computed for a given two-dimensional physical element $\bigtriangleup_{k}$ as

[TABLE]

This generalizes to $d$-dimensions as

[TABLE]

The explicit form of the Jacobian confirms that $\mathbf{J}_{k}$ is constant per element, which shows a characteristic of the isoparametric $\mathbb{P}_{1}$ mapping $\Phi$ being an affine transformation. We further observe, that the columns of $\mathbf{J}_{k}$ are the distance vectors $\mathbf{p}_{i}=\mathbf{x}_{i+1}-\mathbf{x}_{1}$ from node 1 to the other $d$ nodes of the simplex. Clearly, permuting the node numbering changes $\mathbf{J}_{k}$.

3.2 Metric tensor definition

In stabilized finite element flow simulations, the covariant metric tensor $\hat{\mathbf{G}}$ is frequently used to incorporate directional element length information into the stabilization parameter.[3, 4, 5, 17] (We will discuss the details of stabilization in Section 3.5.) The definition of the covariant metric tensor $\hat{\mathbf{G}}$ is based on the inverse Jacobian $\mathbf{J}_{k}^{-1}=\frac{\partial{\boldsymbol{\xi}}}{\partial\mathbf{x}}$

[TABLE]

However, with this definition, $\hat{\mathbf{G}}$ varies under permutations of the physical element’s node numbering for simplex elements. To achieve invariance, we incorporate a regular simplex element $\bigtriangleup_{r}$ in the definition of the metric tensor.[17] The $d$-dimensional regular simplex element $\bigtriangleup_{r}$ is constructed, such that all its edges are of equal length $l$. Requiring further that $\bigtriangleup_{r}$ has the volume of the reference element

[TABLE]

leads to the dimension-dependent edge length of

[TABLE]

The formulas for the volumes in Equation (23) can be proven by induction.[18]

Utilizing the $\mathbb{P}_{1}$ basis functions defined on the reference element $\bigtriangleup_{*}$, we define a mapping $\Psi$ from the reference element $\bigtriangleup_{*}$ to the regular simplex element $\bigtriangleup_{r}$, i.e., $\mathbf{x}_{r}=\Psi({\boldsymbol{\xi}})$. The Jacobian $\mathbf{J}_{r}=\frac{\partial\mathbf{x}_{r}}{\partial{\boldsymbol{\xi}}}$ associated with this mapping can be computed in the same way as $\mathbf{J}_{k}$ using (21). From vol($\bigtriangleup_{r}$) $\overset{!}{=}$ vol($\bigtriangleup_{*}$), follows det$(\mathbf{J}_{r})=1$. The composition of $\Phi$ and $\Psi^{-1}$ we call $\Pi=\Phi\circ\Psi^{-1}$, which maps from the regular simplex element $\bigtriangleup_{r}$ to the physical element $\bigtriangleup_{k}$. A visualization can be found in Figure 2. The Jacobian $\mathbf{L}$ associated with the mapping $\Pi$ is

[TABLE]

If we base the definition of the contravariant metric tensor $\mathbf{G}^{-1}$ on $\mathbf{L}$, then

[TABLE]

is invariant under cyclic permutations of the physical element’s node numbering, and so is the covariant metric tensor $\mathbf{G}$ defined as its inverse

[TABLE]

The mapping $\frac{\partial\mathbf{x}_{r}}{\partial{\boldsymbol{\xi}}}$ is known for each $d$-simplex and can be precomputed, such that the computation of $\mathbf{G}$ simplifies to

[TABLE]

$\mathbf{M}=\mathbf{J}_{r}^{\top}\mathbf{J}_{r}$ prompts the observation that the entries of $\mathbf{M}$ are in fact scalar products of the distance vectors between the nodes of the regular simplex (compare Equation (21))

[TABLE]

By construction, all edges of the regular simplex are of length $l$ (given in Equation (24)) and all angles between two edges are 60 degree. Hence,

[TABLE]

hold for all choices of the base node for the distance vectors, i.e., $\mathbf{x}_{1}$ in Equation (21). It follows that $\mathbf{M}$, characterized by

[TABLE]

exclusively measures edge lengths and angles of the regular simplex and therefore is invariant under permutations of the node numbering. For $d=2,3,4$, $\mathbf{M}$ reads, respectively

[TABLE]

The explanation why $\mathbf{G}$ is invariant under node permutations is now straightforward, we first observe that a permutation of the node numbering of the physical element is equivalent to a permutation of the node numbering of the reference element or the regular simplex element. Applying the node permutation to regular simplex element one obtains

[TABLE]

which is associated with the mapping from the regular simplex element with permuted node numbering to the physical element. The corresponding covariant metric tensor with the permutation applied reads

[TABLE]

The last equality follows from the characterization of $\mathbf{M}$ in Equation (31) and recovers the original definition of $\mathbf{G}$ in Equation (28), showing that $\mathbf{G}$ is indeed invariant under permutations of the physical element’s node numbering. We will use the invariant metric tensor $\mathbf{G}$ to construct a stabilization matrix $\boldsymbol{\tau}$ that is likewise invariant under permutations of the element’s node numbering (see Section 3.5).

3.3 Space-time discretization methods

With the mappings and metric tensors for simplex elements defined, we can now proceed to discretize the subset of the space-time continuum $Q$ in which we wish to solve the compressible Navier–Stokes equations. Thinking of $Q$ being generated by the time dependent spatial domain $\Omega_{t}\subset\mathbb{R}^{{n_{\mathrm{sd}}}}$ evolving over the time interval $I\subset\mathbb{R}$ shows that $Q\subset\mathbb{R}^{{n_{\mathrm{sd}}}+1}$. Based on the discontinuous Galerkin method in time, $Q$ is sliced into $N$ space-time slabs $Q_{n}$. The extent of the space-time slabs in temporal direction is denoted by $\Delta t$. An exemplary slicing of $Q=[x_{0},x_{3}]\times]t_{0},t_{3}]$ is shown in Figure 3a. Each space-time slab is bounded by the spatial domain at the lower and upper time level $\Omega_{\text{l}}$, $\Omega_{\text{u}}$, respectively, as well as by $P\subset\mathbb{R}^{{n_{\mathrm{sd}}}}$, which is the temporal evolution of the spatial domain boundary $\Gamma\subset\mathbb{R}^{{n_{\mathrm{sd}}}-1}$.

In this paper, three methods to discretize $Q_{n}$ are considered. The flat space-time (FST) method extrudes a spatial discretization of $\Omega$ in time generating a discretization of $Q_{n}$ with prismatic space-time elements $Q^{e}_{n}$ (Figure 3a). In the simplex space-time (SST) method, these prismatic elements are further subdivided into simplex elements $Q^{e}_{n}$ (Figure 3b), allowing for local temporal refinement.[9] A third option is the unstructured space-time (UST) method, which generates an unstructured tesselation of $Q_{n}$ with simplex elements $Q^{e}_{n}$, e.g., by Delaunay triangulation (Figure 3c).

In Figure 3, the spatial domain is restricted to one dimension for the sake of a clear visualization. However, the FST and SST discretizations are well established for complex three-dimensional spatial geometries, i.e., $Q\subset\mathbb{R}^{4}$.[14] In case of UST, simplex discretizations of arbitrary $Q\subset\mathbb{R}^{3}$ can be obtained with a variety of commercially or freely available meshing tools. UST discretizations of $Q\subset\mathbb{R}^{4}$ of domains of engineering interest are to the authors’ knowledge an open research problem.

3.4 Weak form

Based on one of the above discretizations of the space-time slabs $Q_{n}$, trial and test functions are generated in each element $Q_{n}^{e}$ by the linear interpolation of the basis functions in the reference element and the mappings $(\Phi_{n}^{e})^{-1}$ from physical elements to the reference element. Further, these functions are constructed to be $C^{0}$-continuous on the closure of the space-time slab $\overline{Q_{n}}$. Hence, we define the $H^{1}$-conformal finite element approximation space as

[TABLE]

In case of an FST discretization, $\mathbb{V}_{1}$ refers to the $\mathbb{PR}_{1}$ basis functions of the prismatic Lagrange finite element. In case of an SST or UST discretization, $\mathbb{V}_{1}$ refers to the $\mathbb{P}_{1}$ basis functions (Equation (18)) of the simplical Lagrange finite element.

Further, we define a projection $\mathbf{P}$ that selects from $\mathbf{Y}^{h}\in[H^{1}_{h,n}]^{{n_{\mathrm{sd}}}+2}$ the degrees of freedom, where Dirichlet boundary data ${\boldsymbol{\mathit{g}}}^{h}$ is prescribed. The parts of the boundary $P_{n}$, where Dirichlet boundary conditions are prescribed, may vary for each of the primitive variables $\left(p,u_{1},u_{2},u_{3},T\right)$. $\mathbf{P}$ allows the definition of the trial function space

[TABLE]

that fulfills the Dirichlet boundary conditions and the test function space

[TABLE]

with functions that vanish where Dirichlet boundary conditions are prescribed.

Considering that the finite element functions are discontinuous at the space-time slab boundaries $\Omega_{l}$ and $\Omega_{u}$, let $(\mathbf{Y}^{h})^{\pm}_{n}$ abbreviate $\lim_{\varepsilon\to 0}\mathbf{Y}^{h}(t_{n}\pm\varepsilon)$. Hence, the weak form of Equation (13) can be stated as: For given initial conditions $(\mathbf{Y}^{h})^{-}_{0}=(\mathbf{Y}^{h})_{0}$, find $\mathbf{Y}^{h}\in S_{h,n}$, such that on each time slab $Q_{n}$ and for all $\mathbf{W}^{h}\in V_{h,n}$

[TABLE]

Therein, the bilinear form reads

[TABLE]

and the linear functional reads

[TABLE]

Here we used

[TABLE]

where $n_{i}$ are the components of the outwards pointing surface normal $\mathbf{n}$, not to be confused with plain $n$ that indexes the space-time slabs form 0 to $N-1$. For test cases with pentatope discretization, the domain boundary consists of tetrahedra in four dimensions. The normals to the boundary tetrahedra are computed using the generalized cross product definition.[13] The diffusive and pressure fluxes in the second integral of Equation (39) are integrated by parts and give rise to the boundary integral over $P_{n}$, while ${\boldsymbol{\mathit{h}}}^{h}$ is evaluated with the solution of the previous Newton step. On the boundary between subsequent time slabs, i.e., $\Omega_{l}$, the third integral in Equation (39) weakly enforces the coupling between the space-time slabs.

3.5 SUPG operator

It is well-known that the Galerkin formulation for the convection-dominated compressible Navier–Stokes equations (Equation (38)) suffers from instabilities. To counteract these instabilities, the weak form is augmented by the following SUPG operator

[TABLE]

[TABLE]

where $m=1,\dots,n_{sd}+1$ runs over $n_{sd}$ spatial dimensions and time. $\hat{\mathbf{A}}_{i}$ are the combined advection matrices for conservation variables

[TABLE]

The advection matrix for conservation variables in time direction is the identity

[TABLE]

For the element-wise stabilization matrix $\boldsymbol{\tau}^{e}$, several definitions are proposed in literature.[3, 5, 19] We select the direct approach [5] and adopt it for unstructured space-time discretizations

[TABLE]

The diffusion matrices for the conservation variables are obtained as $\hat{\mathbf{K}}_{ij}=\mathbf{K}_{ij}\mathbf{A}_{0}^{-1}$. Furthermore, $\boldsymbol{\tau}^{e}$ relies on the covariant metric tensor $\mathbf{G}$ (constructed in Section 3.2). Using $\mathbf{G}$, which is invariant under permutations of the element nodes, ensures that $\boldsymbol{\tau}^{e}$ is invariant as well. In the FST case on non-moving domains, the full metric tensor $\mathbf{G}$ can be decomposed into a spatial part $\left[\mathbf{G}\right]_{{n_{\mathrm{sd}}}\times{n_{\mathrm{sd}}}}$ and a temporal part $\left(\frac{\partial\theta}{\partial t}\right)^{2}=\frac{4}{\Delta t^{2}}$ as

[TABLE]

For SST and UST discretizations, the full metric tensor $\mathbf{G}$ is computed directly from the simplex space-time element with Equation (28). The constant $C_{inv}$ scales the diffusive contribution to $\boldsymbol{\tau}^{e}$. The choice of

[TABLE]

is motivated by an inverse estimate inequality for $\mathbb{P}_{1}$ simplex finite elements.[20] The reference element considered therein is the simplex reference element of the spatial discretization common for FST, SST, and UST discretizations. In case of steady simulations, the indices $m$ and $p$ in Equations (43) and  (46) run only over the spatial dimensions.

The principal square root of the $4\times 4$ or $5\times 5$ matrix in Equation (46) is computed using the direct Schur decomposition method [21] provided through the NAG library.[22] The nonlinear equation system (42) is linearized using the Newton–Raphson technique. The resulting linear equation system is solved with a parallel, preconditioned GMRES implementation.[23]

4 Numerical examples

In the following numerical examples we will use the specific gas constant $R=287\frac{J}{kgK}$ and the ratio of specific heats $\gamma=1.4$. These are the values for air at standard conditions.[24] In the viscous computations in Sections 4.2, 4.3, and 4.4, the temperature dependence of the dynamic viscosity $\mu(T)$ is modeled using Sutherland’s relation

[TABLE]

with $\mu_{ref}=21.7\cdot 10^{-6}$ Pa s, $T_{ref}=373.15$ K, $C=120$ K. The coefficient of thermal conductivity $\kappa$ is coupled to the dynamic viscosity $\mu$ using the Prandtl number $Pr=\frac{\mu}{\kappa}\frac{\gamma R}{\gamma-1}=0.71$. For air and temperatures up to the order of 1000 K, it is sufficient to treat the Prandtl number as constant.[24]

4.1 Supersonic flow over flat plate

The two-dimensional supersonic viscous flow over a flat plate at free-stream $Re_{\infty}=1000$ and $Ma_{\infty}=3$ constitutes a standard test case for compressible flow simulations. The flow develops a curved shock and boundary layer along the flat plate (Figure 5). Already in 1972 the test case was solved numerically.[25] Later, the test case was used to study the stability and spatial accuracy of a space-time finite element formulation for compressilbe flows.[3] More recently, it has been calculated with the following specifications:[4, 5]

The compressible Navier-Stokes equations are solved on a rectangle with $-0.2\leq x\leq 1.2$ and $0\leq y\leq 0.8$. A modified Sutherland’s law

[TABLE]

models the temperature dependence of the viscosity. On the left ($x=0.2$) and top ($y=0.8$) boundary, the free-stream conditions are prescribed for all degrees of freedom

[TABLE]

On the symmetry line ($y=0$ and $x<0$), the symmetry conditions $u_{2}=\tau_{12}=q_{2}=0$ are applied, rendering the viscous terms in the Neumann boundary integral $\int_{P_{n}}\mathbf{W}^{h}\cdot{\boldsymbol{\mathit{h}}}^{h}\,dP$ null. On the solid wall ($y=0$ and $x\geq 0$), the no-slip condition $\mathbf{u}=\mathbf{0}$ and the stagnation temperature of $T_{W}=T_{\infty}\left(1+\frac{\gamma-1}{2}{Ma_{\infty}}^{2}\right)=7.754\times 10^{-4}$ are prescribed as Dirichlet boundary condition. Finally, along the outflow boundary on the right ($x=1.2$) no Dirichlet boundary condition is prescribed. However, we do include the viscous terms ($\tau,q\neq 0$) in the boundary integral $\int_{P_{n}}\mathbf{W}^{h}\cdot{\boldsymbol{\mathit{h}}}^{h}\,dP$ along the outflow boundary, which we find to stabilize the flow in the subsonic part of the outflow. Figure 4 shows the velocity vectors in the bottom right corner of the computational domain, where the subsonic boundary layer touches the outflow boundary. Neglecting the viscous contributions to the boundary integral, we observe a small flow recirculation (see Figure 4a). Including the viscous contributions, a flow field with strictly aligned velocity vectors is obtained (see Figure 4b). Similar findings were reported for the compressible flow simulation with entropy variables.[3]

To compare the results with the findings of Xu et al.,[5] the computations are performed on three uniform meshes with 22400, 89600, and 358400 triangular elements, named ’Coarse’, ’Medium’ and ’Fine’. For the temporal discretization, the flat space-time method (FST) is employed. Figure 5 shows the isocontours of pressure, velocity magnitude, and temperature on the fine mesh. The curved shock and boundary layer are accurately resolved in all solution variables. The pressure coefficient $C_{p}=\frac{2(p-p_{\infty})}{\rho_{\infty}\|\mathbf{u}_{\infty}\|^{2}}$ along the solid wall is shown for the three meshes in Figure 5d. Except for a small kink close to the singularity at the leading edge of the plate, which vanishes with mesh refinement, the solutions of all three meshes are nearly identical and agree well with the results in the literature.[4, 5]

4.2 Pressure pulse

To evaluate the temporal accuracy of the three discretizations methods (FST, SST, UST), we derive the following pressure pulse test case in one space dimension. The test case takes into account the nonlinear characteristics of a sound wave[26] traveling in air at standard temperature and pressure (STP)

[TABLE]

As shown in Figure 6 the computational domain spans four wave lengths in $x$-direction, $0\leq x\leq 4\lambda$. By means of the Heaviside function, $\mathcal{H}(x)$, the perturbation function

[TABLE]

defines the initial distribution of the speed of sound as

[TABLE]

Using isentropic relations, initial conditions for the pressure-primitive variables are defined as

[TABLE]

After the time period $t_{f}=\frac{\lambda}{c_{0}}$, the head and tail of the wave have traveled exactly one wave length $\lambda$. Using the theory on shock formation in one-dimensional unsteady flows,[26] we compute the time $t^{*}$ after which a shock first forms as

[TABLE]

Therein, $\min\left(\frac{\partial c}{\partial x}\right)$ corresponds to the steepest descent along the wave. For our test case, we further request $t_{f}=\frac{1}{2}t^{*}$ to ensure that the nonlinear behavior is clearly visible, but the gradients of the solution remain finite. This condition yields $\alpha=\frac{1}{24\pi}$ in the perturbation function $S(x)$ in Equation (53).

At $x=0$ and $x=4\lambda$, the conditions of the fluid at rest are prescribed as Dirichlet boundary conditions

[TABLE]

For an inviscid fluid, the nonlinear acoustics theory can be used to obtain a reference solution.[26] The reference solution is constructed by advancing the all points of the wave by $(u+c)t_{f}$. Note that $u$ and $c$ vary along $x$. Figure 6 displays the initial condition for the pressure $p(x,t=0)$ and compares the nonlinear reference solution $p(x,t=t_{f})$ to the linear case.

To give a first example of a solution on an unstructured space-time (UST) mesh, Figure 7 shows the solution obtained on an extremely coarse UST mesh. The space dimension is on the horizontal axis and the time dimension on the vertical axis.

Lacking an implementation of the space-time methods for a single space dimension, we discretize the quasi one-dimensional spatial computational domain ($0\leq x\leq 4\lambda$, $0\leq y\leq\frac{\lambda}{100}$) by splitting $400\times 1$ square elements into right triangles. At $y=0$ and $y=\frac{\lambda}{100}$, the discretization has 100 elements per wave length. For all three space-time methods, the spatial discretization of the initial and final time level are identical. The temporal discretization size $\Delta t$ is chosen such that a CFL-number of $C_{\Delta t}=\frac{\Delta tc_{0}}{\Delta x}=2$ is obtained and the temporal discretization error dominates the spatial discretization error.

Looking at the complete wave after travelling one wave length (Figure 8 on the left), the solutions on all three space-time discretizations coincide with the inviscid nonlinear reference solution. Zooming in on the head of the wave (Figure 8 on the right), a small undershoot of approximately 0.2% is visible for the FST and SST solution. The undershoot of the UST solution is significantly smaller (about 0.04%). These results illustrate the good temporal accuracy of all three space-time discretization methods.

4.3 Valve

The UST method allows to solve the compressible Navier–Stokes equations on spatial computational domains that undergo topological changes over time. Namely, with space-time finite elements one can discretize the space-time continuum of a spatial computational domain, which splits into two and reconnects over time, conformingly. This is demonstrated here using a two-dimensional valve test case in which the fluid domains are separated by a valve member. The dimensions and conditions of this test case are inspired by the gas flow through the narrow gaps in the piston ring pack of an internal combustion engine (see Section 4.4). The problem setup is displayed in Figure 9. The fluid domain of the valve consist of a two-dimensional channel of 15 $\mu$m length and 4 $\mu$m height and a valve member that creates a 2.5 $\mu$m gap with the valve bottom in the open configuration. In the closed configuration, the gap between the valve member and the channel bottom vanishes completely ($h_{g}=0$).

The boundary conditions are as follows. On the left domain boundary, the pressure $p_{in}=1.2\times 10^{5}Pa$ and the temperature $T_{in}=373.15K$ is prescribed; the tangential component of the velocity is set to zero $u_{2}=0\frac{m}{s}$. On the right domain boundary, the pressure $p_{out}=1\times 10^{5}Pa$ is prescribed and the tangential component of the velocity is set to zero $u_{2}=0\frac{m}{s}$. On these two open boundaries, the viscous contributions ($\tau,q\neq 0$) are included in the boundary integral $\int_{P_{h}}\mathbf{W}\cdot{\boldsymbol{\mathit{h}}}\,dP$. On the solid walls of the channel bottom, top and the valve member, no-slip boundary conditions are applied and the temperatures $T_{B}=393.15K$, $T_{T}=383.15K$, and $T_{VM}=403.15K$ are prescribed, respectively.

The temporal evolution of the valve cycle is characterized by the gap height $h_{g}(t)$ between valve member and channel bottom plotted in Figure 10. Initially, the valve is held open for 1 $\mu$s, then closed over the following 4 $\mu$s. Subsequently, it is completely shut for 2 $\mu$s, before it is again opened over the next 4 $\mu$s, and finally held open for another 1 $\mu$s. The resulting space-time domain is discretized with the UST method, yielding the mesh shown in Figure 11.

In the open valve configuration, the pressure gradient of 0.2 bar leads to a laminar flow with a maximum Reynolds number $Re=\frac{\rho\,u_{max}\,h_{g}(12\,\mu s)}{\mu_{ref}}\approx 10$ and Mach number $Ma=\frac{u_{max}}{\sqrt{\gamma RT}}\approx 0.21$. The fluid density varies between $0.88\frac{kg}{m^{3}}$ and $1.12\frac{kg}{m^{3}}$. Figure 12 shows the pressure, velocity, and temperature distribution in the closed, half open, and fully opened valve. The renderings are obtained by evaluating the space-time data on the $x$-$y$-planes at $t=6,\,9,\,12\,\mu$s. In the closed configuration, the valve member completely separates the fluids to its left and right. The gas at the left attains the high pressure value prescribed at the inlet, the gas at the right the low pressure value prescribed at the outlet (Figure 12a). As the fluid is at rest (Figure 12b), the temperature distribution is governed by the diffusion between the prescribed values at the inlet and the solid walls (Figure 12c). In the open configurations, the pressure expands from the inlet on the left through the valve to the outlet on the right (Figure 12d and 12g). The fluid velocities for the fully opened configuration (Figure 12h) are significantly higher than the fluid velocities in the half opened configuration (Figure 12e). Figure 12i clearly shows the influence of the fluid flow on the temperature distribution in the valve. In summary, our valve simulation shows, that the UST method is capable to simulate flows on spatial domains with time-varying topology.

4.4 Blow-by past piston rings

This test case derives from a problem of engineering interest and demonstrates the potential to save computation time by local temporal refinement in simplex space-time meshes. As the fuel efficiency of an internal combustion engine is directly related to the performance of the piston ring pack, simulating the blow-by past piston rings is an active field of research.[27] Figure 13a relates the considered computational domain to the piston ring pack of an exemplaric piston. The domain for the gas flow computation consists of the voids in the piston ring pack. The geometry of the computational domain is simplified, such that the three piston lands are only connected through the ring end gaps. Still, scales that need to be resolved range from $\approx 0.1$ mm in the ring end gaps to $\approx 10$ cm (piston diameter) making the mesh generation a intricate task. The tetrahedral discretization of the spatial computaitonal domain (Figure 13a) is strongly adapted to the expected flow field, which leads to highly stretched elements. The mesh was generated with GMSH.[28]

The following boundary conditions characterize this test case. On the uppermost plane of the computational domain, the combustion chamber pressure during intake of $p_{in}=2.174\times 10^{5}Pa$ and a temperature of $T_{in}=352.3K$ are applied. The crank case conditions of $p_{out}=1\times 10^{5}Pa$ and $T_{out}=400K$ are prescribed on the lowest plane of the computational domain. On all solid walls no-slip conditions and a wall temperature of $T_{W}=400K$ are prescribed. A viable initial condition is obtained by gradually increasing the pressure at the inflow until the desired value of $p_{in}$ is reached.

At an engine speed of 6000 rpm, one degree crank angle is simulated. We compare three different time discretizations. First, we consider an FST discretization which has 100 time steps per degree crank angle, opposed to an FST discretization with 400 time steps. Additionally, we apply the time refinement strategy of Behr[9] to generate a four-dimensional SST discretization. The SST discretization has a temporal resolution that corresponds to 400 time steps per degree crank angle in the expected peak velocity region, while the resolution in the remaining spatial domain corresponds to 100 time steps per degree crank angle (see Figure 13d). Note that the spatial discretization for all three simulations is identical.

The flow field after one degree crank angle computed with 400 time steps on the FST discretization is shown in Figures 13b,13c, and 13e. Vortical structures and peak velocities are restricted to the vicinity of the ring end gaps. Away from the ring end gaps, the gas flow along the second piston land is approximately uniform in circumferential direction.

Figure 14 compares the pressure distribution on the cylinder wall below the upper ring’s end gap for the three different temporal discretizations. Despite the identical spatial discretization, flow features along the upper half of the piston land are resolved by the SST and fine FST discretization, but missing in case of the coarse FST discretization. In that sense, a sufficient temporal accuracy can be obtained with 400 time steps on the FST discretization or 100 time steps on the SST discretization. The four-dimensional time refinement allows us to reduce the number of time steps per degree crank angle (${n_{\mathrm{ts}}}$), offers an opportunity to save computational time, and motivates a closer look at the simulations’ timings.

All three simulations have been run in parallel on 120 cores of the Intel Xeon based RWTH cluster using an MPI parallelization. Table LABEL:tab:fall summarizes the mesh data, as well as timings for the three simulations. Comparing the number of the extruded tetrahedral elements (${n_{\mathrm{el}}}$) in the FST meshes to the number of pentatope elements in the SST mesh, one observes a drastic increase. The number of degrees of freedom (${n_{\mathrm{dof}}}$), in contrast, is increased only moderately by four-dimensional time refinement. The two main contributions to the total computation time $t_{T}$, are $t_{A}$, the time to assemble the linear equation system, and $t_{S}$, the time to solve the equation system. In sum, we observe a doubling of the computation time $t_{T}$ comparing the SST simulation to the coarse FST case. Comparing the SST simulation to the fine FST simulation, the computation time is reduced by 50%. This test case demonstrates how temporally refined simplex space-time finite element meshes can help to reduce the overall simulation time.

5 Conclusions

In this paper, we adopted an SUPG stabilized finite element formulation for compressible flows[5] to space-time finite elements. By incorporating a regular simplex element in the definition of the metric tensor $\mathbf{G}$, we obtained a stabilization matrix $\boldsymbol{\tau}$ that is invariant under permutations of the node numbering of the finite elements.

We performed compressible flow simulations on several simplex space-time meshes. Our formulation was validated with super- and subsonic test cases with Reynolds numbers ranging from 10 to 1000. SST and UST discretizations proved to be of comparable temporal accuracy as the well-established FST method.

On a valve test case, we demonstrated the capability of the UST method to simulate compressible flows on spatial computational domains undergoing topological changes. Four-dimensional temporal refinement in a pentatope discretization was successfully applied to the transient simulation of blow-by past piston rings of an internal combustion engine.

6 Acknowledgment

This work was supported by the German Federal Ministry for Economic Affairs and Energy through the Central Innovation Programme for Small and Medium Enterprises in the cooperation project ”Simulationstechnik Kolben-Kolbenring-Zylinder” under grant number KF3462301PO4. Computing resources were provided by JARA–Jülich Aachen Research Alliance and RWTH Aachen University IT Center. Furthermore, we want to thank Philipp Knechtges for his helpful remarks during the preparation of the manuscript.