Numerical Dissipation Rate Analysis of Finite-Volume and Continuous-Galerkin Methods for LES of Combustor Flow-Field

A growing body of literature indicates that element-based high-order methods can exhibit considerable accuracy/cost benefit over conventional second-order finite-volume (FV) methods for large-eddy simulations (LES). This may even hold true for complex configurations relevant to industry that involve under-resolving unstructured grids. However, it is not often clear whether the accuracy/cost benefit stems from the low-dissipative nature of the high-order numerical schemes or from using a different LES approach (implicit/explicit), or a combination of the two. The present paper employs a numerical dissipation rate analysis technique due to Schranner et al. (Comput Fluids 114:84–97, 2015) to better understand the reasons for the high-order benefit seen previously on a complex LES test case related to gas-turbine combustors. It is established that a high(fifth)-order LES run provides better accuracy than its second-order FV counterpart at the same computational cost primarily because of lower numerical dissipation and the LES model dissipation has a secondary role to play. The numerical dissipation is found to contribute 60–90% of the total (numerical and LES model) dissipation.


Introduction
Scale-resolving simulations, and large-eddy simulation (LES) in particular, are becoming increasingly essential in the industry as well as academia for simulating highly unsteady eddy-dominated turbulent flows. This type of flow occurs in various scenarios, such as high-lift wing configurations, aircraft landing gear, and aero-engine gas turbine combustors. The present study is concerned with analysing two different methods for performing LES in a combustor flow scenario. The first method is the well-established second-order accurate finite-volume (FV) method, and the second is a newer high-order continuous-Galerkin (CG) method. More specifically, the former has been widely used for combustion applications, whilst the use of the latter is rare. One of the challenges for LES is to handle complex geometries while introducing minimal numerical error. High numerical accuracy is desired in a LES solver to correctly represent the propagation and interaction of the turbulent eddies over long distances and time. The current mainstream computational fluid dynamics (CFD) packages usually handle complex geometries using unstructured hybrid meshes (mix of hexahedra, tetrahedra, pyramids and prisms) and spatially low-order (second-order or less) accurate FV technique. Some examples of such packages are ANSYS Fluent, Siemens Simcenter StarCCM+, Rolls-Royce PRECISE-UNS and OpenFoam. However, most of these mainstream CFD packages incur significant discretization error due to low-order accuracy, mainly in the form of dissipation (Wang et al. 2013;Vermeire et al. 2017). Thus, solvers based on loworder schemes may not provide the required accuracy for a given cost, or conversely, might be too expensive for the required accuracy. The high-order (third-order or more) accurate schemes are emerging as an alternative to providing the desired cost and/or accuracy benefits (Wang et al. 2013;Slotnick et al. 2014) due to their significantly lower dissipation (Karniadakis and Sherwin 2005) and requirement of lower points-per-wavelength 1 3 for a given accuracy (Jameson 2001;Moratilla-Vega et al. 2020). Several variants of highorder methods exist, but an approach that provides geometrical flexibility and a compact numerical stencil is based on finite-element (FE) type polynomial approximations. The FE approach has resulted in methods such as spectral-element (Patera 1984), (dis)continous-Galerkin (Karniadakis and Sherwin 2005), and flux-reconstruction (FR) (Huynh 2007). Due to their compact nature, these methods can potentially better leverage modern multicore computers because of the enhanced compute-to-communicate ratio (Karniadakis and Sherwin 2005;Wang 2014). Note that the high-order FV method also provides flexibility but uses extended stencil, which results in a lower compute-to-communicate ratio.
The high-order schemes have been greatly successful on simple geometries and highly resolved structured meshes (pertinent to Direct Numerical Simulations or DNS); however, their benefit on complex geometries and under-resolved unstructured hybrid meshes, pertinent to industrial LES, is still a topic of research. The lower quality unstructured hybrid meshes often found in complex geometry cases usually challenge the stability of high-order solvers. In addition, the under-resolution results in non-smooth solutions, and as one moves away from the asymptotic regime (resolution → 0) the high-order schemes may lose their rate of convergence (or simply the accuracy) (Karniadakis and Sherwin 2005). Therefore, it has not been entirely clear if the high-order solvers can provide the promised cost/accuracy benefit while being robust in practical conditions. This requires objectively comparing the traditional low-order FV solvers to upcoming high-order solvers with improved stability techniques.
The international high-order CFD workshops (Wang et al. 2013;Cenaero 2017) have organised studies on accuracy vs cost for high-order schemes on various cases. It has provided a common platform for the community to cross-compare, however, conclusive results on unsteady flows are found only on geometrically simple cases such as Taylor-Green (TG) vortex. Other studies have been performed in a similar spirit, but the high-order benefit has been unclear. For example, two separate studies on plane channel flow (Sprague et al. 2010) and oscillatory plane jet flow (Capuano et al. 2019) compared the second-order OpenFoam solver with a high-order spectral-element solver Nek5000 and showed that comparable accuracy can be obtained by the high-order solver using 7-8 times less solution points. However, the cost incurred was less clear. Also, the meshes were simplistic (hexahedral only) and highly resolving. Vermeire et al. (2017) used unstructured, but still highly resolving, meshes for the comparison of the high-order FR package PyFR to a commercial FV second-order solver on two slightly more complex compressible flow configurations: cylinder-in-cross-flow at Reynolds number Re = 3900 and SD7003 aerofoil at Re = 60, 000 . In both cases, fifth-order accurate PyFR provided a better agreement to the reference but was almost 2 and 6 times, respectively, more expensive than the second-order solver. These studies took a first step to indicate the benefit of high-order solvers. However, neither better accuracy for a given cost nor lower cost for a given accuracy is demonstrated. Furthermore, these studies either involve simple geometries, structured meshes, well-resolved meshes or a combination of all of these.
The above gap is addressed by some recent studies that utilise more complex configurations and unstructured hybrid under-resolving meshes. Jia et al. (2019) compared a high-order FR solver hpMusic and a commercial FV second-order solver for compressible flow over a high-pressure turbine vane at Re = 1.16 × 10 6 . It was shown that a third-order FR simulation costing around one-third of a second-order FV simulation with a 64 times larger mesh provided more accurate results. This configuration, however, consisted largely of relatively thin wall-attached shear-layers and a small wake region. These characteristics do not quite represent the main features of combustor flow, 1 3 e.g., separated free shear-layers and strong re-circulations. A combustor-relevant study was performed by Saini et al. (2022), Saini (2022), where an isothermal configuration mimicking combustor port-flows (as depicted in Fig. 2) was used. A more fundamental TG vortex case, was also included to form a baseline of the study. The cost vs accuracy of second-order OpenFoam runs was compared to fifth-order CG runs of Nek-tar++ package (Cantwell et al. 2015), with experimental data (Hollis 2004) used as reference. The high-order solver was shown to resolve a broader range of fine scales and showed improved unsteady and mean statistics for the same cost. Further, the high-order solver was estimated as 3-8 times cheaper for a fixed accuracy. Even though the results favoured the high-order solver, the complete reasoning behind is unclear. An explicit subgrid scale (SGS) model Wall-adaptive local eddy-viscosity (WALE) (Nicoud and Ducros 1999) was used within the OpenFoam solver, whereas an artificial viscosity approach known as spectral vanishing viscosity (SVV) (Tadmor 1989;Karamanos and Karniadakis 2000) was used in the Nektar++ solver. These approaches were found suitable to the corresponding solvers (see Sect. II.A in Saini et al. 2022). Since different SGS models are used in the two solvers, it is not clear if the accuracy is more affected by the SGS models' dissipation or the inherent numerical dissipation of the solvers.
In the present paper, a technique by Schranner et al. (2015) is extended by including the SGS terms and employed to estimate the dissipation emanating individually from the numerics and the SGS models. This is applicable since in both the SGS approaches (WALE and SVV) model the effect of subgrid scales is purely dissipative (Nicoud and Ducros 1999;Tadmor 1989;Karamanos and Karniadakis 2000). Further, it is indicated in Schranner et al. (2015) that the numerical error generated in a sufficiently large subdomain is dissipative. In general LES practice, it is usually assumed that the SGS dissipation is dominant at subgrid scales. However, there are some studies (e.g., Ghosal 1996; indicating that numerics may have an important role to play as well, especially with low-order solvers. The current paper determines which source of dissipation is dominant in the low-and high-order solvers considered in Saini et al. (2022) and how the accuracy differences observed therein can be explained.
The dissipation estimation technique employed here is a post-processing method that allows the calculation of the numerical dissipation rate in arbitrary subdomains within the full CFD domain. It uses only the solution field provided by the solver under consideration. The technique is inspired by Domaradzki et al. (2003), who used a spectral method as a reference, but this made it difficult to generalise their method to non-periodic domains. The current method (from Schranner et al. 2015) works in the physical space. The procedure has been previously tested on simple and complex flow configurations such as the TG vortex , NACA-0012 airfoil , laminar separation bubble and flow past a sphere (Cadieux et al. 2017). In the present work, the procedure is extended to include the SGS terms, i.e., the WALE and SVV terms used in the FV and CG solvers respectively.
The objective of the paper is to quantify the contribution of the model dissipation (SGS or SVV) and numerical dissipation with respect to the resolved viscous dissipation in the low-and high-order combustor port-flow LES simulations of Saini et al. (2022). The TG vortex runs previously performed Saini 2022) are also considered for validation and to provide a baseline. To the authors' knowledge, this is the first time that the dissipation estimation technique of Schranner et al. (2015) is applied to compare two completely different discretization schemes (low-order FV and high-order CG) and corresponding LES solvers.

