1.1. State of the art and limitations in tracing finite-time singularities

A strategy to identify potential singularities directly is that of time series expansions (Taylor & Green Reference Taylor and Green1937; Morf, Orszag & Frisch Reference Morf, Orszag and Frisch1980; Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983; Pelz & Gulak Reference Pelz and Gulak1997), but it has been found that numerical inaccuracies prevent a definite answer when using this technique. Numerical investigations by means of partial differential equation (PDE) solvers have therefore played the most dominant role in the exploration of finite-time singularities, as detailed below.

A difficulty in identifying singularities with the direct approach by numerical simulations is the inherent conflict that arbitrarily small structures cannot be resolved with a numerical simulation of finite resolution, which renders this problem one of the most challenging topics in computational fluid dynamics. Much work has been done in this field. In the 1980s and 1990s, several early works on direct numerical simulation of both inviscid and high-Reynolds-number viscous incompressible flows reported indications of finite-time singularities for the incompressible Euler equations (Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983; Kerr & Hussain Reference Kerr and Hussain1989; Brachet et al. Reference Brachet, Meneguzzi, Vincent, Politano and Sulem1992; Kerr Reference Kerr1993; Boratav & Pelz Reference Boratav and Pelz1994). Symmetry in the initial conditions plays an important role, as specifically mentioned and addressed in some works (Boratav & Pelz Reference Boratav and Pelz1994; Pelz & Gulak Reference Pelz and Gulak1997; Pelz Reference Pelz2001), raising the question of whether singularities are possible for problems that are not perfectly symmetric. However, these works do not allow a definite answer to the question of finite-time singularities; see also the review articles on this topic (Gibbon Reference Gibbon2008; Hou & Li Reference Hou and Li2008). One of the main reasons why the results of these studies have been inconclusive is that the spatial resolution has been limited due to the computational power and computational approaches available at the time. Numerical results shown in Hou & Li (Reference Hou and Li2008) suggest that dynamic depletion of vortex-stretching could be a mechanism that prevents a finite-time blowup, but the same authors report evidence for a finite-time singularity for a different flow configuration with solid boundaries in a later work (Luo & Hou Reference Luo and Hou2014).

In terms of the flow configuration being studied, numerical investigations on finite-time singularities can be categorised as follows. The Taylor–Green vortex has been analysed in Brachet et al. (Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983), Brachet et al. (Reference Brachet, Meneguzzi, Vincent, Politano and Sulem1992), Brachet (Reference Brachet1991), Shu et al. (Reference Shu, Don, Gottlieb, Schilling and Jameson2005), Cichowlas & Brachet (Reference Cichowlas and Brachet2005) and Bustamante & Brachet (Reference Bustamante and Brachet2012), and for a regularised problem considering the Euler–Voigt equations in Larios et al. (Reference Larios, Petersen, Titi and Wingate2018); the high-symmetry Kida–Pelz initial condition in Hou & Li (Reference Hou and Li2008), Cichowlas & Brachet (Reference Cichowlas and Brachet2005) and Grafke et al. (Reference Grafke, Homann, Dreher and Grauer2008); colliding Lamb dipoles in Orlandi, Pirozzoli & Carnevale (Reference Orlandi, Pirozzoli and Carnevale2012); and other perturbed cylindrical vortex tubes in Kerr & Hussain (Reference Kerr and Hussain1989), Kerr (Reference Kerr1993), Hou & Li (Reference Hou and Li2008), Grauer, Marliani & Germaschewski (Reference Grauer, Marliani and Germaschewski1998) and Kerr (Reference Kerr2013). Most studies use spectral methods as discretisation schemes.

For these direct numerical simulations, common approaches to trace singularities are monitoring the maximum vorticity $\Vert \boldsymbol {\omega } \Vert _{\infty }$ over time (see the Beale–Kato–Majda theorem (Beale et al. Reference Beale, Kato and Majda1984)) and the ‘analyticity strip’ method (see Sulem, Sulem & Frisch Reference Sulem, Sulem and Frisch1983), which aims to capture the smallest scales of the flow. The width of the analyticity strip $\delta (t)$, obtained from fitting the energy spectrum to $E(k,t) = C(t) k^{-n(t)}\exp (-2 k \delta (t))$, is monitored over time for successively finer spatial resolutions up to a resolution for which extrapolations of $\delta (t)$ allow one to conclude whether $\delta (t)$ reaches $0$ in finite time (finite-time singularity) or decreases only exponentially in time (regularity at all times).

Numerical results for the three-dimensional inviscid Taylor–Green vortex shown in Brachet et al. (Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983) and Cichowlas & Brachet (Reference Cichowlas and Brachet2005) indicate only an exponential decay, but this might be due to the limited spatial resolution and also to the fact that only small times of the Taylor–Green vortex flow have been considered, so that a finite-time singularity cannot be excluded by these results. In a later work by Bustamante & Brachet (Reference Bustamante and Brachet2012), a change in regime indicating potentially faster-than-exponential decay is reported, and the results are ‘not inconsistent with the occurrence of a singularity’, but again resolutions higher than the maximum one of $4096^3$ would be required for definite answers. In Cichowlas & Brachet (Reference Cichowlas and Brachet2005) it is estimated that conclusions regarding finite-time singularities using the analyticity strip method would require spatial resolutions of $(16\,\mathrm {k})^3$ to $(32\,\mathrm {k})^3$ for the Kida–Pelz initial data. A recent study by Campolina & Mailybaev (Reference Campolina and Mailybaev2018) suggests that the resolution available via classical direct numerical simulation is not sufficient to investigate blowup.

The development of pancake-like structures with exponentially growing vorticity during the early development of turbulence from smooth initial data is studied in Agafontsev, Kuznetsov & Mailybaev (Reference Agafontsev, Kuznetsov and Mailybaev2015). In Kerr (Reference Kerr2013), a new kind of analysis based on rescaled vorticity moments is proposed (studying anti-parallel vortex tubes), and only double-exponential growth in vorticity is observed, as opposed to the singular behaviour suspected in a previous work by Kerr (Reference Kerr1993). A model describing a cascade of transformations between vortex filaments and sheets potentially explaining the mechanism of singularity formation in the Euler equations is proposed in Brenner, Hormoz & Pumir (Reference Brenner, Hormoz and Pumir2016). Another model has been described recently in Moffatt (Reference Moffatt2019). In McKeown et al. (Reference McKeown, Ostilla-Mónico, Pumir, Brenner and Rubinstein2018), an iterative cascade of instabilities for head-on collisions of vortex rings is investigated both experimentally and numerically.

