Split form ALE discontinuous Galerkin methods with applications to under-resolved turbulent low-Mach number flows

https://doi.org/10.1016/j.jcp.2020.109726Get rights and content

Highlights

  • Construction of entropy stable or kinetic energy dissipative moving mesh DG methods.

  • Analysis of discrete internal and kinetic energy balance for split forms.

  • Split form DG methods efficiently reduce aliasing in under-resolved simulations.

  • Successful application to challenging implicit large eddy simulation test case.

Abstract

The construction of discontinuous Galerkin (DG) methods for the compressible Euler or Navier-Stokes equations (NSE) includes the approximation of non-linear flux terms in the volume integrals. The terms can lead to aliasing and stability issues in turbulence simulations with moderate Mach numbers (Ma0.3), e.g. due to under-resolution of vortical dominated structures typical in large eddy simulations (LES). The kinetic energy or entropy is elevated in smooth, but under-resolved parts of the solution which are affected by aliasing. It is known that the kinetic energy is not a conserved quantity for compressible flows, but for small Mach numbers minor deviations from a conserved evolution can be expected. While it is formally possible to construct kinetic energy preserving (KEP) and entropy conserving (EC) DG methods for the Euler equations, due to the viscous terms in case of the NSE, we aim to construct kinetic energy dissipative (KED) or entropy stable (ES) DG methods on moving curved hexahedral meshes. The Arbitrary Lagrangian-Eulerian (ALE) approach is used to include the effect of mesh motion in the split form DG methods. First, we use the three dimensional Taylor-Green vortex to investigate and analyze our theoretical findings and the behavior of the novel split form ALE DG schemes for a turbulent vortical dominated flow. Second, we apply the framework to a complex aerodynamics application. An implicit LES split form ALE DG approach is used to simulate the transitional flow around a plunging SD7003 airfoil at Reynolds number Re=40,000 and Mach number Ma=0.1. We compare the standard nodal ALE DG scheme, the ALE DG variant with consistent overintegration of the non-linear terms and the novel KED and ES split form ALE DG methods in terms of robustness, accuracy and computational efficiency.

Introduction

The numerical simulation of under-resolved turbulent flows in the regime of low or moderate Mach numbers (the local Mach number is Ma0.3) requires accurate dispersion and dissipation behavior [19], [48], [51]. It is desirable that the dissipation errors are very low for well resolved scales and are very high for scales close to the Nyquist cutoff, to get rid of small scale noise. In addition, the dispersion error should be small for a wide range of scales. This motivates the application of high order methods, e.g. discontinuous Galerkin (DG) methods [8], [24], [34], compact finite difference methods [41] or summation-by-parts finite difference methods [52], to simulate the transition to turbulence and fully turbulent compressible flows.

DG methods are a class of finite element methods using piecewise polynomials as basis functions and Riemann solver based interface numerical flux functions along the element interfaces. The numerical dissipation caused by the interface Riemann solver acts as a filter of high frequency solution components. Furthermore, high order DG methods have excellent dispersion behavior [20], [24]. These observations motivate the application of DG methods in implicit large eddy simulations (iLES). For instance in [6], [14], [15], [49], [62], [66], DG iLES were used to simulate flows at moderate Reynolds numbers. In all these simulations no additional explicit sub-grid-scale models were used, since the high frequency type dissipation of the numerical interface flux functions is interpreted as a dissipative implicit sub-grid-scale model.

In many cases, the iLES DG simulations are negatively affected by aliasing errors due to the strong non-linearity of the flux functions [4], [31]. These errors are generated in the volume integrals and cannot be sufficiently controlled with the dissipation of the numerical surface fluxes alone. A careful treatment of aliasing is however necessary, as these errors may even cause fatal failure of the simulation. Standard de-aliasing techniques in the DG framework are based on either projection of the non-linear flux functions, e.g., [14], [69] or exact evaluations of the integrals in the variational formulation, sometimes termed “overintegration” or “consistent integration”, e.g., [19], [31], [35], [44].

The nodal DG spectral element method (DGSEM), e.g., [34], is constructed with local tensor-product Lagrange polynomial basis functions computed from Legendre-Gauss-Lobatto (LGL) points. The collocation of interpolation and quadrature nodes is used in the spatial discretization. This approach ensures that the derivative matrix in the DGSEM provides a summation-by-parts (SBP) operator [17]. A SBP operator gives a discrete analogue of the integration-by-parts formula [12], [17], [38] and thus ideas from the continuous stability analysis can be mimicked at the discrete level. The SBP property is a powerful tool to construct the numerical approximation in a way that aliasing issues are avoided.