3
The details of the numerical dissipation rate calculation are covered in Sect. 2. A minimal description of LES previously undertaken, along with some results, is provided in Sect. 3. The numerical dissipation calculations are validated in Sect. 4 and the analysis results are presented in Sect. 5. The paper closes with a summary and conclusion in Sect. 6.

Method
The representative incompressible flow solver for the second-order FV technique was pimpleFoam from the OpenFoam package, whereas the high-order solver was Inc-NavierStokesSolver based on the spectral-hp Nektar++ framework (Cantwell et al. 2015). The following Subsections provide descriptions of the solvers' numerical methods, dissipation rate equations and their computations.

Numerical Methods in Solvers
The main features of the numerical methods employed in the two solvers, OpenFoam and Nektar++ are briefly mentioned in Table 1. For a more detailed description the reader is referred to Sect. II.B in Saini et al. (2022).
The contrasting difference between the two solvers is the order of accuracy, which is enabled by the corresponding spatial discretization. Both solvers use unstructured hybrid element meshes and can easily handle complex geometries. Another notable difference is the SGS models. The time discretization is similar, i.e., second-order accurate based on backward scheme. Note that given the small time step size used in LES, the error due to time discretization is considered small. The discretization of the convective term is important for the resulting numerical dissipation in high Re cases. For the second-order FV solver, the TG vortex LES results analyzed here were run with the second-order central scheme. However, the combustor case runs required a 60%-40% blending of second-order central scheme and second-order upwind scheme for stability. This blending ratio is constant across the whole domain. For the viscous term, central scheme was used in all runs. The high-order solver employed Jacobi polynomials over each element to represent the solution. The CG method is used to describe solution  (Nicoud and Ducros 1999) SVV (Karamanos and Karniadakis 2000) over element boundaries. The calculations are run at polynomial order 4 (P4) for velocity and one order less for pressure, to satisfy the inf-sup condition. The convection term uses an upwind scheme. For detailed discussions on discretization aspects, see Sect. II.B in Saini et al. (2022).