1.2. Energy dissipation anomaly

We now focus on the evolution of kinetic energy in incompressible Euler flows. Of particular interest is the question of whether inviscid flows are able to dissipate energy, and if so, by which mechanism such a behaviour can be explained, given that no viscous effects are present. The kinetic energy dissipation equation – valid for incompressible viscous ($\nu > 0$) flows with continuously differentiable solution on a domain with periodic boundaries – reads (Onsager Reference Onsager1949; Eyink & Sreenivasan Reference Eyink and Sreenivasan2006)

(1.1)\begin{equation} \frac{\partial E(t,\nu)}{\partial t} ={-} \int_{\varOmega} \nu \boldsymbol{\nabla} \boldsymbol{u}^{\nu}: \boldsymbol{\nabla} \boldsymbol{u}^{\nu} \,\mathrm{d}\varOmega, \end{equation}

which implies conservation of energy in the inviscid limit $\nu = 0$ provided that the solution is sufficiently regular. However, from phenomenological descriptions of turbulence, there is empirical evidence that the dissipation rate does not tend to zero in the limit $Re\rightarrow \infty$ or $\nu \rightarrow 0$ but takes a positive value independent of $\nu$, which is known as the dissipation anomaly or the zeroth law of turbulence (Eyink Reference Eyink2008; Dubrulle Reference Dubrulle2019). As noted in Eyink (Reference Eyink2008), this was first observed by Taylor (Reference Taylor1935), and also Kolmogorov's similarity theory of turbulence (Kolmogorov Reference Kolmogorov1991) is based on the assumption of a non-vanishing energy dissipation rate in the inviscid limit. Numerical evidence that the dissipation rate is independent of $\nu$ for large $Re$ is for example given in Sreenivasan (Reference Sreenivasan1998), Kaneda et al. (Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003) and Orlandi et al. (Reference Orlandi, Pirozzoli and Carnevale2012), and experimental evidence for example in Pearson, Krogstad & van de Water (Reference Pearson, Krogstad and van de Water2002) and Dubrulle (Reference Dubrulle2019). Under certain regularity or smoothness assumptions (existence of a strong $L^3$ limit; we refer to Duchon & Robert (Reference Duchon and Robert2000), Drivas & Nguyen (Reference Drivas and Nguyen2019) and Drivas & Eyink (Reference Drivas and Eyink2019) for a precise discussion), weak Euler solutions are the $\nu \to 0$ limit of Leray–Hopf weak solutions $\boldsymbol {u}^{\nu }$ of the Navier–Stokes equations, so that the dissipation rate in the inviscid limit equals the viscous dissipation rate in the limit $\nu \rightarrow 0$:

(1.2)\begin{equation} \frac{\partial E_{\nu=0}(t)}{\partial t} = \lim_{\nu \rightarrow 0} \frac{\partial E(t,\nu)}{\partial t} = \lim_{\nu \rightarrow 0} - \nu \int_{\varOmega} \boldsymbol{\nabla} \boldsymbol{u}^{\nu} : \boldsymbol{\nabla} \boldsymbol{u}^{\nu} \,\mathrm{d}\varOmega ={-} D(t) \leq 0, \end{equation}

where anomalous energy dissipation means that $D(t) > 0$ for some (or all) $t > t_*$. (As a consequence, the enstrophy is inversely proportional to the viscosity for large Reynolds numbers in the case of anomalous dissipation with $D(t) > 0$.) In general, weak Euler solutions may be neither unique nor the zero-viscosity limit of weak Navier–Stokes solutions, which might themselves be non-unique (Isett Reference Isett2017; Buckmaster & Vicol Reference Buckmaster and Vicol2019, Reference Buckmaster and Vicol2020; Buckmaster et al. Reference Buckmaster, Masmoudi, Novack and Vicol2021; Daneri, Runa & Székelyhidi Reference Daneri, Runa and Székelyhidi2021). In this sense, we note that the first equality in (1.2) is a conditional one. Let us also refer to Brenier, De Lellis & Székelyhidi (Reference Brenier, De Lellis and Székelyhidi2011) and Wiedemann (Reference Wiedemann2017) regarding the topic of weak–strong uniqueness of Euler solutions.

The above argument already indicates that the theory explaining dissipation of energy in the absence of viscosity is related to the spatial regularity of the solution. According to Onsager (Reference Onsager1949), energy dissipation in three-dimensional incompressible flows can take place in the absence of viscosity through the formation of singularities, with the cascade from large to (arbitrarily) small scales taking place in finite time; see also  Eyink & Sreenivasan (Reference Eyink and Sreenivasan2006). (Hence, the mechanism explaining the occurrence of kinetic energy dissipation in the limit $\nu \rightarrow 0$ is that the velocity gradient might tend to infinity in this limit. Specifically, Onsager (Reference Onsager1949) wrote: ‘In the absence of viscosity, the standard proof of the conservation of energy does not apply, because the velocity field does not remain differentiable!’ Interestingly, Onsager did not consider the energy-dissipating behaviour of inviscid flows (‘ideal turbulence’) an anomalous behaviour but rather a matter of fact.) According to Onsager's conjecture, dissipation of energy may occur if the velocity is Hölder continuous with exponent ${\leq }1/3$ (while Onsager's assertion says that energy is conserved for exponents $> 1/3$; see Eyink (Reference Eyink1994), Constantin, Weinan & Titi (Reference Constantin, Weinan and Titi1994), Duchon & Robert (Reference Duchon and Robert2000) and Cheskidov et al. (Reference Cheskidov, Constantin, Friedlander and Shvydkoy2008) for proofs). For mathematical literature dealing with proofs of Onsager's conjecture we refer to De Lellis & Székelyhidi (Reference De Lellis and Székelyhidi2013), De Lellis & Székelyhidi (Reference De Lellis and Székelyhidi2014), Buckmaster, De Lellis & Székelyhidi (Reference Buckmaster, De Lellis and Székelyhidi2016), Buckmaster et al. (Reference Buckmaster, De Lellis, Székelyhidi and Vicol2018), Isett (Reference Isett2018) and references therein, where dissipative weak Euler solutions up to Onsager's critical regularity have been constructed using convex integration techniques. New insights from these works thus confirm that the Hölder exponent of $1/3$ is indeed the critical one in terms of energy dissipation. As noted in Dubrulle (Reference Dubrulle2019), the original Kolmogorov cascade picture implies irregularities of the velocity field (at least locally) with Hölder exponent $\leq 1/3$, but it was Onsager who established the link between energy dissipation and irregularities of the velocity field for the Euler equations.