The equations of gas dynamics, e.g. the Euler and Navier-Stokes equations (NSE), are equipped with a mathematical entropy equal to the scaled negative thermodynamic entropy. Hence, the mathematical model correctly captures the second law of thermodynamics [1]. A numerical scheme for these equations should reflect the properties of the entropy on the discrete level. In [64] Tadmor gave a discrete entropy criterion to construct entropy conservative (EC) Riemann solver based interface two point flux functions for first oder finite volume (FV) and finite difference (FD) methods to solve systems of conservation laws. This criterion was used to construct EC two point flux functions for the Euler equations in [7], [26], [58]. These low order EC methods can be modified by adding dissipation to the numerical fluxes such that the entropy is decreasing for all times. Then the method becomes entropy stable (ES) such that the entropy of the solution is bounded. The kinetic energy is another important physical quantity for the dynamics of turbulent flows. Therefore, a numerical scheme should also accurately capture its evolution and the transfer of kinetic and inner energy. Jameson [27] gave a discrete criterion to construct kinetic energy preserving (KEP) two point flux functions for first oder FV/FD methods to solve the Euler equations. However, in under-resolved turbulent parts of the simulation, it might be necessary to increase the robustness of the discretization by adding dissipation to the KEP numerical fluxes. Here a suitable dissipation term should be chosen such that the kinetic energy is guaranteed to be dissipated by the scheme (no pile up of energy in higher frequencies). In the current work, methods with this property will be denoted kinetic energy dissipative (KED) methods. It should be noted that a method for the compressible NSE can only be KDE or ES and not KEP or EC, since the model contains viscous and dissipative terms in the momentum and energy evolution.

Gassner [18] used the SBP property to construct a KEP DGSEM for the one dimensional Euler equations. This method is a discrete analogue of a quasilinear skew-symmetric formulation for the Euler equations and can be written as the standard DGSEM with an additional source term that acts as a de-aliasing mechanism. It is important to note that the skew-symmetric formulation is still fully conservative and satisfies the Lax- Wendroff theorem. For under-resolved computations, the source term is active and accounts for the error of the discrete product rule. On the other hand in smooth parts of the solution without large variation the term almost vanishes to zero, since the approximation is consistent. The stabilizing effect of high order conservative skew-symmetric or skew-symmetric-like schemes in under-resolved computations was also observed in other publications, e.g. [11], [29], [39], [47], [56]. However, since the thermodynamic entropy is a logarithmic function depending on density and pressure for the Euler equations and NSE, it is difficult to find an explicit discrete skew-symmetric or skew-symmetric-like formulation.

A framework to construct high order EC schemes in periodic domains has been given by LeFloch et al. [40]. Fisher and Carpenter [13] combined this approach with SBP operators and proved that two-point EC fluxes can be used to construct high order schemes when the derivative approximations in space are SBP operators. Gassner et al. [21], [22] showed that skew-symmetric-like (split form) DGSEM formulations for the Euler equations can be discretely recovered when specific numerical volume fluxes in the flux form volume integral of Fisher and Carpenter are chosen. In particular, the three dimensional version of the skew-symmetric KEP DGSEM in [18] can be recovered by using the two point flux from Morinishi et al. [47] as volume flux. Hence, ES DGSEM can be constructed when the derivative matrix satisfies the SBP property and a two-point EC flux function is used as volume flux in the flux form volume integral of Fisher and Carpenter. The construction of KED DGSEM requires more attention, since it has been shown numerically that there are two point flux functions, e.g. [7], which are KEP in the sense of Jameson [27], but when they are used as volume flux in the DGSEM the kinetic energy is necessarily not preserved (see [21], [58] or Fig. 4.2 in Section 4.2 of this work).

It should be mentioned that the described approach and the associated analysis is for semi-discrete ES or KED high order DGSEM. Moreover, the construction of these methods relies on the assumption that the discrete density and pressure are positive. In general, this might not always be the case for high order methods, since the positivity of the discrete density and pressure can be affected by spurious oscillations in the numerical solution. This phenomenon is for instance observed for strong shock waves which appear in the regime of high Mach numbers, or severe under-resolution of strong shear-layers.

Another important component for the simulation of turbulent flows is an adaptive discretizations, where resolution is increased in regions with large spatial variations. The r-adaptive method involves the re-distribution of the mesh nodes in regions of rapid variation of the solution [65]. In comparison with h-adaptive discretizations, where the mesh is refined and coarsened by changing the number of elements in the tessellation, the r-adaptive method has some advantages, e.g. no hanging nodes appear and the number of elements does not change. On the other hand a r-adaptive method can be only used when the effect of mesh movement is appropriately accounted for the discretization of the system of conservation laws. This can be done by an Arbitrary Lagrangian-Eulerian (ALE) approach [10]. In the last decades several ALE DG methods have been developed, e.g. [36], [42], [46], [50], [55], [68]. Furthermore, in [60] a provably ES moving mesh ALE DGSEM for the three dimensional Euler equations on curved hexahedral elements and in [70] a moving mesh ES spectral collocation scheme for the three dimensional NSE were constructed. The current work is a natural extension of the authors work presented in [60]. The focus herein is on (i) extending the approach to KEP/KED ALE DG methods based on general split forms of the compressible Euler equations, (ii) extending the approach to the viscous terms of the Navier-Stokes equations based on the Bassi and Rebay ansatz [2], and (iii) demonstrating the efficiency of the resulting novel KED ALE DG approach for a complex large scale application.