Dissipation Rate Equations
Gas turbine combustors operate at low speed and are usually modelled as an incompressible flow (McGuirk 2014). Starting from the incompressible flow equations, the derivation of the kinetic energy equations used to calculate the dissipation rate is shown in the following.
In the context of LES, the filtered version of the Navier-Stokes (NS) equations can be written as is the filtered velocity, (x, t) the filtered pressure field with a constant mass density , the kinematic viscosity of the fluid (assumed constant). The SGS stress tensor is denoted by s ij and the filtered strain rate tensor using S ij . No special symbol is used for filtered quantities for the sake of simplicity.
The can be absorbed in the pressure to give another field, say . Multiplying the momentum Eq. (1) with u i , recognising that filtered kinetic energy k = 1 2 u i u i , and using the continuity statement one obtains The second and third terms on the rhs of (2) can be recast as Inserting these into Eq. (2) and grouping the volume and flux terms, one obtains (3) where the repeated index i in u i ∕ x i is replaced by j for consistency. Integration of the kinetic energy transport Eq. (4) over an arbitrary subdomain of volume V leads to where Equation (5) applies to the second-order OpenFoam solver that uses the explicit SGS model. Note that similar balance equations have been obtained in Sagaut (2006) and Cadieux et al. (2017) for a second-order solver. For the high-order Nektar++ simulations, SVV is used instead of an explicit SGS model, so E sgs and F sgs in Eq. (5) are not applicable. The contribution of SVV term to the kinetic energy equation is now derived.
The SVV term can be written as svv x j Q ⋆ u i x j (Karamanos and Karniadakis 2000), where svv is the maximum value of the artificial viscosity, Q is a viscosity operator and ⋆ denotes convolution. (Further details on svv and Q can be found in Sect. 3.) Multiplying this SVV term by u i and denoting Q ⋆ The first term on the rhs resembles the integrand of the SGS flux term F sgs in the sense that there is contribution from the gradient of the product of u i , dV is viscous flux (diffusion by viscous effects), and Note that the svv is directly proportional to the element size in this paper (see Sect. 3), so svv may vary in space for spatially changing element sizes inside a subdomain. It means that svv ∕ x j might be non-zero. However, relatively speaking u i Q u i x j svv ∕ x j would be smaller than svv Q u i x j u i ∕ x j for a large velocity gradient and a slowly varying mesh size. Therefore, one could expect that the second term on the rhs of (7), or H svv in (8) would be negligible. As the subdomains chosen in this study consist of highly turbulent flow (large velocity gradient) and similar mesh size, this is found to be the case for each subdomain studied. Thus, H svv is omitted in the plots for clarity.
Considering the integral kinetic energy balance Eqs. (5) and (8) in their time and space discretized form, the residual on the rhs ceases to be zero because of the discretization errors. The residual is instead a finite value, say −E num , representing the numerical error over the considered subdomain. For the OpenFoam solver, this can be written as where subscript Δ represents the discretized form. Similarly, the Nektar++ formulation becomes The difference in Eqs. 9 and 10 lies in the SGS and SVV terms.
Both dissipative and dispersive types of numerical errors are included in the residual E num . However, for a large enough volume V of a subdomain, the value of E num has been previously seen to be positive Cadieux et al. 2017) and therefore can act as a measure for numerical dissipation rate. For inter-solver comparisons, the ratio of the dissipation from subgrid scale term and numerical term, i.e., E sgs ∕E or E svv ∕E is a useful quantity. It is equivalent to the ratio of corresponding viscosities as follows

Computations
The calculations of the kinetic energy balance terms in Eqs. (9) and (10) must be performed using discretization schemes of the same or higher accuracy as used within the LES solver for consistency . For example, the space derivatives need is the SVV dissipation, x j dV is attributed to spatial variation of svv .
to be calculated with at least second-order accuracy for OpenFoam and fifth-order accuracy for Nektar++ solutions. The subdomain(s) over which to integrate can be arbitrarily chosen. The subdomains can overlap, and a subdomain can even be a single cell/element. However, there are a few points to bear in mind: 1. A subdomain should be large enough for the numerical error to be dissipative in nature. This is because the point-wise values are very noisy, so averaging over a sufficiently large domain is necessary to recover reasonable approximations of the dissipation rate (Sect. 6.2 in Schranner et al. 2015). In the present work, the subdomains are of the order of the size of important geometrical features in the flow configuration. This is similar to other works; for example, in Cadieux et al. (2017) the subdomain size is of the order of sphere diameter in a direct numerical simulation (DNS) of flow around a sphere ( ∼ 60k-500k cells) and in  the size is of the order of maximum aerofoil thickness ( ∼ 20k cells). 2. The location of the subdomain should capture the flow regions of interest Cadieux et al. 2017). One caution is pointed out in  that those regions should be avoided where the flow is wellresolved (and so the numerical error is expectedly low), such as laminar flow regions. This is because the velocity gradients may be small, resulting in erroneous residual estimates. It is noted in  that in the TG vortex case, even though the initial condition is smooth, the gradients are still high enough for sufficient dissipation.
In the calculation process, the time derivative term k∕ t is calculated using a central difference formula. For the nth time snapshot where the kinetic energy k is calculated from the solutions fields at time t n+1 and t n−1 . This scheme is second-order accurate as in the solvers and is also used in Schranner et al. (2015) and Cadieux et al. (2017). The time-step between the snapshots Δt is the same as the simulation's time-step. Hence, all time-local relevant phenomena, e.g., passing of large vortical structures, are captured with high accuracy. The volume integrals can be calculated either by direct volume integration in the discretized space or via calculating the flux integrals ( F k , F ac etc.) by integrating over the surface encapsulating the subdomain (Gauss theorem). It has been shown in Schranner et al. (2015) that both approaches (volume and surface integrals) yield very close results. In this work, all the integrals are calculated using the volume integration approach.
The spatial derivatives are calculated using the numerics available in the respective solvers. For the OpenFoam solver all the terms in Eq. (4) are calculated on-the-fly in the entire domain. The solver is modified to include extra calculations following Cadieux et al. (2017). The SGS term calculation is added in the current work. A snippet can be found in Sect. A.2 of Saini (2022). The full domain snapshots containing all the calculated terms are saved to the disk. The integration over arbitrary subdomains is performed as a postprocessing step in Paraview software. As the representation of the data in Paraview is the same as the solver, this process introduces negligible error in the calculation of balance terms (9) for the OpenFoam solver.