The above considerations might explain why this phenomenon is known as the ‘kinetic energy dissipation anomaly’; an alternative term used e.g. in Dubrulle (Reference Dubrulle2019) is ‘inertial dissipation’ (as opposed to viscous dissipation). The one-dimensional inviscid Burgers equation (Burgers Reference Burgers1948) with formation of a shock and the associated dissipation of energy serves as a prominent and well-understood example; it is discussed, for example, in Sulem et al. (Reference Sulem, Sulem and Frisch1983) in the context of finite-time singularities and in Dubrulle (Reference Dubrulle2019) in the context of inertial energy dissipation. In Josserand, Pomeau & Rica (Reference Josserand, Pomeau and Rica2020), the phenomenon of energy dissipation through finite-time singularities is illustrated for another one-dimensional model problem, the nonlinear Schrödinger equation. For the two-dimensional incompressible Euler equations, it is known that singularities cannot develop in finite time from smooth initial data (Eyink & Sreenivasan Reference Eyink and Sreenivasan2006). Regarding three-dimensional turbulent flows, Onsager's conjecture appears to be widely accepted by now, with the occurrence of singularities representing a building block of modern understandings of turbulence (Dubrulle Reference Dubrulle2019).

1.3. Interplay between physics and numerics

The present work focuses on the numerical solution of turbulent flows and particularly the inviscid limit. From a closer look at numerical simulations of the inviscid Taylor–Green vortex, it can be observed that many of these simulations have been performed mainly for small times up to $t\approx 5$ (up to $t=4$ in  Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983, Brachet et al. Reference Brachet, Meneguzzi, Vincent, Politano and Sulem1992, Cichowlas & Brachet Reference Cichowlas and Brachet2005 and Bustamante & Brachet Reference Bustamante and Brachet2012, and up to $t\leq 6$ in Shu et al. Reference Shu, Don, Gottlieb, Schilling and Jameson2005 and Chapelier, De La Llave Plata & Renac Reference Chapelier, De La Llave Plata and Renac2012), but not beyond the time at which finite-time singularities have been suspected, especially not up to the time at which the transition to a fully turbulent state takes place, with maximum kinetic energy dissipation rate at time $t = 8 - 9$ (expected from high-Reynolds-number viscous simulations (see Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983)) and subsequent decaying turbulence. As mentioned by some of these works, one reason for this is that the results for a specific resolution are no longer reliable at later times once the flow becomes under-resolved (note that the spatial resolution is severely limited by computational resources, even with the large supercomputers available today) . The key aspect, however, is that numerical simulations of the incompressible Euler equations are very challenging in terms of energy stability and numerical blowup of the discretisation scheme. Often, a lack of robustness of the numerical discretisation scheme is reported for this challenging inviscid Taylor–Green vortex problem (see for example Chapelier et al. Reference Chapelier, De La Llave Plata and Renac2012 and Winters et al. Reference Winters, Moura, Mengaldo, Gassner, Walch, Peiro and Sherwin2018).

The lack of understanding of what is to be expected in terms of kinetic energy dissipation from a physical perspective (the energy dissipation anomaly discussed in § 1.2) is manifested in an uncertainty regarding the optimal design of discretisation schemes from a numerical perspective. In Moura et al. (Reference Moura, Mengaldo, Peiró and Sherwin2017b), Moura et al. (Reference Moura, Mengaldo, Peiró and Sherwin2017a) and Piatkowski (Reference Piatkowski2019), it is argued that a fundamentally different behaviour in terms of energy dissipation and time reversibility is expected between viscous flows in the limit $Re\rightarrow \infty$ and inviscid flows at $Re = \infty$. Especially in numerical studies, it is often assumed that energy conservation holds for an exact solution of the Euler equations not only in two space dimensions but also in three space dimensions; see for example Shu et al. (Reference Shu, Don, Gottlieb, Schilling and Jameson2005), Bustamante & Brachet (Reference Bustamante and Brachet2012), Grauer et al. (Reference Grauer, Marliani and Germaschewski1998), Chapelier et al. (Reference Chapelier, De La Llave Plata and Renac2012), Winters et al. (Reference Winters, Moura, Mengaldo, Gassner, Walch, Peiro and Sherwin2018), Schroeder (Reference Schroeder2019), Krais et al. (Reference Krais, Schnücke, Bolemann and Gassner2020) and the recent review article by Coppola, Capuano & de Luca (Reference Coppola, Capuano and de Luca2019), to mention just a few. Numerical schemes that are exactly energy-conserving can indeed be constructed and have the advantage that nonlinear blowup of the numerical discretisation scheme can be avoided in the challenging inviscid limit. For these reasons, energy-conserving schemes appear to be the current gold standard for the simulation of this type of problem. Inviscid Taylor–Green vortex simulations performed in Schroeder (Reference Schroeder2019) using exactly divergence-free, energy-conserving discretisation methods result in an exact conservation of energy, and the results are considered superior to those of simulations with upwind fluxes that show a dissipative behaviour.

