Space-Time NURBS-Enhanced Finite Elements for Solving the Compressible Navier-Stokes Equations

TL;DR

This paper introduces NEFEM, a finite element method that uses NURBS for boundary representation, improving the simulation of compressible Navier-Stokes equations by integrating CAD data directly into the analysis.

Contribution

The paper presents NEFEM, a novel finite element approach that incorporates NURBS boundary representations from CAD tools without requiring volume splines, enhancing compressible flow simulations.

Findings

NEFEM accurately models boundary conditions using CAD-based NURBS.

NEFEM shows improved accuracy over conventional FEM in compressible flow problems.

Numerical examples demonstrate NEFEM's efficiency and precision.

Abstract

This article considers the NURBS-Enhanced Finite Element Method (NEFEM) applied to the compressible Navier-Stokes equations. NEFEM, in contrast to conventional finite element formulations, utilizes a NURBS-based computational domain representation. Such representations are typically available from Computer-Aided-Design tools. Within the NEFEM, the NURBS boundary definition is utilized only for elements that are touching the domain boundaries. The remaining interior of the domain is discretized using standard finite elements. Contrary to isogeometric analysis, no volume splines are necessary. The key technical features of NEFEM will be discussed in detail, followed by a set of numerical examples that are used to compare NEFEM against conventional finite element methods with the focus on compressible flow.

Tables

Table 1. Table 1: Flow conditions for the supersonic flow around the cylinder.

Flow conditions
\svhline Mach	1.7
Re	$2.0\times {10}^5$
${\rho }_{in}$	1.0
$u_{in}$	1.0
$v_{in}$	0.0
$e_{in}$	1.1179

Table 2. Table 2: Grids used for the grid refinement study.

Grid #	$n_{en}$	$n_{e{n}_{wall}}$
\svhline 0	6.72K	64
1	26.88K	128
2	107.52K	256
3	430.08K	512
4	1.72M	1028

Table 3. Table 3: Flow conditions for the transonic flow around the NACA0012 airfoil.

Flow conditions
\svhline Mach	0.8
$\rho$	1.0
$u$	1.0
$v$	0.0
$e$	3.29

Equations

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial {\mathbf{F}}_i}{\partial x_i}-\frac{\partial {\mathbf{E}}_i}{\partial x_i}$$

$$\begin{gathered} {\mathbf{F}}_i=\left(\begin{array}{c}{u}_i\rho \\ u_i\rho u_1+{\delta }_{i1}p\\ u_i\rho u_2+{\delta }_{i2}p\\ u_i\rho u_3+{\delta }_{i3}p\\ u_i(\rho e+p)\end{array}\right), \\ {\mathbf{E}}_i=\left(\begin{array}{c}0\\ {\tau }_{i1}\\ {\tau }_{i2}\\ {\tau }_{i3}\\ -q_i+{\tau }_{ik}{u}_k\end{array}\right), \end{gathered}$$

$$\mathbf{U}\cdot {\mathbf{e}}_d$$

$$(n_i{\mathbf{E}}_i)\cdot {\mathbf{e}}_d$$

$$\mathbf{U}(\mathbf{x},0)$$

$$\displaystyle\frac{\partial\mathbf{U}}{\partial t}+\mathbf{A}_{i}\frac{\partial\mathbf{U}}{\partial x_{i}}-\frac{\partial}{\partial x_{i}}\left(\mathbf{K}^{h}_{ij}\frac{\partial\mathbf{U}}{\partial x_{j}}\right)=\mathbf{0}\quad\mbox{on $\Omega_{t}\quad\forall\>t\in(0,T)$},$$

$$({\mathcal{S}}_U^h{)}_n$$

$$({\mathcal{V}}_W^h{)}_n$$

$${\int }_{Q_n}{\mathbf{W}}^h$$

$$\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\mathbf{\tau }}_{mom}\left[({\mathbf{A}}_k^h{)}^T\frac{\partial {\mathbf{W}}^h}{\partial x_k}\right]\cdot \left[\frac{\partial {\mathbf{U}}^h}{\partial t}+{\mathbf{A}}_i^h\frac{\partial {\mathbf{U}}^h}{\partial x_i}-\frac{\partial }{\partial x_i}\left({\mathbf{K}}_{ij}^h\frac{\partial {\mathbf{U}}^h}{\partial x_j}\right)\right]dQ+$$

$$\sum _{e=1}^{(n_{el}{)}_n}{\int }_{Q_n^e}{\mathbf{\tau }}_{DC}\left(\frac{\partial {\mathbf{W}}^h}{\partial x_i}\right)\cdot \left(\frac{\partial {\mathbf{U}}^h}{\partial x_i}\right)dQ+$$

$${\int }_{{\Omega }_n}({\mathbf{W}}^h{)}_n^{+}\cdot \left(({\mathbf{U}}^h{)}_n^{+}-({\mathbf{U}}^h{)}_n^{-}\right)d\Omega ={\int }_{(P_n{)}_h}\mathbf{W}\cdot {\mathbf{h}}^{h}dP.$$