3
The essence of the process is similar for the Nektar++ solver. The snapshots of the entire domain containing velocity and pressure fields are saved to the disk, and all terms are calculated thereafter. The calculations (including the spatial derivatives) are performed using Nektar++ package's post-processing utility FieldConvert. This ensures that the derivatives are of the same order of accuracy as the simulation itself. However, the convolution part Q ⋆ u i x j involved in the SVV term was not straightforward to calculate using the pre-existing utilities. So, the solution field was interpolated onto a uniform grid, and the convolution calculated via the Fourier transform. The details are as follows.
The velocity field in each subdomain is interpolated onto a three-times finer directionally-uniform Cartesian grid to sufficiently capture the polynomial variations (Appendix B in de Wiart et al. 2014). The interpolation is performed using the utility FieldConvert such that the order of interpolation is the same as that of the solution. Then, using a Python script, the convolution is calculated as a product of the Fourier coefficients of the velocity field û and the SVV kernel Q . Mathematically, where k is the wavenumber (Moura et al. 2016). The velocity fields may not be periodic, so they are pre-multiplied by a window function (Tukey window Harris 1978 with = 0.25 ) prior to applying the Fourier transform. The windowing reduces the 'leakage' of the signal energy to high-frequencies. The current method to calculate the convolution term is not exact near the boundaries. However, it sufficiently captures the energy as seen in the validation exercise (Sect. 4). The element size h is required for svv and is calculated as the cube root of the element volume. This method to calculate h is sufficiently valid for nearly isotropic meshes. Considering this constraint, the subdomain (SDs) are kept away from large aspect ratio elements (e.g. near walls). Overall, the procedure for calculation of balance terms for Nektar++ solver (10) is less general and non-exact but care has been taken to provide sufficiently accurate results.

Large-Eddy Simulations
The LES runs on the two configurations have been performed and reported previously Saini 2022). The detailed mean velocity and turbulent kinetic energy along with unsteady features such as flow scales, spectra and probability density function have been presented therein. Thus only a brief overview of these LES is provided here. The cost is defined as the total core hours for the whole simulation to finish. This is applicable because a single type of machine was employed for all the simulations. In the following, the Nektar++ and OpenFoam solvers are frequently referred to as Npp and OF, respectively.
The OF solver used WALE model (Nicoud and Ducros 1999) with model constant as C w = 0.325 . For the reasons of choice of the WALE model and the constant value, see Saini et al. (2022). The Npp solver used a recently formulated SVV operator with a power-kernel (Moura et al. 2016). In this operator, the SVV magnitude is proportional to a measure of the local advection speed U l and local mesh spacing h∕P , which gives where 0 is a fixed parameter, h is the local elemental mesh size and P is the polynomial order. Further, the spectral variation ( ⋅ ) of the kernel is defined as where m is the mode index and P SVV relates to the activation threshold (a low P SVV value generates a stronger damping of the higher modes). For complete details, see Sect. II.A in Saini et al. (2022).

TG Vortex
In this case, the flow transitions from well-defined laminar conditions to a fully developed turbulent regime in a triply periodic cubic domain. The flow is well-resolved to start with, and sub-grid scales begin to appear as time progresses. The coarser the mesh, the earlier in time the sub-grid scales appear. Due to simple geometry and initial conditions, this configuration is a commonly used test of numerical schemes for solving unsteady NS Eqs. (Wang et al. 2013;Vermeire et al. 2017). Some specifics are provided briefly before presenting the results. The computational domain extents are − L ≤ x i ≤ L . The Reynolds number is Re = 1600 , based on the initial velocity scale V 0 and length scale L of the vortices. The exact initial conditions are omitted here but can be found in Saini et al. (2022). The non-dimensional time is defined as t * = tV 0 ∕L , and the simulations were run for t * = 10.
Several mesh-resolution runs were performed in the reference LES studies Saini 2022). Two second-order OpenFoam and two fifth-order Nektar++ runs are considered here for the dissipation rate analysis. These represent a wide range of cost vs accuracy, as explained later. The OpenFoam runs consist of a uniformly distributed mesh with 64 and 128 points per direction, and are labelled as 64 × 1 OF and 128 × 1 OF respectively. In these labels, 1 refers to the polynomial order P. The Nektar++ runs have 16 and 8 mesh points and are labelled as 16 × 4 Npp and 8 × 4 Npp, 4 being the polynomial order. The cost of the baseline run 64 × 1 OF was 285 s on 56 cores .
A quantity of interest is the integral kinetic energy dissipation rate ( D t ) (14) SVV = 0 U l h∕P, The temporal evolution of D t is depicted in Fig. 1 for all four runs along with a reference result obtained from a spectral solver (van Rees et al. 2011). The 16 × 4 Npp run is just 1.3 times more computationally expensive than the baseline 64 × 1 OF, but it is about ten times more accurate (in the L 2 sense, see Saini et al. 2022). The 128 × 1 OF provides better accuracy, but is eight times more costly than the baseline 64 × 1 OF. The 8 × 4 Npp provides similar accuracy as 128 × 1 OF until t * = 5 (thereafter they diverge from reference at different times), for a fraction of the cost, i.e., 0.15 times the baseline cost. Overall, the results favour the high-order solver, but what is less clear is whether numerical dissipation or SGS model dissipation is affecting the accuracy more. This is where the dissipation quantification from different sources becomes interesting.
Another relevant quantity is the Enstrophy based dissipation rate ( D ): where Ω is the flow domain volume, = ∕ 0 is the kinematic viscosity and is the enstrophy, defined as where is the vorticity. Note that the expression (18) is valid only for incompressible flows (DeBonis 2013). In the case of incompressible flow simulations, the kinetic energy dissipation mainly comes from the viscous action, the numerics (discretization error) and the SGS model, if used. The dissipation rate calculated using the enstrophy integral (18), referred to as D , relates to the viscous/turbulent contribution (Bull and Jameson 2015). On the other hand, the dissipation rate D t calculated using (16) gives the overall (or total) dissipation. The difference D t − D v (= D n ) therefore represents the dissipation caused collectively by the numerics and the SGS treatment (including any filtering) (DeBonis 2013; Bull and Jameson 2015). This means that for a well-resolved or DNS simulation, D t = D v . However, in under-resolved scenarios the difference D n grows as sub-grid scales appear in time. This observation will be useful in the later Sections. (17)

