Split form ALE discontinuous Galerkin methods with applications to under-resolved turbulent low-Mach number flows
Introduction
The numerical simulation of under-resolved turbulent flows in the regime of low or moderate Mach numbers (the local Mach number is ) 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 and Mach number 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 by These equations are considered on a time-dependent domain with suitable initial data and boundary conditions. The state vector is given by where ρ is the mass density, 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 constructed from Legendre Gauss Lobatto (LGL) points . We note that and . The Lagrange basis functions satisfy the cardinal property On the reference element 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.,
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)
- et al.
A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations
J. Comput. Phys.
(1997) - et al.
The effect of the formulation of nonlinear terms on aliasing errors in spectral methods
Appl. Numer. Math.
(1996) - et al.
Conservative and entropy stable solid wall boundary conditions for the compressible Navier–Stokes equations: adiabatic wall and heat entropy transfer
J. Comput. Phys.
(2019) - et al.
High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows
J. Comput. Phys.
(2000) - et al.
Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations
Comput. Fluids
(2014) - et al.
High-order entropy stable finite difference schemes for nonlinear conservation laws: finite domains
J. Comput. Phys.
(2013) - et al.
On the use of kinetic energy preserving DG-schemes for large eddy simulation
J. Comput. Phys.
(2017) - et al.
Simulation of underresolved turbulent flows by adaptive filtering using the high order discontinuous Galerkin spectral element method
J. Comput. Phys.
(2016) - et al.
Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations
J. Comput. Phys.
(2016) On the symmetric form of systems of conservation laws with entropy
J. Comput. Phys.
(1983)
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.
Affordable, entropy-consistent Euler flux functions ii: entropy production at shocks
J. Comput. Phys.
Reduced aliasing formulations of the convective terms within the Navier–Stokes equations for a compressible fluid
J. Comput. Phys.
Low-storage, explicit Runge–Kutta schemes for the compressible Navier–Stokes equations
Appl. Numer. Math.
De-aliasing on non-uniform grids: algorithms and applications
J. Comput. Phys.
De-aliasing on non-uniform grids: algorithms and applications
J. Comput. Phys.
A provably stable discontinuous Galerkin spectral element approximation for moving hexahedral meshes
Comput. Fluids
Kinetic energy and entropy preserving schemes for compressible flows by split convective forms
J. Comput. Phys.
Compact finite difference schemes with spectral-like resolution
J. Comput. Phys.
A discontinuous Galerkin ALE method for compressible viscous flows in moving domains
J. Comput. Phys.
Dealiasing techniques for high-order spectral element methods on regular and irregular grids
J. Comput. Phys.
Discontinuous Galerkin spectral element approximations on moving meshes
J. Comput. Phys.
Skew-symmetric form of convective terms and fully conservative finite difference schemes for variable density low-Mach number flows
J. Comput. Phys.
Linear dispersion–diffusion analysis and its application to under-resolved turbulence simulations using discontinuous Galerkin spectral/hp methods
J. Comput. Phys.
On the eddy-resolving capability of high-order discontinuous Galerkin approaches to implicit LES/under-resolved DNS of Euler turbulence
J. Comput. Phys.
An arbitrary Lagrangian–Eulerian discontinuous Galerkin method for simulations of flows over variable geometries
J. Fluids Struct.
A low-dissipation DG method for the under-resolved simulation of low Mach number turbulent flows
Comput. Math. Appl.
Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier–Stokes equations
J. Comput. Phys.
Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations
J. Comput. Phys.
Discontinuous Galerkin solution of the Navier–Stokes equations on deformable domains
Comput. Methods Appl. Mech. Eng.
Generalized conservative approximations of split convective derivative operators
J. Comput. Phys.
ALE–DGSEM approximation of wave reflection and transmission from a moving medium
J. Comput. Phys.
A comparative study on polynomial dealiasing and split form discontinuous Galerkin schemes for under-resolved turbulence computations
J. Comput. Phys.
Entropy stable spectral collocation schemes for the 3-D Navier–Stokes equations on dynamic unstructured grids
J. Comput. Phys.
Numerical methods for gasdynamic systems on unstructured meshes
Cited by (12)
A high-order fluid–structure interaction framework with application to shock-wave/turbulent boundary-layer interaction over an elastic panel
2023, Journal of Fluids and StructuresStudy on the impact characteristics of underwater explosion bubble jets induced by plate structure
2022, Ocean EngineeringCitation 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.
On the entropy conserving/stable implicit DG discretization of the Euler equations in entropy variables
2022, Computers and FluidsAn efficient sliding mesh interface method for high-order discontinuous Galerkin schemes
2021, Computers and Fluids