However, the use of energy-conserving schemes is accompanied by a major limitation, namely that it excludes – by construction – the occurrence of anomalous energy dissipation. Onsager's conjecture forces one to rethink whether it is really a desirable quality criterion that a numerical method preserves the kinetic energy exactly in the inviscid limit $\nu = 0$. If Onsager's hypothesis is true, there is an inconsistency between the physical dissipation behaviour and the numerical dissipation behaviour of energy-conserving discretisation methods. Thus, energy-conserving numerical methods would result in an $O(1)$ error in the case of inviscid flows with anomalous/inertial energy dissipation. Since no energy can leave the system for such a discretisation scheme, energy accumulates in small scales, a well-known phenomenon called thermalisation. In terms of the kinetic energy spectrum, an energy-conserving numerical scheme typically leads to an unphysical equipartitioning of energy when simulating beyond the time of the finite-time singularity (Orlandi Reference Orlandi2009; Ray et al. Reference Ray, Frisch, Nazarenko and Matsumoto2011; Orlandi et al. Reference Orlandi, Pirozzoli and Carnevale2012). The energy-conserving results of Schroeder (Reference Schroeder2019) indeed show such a behaviour. One may conclude that the application of energy-conserving numerical methods is only reasonable for times $t < t_*$ before a potential singularity forms, since anomalous dissipation might occur afterwards. Further, one may formulate that a numerical scheme must contain mechanisms of dissipation as a minimal requirement in order to address the topic of anomalous energy dissipation. To describe discretisation schemes suitable for investigating anomalous dissipation more precisely, we consider it a prerequisite to use consistent and stable discretisation schemes whose dissipation mechanisms are coupled to under-resolution effects in the numerical approximation of the solution, i.e. the dissipation mechanisms act on the finest resolved scales and shift to smaller scales under mesh refinement. As we discuss in more detail below, it appears to be unclear mathematically whether such a scheme is able to find a dissipative weak solution of the Euler equations.

By studying Galerkin-truncated, energy-conserving simulations of the one-dimensional Burgers equation, the work of Ray et al. (Reference Ray, Frisch, Nazarenko and Matsumoto2011) describes an interesting phenomenon, called ‘tyger phenomenon’ in that work, where short-wavelength oscillations occur out of the blue in the presence of singularities, that finally lead to thermalisation. The importance of numerical dissipation to avoid the effect of thermalisation for this one-dimensional Burgers problem is emphasised in a recent study by Murugan et al. (Reference Murugan, Frisch, Nazarenko, Besse and Ray2020), which can therefore be seen in close analogy to the present work focusing on three-dimensional Euler problems.

The situation is less complicated for two-dimensional Euler flows that are non-dissipative. In that case, it can be expected that the kinetic energy dissipation rate converges to zero under mesh refinement for a consistent and energy-stable discretisation scheme, and that there is per se no conflict with physics if an energy-conserving scheme is applied.

1.4. An indirect approach to identify finite-time singularities through energy considerations

The indirect approach to identify finite-time singularities relies on the physical intuition that the appearance of anomalous energy dissipation in free decay from smooth initial data requires a finite-time singularity. The basic idea is to capture the temporal evolution of the kinetic energy by a numerical method with appropriate inbuilt dissipation mechanisms as described in § 1.3. If grid-convergence to a dissipative solution with non-zero kinetic energy dissipation rate can be demonstrated numerically, indirect evidence of a finite-time singularity is provided by the following line of argument:

1. (i) Assume convergence of a sequence of numerical solutions to a dissipative weak Euler solution for $h \to 0$.

2. (ii) Weak–strong uniqueness holds for dissipative weak Euler solutions (Lions Reference Lions1996; Brenier et al. Reference Brenier, De Lellis and Székelyhidi2011; Wiedemann Reference Wiedemann2017).

3. (iii) Supposing that an energy dissipation anomaly with non-zero kinetic energy dissipation rate is observed, it follows from items (i) and (ii) that a strong solution cannot exist but must have become singular.

The conclusion in item (iii) is based on an indirect proof. Assume that a strong solution exists. This strong solution is energy-conserving. By weak–strong uniqueness, the weak solution must be identical and, therefore, energy-conserving. By contradiction it follows that a strong solution cannot exist.

Note that assumption (i) could be weakened. Rather than (i), only convergence to a dissipative generalised weak solution in the sense of Lions (Reference Lions1996) or DiPerna & Majda (Reference DiPerna and Majda1987) is required. Interestingly, this weaker assumption might be provable for the limit $h \to 0$, as it has in fact been proven for the limit $\nu \to 0$ (at least along subsequences). Since the existence of generalised weak Euler solutions as limits along subsequences relies on very general compactness arguments that require only $L^2$ (kinetic energy) bounds for both the Lions and DiPerna & Majda theories, a proof for the limit $h \to 0$ for a discontinuous Galerkin (DG) Euler discretisation scheme as used here might be conceivable in analogy to what has been shown for the limit $\nu \to 0$. To the best of our knowledge, such a proof is still outstanding.

Let us note that a similar idea to identify singularities experimentally based on energy arguments has been used by Saw et al. (Reference Saw, Kuzzay, Faranda, Guittonneau, Daviaud, Wiertel-Gasquet, Padilla and Dubrulle2016) and Kuzzay et al. (Reference Kuzzay, Saw, Martins, Faranda, Foucaut, Daviaud and Dubrulle2017), who calculate the inertial dissipation at scale $l$ from particle image velocimetry measurements. To the best of the authors’ knowledge, the present study makes the first attempt to use energy arguments for singularity detection in numerical simulations of three-dimensional Euler flows.

From numerical simulations of viscous problems at finite Reynolds number, there are indications that the zeroth law of turbulence holds for the Taylor–Green vortex problem. Numerical results for the kinetic energy dissipation rate for increasing Reynolds numbers up to $Re = 3000$ in Brachet et al. (Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983) and additional results for higher Reynolds numbers of $Re = 5000$ in Brachet (Reference Brachet1991), $Re=10\,000$ in Arndt et al. (Reference Arndt, Fehn, Kanschat, Kormann, Kronbichler, Munch, Wall and Witte2020) and $Re=20\,000$ in Lamballais et al. (Reference Lamballais, Dairay, Laizet and Vassilicos2019) strongly suggest that the function $D(t)$ does not tend to zero in the limit $\nu \rightarrow 0$. This argument is summarised in figure 1, which shows results for viscous and inviscid simulations of the Taylor–Green vortex obtained with the present discretisation approach. For each Reynolds number, results are shown for two resolutions of the numerical discretisation approach to judge whether the results are mesh-independent. We achieve grid-converged results for all finite Reynolds numbers shown in figure 1. In the inviscid limit, the temporal evolution of the kinetic energy is almost indistinguishable for the two finest resolutions, while small differences in the energy dissipation rate are still visible between the two resolutions $4096^3$ and $8192^3$ at later times around the dissipation maximum and beyond. However, the onset of dissipation around $t\approx 6$ appears to be grid-converged also for this challenging inviscid simulation. These results are consistent with grid-convergence to a dissipative solution of the incompressible Euler equations for the three-dimensional Taylor–Green problem. Figure 1 therefore summarises the main result of the present work.