Combustor Port Flow
The port-flow case represents the dilution jets employed in rich-burn type gas turbine combustors. These cross-flow jets provide the oxidant, support the recirculation zone, enhance the turbulent mixing and thus, improve overall combustion efficiency (McGuirk 2014). The present configuration is an isothermal flow in a representative geometry, designed and studied experimentally by Spencer (1998), Spencer and McGuirk (2001), and later by Hollis (2004) using particle image velocimetry (PIV). The latest PIV data (Hollis 2004) was used as a reference for accuracy in Saini et al. (2022). The geometry consists of two concentric pipes. The inner 'core' pipe is fed from the outer 'annulus' through six equally spaced circular ports, thus creating six jets impinging onto each other (see Fig. 2 for a model). Three parameters define the flow condition: (1) ratio of mass flow rate through ports and annulus, or bleed ratio, B, (2) ratio of bulk velocity through the ports and the core inlet, , and (3)  The modelled geometry is shown in Fig. 2. The origin is at the intersection of the axes of the core and the ports. Three runs, conducted previously , are under consideration. These are listed in Table 2. The first run is a baseline 16 million cells secondorder OpenFoam run, whereas the second is a 130,000 element polynomial order 4 (total 8.4 million solution points) Nektar++ run. The mesh of the second run was designed to incur an almost identical computational cost as the first baseline run. The third run is again a second-order OpenFoam run but with higher mesh resolution (48 million cells), targeted to match the turbulence spectra of the high-order Nektar++ run. This third OF run costs around 5 times the previous two runs. The cost of the baseline 16 million OF run was 178 hours on 336 cores ). The mesh topology, element-type and size distribution are equivalent in the three runs. For a first idea of the mesh design, see Fig. 9, which displays the elemental mesh of the Nektar++ run. The reader is referred to Saini et al. (2022), Saini (2022) for the detailed rationale behind the model geometry, the mesh design and time step size.
There are particular flow features worth revisiting. Due to the strong interaction of the jets, a prominent central recirculation zone (CRZ) (a 'bubble' with negative x-velocity) is formed, as labelled in Fig. 3a. This CRZ is observed to grow and shrink periodically with a low frequency (Fig. 9a in Saini et al. 2022). Further, a high level of turbulence is present in the domain that peaks in the vicinity of the impingement point ( x ≈ 15 mm, r = 0 ). See Fig. 4a. Note that the overbar u denotes the ensemble mean and u ′ u ′ the variance. In Saini et al. (2022), the high-order results were shown to be more accurate for a given cost. For instance, it may be noticed from x-velocity contours in Fig. 3 that the upstream reach of the CRZ is better captured by the high-order 8.4M-Npp than the second-order 16M-OF. The CRZ is reproduced slightly more accurately by the 48M-OF but at a five times higher cost. Similarly the low-frequency movement of the CRZ was better captured by the high-order run for a given cost (not shown here). This was attributed to excessive dissipation in 16M-OF that, if estimated correctly, would have allowed the fluctuations to reach further in the upstream direction. Further, a comparison of the resolved mean turbulence kinetic energy (k) was performed in Saini et al. (2022). An estimate of the turbulence kinetic energy, as available from the two-component PIV data (Hollis 2004), was used. The contours of this quantity are reproduced in Fig. 4. The three runs exhibit a similar overall distribution, however the second-order OF runs significantly under-predict the peak value in the impingement region ( x ≈ 15 mm, r = 0 ). This under-prediction can be directly attributed to the larger dissipation rate in the second-order runs. The picture becomes clearer when the breakdown of the dissipation is presented later in Sect. 5.2.

Validation
The numerical error calculation procedure outlined in Subsect. 2.3 is validated and applied to the selected TG vortex runs in this Section. Following Schranner et al. (2015) and Castiglioni and Domaradzki (2015), a test case with Re = 3000 is chosen for validation. The mesh resolution is kept at 64 solution points per direction, i.e., 64 cells for OpenFoam and 16 elements with polynomial order 4 for Nektar++ per direction. Depending on the extents of the SD considered, the following three validation tests are used-1. Full domain: Given the periodic nature of the TG vortex domain, the flux terms vanish in the volume integral sense and the contribution to E num comes only from the timederivative ( K t ) and the dissipation terms . So, the terms −K t and E calculated with the presented method (Sect. 2.3) must match the previously (independently) calculated D t and D terms respectively. Further, as explained in Sect. 3.1, the difference D t − D should be equal to the total dissipation rate, i.e., numerical plus the SGS (E num + E sgs ). 2. Octant subdomain: Due to flow symmetries, the flow in the full domain can be represented by any one octant, say 0 ≤ x, y, z ≤ L . This means that the sum of the fluxes over each octant's boundary cancels and provides the same average dissipation rates as the full domain. 3. Arbitrary internal subdomain: The different terms are calculated over a particular cubical internal subdomain and compared to Schranner et al. (2015) and Castiglioni and  . The results are not expected to match exactly because of the different solvers used (especially after the time range when the subgrid scales begin to appear). Nevertheless, similar trends would increase confidence in the current implementation.
The various terms of Eqs. (9) and (10) integrated over the full TG vortex domain are plotted in Fig. 5. The markers in the plots represent the individual instantaneous fields. First, the fluxes (denoted by F terms) are individually close to zero. Further, for both solvers, there is an excellent match of −K t and E to D t and D respectively. This verifies the volume integration and time derivative computations. In addition, the difference term D t − D matches the SGS model plus the numerical dissipation rate (square marker) very well. An OpenFoam run that does not use the SGS model is also plotted (Fig. 5b) to add confidence in the implementation. In this case E sgs = 0 is confirmed. Finally, due to wellresolvedness during the initial period of the simulations, the dissipation from the numerics, E num , as well as the model, E sgs (or E svv ), is nearly zero. As seen from the Nektar++ run (Fig. 5c), the accuracy of the SVV terms calculated using the Fast Fourier transform (FFT) convolution procedure is also satisfactory. Similar observations are made for the octant subdomain. For brevity, the results of the octant SD are moved to Appendix 1 (Fig. 17).
As the third and final test, an arbitrary cubical box consisting of cells from 10 to 56, i.e., is chosen following Schranner et al. (2015) and . The reference data is taken from  where a commercial second-order FV solver is used without any SGS model. For a fair comparison, the present OpenFoam solver is also run without the SGS model. The results are shown in Fig. 6a. As expected, the match is excellent in the initial well-resolved time range ( t * < 3 ) but deteriorates later. Anyhow, the trends are reproduced considerably well. The Nektar++ simulation cannot exactly represent this subdomain due to the discretization over 16 elements per direction. Thus, a slightly larger domain of 2 8 64 L − L ≤ x, y, z ≤ 2 56 64 L − L had to be chosen. This results in a slight mismatch, but the trends agree with the reference values, as plotted in Fig. 6b.
The results here show that the calculation procedure for the terms in the kinetic energy balance Eqs. (9) and (10) functions satisfactorily. The dissipation rate due to the numerical and subgrid scale sources for the low-and high-order solvers can now be compared.

