Viscous Diffusion Effects in the Eigenanalysis of (Hybridisable) DG Methods

We present the first eigenanalysis of hybridisable discontinuous Galerkin (HDG) schemes for the advection-diffusion equation in one dimension. This study is also one of the first to include viscous diffusion effects in the eigenanalysis of discontinuous spectral element methods. The interplay between upwind dissipation and viscous diffusion is discussed and preliminary insights deemed relevant to (under-resolved) turbulence computation approaches are presented.


Introduction
When numerically solving partial differential equations, numerical errors are likely to impact not only solution accuracy, but also the stability/robustness of the computation. This is particularly the case in eddy-resolving approaches to turbulent flows, such as large-eddy simulation (LES) and direct numerical simulation (DNS). Also, in the so-called implicit LES / under-resolved DNS strategies [1], where numerical error (specifically dissipation) provides small-scale regularisation in lieu of a turbulence model, understanding the nature of numerical errors is crucial. These typically appear in the form of dispersion and diffusion errors, where the former distorts the solution, while the latter is responsible for its damping. A useful framework for the assessment of such numerical errors is the eigensolution analysis technique [2,3].
We present the first eigenanalysis of hybridisable discontinuous Galerkin (HDG) methods. This is also one of the first studies to consider viscous diffusion effects in the eigenanalysis of discontinuous SEM (spectral element methods), as it addresses the advection-diffusion equation in one dimension. Focus is given to the temporal analysis approach [2,5], which is suited for problems with periodic boundary conditions. The spatial analysis [3,4], suited for inflow-outflow problems, will be considered in subsequent studies. Here, we offer preliminary results on (i) the effects of the Peclét number (a cell-based Reynolds number), and (ii) the interplay between upwind (numerical) dissipation and viscous (physical) diffusion. We highlight how these results improve upon our understanding and practice of implicit LES / under-resolved DNS approaches.
We note that, although a non-modal eigenanalysis strategy better suited for turbulence computations has been recently proposed [6], the present work will focus on more fundamental aspects and follow therefore the classical eigenanalysis. Finally, the results presented here are representative of a broader class of discontinuous SEM, given the well established connections within this class-see e.g. [7]. This paper is organized as follows. Section 2 introduces the HDG discretisation as applied to the linear advection-diffusion equation in one dimension. Section 3 details the temporal eigenanalysis framework and presents our preliminary results. Finally, in Sect. 4, our conclusions are summarised and future research topics are outlined.

HDG Discretisation
In one dimension, the linear advection-diffusion equation is given by where the advection velocity a and the viscosity μ are positive constants. This equation can be written in conservation form through the flux function f (u, g) = au − μg, as the system where g is the auxiliary gradient variable. The discretisation procedure is similar to that of traditional DG methods.
After the (1D) physical domain is partitioned into non-overlapping elemental regions of size h, the numerical solution and its gradient are locally approximated by polynomial expansions in the form where φ j are polynomial basis functions of degree up to P , defined in the standard domain st = [−1, 1]. A linear mapping relation is assumed between the physical coordinate x of element and the coordinate ξ ∈ st . Multiplying Eqs.
(2)-(3) by φ i , integrating over element and applying integration by parts leads respectively to where and ⊕ denote the left and right boundaries of element , in that order. As typical, expansions in (4) are to be inserted into (5)- (6), which are then required to hold for i = 0, . . . , P . Note that the integrals above have been moved to st and interface quantities u and f have been introduced. The state average u is peculiar to HDG in that it represents a uniquely defined interface variable whose value stems indirectly from the enforced continuity of the numerical flux f . This continuity ensures local conservation for HDG methods, regardless of the chosen flux formula. For the advection-diffusion problem at hand, the interface fluxes on either side of a given element (cf. Fig. 1, left diagram) can be taken in the form Also, τ = β|a| + σ is a stabilisation constant combining an upwinding parameter β and a penalty term σ that accounts for the partially diffusive character of the model equation considered. This work however assumes σ = 0 as it focuses on advectiondominated cases, which are typically stable without the penalty term σ , even within the context of turbulence simulations [8].
Flux formulas (7)- (8) are inspired in Ref. [9]. In the case of pure advection (with σ = 0), the interface solution variable becomes the simple average u = u L ⊕ + u R of the adjacent states from the left (L) and right (R) elements sharing the considered interface. Under this case, it is also easy to show that the fluxes in (7)-(8) recover those used in traditional DG methods, whereby HDG exactly reproduces DG. This does not hold, however, when diffusion is taken into account, in which case u is only implicitly defined from the flux continuity condition enforced at interfaces, where g L ⊕ and g R depend on values of u at two other interfaces via (6). The diagram on the right-hand-side of Fig. 1 should help clarify the notation adopted.
Using vectorsû = {û 0 , . . . ,û P } T andĝ = {ĝ 0 , . . . ,ĝ P } T , the flux continuity condition (11) becomes (6) can be written as in which matrices M and D have been introduced, namely Finally, (5) becomes with Note that (12) is a scalar equation written from the point of view of a given interface, whereas (13) and (15) are vector equations written from the viewpoint of an arbitrary element of size h.
It is now convenient to eliminateĝ and work with variablesû and u alone. This can be done by solving (13) forĝ and substituting the resulting expression in both (12) and (15). The former substitution leads, after some algebra, to where Pe = |a|h/μ denotes the Péclet number, for which a uniform mesh spacing is assumed. Moreover, four scalar constants 'm' have been introduced, defined as In addition, the following matrices appear in (18) Note that (18) relates the solution vectorsû of two adjacent elements ( L and R ) with the three interface states u associated to the boundaries of these elements. The second step consists in usingĝ from (13) into (15), not forgetting to take the fluxes (16)-(17) into account. After some more algebra, one arrives at whose matrices now introduced are given by Note that (21) links the solution vectorû and its time derivative to the two interface variables u at the boundaries of the considered element.
In the actual context of simulations, (21) would be first solved (analytically) for u after an implicit time-stepping scheme is chosen. This is possible since it entails expressing dû/dt in terms ofû at the current as well as previous time levels. The next step would be to insert the resulting expression forû into (18), from which a scalar equation whose only unknowns are u at various interfaces is obtained. This equation is finally used for the assembly of a global system given suitable boundary conditions, which can be solved via direct or iterative techniques. Since the system's solution grants u for all interfaces,û can be obtained locally for each element from the time-discrete version of (21). The reader is referred to [9] for the details of this procedure. In this work, however, as we are interested in the eigenanalysis of HDG, a different strategy is adopted, as outlined next.