Figure 1.

Temporal evolution of kinetic energy (a) and kinetic energy dissipation rate (b) for the three-dimensional Taylor–Green vortex problem for increasing Reynolds numbers of $Re=100,200,400,800,1600,3200,\infty$. For each $Re$, results are shown for two mesh resolutions (fine mesh as solid line, coarse mesh as dashed-dotted line). The effective resolutions (see § 3 for a definition) are $64^3, 128^3$ for $Re=100$, $128^3, 256^3$ for $Re=200, 400$, $256^3, 512^3$ for $Re=800$, $1024^3, 2048^3$ for $Re=1600$, $2048^3, 4096^3$ for $Re=3200$, and $4096^3, 8192^3$ for $Re=\infty$. The results suggest that the kinetic energy reduces to a value as low as approximately $0.02$ at time $t=20$ for large Reynolds numbers, and that a similar amount of energy dissipation also takes place in the inviscid limit.

While the accumulation of energy in small scales in case of energy-preserving schemes is unphysical, it can be exploited under certain circumstances in order to gain insights into the physical dissipation behaviour. In Cichowlas et al. (Reference Cichowlas, Bonaïti, Debbasch and Brachet2005), an effective dissipation is estimated from the small-scale thermalised energy of energy-conserving, spectrally truncated Euler simulations, and it is found that the large-scale Euler dynamics are similar to high-Reynolds-number Navier–Stokes dynamics. Although the underlying numerical methods in that work are different from those of the present study (the energy dissipation rate is derived by a postprocessing of results in Cichowlas et al. (Reference Cichowlas, Bonaïti, Debbasch and Brachet2005), whereas it is simulated directly in the present work), there are interesting parallels. The onset of dissipation around $t=5 - 6$ and the dissipation maximum around $t = 8 - 9$ appear to be very similar to the present results, so that the two studies can be seen to complement and verify each other. From this perspective, the study by Cichowlas et al. (Reference Cichowlas, Bonaïti, Debbasch and Brachet2005) supports the main conclusions of our work.

Instead of considering a two-parameter limit problem $h \to 0, \nu \to 0$ as illustrated in figure 1, the remainder of this work focuses on the one-parameter limit $h \to 0$ for the inviscid limit $\nu = 0$. Let us explain this decision in more detail. A two-parameter study $h \to 0, \nu \to 0$ would technically not be realisable due to the large amount of computational costs required for such simulations: already for moderate Reynolds numbers of $Re={O}(10^4)$, as considered for example in Arndt et al. (Reference Arndt, Fehn, Kanschat, Kormann, Kronbichler, Munch, Wall and Witte2020) and Lamballais et al. (Reference Lamballais, Dairay, Laizet and Vassilicos2019), the spatial resolutions required for grid-convergence are comparable to the highest-resolution simulations that we are able to realise in the present work for the inviscid limit. These highest-resolution simulations require computational costs of tens of millions of CPUh in the present study, despite the fact that we use a highly efficient implementation that is well optimised for the hardware under consideration. With the goal of realising spatial resolutions as high as those shown in the present study, computational costs would allow us to consider only a single finite Reynolds number beyond what is shown in figure 1, which explains why we immediately address the inviscid limit $\nu = 0$.

1.5. Outline

The rest of this article is organised as follows. We describe the mathematical model of the incompressible Navier–Stokes or Euler equations and its numerical discretisation in space and time in § 2. Section 3 shows results for the one-dimensional Burgers equation with formation of a shock, a two-dimensional shear layer problem with a numerical investigation of the kinetic energy dissipation in the limit $\nu \rightarrow 0$, and finally the three-dimensional Taylor–Green vortex problem that has been suspected to exhibit finite-time singularities in the inviscid limit. In § 4, we summarise our results, draw conclusions and raise questions based on the present observations.

2.1. Temporal discretisation – high-order projection method

Discretisation in time is based on projection methods which aim at obtaining computationally efficient incompressible flow solvers by decoupling the velocity and pressure unknowns (Karniadakis & Sherwin Reference Karniadakis and Sherwin2005). A convection–diffusion-type problem is solved for the velocity and a Poisson equation for the pressure, with a subsequent projection of the velocity onto the space of solenoidal vector fields according to the Helmholtz decomposition. We use the high-order dual splitting scheme proposed in Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991), which consists of the following four sub-steps:

(2.4)$$\begin{gather} \frac{\displaystyle\gamma_0^n\hat{\boldsymbol{u}}-\sum_{i=0}^{J-1}\alpha_i^n\boldsymbol{u}^{n-i}}{\Delta t_n} ={-} \sum_{i=0}^{J-1}\beta_i^n \boldsymbol{\nabla} \boldsymbol{\cdot} {(\boldsymbol{u}^{n-i} \otimes \boldsymbol{u}^{n-i})}, \end{gather}$$

(2.5)$$\begin{gather}-\nabla^2 p^{n+1} ={-}\frac{\gamma_0^n }{\Delta t_n}\boldsymbol{\nabla} \boldsymbol{\cdot} {\hat{\boldsymbol{u}}}, \end{gather}$$

(2.6)$$\begin{gather}\hat{\hat{\boldsymbol{u}}} = \hat{\boldsymbol{u}} - \frac{\Delta t_n}{\gamma_0^n} \boldsymbol{\nabla}{p^{n+1}}, \end{gather}$$

(2.7)$$\begin{gather}\frac{\gamma_0^n}{\Delta t_n} \boldsymbol{u}^{n+1} - \nu \nabla^2 \boldsymbol{u}^{n+1} = \frac{\gamma_0^n}{\Delta t_n}\hat{\hat{\boldsymbol{u}}}, \end{gather}$$