$$\mathbf{C}(\xi )=\sum _{i=1}^n{\mathbf{R}}_i^p(\xi ){\mathbf{B}}_i,$$

$$\mathbf{\Phi }(s,r)=(1-s-r)x_3+(s+r)\mathbf{C}\left(\frac{s{\xi }_1+r{\xi }_2}{s+r}\right).$$

Full text

The key technical features of NEFEM will be discussed in detail, followed by a set of numerical examples that are used to compare NEFEM against conventional finite element methods with the focus on compressible flow.

1 Introduction

Geometries in engineering applications are commonly designed with the use of Computer-Aided-Design (CAD) tools. In general, these tools utilize Non-Uniform Rational B-Splines (NURBS) to accurately represent complex geometric domains by means of surface splines. When performing numerical analysis on such domains, it is common practice to first discretize the domain into finite sub-domains or elements. This discretization process, typically results in loss of the exact geometry.

An alternative approach, known as Isogeometric Analysis (IGA), was proposed in hughes2005 . The key idea of IGA is to use the NURBS basis functions not only for the geometric representation, but also for the numerical solution itself. By doing so, numerical analysis can be applied to the CAD model directly without the loss of geometric accuracy caused by discretizing the computational domain. Numerical analysis of fluid flow problems, however, commonly involves complex three-dimensional volume domains. Parametrizing such domains using closed volume splines can be challenging.

An alternative was proposed in sevilla2008 , and further extended for space-time finite elements and free-surface flows in stavrev2016 . This method was then modified for interface-coupled problems in hosters2018 . This approach suggests to use standard finite elements in the interior of the computational domain supplemented with so-called NURBS-enhanced finite elements along domain boundaries. These elements make use of NURBS to accurately represent complex geometries. The NURBS-Enhanced Finite Element Method (NEFEM) allows for maintaining as much as possible the proven computational efficiency of standard finite element methods, while utilizing the accurate geometric representation provided by the NURBS.

In this work, we apply NEFEM to supersonic flow problems. For this type of problems, accurate geometry representation can be important, especially due to the presence of shock waves and their interaction with solid walls.

2 Quasi-Linear Form of the Navier-Stokes equations

Before presenting the NEFEM concept, first the governing Navier-Stokes equations are presented. For this, let $\Omega_{t}\subset\mathbb{R}^{n_{sd}}$ and $t\in(0,T)$ be the spatial and temporal domains respectively, and let $\Gamma_{t}$ denote the boundary of $\Omega_{t}$. Then the model problem, written as a generalized advective-diffusive system, is given by:

[TABLE]

where $\mathbf{U}=(\rho,\rho u_{1},\rho u_{2},\rho u_{3},\rho e)^{T}$ is the solution vector. $\rho$, $u_{i}$, and $e$ represent the density, velocity components, and total energy per unit mass respectively. For the three-dimensional case, the Euler and viscous flux vectors $\mathbf{F}_{i}$ and $\mathbf{E}_{i}$ are defined as:

[TABLE]

where $q_{i}$, $\tau_{ik}$, $p$, and $\delta$ represent the heat flux, viscous stress tensor, pressure, and the Kronecker delta respectively. The boundary and initial conditions are given by:

[TABLE]

where $(\Gamma_{t})_{g_{d}}$ and $(\Gamma_{t})_{h_{d}}$ are the subsets of $\Gamma_{t}$, $\mathbf{e}_{d}$ is a basis in $\mathbb{R}^{n_{dof}}$, and $n_{dof}$ is the number of degrees of freedom. The quasi-linear form of Equation (1) is written as:

[TABLE]

where $\mathbf{A}_{i}=\frac{\partial\mathbf{F}_{i}}{\partial\mathbf{U}}$ represent the Euler Jacobians, and $\mathbf{K}^{h}_{ij}\frac{\partial\mathbf{U}}{\partial x_{j}}=\mathbf{E}_{i}$ the diffusivity matrices. $\mathbf{A}_{i}$ and $\mathbf{K}^{h}_{ij}$ are defined according to the set of solution variables (conservation variables in this case). For a detailed discussion on the various variable sets and corresponding $\mathbf{A}_{i}$ and $\mathbf{K}^{h}_{ij}$ matrices see, e.g., hauke1998 .

3 Stabilized Space-Time Finite Element Formulation

Following aliabadi1995 , the deformable spatial domain/stabilized space-time (DSD/SST) finite element formulation is derived for the quasi-linear form in equation (6).

