A natural decomposition of viscous dissipation in DG methods for turbulent incompressible flows

In this note we aim at a characterisation of the discretisation of viscous dissipation which allows to distinguish `physical' (also frequently called `molecular', or `resolved') from `numerical' dissipation in DG-discretised incompressible flow simulations.


Introduction
Let us consider an incompressible Navier-Stokes problem without acting outer forcing terms and periodic or no-stress boundary conditions. Given a physical domain Ω ⊂ R d , the strong form of such a problem, equipped with a suitable initial condition u 0 : Ω → R d , reads ∂ t u − ν∆u + (u · ∇)u + ∇p = 0, subject to ∇ · u = 0. Here, u : (0, T ) × Ω → R d indicates the velocity field, p : (0, T ) × Ω → R is the (zero-mean) kinematic pressure, and the underlying fluid is assumed to be Newtonian with kinematic viscosity 0 < ν 1. We are especially interested in the situation where the corresponding Reynolds number is large enough such that a turbulent flow is expected and its approximation is performed in a strongly under-resolved setting.
For a finite element pair V h /Q h for velocity/pressure, and assuming that the simulation is performed up to the time instance T > 0, a typical DG scheme (in primal form) for discretising the Navier-Stokes equations reads as follows: The form a h treats the viscosity effects, c h is the nonlinear convection term, b h connects pressure and incompressibility condition and j h is a possible additional stabilisation and/or turbulence model [4,6]. In this note we focus on the viscous term a h .
Testing (1) symmetrically with (v h , q h ) = (u h , p h ) leads to the discrete kinetic energy balance The only physical dissipation process present in the original Navier-Stokes model is due to viscosity. Therefore, in our opinion, every additional energy-dissipating (or even energy-producing) mechanism, which is frequently incorporated in c h and j h , has to be characterised as an artificial (numerical) contribution.
To introduce a mathematically rigorous notion of viscous dissipation processes, let a phy h denote the non-negative part in the discretisation of the viscous term that represents physical dissipation, which is supposed to fulfil a phy h (u, u) = ∇u 2 L 2 for the exact solution u. We assume that the remainder of a h is a non-negative bilinear form a num h , which describes numerical dissipation in the discretisation of the viscous term. The decomposition is required to be consistent in the sense that a num h (u h , u h ) vanishes for h/k → 0, with u h being a discrete solution converging to u as h/k → 0. Here, h denotes the underlying mesh size and k is the polynomial order of discrete velocities belonging to V h .
Let us emphasise that the requirement that both parts of the decomposition be non-negative is a restriction and disallows some choices for a phy h which may seem intuitive at first glance; cf. [5,Sec. 2]. Being able to identify the physical dissipation, the total numerical dissipation ε tot h of the scheme can be defined as ε tot . This ε tot h then fulfils the reasonable and widely accepted expectation (for a meaningful discretisation) that it is non-negative; that is, ε tot h 0. Especially, the frequently used characterisation a phy can be misleading in DG methods as will be shown below. However, let us mention that the difference between different notions of physical dissipation in DG methods is only relevant in the under-resolved case. We will demonstrate that a lifting technique can be used to define a suitable/natural decomposition of the total viscous dissipation into a physical and a numerical contribution. In doing so, we restrict ourselves to the symmetric interior penalty (SIP) method as a very frequently used DG method.

A natural decomposition of viscous dissipation for DG methods
The SIP bilinear form is given by [4] where λ > 0 is a sufficiently large (due to a discrete inverse inequality) penalty parameter.
We can interpret the DG formulation (3) in a mixed setting, cf. [1], which gives a natural definition of a discrete diffusive flux (scaled with ν −1 ) σ h = σ h (u h ), which is defined element-wise for all τ h ∈ ∇ h V h K by the following operation on any The second equality is due to integration by parts and u h denotes a 'numerical trace' which characterises different DG methods, see [1, Table 3.1]. In the following, for simplicity, we only want to consider the SIP method where With the notion (5) of the lifting operator L L L, (4) can finally be used to obtain the characterisation σ h (u h ) = ∇ h u h −L L L( u h ). Using this definition of σ h , one can rewrite the symmetrically tested bilinear form a h from (3) as follows: We notice that the usual assumption on the parameter λ guarantees that both parts a phy h (u h , u h ) and a num h (u h , u h ) are nonnegative for any discrete function u h ∈ V h ; a detailed explanation for this statement can be found in [5,App. B]. Further, note that a num As shown in detail in [5, App. A] the bilinear form a phy h in (6) corresponds to the DG method by Bassi and Rebay [2]. It can be seen as a central flux approximation to diffusion/viscosity (to the corresponding first order system). Moreover, note that for SIP a piecewise constant function will not induce physical dissipation if exclusively the broken gradient is used for the definition of a phy h . In contrast, using the definition (6) [3] at Re = 1600. We consider four different choices for the SIP stabilization parameter, λ ∈ {λ * , 1.25λ * , 1.5λ * , 2λ * }. Here, λ * is the smallest parameter that guarantees stability based on usual arguments following trace inverse estimates. We refer to [5,Sec. 4] for details on the setup. Using the broken gradient ∇ h u h in the definition of the physical dissipation (left) leads to negative numerical viscous dissipation, while using σ h (right) as in (4) does not.
In Fig. 1 the difference between two different interpretations of physical and numerical dissipation is demonstrated for a highly under-resolved homogeneous decaying turbulence example. We observe that the broken gradient is not a good metric for characterising physical dissipation in that example as it -misleadingly -suggests that the numerical method contributes with negative numerical dissipation. Using σ h as in (4), on the other hand, leads to a physically meaningful interpretation.