where $n$ denotes the current time step, in which the equations are integrated from time $t_n$ to time $t_{n+1} = t_n + \Delta t_n$. The temporal discretisation is based on a backward differentiation formula of order $J$ with coefficients $\gamma _0^n$ and $\alpha _i^n,\ i=0,\ldots,J-1$. In the first sub-step (2.4), the convective term is treated explicitly in time by using a high-order extrapolation scheme of order $J$ with coefficients $\beta _i^n, i=0, \ldots ,J-1$. In the next two sub-steps, a pressure Poisson equation (2.5) is solved and a divergence-free velocity is obtained in the projection step (2.6). Finally, the viscous term is taken into account in the last sub-step (2.7), which can be omitted in the case of the Euler equations, $\boldsymbol {u}^{n+1} = \hat {\hat {\boldsymbol {u}}}$. The explicit treatment of the convective term implies a restriction of the time step size according to the Courant–Friedrichs–Lewy (CFL) condition. For reasons of computational efficiency, we use adaptive time stepping with variable time step sizes $\Delta t_n$, readjusting the time step size after each time step in such a way as to maximise the time step size and operate close to the CFL stability limit. We introduce the CFL condition below since it also depends on the spatial discretisation scheme. A second-order accurate time integration scheme with $J=2$ is used in the present work.

2.2. Spatial discretisation – high-order discontinuous Galerkin method

A key element of the present study is the use of a novel high-order DG discretisation for incompressible flows. This discretisation approach has been developed recently and is documented in a series of publications (Fehn, Wall & Kronbichler Reference Fehn, Wall and Kronbichler2017; Krank et al. Reference Krank, Fehn, Wall and Kronbichler2017; Fehn, Wall & Kronbichler Reference Fehn, Wall and Kronbichler2018b; Fehn et al. Reference Fehn, Kronbichler, Lehrenfeld, Lube and Schroeder2019), with a focus on the stability of projection methods in Fehn et al. (Reference Fehn, Wall and Kronbichler2017), and a focus on the stability and dissipation characteristics of DG discretisations for under-resolved turbulence in Fehn et al. (Reference Fehn, Wall and Kronbichler2018b, Reference Fehn, Kronbichler, Lehrenfeld, Lube and Schroeder2019). Of particular importance are the inbuilt dissipation mechanisms of DG methods acting on the finest resolved scales of the flow. In order to investigate the phenomenon of anomalous energy dissipation, we consider it a prerequisite to use a discretisation scheme that provides the flexibility to find both dissipative and non-dissipative solutions. Let us point out that all simulations shown in this work are performed as direct numerical simulations, i.e. without any form of explicit turbulence model or numerical viscosity. More precisely, upwind-type numerical fluxes represent the main dissipation mechanism of the present scheme; see also Ainsworth (Reference Ainsworth2004) for a general dispersion and dissipation analysis and distribution of dissipation among the resolved scales. Note that the $H(\text {div})$-stabilisation presented below also contributes with a certain amount of dissipation. For reasons of brevity and to focus on the main aspects, we avoid technical aspects related to the imposition of boundary conditions in the following, but refer to the original publications.

We assume a computational domain $\varOmega _h= \bigcup _{e=1}^{N_{{el}}} \varOmega _{e} \in \mathbb {R}^d$ consisting of conforming quadrilateral or hexahedral elements $\varOmega _e$, $e=1,\ldots, N_{{el}}$. A common abstraction of finite element methods is to define a mapping $\boldsymbol {x}^e(\boldsymbol {\xi })$ from a reference element $\tilde {\varOmega }_e = [0,1]^d$ with Cartesian coordinates $\boldsymbol {\xi }$ to element $\varOmega _e$ in physical space with coordinates $\boldsymbol {x}$, and approximate the solution within each element by polynomials defined in reference coordinates. Here, the numerical solution in $d=2,3$ is represented by a tensor product of one-dimensional Lagrange polynomials with a Legendre–Gauss–Lobatto point distribution, and is allowed to exhibit discontinuities between elements in an $L^2$-conforming sense. The spaces of shape functions are then given as

(2.8)\begin{align} \mathcal{V}^{u}_{h} &= \lbrace\boldsymbol{u}_h\in [L^2(\varOmega_h)]^d: \boldsymbol{u}_h\left(\boldsymbol{x}^e(\boldsymbol{\xi})\right)\vert_{\varOmega_{e}}= \tilde{\boldsymbol{u}}_h^e(\boldsymbol{\xi})\vert_{\tilde{\varOmega}_{e}}\in \mathcal{V}^{u}_{h,e}=[\mathcal{Q}_{k}(\tilde{\varOmega}_{e})]^d,\nonumber\\ &\quad \vphantom{\left[L^2(\varOmega_h)\right]^d} \forall\ e=1,\ldots,N_{{el}}\rbrace, \end{align}

(2.9)\begin{align} \mathcal{V}^{p}_{h} &= \lbrace p_h\in L^2(\varOmega_h): p_h(\boldsymbol{x}^e(\boldsymbol{\xi}))\vert_{\varOmega_{e}} = \tilde{p}_h^e(\boldsymbol{\xi})\vert_{\tilde{\varOmega}_{e}}\in \mathcal{V}^{p}_{h,e}=\mathcal{Q}_{k-1}(\tilde{\varOmega}_{e}), \nonumber\\ &\quad \vphantom{\left(\boldsymbol{x}^e(\boldsymbol{\xi})\right)\vert_{\varOmega_{e}} = \tilde{p}_h^e(\boldsymbol{\xi})\vert_{\tilde{\varOmega}_{e}}\in \mathcal{V}^{p}_{h,e}} \forall\ e=1,\ldots,N_{{el}}\rbrace. \end{align}