Dissipation Rate Comparison
The dissipation rate emanating from the viscous action, SGS models and the numerics are compared for the TG vortex and port-flow LES runs mentioned previously in Sect. 3. The related differences between the low-and high-order solvers are highlighted.

TG Vortex
The analysis here is performed over the full TG vortex domain so that the flux terms simplify to zero. Recall that simulations with Re = 1600 are considered. Figure 7 plots the non-zero terms of the kinetic energy balance integrated over the full domain from two meshes of the two solvers. The normalised cost of the runs is included. The sum (E model + E num ) , called from here on as total non-viscous dissipation (TND), is also plotted. Note that E sgs and E svv are collectively referred to as E model . For 64 3 OF simulation (Fig. 7a), a prominent local peak of the TND develops early at t * ≈ 5 with a dominant contribution ( ∼75%) from E num . This peak represents the excessive dissipation due to which the D t peak shifts to an early t * compared to the reference, as seen previously in Fig. 1. The early peak of (E sgs + E num ) is lowered with a higher 128 3 resolution (and eight times cost), as depicted in Fig. 7b. This improves the resulting D t curve as noticed in Fig. 1. The numerical contribution in TND is still dominant (65-70% for t * > 5 ). As an aside, the numerical dissipation in the second-order solver is seen to 'adjust' to a higher value when the SGS model is absent (compare Fig. 5a and b).
Turning attention to the high-order runs, the 16 3 P4 run in Fig. 7c displays negligible TND (square symbol) until t * = 3 . Whereas, the second-order 64 3 run in Fig. 7a display non-zero TND even in the well-resolved period ( t * < 3 ). Further in time, the TND at t * = 5 is only a quarter of 64 3 OF run. Thus, the increment in TND is more controlled and the agreement of D t with the reference data ( Fig. 1) confirms that a correct amount of dissipation is imparted. The TND consists mainly of the numerical part E num (90-95%) and contribution from the model part E svv is very small. For the coarser mesh of 8 3 P4 in Fig. 7d, the TND is again tiny up until t * = 3 which is comparable to the four times finer 128 3 second-order run. However, the escalation thereafter in numerical dissipation is quite high, with the relative contribution from SVV dissipation still significantly lower (10-20%). The absolute values of E svv are even smaller than 64 3 secondorder E sgs . The TND peak values in high-order runs at later t * are similar to the secondorder run, but this is viewed as a consequence of the preservation of energy at early t * . In this non-stationary flow case, the flow history is important. The correct representation at early times leads to accurate results at later time range. This is a notable difference between this canonical case and the practical cases.

3
A quantification for the 'resolvedness' is provided by the ratio of numerical or model dissipation to the viscous dissipation as defined in Eq. (11). Figure 8 plots these quantities for the runs discussed above. The peak at t * = 5 in the second-order runs is more prominent here. The relative TND (E model + E num )∕E is approaching 4 in 64 3 , which indicates that the simulation actually runs at a much lower effective Re than intended. Further, it is more apparent here that the SGS dissipation is significant even in the initial 'resolved' period. These factors lead to excessively dissipative, undesirable results. The situation dramatically improves with the 16 3 P4 run for a similar cost with timely intervention of numerics and SVV. Here, the relative TND is at maximum 1.2 times the viscous dissipation. It is noted that for high-order runs (E num ∕E ) is slightly negative in the laminar flow regime, which is attributed to the division of two small quantities, as observed previously by .
In summary, there are two main findings. First, the numerics are revealed as the main contributor to the dissipation rate of kinetic energy (65-75% of TND in OF second-order runs and 85-95% in Npp high-order runs). In other words, numerical dissipation is the dynamics controlling term. This is important as it is usually assumed that the subgrid model provides dissipation at the smallest scales. Second, the TND in high-order runs is significantly lower (1/4th) than second-order runs for a similar cost. As the numerical dissipation E num (identified as controlling term) from the high-order solver exercises lower