In order to construct the finite element function spaces for the DSD/SST formulation, the time interval $(0,T)$ is decomposed into subintervals $I_{n}=(t_{n},t_{n+1})$, where $t_{n}$ and $t_{n+1}$ are part of the ordered series: $0=t_{0}<t_{1}<\cdots<t_{N}=T$. Additionally, we define $\Omega_{n}=\Omega_{t_{n}}$, and $\Gamma_{n}=\Gamma_{t_{n}}$. $Q_{n}$ now represents a so-called space-time slab, which is the domain enclosed by $\Omega_{n}$, $\Omega_{n+1}$ and $P_{n}$. Here, $P_{n}$ is the surface described by the boundary $\Gamma_{t}$ along interval $I_{n}$ (see Figure 1).

Similar to $\Gamma_{t}$ given in Section 2, $P_{n}$ can be decomposed into $(P_{n})_{g_{d}}$ and $(P_{n})_{h_{d}}$. The discrete finite element space-time function spaces for the trial and weighting functions $\mathbf{U}$ and $\mathbf{W}$ are given by:

[TABLE]

Using the Streamline-Upwind Petrov-Galerkin (SUPG) formulation, the weak form of Equation (6) states: given $(\mathbf{U}^{h})_{n}^{-}$, find $\mathbf{U}^{h}\in(\mathcal{S}_{U}^{h})_{n}$, such that $\forall\>\mathbf{W}^{h}\in(\mathcal{V}_{W}^{h})_{n}$:

[TABLE]

Here, the first two integrals on the left-hand side and the integral on the right-hand side represent the standard Galerkin form. The third and fourth left-hand side integrals represent the SUPG stabilisation and shock-capturing terms respectively. Continuity over the time-slab interface $\Omega_{n}$ is weakly imposed by the jump term, i.e., the fifth left-hand side integral in Equation (9). Here, the $\pm$ subscripts refer to the upper and lower time-slab solutions at time $n$. Equation (9) is solved sequentially for all space-time slabs $Q_{1},Q_{2},\dotsc,Q_{N-1}$ with initial condition $(\mathbf{U}^{h})_{0}^{-}=\mathbf{U}_{0}$.

For $\boldsymbol{\tau}_{mom}$ in the SUPG stabilization term in Equation (9), the formulation proposed in hughes1986e was used. The shock-capturing parameter $\boldsymbol{\tau}_{DC}$ used in this work, is similar to that given in kirk2009 , which is a modification of the original definition presented in hughes1986e . For brevity, a detailed discussion on the stabilization and shock-capturing formulations is omitted. For an extensive discussion on DSD/SST finite elements for compressible flow problems including SUPG-stabilization and shock-capturing, please refer to hughes2010 .

4 NURBS-Enhanced Finite Elements

In this section, the NURBS-enhanced finite element method as proposed in hosters2018 will be presented. The key idea of this method, is to use a NURBS definition of the computational domain to enhance the finite elements along the domain boundary. On all remaining elements in the interior of the domain a standard finite element formulation is used (cf. Figure 2).

Before discussing how the boundary elements make use of the NURBS boundary, let us first define a NURBS-curve of degree $p$. Such a curve is composed of piecewise rational basis functions $\mathbf{R}_{i}^{p}(\xi)$, and control points $\mathbf{B}_{i}$. The curve is then expressed by means of parametric coordinate $\xi\in(0,1)$ as follows:

[TABLE]

where $n$ denotes the total number of control points.

The elements that touch the NURBS domain boundary make use of a non-linear mapping between a reference element and the element in physical coordinates. This mapping, was proposed in hosters2018 as Triangle-Rectangle-Triangle (TRT) mapping. The mapping $\boldsymbol{\Phi}(s,r)$ is given by:

[TABLE]

Here, $s$ and $r$ are the parametric coordinates of the triangular reference element, $x_{3}$ denotes the physical coordinate of the interior node, and $\xi_{1}$ and $\xi_{2}$ are the parametric coordinates of the NURBS curve at which the element boundary nodes are located. A graphical representation of this mapping is shown in Figure 3.

By using the TRT mapping, the NURBS definition can be incorporated into the numerical analysis. As a result, the distribution of the integration points is determined from the exact geometry and not the erroneous discretized geometry (cf. Figure 4).

Furthermore, the shape functions corresponding to the interior nodes of the boundary elements remain zero along the NURBS (cf. Figure 4). This has the advantage that there is no contribution of the interior nodes when considering Dirichlet boundaries or boundary integrals. Especially for interface-coupled problems this can be important, where Dirichlet boundaries and boundary integrals are used to compute the coupling conditions hosters2018 .

5 Numerical Examples

To demonstrate the performance of the NEFEM in comparison to standard finite elements (SFEM), two test cases are considered next: 1) 2D supersonic viscous flow around a cylinder; 2) 2D transonic inviscid flow around a NACA0012 airfoil.

5.1 Cylinder Flow

The supersonic flow around a 2D cylinder is computed using the NEFEM and the SFEM. The flow conditions and the computational domain are shown in Table 5.1 and Figure 5 respectively. For the NEFEM computations, the cylinder is represented by a second order NURBS-curve (cf. Figure 5).