Here, $\mathcal {Q}_{k}$ denotes the space of tensor-product polynomials of degree $k$, and we highlight the ambiguity in notation related to the spatial wavenumber $k$ in the context of energy spectra. The polynomial degree of the velocity shape functions is $k$, while it is $k-1$ for the pressure for reasons of inf–sup stability. In this work, we only consider problems with Cartesian meshes so that the mapping $\boldsymbol {x}^e(\boldsymbol {\xi })$ is an affine transformation. At the interface between two elements $e^-$ and $e^+$, the numerical solution is not unique, and we use the superscripts $-$ and $+$ to denote the solution from the two sides of an interface. When referring to the element $e$, we denote interior information by $-$, and the outward-pointing unit normal vector by $\boldsymbol {n}=\boldsymbol {n}^-$. To define numerical fluxes, we introduce the average operator $\{\!\{u\}\!\}=(u^- + u^+)/2$ and the jump operator $\left[\kern-0.15em\left[ {u} \right]\kern-0.15em\right] =u^- \otimes \boldsymbol {n}^- + u^+ \otimes \boldsymbol {n}^+$, where $\boldsymbol {n}^+ = - \boldsymbol {n}^-$. We apply the usual abbreviations of integrals, $\left ( {v},{u} \right )_{\varOmega _{e}} = \int _{\varOmega _e} v \odot u\, \mathrm {d}\varOmega$ for volume integrals and $\left ( {v},{u} \right )_{\partial \varOmega _{e}} = \int _{\partial \varOmega _e} v \odot u \,\mathrm {d} \varGamma$ for surface integrals with $\odot$ denoting an inner product. We now state the weak formulation for all sub-steps of the projection scheme. We obtain the numerical solutions $\hat {\boldsymbol {u}}_h,\hat {\hat {\boldsymbol {u}}}_h, \hat {\hat {\hat {\boldsymbol {u}}}}_h ,\boldsymbol {u}_h^{n+1}\in \mathcal {V}^u_h$ and $p_h^{n+1}\in \mathcal {V}^p_h$ by testing with all test functions $\boldsymbol {v}_h \in \mathcal {V}^{u}_{h,e}$, $q_h \in \mathcal {V}^{p}_{h,e}$ for all elements $e=1,\ldots,N_{{el}}$.

In the first step, the convective term is discretised by the local Lax–Friedrichs flux

(2.10)\begin{align} & \left( {\boldsymbol{v}_h},{\frac{\gamma_0^n \hat{\boldsymbol{u}}_h-\sum_{i=0}^{J-1}\alpha_i^n\boldsymbol{u}^{n-i}_h}{\Delta t_n}} \right)_{\varOmega_{e}} ={-} \sum_{i=0}^{J-1} \beta_i^n \left(- \left( {\boldsymbol{\nabla}{\boldsymbol{v}_h}},{\left(\boldsymbol{u}_h \otimes \boldsymbol{u}_h\right)^{n-i}} \right)_{\varOmega_{e}} \vphantom{\left(\{\!\{\boldsymbol{u}_h \otimes \boldsymbol{u}_h\}\!\} +\frac{\varLambda}{2}\left[\kern-0.15em\left[ {\boldsymbol{u}_h} \right]\kern-0.15em\right] \right)^{n-i}}\right.\nonumber\\ &\quad \left.+ \left( {\boldsymbol{v}_h},{\left(\{\!\{\boldsymbol{u}_h \otimes \boldsymbol{u}_h\}\!\} + \frac{\varLambda}{2}\left[\kern-0.15em\left[ {\boldsymbol{u}_h} \right]\kern-0.15em\right] \right)^{n-i}\boldsymbol{\cdot} \boldsymbol{n}} \right)_{\partial\varOmega_{e}} \right), \end{align}

where the stabilisation parameter is chosen as $\varLambda = \max ( \vert \boldsymbol {u}_h^{-} \boldsymbol {\cdot } \boldsymbol {n}\vert, \vert \boldsymbol {u}_h^{+} \boldsymbol {\cdot } \boldsymbol {n}\vert )$. The pressure Poisson operator is discretised by the symmetric interior penalty Galerkin (SIPG) method

(2.11) \begin{align} & ({\boldsymbol{\nabla}{q_h}},{\boldsymbol{\nabla}{p^{n+1}_h}})_{\varOmega_{e}} - \left( {\boldsymbol{\nabla}{q_h}},{\frac{1}{2}\left[\kern-0.15em\left[ {p^{n+1}_h} \right]\kern-0.15em\right] } \right)_{\partial\varOmega_{e}} - ({q_h},{\{\!\{\boldsymbol{\nabla}{p^{n+1}_h}\}\!\}\boldsymbol{\cdot}\boldsymbol{n}})_{\partial\varOmega_{e}} \nonumber\\ &\quad + ( {q_h},{\tau\left[\kern-0.15em\left[ {p^{n+1}_h} \right]\kern-0.15em\right] \boldsymbol{\cdot}\boldsymbol{n}} )_{\partial\varOmega_{e}} ={-} \frac{\gamma_0^n}{\Delta t_n} ( - ( {\boldsymbol{\nabla}{q_h}},{\hat{\boldsymbol{u}}_h} )_{\varOmega_{e}} + ( {q_h},{\{\!\{\hat{\boldsymbol{u}}_h\}\!\}\boldsymbol{\cdot}\boldsymbol{n}})_{\partial\varOmega_{e}}), \end{align}

where a central flux is used for the velocity divergence operator on the right-hand side of the pressure Poisson equation. A central flux is also used for the pressure gradient term in the projection step

(2.12) \begin{equation} ( {\boldsymbol{v}_h},{\hat{\hat{\boldsymbol{u}}}_h} )_{\varOmega_{e}} = ( {\boldsymbol{v}_h},{\hat{\boldsymbol{u}}_h} )_{\varOmega_{e}} - \frac{\Delta t_n}{\gamma_0^n} ( - ( {\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}_h}},{p^{n+1}_h})_{\varOmega_{e}} + ( {\boldsymbol{v}_h},{\{\!\{p^{n+1}_h\}\!\}\boldsymbol{n}} )_{\partial\varOmega_{e}}). \end{equation}

The viscous term is also discretised by the SIPG method,