Combustor Port-Flow
The dissipation rate of kinetic energy in the low-and high-order runs, as listed in Table 2, is examined. The balance terms are calculated from Eqs. (9) and (10) over several SDs, where the SDs are chosen to represent the key flow regions.
The investigated SDs lie primarily on the central axis and cover the impingement region, CRZ and downstream. An additional SD covers the port-jet shear layer. These are visualised in Fig. 9 and details are documented in Table 3. As explained in Subsect. 2.3, the SDs are kept away from large aspect ratio elements. The extents y, z ∈ [−12, 12] mm keep the 'central-line' SDs inside the tetrahedral mesh block and provide a volume large enough to result in a smooth numerical dissipation ( E num ) time history. The analysis is performed over a physical time of 0.5 s. In non-dimensional terms, this is equivalent to 25D p ∕U p , where D p is port diameter and U p is port flow velocity. This sampling duration is restricted by the available disk size, as the flow fields need to be stored. As we are considering integral quantities over considerably large SDs, the sampling period for time-averaging is less restrictive than for single-point statistics. For example, standard error on E num varies  between 2.5-5.0% for all SDs except the jet stream SD. For jet SD, the standard error is about 20%, which is considerably large. For the Npp run, it is verified that the H svv term is consistently smaller than 5% of E svv for the subdomains considered. An exception is the jet stream SD where it reaches around 10% of E svv due to lower turbulence levels (and thus smaller velocity gradients). The small value of the H svv term is as anticipated in Sect. 2.2, and is therefore omitted in the following plots. All four subdomains are discussed in the following. Impinge: Figure 10 plots the balance terms integrated over the Impinge SD for the three LES runs. This SD is important as it contains the peak of turbulence kinetic energy. It also plays a key role in setting the fuel-air ratio in the CRZ, which is critical in . The solid circles of E num indicate the instants of the calculations. The typical time interval between data points is 0.01 s (0.0125 s for 8.4M-Npp, 0.01 s for 16M-OF, 0.0132 s for 48M-OF). As a reference, the port flow time scale is D p ∕U p = 0.02 s. Recall that the time-derivative K t is calculated using a local time-step that is the same as the simulation's time-step. In every run, the surface transport due to viscous effects ( F term) is negligibly small, as can be expected due to high Re. A similar observation is made by Cadieux et al. (2017). The subgrid model flux ( F svv or F sgs ) is also small. The dominant contributing terms are the time rate of change K t , advective flux F k , pressure work F ac and the dissipation terms E . The first three terms fluctuate more than the dissipation terms. The numerical dissipation E num is significantly larger than the viscous dissipation E and subgrid model dissipation E model (the E sgs and E svv are collectively referred to as E model ). However, the relative contribution from viscous, subgrid and numerical sources varies from case to case. It is useful to calculate the ratio of numerical or model dissipation to the resolved viscous dissipation as defined in Eq. (11).
The sum (E model + E num )∕E indicates the level of 'resolvedness' in each LES run; the higher this ratio, lower the resolvedness. Also, this ratio indicates the 'effective' Reynolds number at which the simulation is run. Figure 11 plots E model ∕E , E num ∕E and their sum for the Impinge SD. To begin with, the time variation is plotted for the three LES runs. The time averages are then collected in Subfigure d. The variance in time histories is small, so the time average has a low statistical error. Note that the fluctuations in E model ∕E (dash-dot lines) are much smaller than E num ∕E (solid lines) indicating that numerical dissipation is more sensitive to time-local flow than the model dissipation.
For further discussion, the sum of the model and numerical dissipation is called TND (as in the previous Subsection). As expected from the resolved k 2C contours in Fig. 4, the second-order 16M-OF has the highest level of TND. This is depicted by the time averages in Fig. 11d (third group of bars). This is almost 30% lower in the finer-mesh 48M-OF run with respect to the 16M-OF run. The 8.4M-Npp provides the lowest TND, which is 40% lower than the same-cost 16M-OF run. Notably, the model dissipation in the 16M-OF and the 8.4M-Npp are approximately the same (first group of bars), and the main difference in the TND comes from the numerical dissipation (second group of bars). This shows that the high-order solver resolves more kinetic energy mainly due to lower numerical dissipation and the model dissipation has a smaller role to play. This is further highlighted by the fact that even though the SVV model dissipates higher in 8.4M-Npp than 48M-OF (first group of bars), the overall dissipation is dominated by the numerics, which is smaller in 8.4M-Npp. Finally, note the relative contribution from the model and numerics within each run. The model contribution is higher in 8.4M-Npp at 20% as compared to OF runs at around 10%.
CRZ: The CRZ subdomain is important as this region represents a typical combustor's primary zone (McGuirk 2014). The balance terms over the CRZ subdomain are plotted in Fig. 12. Again, the numerical dissipation E num is significantly (an order of magnitude) larger than the viscous and subgrid model dissipation rates. Similar to the Impinge SD, the F term and the subgrid model flux ( F svv or F sgs ) are relatively small. The F k , F ac and K t terms are dominant and highly fluctuating. On the contrary, the dissipation rates are relatively smooth in this SD as well. Finally, numerical dissipation dominates the viscous and model ones.
The ratios of numerical and model dissipation to the resolved viscous dissipation are plotted in Fig. 13. The 16M-OF displays the highest TND (third group of bars in Fig. 13d) as before. However, both the Npp and 48M-OF runs stand equally, 25% lower. As in the previous SD the difference in the numerics (second group of bars) is the dominant factor compared to the model. Finally, the relative contribution from the model in 8.4M-Npp is around 15% and in OF runs around 8%. In previous experiments (Hollis 2004) and LES Saini 2022), the recirculation region has been shown to inflate and shrink periodically with a low frequency (see Fig. 9 in Saini et al. 2022). In this context, the current dissipation distribution is not expected to change significantly with time as the CRZ SD does not extend into the region of inflation ( x < −20 mm in Fig. 9). DownStr: Figure 14 plots the quantities for the DownStr SD. Here, a considerably higher viscous dissipation E is resolved by the Npp run as compared to the OF runs (four times of 16M-OF and three times of 48M-OF). As a response, the SVV dissipation is also higher in Npp. It indicates that comparatively more small scales are being resolved in this SD by Npp and thus more of viscous dissipation spectrum is resolved. This is visually  Fig. 7 of Saini et al. (2022). Perhaps, the hexahedral mesh present in this SD (see Fig. 9) results in low dissipation for Npp and hence the difference.
The normalised numerical and model dissipation clarify this further in Fig. 15. Between the OF runs, the normalised SGS and numerical dissipation drop (by 40% and 30% respectively) when the mesh is refined from 16M to 48M. This accumulates to 35% lower normalised TND in 48M run, at a five times higher cost. The normalised TND of Npp is much lower (one-fifth of 16M-OF). Surprisingly, the SVV and numerical contributions from Npp are even smaller than that of the expensive 48M-OF run. These quantitative results are reflected in the downstream region Q-criterion illustrations in Fig. 5.7 of Saini (2022), where only large-scale structures are seen in the 48M-OF. Finally, the relative contribution of SVV in Npp changed from 15-20% in earlier SDs to 44% here. This occurred in the OF runs as well, as SGS contribution changed from 10% in earlier SDs to 20% in DownStr SD. The exact reason for this peculiarity is unclear. Possible reasons could be different mesh-type or the level of 'resolvedness' or a combination of these. The latter refers to the fact that the current SD features much higher (E model + E num )∕E values than the previous SDs. For instance, 16M-OF is approximately 13 for Impinge SD, 11 for CRZ SD but 23 for DownStr SD.