The structure of the present work is as follows: First, the approaches [22], [60] are used to construct a provably ES ALE DGSEM for the three dimensional NSE. Then a provable KEP ALE DGSEM for the Euler equations on curved moving elements is constructed. The discrete energy balance of the discrete kinetic, internal and total energy is analyzed for the split form ALE DGSEM when different well known split forms [7], [29], [39], [56], [58] are used. Afterward, the constructed KEP DGSEM is combined with the approach in [22] to construct a KED ALE DGSEM for the three dimensional NSE. The theoretical findings and the constructed methods are numerically investigated for the Taylor-Green Vortex (TGV) problem [61]. Finally, the constructed ALE DGSEM is applied in an iLES for a complex aerodynamics application. An interesting area of application for problems with moving boundaries in aerodynamics is the moderate Reynolds number flow around pitching airfoils. These flows can serve as a model for the wing motion of natural flyers or the situation in micro unmanned vehicles. As a specific example, the transitional flow around a plunging SD7003 airfoil at Reynolds number Re=40,000 and Mach number Ma=0.1 is simulated. This setup, among others, was investigated by Visbal [67], which gives the possibility for a comparison to available experimental measurements and numerical simulation results from literature [16], [43], [53], [57], [71].

Section snippets

The Navier Stokes equations (NSE) in three dimensions

We apply the block vector notation as in [22] to present the three dimensional NSE. In Appendix A the key elements of this notation are briefly summarized. The three dimensional NSE in compact block vector form are given byut+xf(u)=xfv(u,xu). These equations are considered on a time-dependent domain Ω(t)R3 with suitable initial data and boundary conditions. The state vector is given byu=[ρ,ρu,ρe+12ρ|u|2]T, where ρ is the mass density, u=[u1,u2,u3]T the velocity, p is the

Building blocks for the spectral element approximation

The spectral element approximation is based on a nodal approach with Lagrange basis functions {j()}j=0N constructed from Legendre Gauss Lobatto (LGL) points {ξi}i=0N. We note that ξ0=1 and ξN=1. The Lagrange basis functions satisfy the cardinal propertyi(ξj)=δji. On the reference element E=[1,1]3 the solution and fluxes of the system (2.18a), (2.18b), (2.18c) are approximated by tensor product Lagrange polynomials of degree N, e.g.,J(ξ1,ξ2,ξ3,t)J(ξ1,ξ2,ξ3,t):=i,j,k=0NJijk(t)i(ξ1)j(ξ2)k

Numerical results

The proposed split form ALE DGSEM is implemented in the open source high order DG solver FLEXI3 [37]. It provides the necessary framework for the implementation of different split forms for high order unstructured meshes, was successfully applied to under-resolved simulations in fluid dynamics before [3], [15] and shows excellent scaling properties, making it a suitable choice for large-scale simulations, as presented in the next chapter.

Simulation of transitional flow past a plunging airfoil

In this section, a complex application of our novel split form ALE DGSEM is presented. As our test case, we choose the low Reynolds number flow around an airfoil that undergoes an unsteady plunging motion. The different fluid dynamics processes are important for both the understanding of flapping flight and the flows around micro unmanned vehicles. Such flows have thus attracted considerable attention from engineers and scientists in the past, but efficient simulation of those cases still

Conclusions

High order accurate DG methods might be affected by aliasing errors due to the non-linearity of the flux functions when solving under-resolved turbulent vortex dominated flows. One possibility to avoid aliasing issues in the discretization is the construction of KED or ES high order split form DG methods [14], [21], [22], [69].

In this work, KEP and EC high order ALE DGSEM for the Euler and KED and ES high order ALE DGSEM for the NSE were analyzed. Here the key element in the approximation is

CRediT authorship contribution statement

Nico Krais: Investigation, Software, Validation, Visualization, Writing - original draft, Writing - review & editing. Gero Schnücke: Conceptualization, Formal analysis, Methodology, Writing - original draft, Writing - review & editing. Thomas Bolemann: Investigation, Software. Gregor J. Gassner: Conceptualization, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Gero Schnücke and Gregor Gassner are supported by the European Research Council (ERC) under the European Union's Eights Framework Program Horizon 2020 with the research project Extreme, ERC grant agreement no. 714487. The authors gratefully acknowledge the support and the computing time on “Hazel Hen” provided by the HLRS through the project “hpcdg”.