(2.13) \begin{array}{l} \textrm{ }{\left( {{\boldsymbol{v}_h},\frac{{\gamma _0^n}}{{{\Delta }{t_n}}}{{\hat {\hat {\hat {\boldsymbol{u}}}}}_h}} \right)_{{{\Omega }_e}}} + \textrm{ }{\left( {\nabla {\boldsymbol{v}_h},\nu \nabla {{\hat {\hat {\hat {\boldsymbol{u}}}}}_h}} \right)_{{{\Omega }_e}}} - \textrm{ }{\left( {\nabla {\boldsymbol{v}_h},\frac{\nu }{2}\left[\kern-0.15em\left[{{\hat {\hat {\hat {\boldsymbol{u}}}}}_h}\right]\kern-0.15em\right]} \right)_{\partial {{\Omega }_e}}}\\ \quad - \textrm{ }{\left( {{\boldsymbol{v}_h},\nu \textrm{ \{ \{ }\nabla {{\hat {\hat {\hat {\boldsymbol{u}}}}}_h}\} \} \cdot {\boldsymbol{n}}} \right)_{\partial {{\Omega }_e}}} + \textrm{ }{\left( {{\boldsymbol{v}_h},\nu \tau \left[\kern-0.15em\left[{{\hat {\hat {\hat {\boldsymbol{u}}}}}_h}\right]\kern-0.15em\right] \cdot {\boldsymbol{n}}} \right)_{\partial {{\Omega }_e}}} = \textrm{ }{\left( {{\boldsymbol{v}_h},\frac{{\gamma _0^n}}{{{\Delta }{t_n}}}{{\hat {\hat {\boldsymbol{u}}}}_h}} \right)_{{{\Omega }_e}}}, \end{array}

while this step is skipped, $\hat {\hat {\hat {\boldsymbol {u}}}}_h = \hat {\hat {\boldsymbol {u}}}_h$, in the inviscid limit when solving the incompressible Euler equations. In a final postprocessing step, consistent divergence and continuity penalty terms are applied to weakly enforce the incompressibility constraint and normal continuity of the velocity field (Fehn et al. Reference Fehn, Wall and Kronbichler2018b, Reference Fehn, Heinz, Wall and Kronbichler2021):

(2.14) \begin{align} & ( {\boldsymbol{v}_h},{\boldsymbol{u}_h^{n+1}} )_{\varOmega_{e}} + ( {\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}_h}},{\tau_{D}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}_h^{n+1}}} )_{\varOmega_{e}} \Delta t_n \nonumber\\ &\quad + ( {\boldsymbol{v}_h\boldsymbol{\cdot} \boldsymbol{n}},{\tau_{C} (\boldsymbol{u}_h^{-} - \boldsymbol{u}_h^{+})^{n+1}\boldsymbol{\cdot} \boldsymbol{n}} )_{\partial\varOmega_{e}} \Delta t_n = ( {\boldsymbol{v}_h},{\hat{\hat{\hat{\boldsymbol{u}}}}_h} )_{\varOmega_{e}} . \end{align}

The above postprocessing step is specifically related to the type of spatial discretisation used in this work, i.e. the $L^2$-conforming DG method, and disappears in the continuous case. It was found that the divergence and normal-continuity penalty terms are crucial in terms of mass conservation and energy stability in order to obtain a method that is robust for large Reynolds numbers and coarse spatial resolutions, and we refer to Fehn et al. (Reference Fehn, Wall and Kronbichler2018b, Reference Fehn, Kronbichler, Lehrenfeld, Lube and Schroeder2019) for detailed numerical justification and validation of this approach. In this context, we note that alternatives to this stabilised approach exist, e.g. by using tailored finite element function spaces. For example, function spaces can be used that result in an $H(\text {div})$-conforming (normal-continuous) velocity field and that is divergence-free in every point of the computational domain. The present $L^2$-conforming approach does not fulfil these two properties exactly, but in a weak finite element sense. Definitions of the penalty parameters $\tau, \tau _{{D}}, \tau _{{C}}$ are given in Fehn et al. (Reference Fehn, Wall and Kronbichler2018b), where the default penalty factor $\zeta =1$ is used in the present work unless specified otherwise. In terms of the occurrence of finite-time singularities and anomalous energy dissipation, the present work makes use of the argument that – due to the weighted residual formulation – the present discretisation can be applied to problems which lack regularity and for which the differential form of the equations is no longer an appropriate description. With respect to the implementation of the method, integrals in the weak form are evaluated by means of Gaussian quadrature, with consistent integration according to the $3/2$ rule for the nonlinear convective term (also known as polynomial dealiasing). Since we consider uniform Cartesian meshes (with elements of size $h$ in all coordinate directions) in the present work, integrals are calculated exactly.

To obtain the size of the next time step within the adaptive time stepping scheme, we use the following local CFL condition evaluated at all quadrature points $q$ of element $e$:

(2.15)\begin{equation} \Delta t = \min_{e=1,\ldots,N_{{el}}} \left(\min_{q}\frac{Cr}{k^{1.5}} \left.\frac{h}{\Vert \boldsymbol{u}_h \Vert }\right\vert_{q,e} \right).\end{equation}

where we ensure stability of the time integration scheme by selecting a Courant number ${Cr} < {Cr}_{{crit}}$. Note that – due to the CFL condition – an increase in spatial resolution by reducing $h$ or increasing $k$ automatically implies an increase in temporal resolution through a reduced time step size.

The present incompressible flow solver is implemented in C++ and forms a core module of the ExaDG (High-Order Discontinuous Galerkin for the Exa-Scale) software, which is publicly available on github (see https://github.com/exadg/exadg). The ExaDG software primarily makes use of the open-source finite element library deal.II (Alzetta et al. Reference Alzetta2018; Arndt et al. Reference Arndt, Bangerth, Davydov, Heister, Heltai, Kronbichler, Maier, Pelteret, Turcksin and Wells2021). The computational efficiency of the present discretisation methods in terms of fast implementations and fast iterative solvers is discussed in Fehn, Wall & Kronbichler (Reference Fehn, Wall and Kronbichler2018a), Fehn et al. (Reference Fehn, Munch, Wall and Kronbichler2020), Kronbichler & Kormann (Reference Kronbichler and Kormann2019) and Arndt et al. (Reference Arndt, Fehn, Kanschat, Kormann, Kronbichler, Munch, Wall and Witte2020).