JetStr:
The jet-stream region is important to predict any possible instabilities in combustion (see page 21 in Saini et al. 2022). The terms in the JetStr SD are plotted in Fig. 16.
Here the flux terms are considerably larger than the dissipative contributions (note that the dissipation terms are arbitrarily multiplied by 10 for the sake of readability). This is attributed to fast-moving large coherent vortical structures and negligible stochastic small-scale turbulent fluctuations, as visible in Fig. 9. In fact, this is the kind of situation where the current method might become unreliable, as explained previously in point 2 of Sect. 2.3. Curiously, higher numerical dissipation is observed in the finer 48M-OF case than in the 16M-OF, which points to the unreliability of the analysis. In addition, the variance in the numerical dissipation E num curve is significant (standard error around 20%). Therefore, the rest of the plots are omitted for this SD as strong conclusions cannot be drawn.
The main findings are summarised as follows. In the key regions or SDs with fully developed turbulence (Impinge, CRZ and DownStr), 1. numerical dissipation rate ( E num ) is generally much larger than that of the subgrid model ( E svv or E sgs ) and viscosity ( E ). Numerical dissipation forms the dominant factor for the differences observed between the three LES runs. This can be quantified in several ways: (a) E num is 4-12 times of E sgs for second-order runs, and 1.3-6 times for high-order runs. (b) Between numerical and model dissipation, the former is a dominant factor in determining differences in the (E model + E num )∕E among different LES runs. For instance, the E model ∕E is almost the same for both 16M-OF and 8.4M-Npp in the Impinge SD and all the difference comes from the E num ∕E . Comparing 16M-OF and 48M-OF, 15% of the difference in (E model + E num )∕E is due to the model and the rest 85% is due to numerics. This is also true for the CRZ and DownStr subdomains as well. (c) In terms of the percentage contribution of numerics and model in TND for each individual run, the numerics amount to 80-90% for the OF runs, and a somewhat lower 56-86% for the Npp run. The model contribution in OF may be lower because numerical dissipation in OF already suppresses the (square of) velocity gradients, and thus there is less to do for the SGS model.
2. the normalised TND (E model + E num )∕E or the 'effective viscosity' Eq. (11) varies in range 10-23 for 16M-OF, 8-15 for 48M-OF and 2.5-7 for 8.4M-Npp. The 48M-OF range is expectedly lower than 16M-OF due to mesh refinement. More interestingly, the 8.4M-Npp range is consistently lower than the same-cost 16M-OF. It is also lower as compared to the five times expensive 48M-OF except in the CRZ subdomain, where they are approximately equal. Thus, the high-order run here is at least capable of resolving as much energy as the five times more expensive second-order run. As an aside, for the TG vortex simulations the overall (E model + E num )∕E was lower (e.g., 1-4 for second-order runs), meaning that the resolvedness was higher. This highlights the fact that usually the industrial flow simulations are higher Re simulations than the canonical test cases. 3. the fluctuations in E model ∕E are much smaller than E num ∕E . This could indicate that numerical dissipation might be more responsive to instantaneous flow than the model dissipation.

Conclusion
A kinetic energy dissipation rate analysis of previously performed LES of a canonical configuration and a combustor port-flow case ) was performed. The LES had compared a standard second-order FV solver with a newer high-order CG solver and had demonstrated cost vs accuracy advantage of the latter. However, since different order-of-accuracy and SGS models were used in Saini et al. (2022), the present dissipation rate analysis aimed to clarify the role of the numerics and SGS model in the overall dissipation rate and thus overall accuracy. The various dissipation terms were calculated over physically important subdomains by extending the method of Schranner et al. (2015) to include SGS model terms.
It was found that the numerical dissipation (E num ) is a dominant term, and SGS model dissipation (E model ) is minor in both second-and high-order solvers. In the combustor case, E num was as high as 90% of the total non-viscous dissipation (E num + E model ) for the secondorder solver, and 85% for the high-order solver. This distribution aligns with the analysis of Ghosal (1996), who suggested that numerical dissipation due to second-order discretisation may be similar or even larger than that of the subgrid model over energy containing low wavenumber range. This challenges the usual assumption in LES practice that all or most of the required subgrid dissipation is imparted by the SGS model and that from the numerics is negligibly small.
A major finding was related to the solvers' comparison. The total non-viscous dissipation normalised by the viscous dissipation, (E num + E model )∕E , from the high-order solver was significantly lower (around 1/4th) than the similar-cost second-order runs. The major factor in this difference was the numerical term E num ∕E . For instance, in the impingement region of the combustor case, almost all the difference was due to E num ∕E . This provides a clear explanation to the accuracy vs cost gains observed in the original comparative study . Note that the quantity (E num + E model )∕E is directly proportional to the effective viscosity in LES. Its lower values in the highorder runs mean that one operates at a higher effective Reynolds number for a given cost.
A limitation of the current study is that the grid-effects, i.e., the effect of grid types on dissipation rate, were not explored in the simpler test case. This could have provided an even better insight into the results of combustor port flow.
The findings highlight that numerical dissipation plays a major role in scale-resolving simulations of combustor flow scenarios. Since the high-order solver imparts low dissipation at energy-containing large scales, the energy at low and intermediate wavenumbers is better preserved, and the final LES flow field is more accurate. This is important for turbulent combustion predictions as these rely heavily on the resolution of instantaneous large-and small-scale flow structures. Combined with the ability to handle complex hybrid meshes, the high-order CG solvers with low numerical dissipation are therefore, promising tools for combustor flow simulations. Further studies could focus on swirling flow configurations, another prominent feature relevant to modern combustors.

Conflict of interest
The authors declare that they have no conflict of interest.

Informed Consent
No human or animal participants were involved in the research. All authors have agreed to publish the present article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, 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 article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article'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. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.