References (71)

  • T.J. Hughes et al.

    A new finite element formulation for computational fluid dynamics: I. symmetric forms of the compressible Euler and Navier–Stokes equations and the second law of thermodynamics

    Comput. Methods Appl. Mech. Eng.

    (1986)
  • F. Ismail et al.

    Affordable, entropy-consistent Euler flux functions ii: entropy production at shocks

    J. Comput. Phys.

    (2009)
  • C.A. Kennedy et al.

    Reduced aliasing formulations of the convective terms within the Navier–Stokes equations for a compressible fluid

    J. Comput. Phys.

    (2008)
  • C.A. Kennedy et al.

    Low-storage, explicit Runge–Kutta schemes for the compressible Navier–Stokes equations

    Appl. Numer. Math.

    (2000)
  • R.M. Kirby et al.

    De-aliasing on non-uniform grids: algorithms and applications

    J. Comput. Phys.

    (2003)
  • Robert M. Kirby et al.

    De-aliasing on non-uniform grids: algorithms and applications

    J. Comput. Phys.

    (October 2003)
  • D.A. Kopriva et al.

    A provably stable discontinuous Galerkin spectral element approximation for moving hexahedral meshes

    Comput. Fluids

    (2016)
  • Y. Kuya et al.

    Kinetic energy and entropy preserving schemes for compressible flows by split convective forms

    J. Comput. Phys.

    (2018)
  • S.K. Lele

    Compact finite difference schemes with spectral-like resolution

    J. Comput. Phys.

    (1992)
  • I. Lomtev et al.

    A discontinuous Galerkin ALE method for compressible viscous flows in moving domains

    J. Comput. Phys.

    (1999)
  • G. Mengaldo et al.

    Dealiasing techniques for high-order spectral element methods on regular and irregular grids

    J. Comput. Phys.

    (2015)
  • C.A.A. Minoli et al.

    Discontinuous Galerkin spectral element approximations on moving meshes

    J. Comput. Phys.

    (2011)
  • Y. Morinishi

    Skew-symmetric form of convective terms and fully conservative finite difference schemes for variable density low-Mach number flows

    J. Comput. Phys.

    (2010)
  • R.C. Moura et al.

    Linear dispersion–diffusion analysis and its application to under-resolved turbulence simulations using discontinuous Galerkin spectral/hp methods

    J. Comput. Phys.

    (2015)
  • R.C. Moura et al.

    On the eddy-resolving capability of high-order discontinuous Galerkin approaches to implicit LES/under-resolved DNS of Euler turbulence

    J. Comput. Phys.

    (2017)
  • V.-T. Nguyen

    An arbitrary Lagrangian–Eulerian discontinuous Galerkin method for simulations of flows over variable geometries

    J. Fluids Struct.

    (2010)
  • A. Nigro et al.

    A low-dissipation DG method for the under-resolved simulation of low Mach number turbulent flows

    Comput. Math. Appl.

    (2019)
  • J. Nordström et al.

    Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier–Stokes equations

    J. Comput. Phys.

    (1999)
  • M. Parsani et al.

    Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations

    J. Comput. Phys.

    (2015)
  • P.-O. Persson et al.

    Discontinuous Galerkin solution of the Navier–Stokes equations on deformable domains

    Comput. Methods Appl. Mech. Eng.

    (2009)
  • S. Pirozzoli

    Generalized conservative approximations of split convective derivative operators

    J. Comput. Phys.

    (2010)
  • A.R. Winters et al.

    ALE–DGSEM approximation of wave reflection and transmission from a moving medium

    J. Comput. Phys.

    (2014)
  • A.R. Winters et al.

    A comparative study on polynomial dealiasing and split form discontinuous Galerkin schemes for under-resolved turbulence computations

    J. Comput. Phys.

    (2018)
  • N.K. Yamaleev et al.

    Entropy stable spectral collocation schemes for the 3-D Navier–Stokes equations on dynamic unstructured grids

    J. Comput. Phys.

    (2019)
  • T.J. Barth

    Numerical methods for gasdynamic systems on unstructured meshes

  • Cited by (12)

    • Study on the impact characteristics of underwater explosion bubble jets induced by plate structure

      2022, Ocean Engineering
      Citation Excerpt :

      In order to further obtain the influence of the plate on the water jet shape and angle, and obtain the plate stress and deflection induced by the water jet impact on the plate, the numerical analysis should be adopted. ALE algorithm, as an effective numerical method for studying explosion and impact problems (Cai et al., 2021; Krais et al., 2020; Zheng et al., 2022), can accurately deal with the complex multi-material interface problems in the underwater explosion. Therefore, the underwater explosion bubble jet characteristics and its impact characteristics on the plate structure at different explosion distances and azimuth angles are deeply studied by ALE method.

    View all citing articles on Scopus
    View full text