Temporal Eigenanalysis
In the eigenanalysis of spectral element methods [2,5], it is typical to assume wavelike solutions in the formû ∝ exp[i(κx − ωt)], wherebyû L =û exp(−iκh) and u R =û exp(+iκh). Here,û is the solution vector of a "central" element, whereasû L andû R refer to solution vectors of neighbouring elements from the left (L) and from the right (R), respectively. For the HDG formulation, an additional assumption can be made regarding a wave-like behaviour for u. We assume that u L = u exp(−iκ h) and u R ⊕ = u exp(+iκ h), where now u is the interface variable shared by two adjacent elements, whereas u L and u R ⊕ refer to interface variables at the nearest interfaces from the left/right (L/R). This second assumption is only natural given the connection betweenû and u. Actually, we now show that κ = κ, which is not surprising.
We start from (21) assuming wave-like behaviour forû, obtaining which uniquely definesû from u ⊕ and u . If the above is written for another element, say, the adjacent element from the right (a translation x → x + h), one has which then implies where (25) has been used on the right-hand side. Comparing the left-and rightmost expressions above leads to exp(iκ h) = exp(iκh), which means κ h = κh + 2nπ , for n integer. This phase ambiguity can be sorted out by the evaluation of the x-derivative of (25) at x + h, given by which yields κ = κ. This last step about the phase is, however, not really necessary to the eigenanalysis because only the complex exponential factors appear throughout the relevant equations, hence knowing that exp(iκ h) = exp(iκh) is sufficient.
In the remainder of the study, orthonormal Legendre basis functions are assumed, whereby M = I . We note that numerical dispersion and diffusion eigencurves, which are the focus of the study, do not change depending on the basis functions adopted, provided that exact integrations are used in the spatial discretisation.
In the temporal analysis, an eigenvalue problem is set where, given a realvalued wavenumber κ, multiple (P + 1) eigenvalues of the relevant eigenmatrix are associated to admissible complex-valued numerical frequencies ω = ω(κ). The procedure to obtain this eigenvalue problem is described below.
We begin from (18), assuming u L = u exp(−iκh) and u R ⊕ = u exp(iκh), to find Then, (29) is used into (21), relating the solution vectorû at a given element to the state vectors of its left (û L ) and right (û R ) neighbours. From the wave-like behaviour ofû and the relationsû L =û exp(−iκh) andû R =û exp(+iκh), one can arrive at where = ω/a and matrix Z = Z(κh; Pe, β) is given by in which ⊕ ⊕ and are given by (24), whereas In (31), we have the desired eigenvalue problem of size P + 1, which thus supports this same number of eigenvalues λ j . These are related to the (normalised) numerical frequencies j via Typically, one of the eigenvalues represents the so-called primary eigenmode, while the remaining ones can be regarded as secondary as they simply replicate the behaviour of the primary mode on shifted wavenumber ranges. This formally allows us to focus on the analysis of the primary eigenmode and on its dispersion and diffusion eigencurves. The reader is referred to [2,5] for the concepts relevant to the separation of primary and secondary modes adopted in this work.
Once the primary mode is identified, the scheme's numerical diffusion behaviour can be assessed in wavenumber space through the imaginary part of * h, where the asterisk subscript denotes the primary mode from (34). Note that numerical diffusion is especially relevant to turbulence computations as it impacts not only accuracy, but also stability. Note that eigencurves are entirely defined by the polynomial order P , the upwinding parameter β and, in case viscosity is present, the normalised Péclet number Pe = |a|h/μ, withh = h/(P + 1). Standard upwinding is here assumed. Figure 2 depicts a comparison between HDG's primary dissipation curves for pure advection and for advection-diffusion at Pe = 100 for P = 1, 4 and 7. As explained further below, this is about the lowest value of Pe one achieves (domainwise) in a turbulent flow computation. However, at this Pe , viscous effects are still somewhat weak in regular (linear-scale) plots of ih vs. κh, where i is the absolute value of 's imaginary part. This is especially true for P ≤ 4. Hence, Fig. 2 also shows these plots in log-log scale, highlighting what happens at wellresolved wavenumbers.
The log-log plots in Fig. 2 are revealing. They make clear that HDG's numerical diffusion follows the correct diffusive behaviour up to a certain wavenumber, hereinafter named κ c , beyond which upwind dissipation overcomes viscous diffusion. The exact diffusive behaviour, as derived from our model problem, is given by ih = (κh) 2 / Pe or log 10 ( ih ) = 2 log 10 (κh) − log 10 (Pe ) , showing that, as Pe increases, the reference line of exact diffusive behaviour shifts downwards, reducing the value of κ ch . Also, for a given number of DOFs, i.e. fixed h, increasing the discretisation order increases κ c . This type of analysis reveals how upwind dissipation and viscous diffusion complement each other, allowing also for the estimation of the wavenumber κ c after which upwinding dominates. The latter, though important for small-scale regularisation and stability, is not entirely physical in the sense of subgrid-scale modelling. Hence, k c values could be used as quality criteria for implicit LES / under-resolved DNS approaches based on discontinuous SEM. For transitional flows, where small numerical dissipation is particularly important, this kind of analysis might prove very useful. Although specific estimates would be needed for different schemes, the analysis strategy should be similar. Finally, it is now explained why Pe = 100 is about the lowest Péclet value one may find in a turbulent flow simulation. As candidates for very small Pe , one could think of the near-wall region of turbulent boundary layers, given the low velocity and small mesh spacing in typical wall-resolved LES. For the viscous sublayer, where u + < 5, the streamwise Peclét number can be evaluated using wall quantities: where by definition ν = u τ δ ν , being u τ the friction velocity and δ ν the associated viscous lengthscale. Our argument is then concluded since 50 < x + =h + < 150 in typical wall-resolved LES or under-resolved DNS approaches, cf. e.g. [10].

Concluding Remarks
We presented a preliminary study of the numerical dispersion and diffusion characteristics of HDG methods for linear advection-diffusion problems using the temporal eigenanalysis technique. To the authors' knowledge, this is the first eigenanalysis of HDG methods, and also one of the first of such analyses of a discontinuous SEM to consider viscous diffusion effects, cf. also [11]. It was shown that, for the range of Péclet numbers encountered in underresolved turbulence simulations, upwind (numerical) dissipation dominates viscous (physical) diffusion in the smallest resolved scales. Only in the large scales, the effect of viscous diffusion becomes significant. The wavenumber beyond which upwind dissipation overcomes viscous diffusion, and its dependence on the polynomial order, can be estimated through eigenanalysis, and this can be used as quality criterion for LES and DNS in general, and for implicit LES/under-resolved DNS in particular.
Future work includes further analysing the interplay between viscous and upwind diffusion, investigating other numerical fluxes (e.g. over-upwinding β 1, nearly central fluxes β ≈ 0, non-zero viscous stabilization σ = 0), and testing eigenanalysis against actual turbulence simulations. Finally, the dispersion-diffusion characteristics of HDG methods for spatially developing simulations could be investigated using spatial eigenanalysis techniques.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.