Smoothed Particle Hydrodynamics Physically Reconsidered - The Relation to Explicit Large Eddy Simulation and the Issue of Particle Duality

In this work we will identify a novel relation between Smoothed Particle Hydrodynamics (SPH) and explicit Large Eddy Simulation (LES) using a coarse-graining method from Non-Equilibrium Molecular Dynamics (NEMD). While the current literature points at the conclusion that characteristic SPH issues become restrictive for subsonic turbulent ﬂows, we see the potential to mitigate these SPH issues by explicit subﬁlter stress (SFS) modelling. We verify our theory by various simulations of homogeneous, isotropic turbulence (HIT) at Re = 10 4 and compare the results to a Direct Numerical Simulation (DNS) reported in Ref. 1 Although the simulations substantiate our theory, we see another issue arising, which is conceptually rooted in the particle itself, termed as Particle Duality . Finally, we conclude our work by acknowledging SPH as coarse-graining method for turbulent ﬂows, highlighting its capabilities and limitations.


I. INTRODUCTION
Since its first introduction for astrophysical flow problems by Lucy and Gingold & Monoghan in 1977 2,3 , the success of SPH as a viable Lagrangian method in the Computational Fluid Dynamics (CFD) community is undeniable.In the last decades there was a considerable research effort to increase the fundamental maturity of the method, summarized in several reviews [4][5][6][7][8][9] , which was in parallel accompanied by progress regarding applications of higher complexity, e.g.Ref. [10][11][12] .
One of the most fundamental problems of classical Lagrangian SPH is that it suffers from zeroth order errors, which result in a substantial amount of noise compared to grid based Eulerian methods 13,14 .Physically, this noise causes excessive dissipation [13][14][15][16] by numerically induced small scale vorticity 16 .Although it could already be hypothesized based on the work of Ellero et al. 15 that this might become a severe issue for subsonic turbulence, a rigorous and detailed analysis proving this fact for forced homogeneous, isotropic turbulence (HIT) was presented in the seminal work of Bauer & Springel 13 .Up to date, the shortcomings of SPH for subsonic turbulence as discussed by the authors persist, namely that large scale turbulent structures can be qualitatively captured but at comparably high computational cost taking alternative CFD methods into account.This is quite unsatisfactory given that turbulence is a key aspect of most fluid flows.
However, it might be argued that the results obtained by Bauer & Springel 13 are the consequence of a missing turbulence model and that they are only valid for underresolved Direct Numerical Simulations (uDNS).There were several publications on turbulence modelling in SPH [17][18][19][20][21][22][23] , but most of them either show a marginal improvement or are rather inconclusive for three dimensional subsonic turbulence.The latter can be attributed to the fact that the models are tested only with scarce validation runs, on setups which contain complex boundaries adding other SPH specific uncertainties on top of the actual turbulent flow 24 or are validated for two dimensional turbulence, which behaves qualitatively different 25 and where the use of usual turbulence closure models is unjustified, e.g.Ref. 26,27 .In our opinion, the most promising approach so far was presented in the pioneering work of Di Mascio et al. 22 and only recently extended by Antuono et al. 23 .In these works the authors explore SPH from a Large Eddy Simulation (LES) perspective, which represents a natural option as already noticed three decades ago 28 .Despite the fact that the derived SPH-LES approach with its various additional terms is an important step for SPH towards turbulent flows, both works do not evaluate the foundations of the classical LES subfilter stress (SFS) in a SPH framework, which is the central quantity in LES.In our opinion, the latter objective is vital, because the classical Lagrangian SPH features a turbulent kinetic energy deficit 13 , which questions the intention of introducing mostly dissipative SFS models from the beginning 29 .
The work presented in this paper focuses on connecting explicit LES and SPH with a coarse-graining method from Non-Equilibrium Molecular Dynamics (NEMD).Hence, this study can be viewed as a sequel of our recent publications 30,31 .We will demonstrate that the only additional term emerging from our theory is the SFS term as it is known from Eulerian LES methods.This is contrary to the statement of Di Mascio et al., concluding that a proper LES interpretation of SPH necessitates the consideration of additional SPH exclusive terms, though the authors ascertain that these terms play a minor role 22,23 .Consequently, we will be able to discuss the rationality of SFS models for the SPH simulation of subsonic turbulent flows.We are not considering further heuristic noise-mitigating techniques.Most importantly, our work is motivated by the following central question: Can resolved large scale structures profit from the reduction of SPH typical small scale noise by explicit use of SFS models?
In order to elaborate this hypothesis, this paper is structured as follows: We start with a short review of the main characteristics of subsonic HIT and how the large scale dynamics of such turbulent flows can entirely be described by coarse-graining regularization of the fluid dynamic balance equations 32 .This technique is most commonly known as LES.Then, we will relate this coarse-grained picture of arXiv:2206.10127v4[physics.flu-dyn]22 Nov 2022 subsonic turbulence to SPH.We constitute that SPH can be viewed as a Lagrangian quadrature technique for the governing equations of explicit LES and discuss the significant implications of this approach.To verify our theory, we will subsequently present various results of subsonic HIT simulations at Re = 10 4 and compare the results to a DNS solution reported in Ref. 1 .Finally, we will draw a conclusion on the rationality of SFS models in SPH.

II. LES: COARSE-GRAINED DYNAMICS OF SUBSONIC TURBULENT FLOWS
Despite the omnipresence and extraordinary beauty of turbulent flows, a comprehensive theory is still missing.However, it is agreed by the fluid dynamics community that the concept of the energy cascade is an important cornerstone of turbulence theory 32 .The cascade process was metaphorically described by Richardson 33 in 1922 for the first time, before it was quantified for incompressible HIT by Kolmogorov 34 , Obukhov 35 , Onsager 36 and Heisenberg 37 about 20 years later.For fully developed turbulence it was already demonstrated back then that, in a statistically averaged sense, a range of large scales exists in which the kinetic energy of velocity fluctuations of a specific wavenumber, namely E(k), is transferred from larger to smaller scales in the absence of viscous dissipation effects.This range is known as inertial range 32 and its scaling characteristics of serves as an important benchmark for CFD solvers to prove their capability to reproduce large scale dynamics of strongly subsonic turbulent flows, e.g.Refs. 1,13.In contrast, latest research activities on strongly subsonic HIT focus on the smallest scales of the cascade process, namely the dissipation range and beyond [38][39][40][41][42][43] .Among others, it has been rigorously argued by Sreenivasan & Yakhot 40 , with the aid of a novel anomalous scaling theory, that the actual smallest turbulent length scale with wavenumber k ∞ falls below the Kolmogorov scale with k η .Another example is the importance of thermal fluctuations beyond the classical dissipation range, i.e. k η < k < k m f p with k m f p denoting the wavenumber of the mean free path, eventually leading to a E(k) ∼ k 2 scaling [41][42][43] .Overall, all these insights can be vividly condensed in the turbulent kinetic energy spectrum as depicted in FIG. 1, highlighting the different turbulent regimes.Although all scales are of significant importance for a holistic view of HIT, based on the work of Bauer & Springel 13 we can already conclude that classical Lagrangian SPH struggles to directly resolve turbulent small scale dynamics.Hence, it seems more convenient to combine SPH with a coarsegrained model, which intrinsically focuses on the turbulent large scales with its inertial range characteristics (Eq.( 1)).This approach is more likely to be compatible, as the feedback of the smallest scales only has to be modelled and not resolved.A coarse-grained model of that kind can be derived by means of a technique from the Non-Equilibrium Molecular Dynamics (NEMD) community, as we demonstrated in Smaller structures: lg k 1. Schematic of the turbulent kinetic energy spectrum for fully developed HIT in the strongly subsonic regime, including the thermal range beyond the classic dissipation range according to Ref. 41 .
Ref. 31 .In the following we will shortly summarize the main ideas of this generalization of Hardy's theory from 1982 44 .
The key aspect of the Hardy theory is that it transfers arbitrary Lagrangian particles into coarse-grained particles by appropriate averaging, satisfying axiomatic conservation properties scale independently.The averaging is mathematically accomplished by the introduction of a normalized, symmetrical, positive, and monotonously decaying function W h with compact support supp{W h } ⊂ R 3 , in short kernel.It is assumed that the latter is spherical and its spatial extent is quantified by the scalar index h ∈ R + .Originally, as depicted in FIG. 2, the Hardy theory was used to link the dynamics of discrete molecules and individual fluid elements governed by their continuum balance equations, e.g.Navier-Stokes for Newtonian flows.Even more important for this work is the fact that the Hardy theory can be generalized to a continuous set of Lagrangian particles as well, i.e. fluid elements (FIG.2).This generalization leads to the governing equations of coarse-grained super fluid elements, which for h = const are completely equivalent to the governing equations of LES.The latter represent a deterministic fluid flow model, which inherently focuses on scales above the kernel size and is able to capture the turbulence cascade 32 .Thus, we interpret it as a physical description, which matches well with the SPH image for subsonic turbulence provided by Bauer & Springel 13 .
Although this LES perspective based on Lagrangian particles requires an additional mathematical effort compared to the common commutable filtering operation 45 , it reveals a striking similarity between LES and SPH.A relation between the latter was already noticed 30 years ago by Bicknell 28 , but we take the view that our perspective on LES vividly strengthens this point.As common feature, both methods perform a coarse-graining of Lagrangian particles by means of a kernel and we will use this property in Sec.III to argue that SPH can be reinterpreted as a Lagrangian quadrature technique of the governing equations of LES.
However, before we proceed with this objective, we will first concentrate on a fluid flow model, which is in general capable to describe isothermal, strongly subsonic HIT.There- fore, we select a Newtonian, barotropic fluid flow completely specified by its density, pressure and velocity field, namely ρ, p and v, and apply the Hardy theory to it as detailed in Ref. 31 .With x ∈ R 3 , we denote an individual kernel center position and its corresponding support is abbreviated as The equations Eqs.(2a), (2b) & (2c) represent the averaged continuity equation 31 , the averaged momentum transport equation and a linear barotropic equation of state (EOS) with ρ re f , p re f & K as constants.These describe the reference density of the strongly subsonic flow, the reference pressure and a stiffness constant, which are highly dependent on the problem.Their choice will be specified in Sec.IV.Furthermore, from Eq. (2a), the meaning of the overline notation for an arbitrary field f can be deduced.It describes a spatial average over a superfluid element V x of size ∼ h, namely with dy as volume differential of a Lagrangian fluid element.Moreover, the momentum transport equation Eq. (2b) em-ploys a density-weighted averaged velocity ṽ over V x as indicated by the tilde notation.Generally, this density-weighted average for a field f is termed Favre average 48,49 , although it was already suggested by Reynolds 50 in 1895.It is defined as The use of Favre averages is not mandatory but handy, as it circumvents correlation terms related to the density field 49 .
It is important to distinguish between quantities according to Eq. ( 3) & Eq. ( 4), which refer to super fluid elements in the LES framework, and fluid element quantities, which are simply noted without an overline or tilde.Additionally, as a consequence of the coarse-graining regularization of the balance equations by Eq. ( 3), an extra term div x [τ SFS ] (x,t) appears in Eq. (2b) 31,32 .It is the contribution from scales below V x to the momentum transport of super fluid elements.The SFS tensor τ SFS can be written as covariance tensor of the velocity field 31 τ SFS (x,t) := which is an interesting representation as the discretized version of Eq. ( 5) localizes flow subdomains, where SPH struggles with accurate approximations 30 .Due to its relevance for this study, we will elaborate on this in more detail in Sec.III B. Finally, it might be surprising that the averages on the right hand side of the transport equations in Eq. ( 2) are explicitly noted for each V x and not abbreviated by Eq. (3).By that we intend to emphasize that we follow the philosophy of explicit LES methods, e.g.Refs. 51,52.Contrary to the usual procedure in explicit LES, where the nonlinear convective term is explicitly filtered, the filter is explicitly applied to the right hand side of the transport equations.This is due to the Lagrangian perspective we take, in which the convective term is not directly considered but rather a consequence of the individual forces on the right hand side of Eq. (2b).

III. SPH AS A LAGRANGIAN QUADRATURE OF EXPLICIT LES
As explained in the last section, the governing equations of LES can generally be derived by coarse-graining of Lagrangian fluid elements using the Hardy theory from NEMD.This explains the conceptual similarity of LES and SPH, which becomes also evident from FIG. 2, being a reminder of how SPH is often vividly introduced, e.g. in the work of Price 6 .
The objective of this section is to argue that SPH should be generally viewed as a Lagrangian quadrature technique for the governing equations of explicit LES.This general fluid dynamic framework includes the kernel concept from the beginning and requires identical to SPH a joint limit for formal convergence 53 , in which the ratio of filter width to grid spacing ∆/∆l → ∞ and ∆ → 0 52 .In the following, we will present the resulting SPH model and discuss the implications of the explicit LES perspective.

A. The SPH-LES Model and its Implications
Decomposing the fluid domain into a finite number of Lagrangian SPH particles i ∈ {1, ..., N} that are connected to the kernel center positions, i.e. ∀i ∈ {1, ..., N} : x i = y i , one can derive the final SPH model.The discretization procedure of Eq. ( 2) is detailed in the Appendix A (Eqs. (A3), (A10) & (A16)).It is important to highlight that the SPH particles only have to be Lagrangian representatives of super fluid elements V x with the arbitrary length scale h instead of fluid elements.This is a significant difference to the usual SPH approach because traditionally SPH particles suffer from pseudo-Lagrangian behaviour at finite resolution 30,54 .For an individual particle i with j ∈ {1, ..., N ngb } neighbors the model reads and the particle trajectories follow from the kinematic condition dx i dt = ṽi .
Formally, the emerging system of Eqs. ( 6) & ( 7) is identical to the SPH discretization of the weakly-compressible Navier-Stokes equations (WCSPH) except for the SFS term div [τ SFS ] i in Eq. (6b).The latter is a direct consequence of the coarse-graining at the arbitrary kernel scale h, compensating for subkernel effects.Contrary, in traditional SPH, the choice of the scale h is merely a matter of convergence.We understand this as a physically convincing argument, going beyond empty formalities, to state that SPH should be understood as a Lagrangian quadrature technique intrinsically connected to explicit LES.Then, from this LES perspective, deficits introduced at the kernel scale for a specific choice of h could potentially be compensated by a proper modelling of the SFS tensor τ SFS in Eq. (2b).We believe that empirical evidence for this reconsideration is also given by the fact that already in the pioneering SPH works of Lucy and Gingold & Monoghan 2,3 artificial damping terms were used.These can be interpreted as the first SFS models accounting for subkernel deficiencies.Eventually, the reinterpretation of SPH as an intrinsic Lagrangian quadrature of explicit LES comes with two significant implications: 1. Implication: Ideally, the physical resolution of SPH is limited by the kernel scale of V x or more precisely the kernel diameter D K .Thus, SPH is unsuited as a DNS method.
2. Implication: Deficits introduced below the kernel scale of V x might be resolved by explicit consideration of the SFS term, from which structures above the kernel scale of V x could profit.
While the first implication can be easily understood and there is already empirical evidence proving it, e.g. the work of Bauer & Springel 13 , the second implication should be interpreted as a working hypothesis, which we will test by numerical experiments in Sec.V.However, before this, we will recapitulate the current knowledge about the origin of numerical dissipation in SPH in the next section and explain its relation to the SFS.

B. Numerical Dissipation and the Role of the SFS
In this paragraph, we will summarize how numerical dissipation in kernel-based particle methods like SPH emerges and how it can be localized by usage of Hardy theory again.Exemplary, we explain the coherencies with the aid of a SPH solution of the forced two dimensional turbulent Kolmogorov flow of Rivera et al. 55,56 .This was the main topic of our work in Ref. 30 , which subsequently is summarized.As the discussion will reveal, the current understanding of numerical dissipation in SPH is rather heuristic than rigorous compared to conventional grid-based CFD methods.The most important metrics for the analyzed Kolmogorov flow are depicted in FIG. 3 and non-dimensionalized with the absolute maximum value of the viewed snapshot.
Although the Kolmogorov flow can be generally reproduced by proper calibration with particle discretization methods in terms of the energy characteristics, the analysis shows that the dissipation rate is strongly overpredicted.From current understanding, the excessive numerical dissipation is rooted in increasing particle disorder at higher Reynolds number causing local SPH particle transport perpendicular to the main flow direction 15 .The resulting unphysical momentum transport prevails especially in zones with large velocity gradients and introduces artificial vorticity fluctuations 16 .For the considered problem, exemplary evidence is given by the noisy vorticity field in FIG. 3 (a).Since the turbulent flow is periodic, it can be analytically argued that the vorticity fluctuations ω will quadratically increase the averaged dissipation rate ε according to 57 with • V denoting a volume average.Moreover, based on Eq. ( 8), the vorticity fluctuations, causing numerical dissipation, are evidently linked to the kinetic energy spectrum E(k).An analysis of the extracted spectrum in FIG. 3 (b) for the Kolmogorov flow demonstrates that the scalings of the inverse and direct cascades known from the experiment can be matched above the kernel scale (k < 2π/D K ) 56 .Merely the saturation of E(k) below the kernel scale (k > 2π/D K ) deviates from the ideal form.Since dissipation takes place at small scales due to the ∼ k 2 weighting of E(k) in Eq. ( 8), it seems likely that this saturation represents the spectral signature of numerical dissipation.In accordance with FIG. 1, we will term this bottleneck as artificial thermal range or artificial thermalization.It should be emphasized that this spectral signature was also observed in other works, e.g.Refs. 13,22,23.
The upper paragraph should be understood as a heuristic description of numerical dissipation in kernel-based particle methods as SPH.However, it does not locally explain how numerical dissipation in the flow field takes place.Motivated by this deficit, we applied Hardy theory to a numerical particle set in Ref. 30 deriving a tensor, namely a discrete approximation of the SFS tensor, which was shown to serve this purpose.It reacts to large velocity gradients and particle noise, which are believed to be the root of the numerical dissipation.In Ref. 30 the tensor was termed as molecular stress, since we were not aware of its connection to explicit LES.Interestingly, in kernel-based particle methods this tensor can always be estimated, even if no LES perspective is employed.Considering the fact that two different spatial resolution scales exist, i.e. the kernel diameter D K and the particle size ∆l < D K , a Lagrangian quadrature of Eq. ( 5) gives the following estimate For the considered Kolmogorov flow, an evaluation of the Frobenius norm of the SFS tensor according to Eq. ( 9) clearly demonstrates that the SFS estimate fits into the former dissipation characteristics.A visual comparison between the vorticity (FIG.3 (a)) and the SFS tensor (FIG.3 (c)) reveals an unambiguous relationship, which is supported by the corresponding bivariate probability density function in FIG. 3 (d).
The latter describes a cone like structure indicating that high levels of absolute vorticity and SFS are connected 30 , as well as that the variance of the absolute vorticity increases with the SFS norm up to τ * SFS 0.4.Consequently, we conclude that the quantities ω , E(k) & τ SFS are evidently correlated with each other, containing different levels of information about the numerical dissipation dynamics.
With the former relations, the role of the SFS term div [τ SFS ] i in the discretized explicit LES equations (Eq.( 6)) becomes apparent.Since the SFS term behaves diffusive regarding the velocity field in a statistically averaged sense 32,58 , an explicit consideration of this term will attempt to locally homogenize the velocity field.Hence, a mitigation of the SFS norm according to Eq. ( 9) will be caused, which is expected to reduce the vorticity variance according to FIG. 3 (d).Then, based on Eq. ( 8), the artificial thermalization of the kinetic energy spectrum should be likewise reduced.Considering again that dissipation takes place at small scales due to the ∼ k 2 weighting of E(k) in Eq. ( 8), this could potentially enable a reduction of numerical dissipation from which large scale structures might profit.Thus, the main idea of this work is to replace numerical dissipation by an explicit, dissipative SFS model, which counteracts the numerically induced small scale fluctuations.
To verify these anticipated causalities, we deem the eddy viscosity concept in connection with Boussinesq's hypothesis 45,[58][59][60] for the modelling of the SFS term div [τ SFS ] i in Eq. (6b) as adequate.The following approaches will be utilized: • SMAG: This represents the classical Smagorinsky model discretized according to Eqs. (A20), (A21) and (A24).It is angular momentum conserving in the continuum limit 61 .
• SIGMA: This represents the superior σ -model of Nicoud et al. 62 discretized according to Eqs. (A20), (A21) and (A25).It is also angular momentum conserving in the continuum limit 61 but should overcome severe drawbacks of the Smagorinsky model, e.g.non-vanishing subfilter dissipation in laminar regions 60,62 .
• SMAG-MCG: This represents the classical Smagorinsky model, however, discretized in the Monoghan-Cleary-Gingold (MCG) form 63,64 .It is angular momentum conserving on the particle level as well.
It is of paramount importance that the SFS dissipation is introduced only on subkernel scales in order to guarantee a successful application of the explicit SFS model, eventually reducing the artificial thermalization.However, the SPH discretization requires non-local approximations, which might jeopardize this goal a priori.This can be vividly illustrated by the concept of Particle Duality, which results from the coarsegraining perspective.

C. Particle Duality and Numerical Dispersion
In order to understand the concept of Particle Duality, it is necessary to precisely define the terminology of explicit LES.According to Sec.II of this work, explicit LES is introduced as a general fluid dynamic framework, in which fluid elements are coarse-grained by an explicit kernel to so called super fluid elements.This is illustrated in the left part of the schematic in FIG. 4. From the schematic, an averaging over a fluid el-ement collective (grey particles) exactly determines the properties of a single super fluid element (red particle) with its specific kernel support.This corresponds to a truly explicit LES.However, in a SPH model the fluid element properties are unknown, which is synonymous to the closure problem of turbulence.This issue is resolved in a SPH framework by a direct substitution of the fluid elements (grey particles) by super fluid elements (red particles).Only then an averaging is performed to estimate the properties of a single super fluid element itself.As a consequence the SPH particles must represent super fluid element approximants and fluid element surrogates at the same time, which is what we term as Particle Duality.Practically, this implies that super fluid element approximants interact with each other, which are not direct neighbors but rather separated by some particles in between.This occurs as long as a the particles share the same kernel support and implicitly causes an increase of the effective interaction distance.Physically, however, the considered interaction is inadequate as the governing LES equations are a local model in terms of the super fluid element quantities.Thus, the Particle Duality as a manifestation of the non-locality introduced by the SPH discretization gives an picturesque description why the consideration of an explicit dissipative SFS model might fail to result in an improvement as anticipated in Sec.III B. Conceivably, the SFS model will not solely remove kinetic energy from the problematic artificial thermal range but also affect resolved scales larger than the kernel.
From the phenomenological issue described above, one might not exclusively expect problems in terms of the application of dissipative SFS models.It is also likely that numerical dispersion effects become relevant.Since the Particle Duality results in a decreased effective resolution, it should be anticipated that modes of a certain wavenumber k will be properly resolved only if the kernel wavenumber k kern exceeds k significantly, i.e. k k kern , and N ngb → ∞.Indeed generic SPH dispersion studies for the Euler equations demonstrate that the resulting dispersion errors depend on (k/k kern ), the kernel type, the number of neighbors N ngb inside the kernel, the reference density ρ re f and pressure level p re f , e.g.Ref. 65 .However, from these insights it can be hardly predicted how turbulent flows will be affected.A recent study of Yalla et al. 66 investigates the effect of numerical dispersion effects on the energy cascade in grid-based LES.The authors conclude that numerical dispersion can spoil the Galilean invariance in numerical simulations, finally inhibiting the energy transfer to the small, most dispersive modes, hence causing a degradation of the inertial range in the energy spectrum.Though the one-to-one transfer of these results to Lagrangian kernelbased particle methods as SPH is questionable, we will still test in Sec.V whether Galilean invariance holds.This might give an indication about the relevance of numerical dispersion errors in our study.
Having laid the foundations of SPH as a discretization intrinsically connected to explicit LES, it is now indispensable to answer our central question: Can resolved large scale structures profit from the reduction of SPH typical small scale noise by explicit use of SFS models?Therefore, a thorough investigation of the influence of explicit SFS models for a well- defined HIT problem are vital, which is why we proceed with the description of such in the next section.

IV. THE HIT PROBLEM
The HIT problem that will be subsequently investigated is the Taylor-Green flow 67,68 at Re = 10 4 presented in the work of Dairay et al. 1 (other values of Re are discussed in Appendix B).Their DNS solution will serve as reference for our study, as well as an second order accurate (space & time) solution of a standard Smagorinsky LES from a Finite Volume Methode (FVM).The latter was computed with OpenFOAM 7 on a Cartesian grid with 384 3 cells and closely matches the energy characteristics of the DNS.The initial velocity field in the triperiodic domain Ω := [0, 2π] 3 of the freely decaying flow is specified as from a Cartesian lattice arrangement, the particles were regularized into a stable configuration in corresponding pre-runs following the particle packing scheme of Colagrossi et al. 70 .Only then, the fields given by Eqs. ( 10) & ( 11) were mapped onto the particles.The initial velocity field magnitude for a N = 192 3 Case is depicted in FIG. 5 (a) highlighting two shear flow planes at z = π/2 & z = 3π/2 and the rotational direction at the plane z = 2π.In order to match the pressure field in the initial time step and avoid artificial dynamical effects beyond the one resulting from the initial particle configuration, a consistent mass distribution M 0,i = ρ 0,i ∆l 3 was imposed.It is illustrated in FIG. 5 (b).The density field ρ 0,i is given by the combination of the initial pressure field in Eq. ( 11) and the EOS in Eq. (6c).For the latter a reference density of ρ re f = 1 kg/m 3 , a stiffness constant K = ρ re f c 2 a = 25 Pa and a reference pressure p re f = K/4 = 6.25 Pa were chosen.The stiffness constant implies an artificial speed of sound of c a = 5 m/s, which corresponds to an initial Mach number Ma 0 = 0.2 and justifies to neglect ∇ • v based forces 47 .We want to emphasize that different values of p re f were tested, namely p re f ∈ {K/10, K/4, K/2, K}, however, the value p re f = K/4 yielded the best trade-off between stability and numerical dissipation.The results in the following were all computed using a Wendland C4 kernel and a kernel diameter of D K = 8∆l resulting in N ngb ≈ 250.Other comparative simulations were conducted with D K ∈ {4, 6, 8} ∆l (see Appendix C) and a quintic B-spline kernel but only the chosen configuration was capable to provide reasonable numerical convergence with increasing N avoiding pairing instabilities at the same time 53,65 .If not mentioned otherwise, the before described strategy will be used for all runs in this study.
To facilitate the discussion of explicit SFS models for the SPH method, the results in Sec.V will follow the subsequent argumentation sequence: First, we will demonstrate that the system of Eqs. ( 6) & ( 7), shows a convergent tendency in the numerical sense (TABLE I: Case 1-4).These runs cor- respond to usual WCSPH simulations or, from our coarsegraining perspective in Sec.III, to a Lagrangian quadrature of the explicit LES equations (Eq.( 2)) without an explicit SFS model.Hence, the SFS will only be implicitly considered.The Cases 5-12 will represent runs in which the SFS is explicitly added by means of eddy viscosity approaches (Sec.III B).To proof that our conclusions are generally valid for for eddy viscosity approaches, not only the standard Smagorinsky SMAG (TABLE I: Case 5-8) model will be evaluated but also the σ -model SIGMA (TABLE I: Case 9), which should ensure vanishing subfilter dissipation in the initial laminar phase of the HIT problem.Moreover, the Smagorinsky model will be additionally combined with the angular momentum conserving MCG form in Eq. (A26) called SMAG-MCG (TABLE I: Case 10), demonstrating the robustness of our observations.In order to refute that the obtained results are the consequence of a wrong calibration of the model constant in the eddy viscosity models 58 , e.g.C S for Smagorinsky, we will further present two SMAG runs in which the standard value C S = 0.15 was either halved or doubled (TABLE I: Case 11-12).All eddy viscosity runs were performed with a filter width ∆ = D K /2 = R K being equivalent to the kernel radius R K and in accordance with Rennehen 29 .From our explicit LES perspective in Sec.III the most consistent choice would correspond to ∆ = D K but some tests led to the conclusion that only the overall dissipation is enhanced without any further physical improvements.Interestingly, for the problem considered, the choice ∆ = ∆l had a nearly negligible effect on our solutions.We interpret this as evidence in favor of the intrinsic connection between explicit LES and SPH, in which the particles should approximate LES super fluid elements and not fluid elements itself.This will be confirmed in Sec.V C. The final Case 13 is identical to Case 3 except for the fact that in Eq. ( 10) a constant velocity is superimposed on the third component, namely v 0,z (x, y, z) = π.By that we intend to test, according to Sec.III C, whether numerical dispersion might be of relevance for our investigation.If so, the Galilean invariance of the model could be spoiled, finally altering the shape of the spectrum obtained in Case 3 according to Yalla et al. 66 .We want to note that the choice v 0,z (x, y, z) = π is an order of magnitude larger than the velocity fluctuations at time t = 14 s, on which our analysis mostly focuses on.This implies a turbulence intensity of O(Tu) ≈ 0.1.Moreover, the maximum of the velocity magnitude in Case 13 corresponds to v max = √ 1 + π 2 ≈ 3.3 m/s.To ensure initial Mach similarity with Case 1, a speed of sound c a = 16.5 m/s is required.This also alters the initial mass distribution at the beginning according to the strategy above and the choice of the reference pressure p re f = ρ re f c 2 a /4 ≈ 68 Pa to ensure stability.To evaluate the quality of the results, different metrics will be invoked.On the one hand the assessment of the overall dissipation inside the domain will be based on the density weighted averaged kinetic energy e v (t) := 1 2 This metric is in accordance with the definition of Dairay et al. 1 except for the density weighting.It only has a minor influence as we could verify, but should be included for consistency with the weak compressibilty approach.On the other hand the overall dissipation will be assessed by the corresponding averaged dissipation rate, which can be computed from a finite difference approach for sufficient temporal sampling 1 using the relation Furthermore and most importantly, we will compute the kinetic energy spectra E(k) at time t = 14 s where HIT with the characteristic inertial range scaling of Eq. ( 1) should be present up to a wavenumber of k DNS ≈ 50 1/m, according to Ref. 1 .Therefore, we employ the nearest neighbor sampling technique of Bauer & Springel 13 on a Cartesian grid with ∆l/2 in combination with the method of Durran et al. 71 .This methodology is kinetic energy conserving or in other words satisfies the discrete Parseval relation.Due to the Nyquist criterion, spectra will only be presented up to wavenumbers corresponding to 2∆l.The interpolation method of Shi et al. 72 will not be considered as is unclear whether it might introduce smoothing in the artificial thermal range, which we want to avoid.Since the values e v (t = 14 s) for the different cases in TABLE I can significantly differ, we will normalize the corresponding spectra with the product e v (t = 14 s)L c and L c = 1 m for all runs to enable a relative comparison.Observations in spectral space will further be related to physical space by means of the Frobenius norm of Eq. ( 9), namely ||τ SFS || F , the backward finite-time Lyapunov exponent (FTLE) 73,74 in the time range [11, 14] s and the vorticity component While the first is indicative for small scale structures (see Sec. III B), the other will be used to assess the large scale structures.Additionally, we introduce a signal to noise (SNR) metric for the kinetic energy spectra defined by the ratio of energy above the kernel scale (with kernel wavenumber k kern ) in relation to the overall energy, namely It is important to stress that the SNR metric is only indicative for the reduction of the artificial thermalization but gives no insight about the solution quality above the the kernel scale.All computations in this work were performed for a time range of I T = [0; 15] s with the in-house SPH code turboSPH.The latter was developed for the prediction of primary atomization.For details, please refer to the work of Chaussonnet et al. 12 .

A. Numerical Convergence of the WCSPH Scheme
We will first start with a qualitative discussion of the numerical convergence of the discretized LES equations resulting from the Lagrangian quadrature (Eq.( 6)).For now the SFS term is neglected, formally leading to the standard WCSPH discretization of the Navier-Stokes equations (TABLE I: Case 1-4).Hence, these results will serve as reference to evaluate the effect of explicit SFS models in SPH.
From FIG. 6 (a), which illustrates the temporal evolution of the averaged kinetic energy e v (Eq.( 12)) for increasing particle counts N (darker colors), one can apparently conclude that the metric is numerically converging towards the DNS reference solution of Dairay et al. 1 (solid black line).However, even for N = 512 3 ≈ 130 Mio.particles, a significant gap remains compared to the DNS in the interval t ∈ [2.5; 10] s.The reason for this gap can be understood from the averaged dissipation rate ε t (Eq.( 13)) as depicted in FIG.6 (c).Especially in the initial timeframe t ∈ [0; 2.5] s the vanishing dissipation rate of the DNS solution is strongly overestimated by the WC-SPH scheme.This is probably linked to the strong anisotropic particle rearrangement of the initially laminar vortex configuration, causing numerical dissipation effects.Interestingly, although the highest resolution contains 64 times more particles than the lowest resolution, the curves only slowly approach the vanishing DNS level.The overestimation of the dissipation rate continues until t ≈ 7.5 s, where the blue WCSPH lines cross the black DNS line, and passes over to a systematic underestimation until the end of the simulation, except for the N = 512 3 case, which overestimates the dissipation rates again for t 12.5 s.Nevertheless, this global dissipation characteristic leads to the consequence that the energy levels in FIG.6 (a) approach the DNS solution again for t 7.5 s.Despite these quantitative deviations, it should be highlighted that the qualitative agreement of the dissipation rates ε t in FIG.6 (c) is reasonable.For all N the temporal occurrence of the dissipation peak at t ≈ 9 s is matched and for increasing N the formation of the second local dissipation peak at t ≈ 11 s is also evident.Moreover, it must be positively emphasized that after the first dissipation peak the dissipation rates are consistent with theoretical predictions 75 , which state that e v (t) ∼ t −1.2 .This is illustrated in FIG.6 (a) by the two black dashed lines, which only differ in the time shift parameter t 0 .The latter can be understood as an indicator for the begin of the decay process and is consequently smaller for the WCSPH runs than for the DNS.
Having verified that the global dynamics of the kinetic energy is reasonably well approximated by the WCSPH scheme, it remains to clarify whether HIT with the characteristic inertial range scaling in Eq. ( 1) develops after the dissipation peak at t ≈ 9 s.Therefore, we consider the normalized kinetic energy spectra as explained in Sec.IV for t = 14 s, which is in accordance with the work of Dairay et al. 1 .From their work an inertial range scaling should prevail up to a wavenumber of k DNS ≈ 50 1/m.The computed spectra (TABLE I: Case 1-4) are visualized in FIG.6 (e).For better orientation, the latter also contains a solid black line representing the k −5/3 scaling and a dashed black line with a stronger k −4 scaling.Moreover, the diagram includes the integral scale of the HIT problem (fawn dashed line) and the kernel scale of the simulation run with the highest particle count N (fawn solid line).The corresponding SNR values (Eq.( 14)) are also listed.Indeed it can be observed that the utilized WCSPH scheme without SFS model is able to recover a significant amount of the inertial range scaling for increasing N.While for the lowest N the scaling is only evident in the range k ∈ [4; 7] 1/m, the scaling range for the highest N covers nearly an order of magnitude in wavenumber, namely k ∈ [4; 32] 1/m.Hence, we can infer that WCSPH without SFS model is generally able to capture subsonic HIT, though at significant expense of N = 512 2 and N ngb ≈ 250.This is consistent with the convergence properties described by Zhu et al. 53 .The artificial thermalization as an anticipated origin of numerical dissipation is simultaneously reduced with increasing N, which is confirmed by an increase of the SNR values in FIG.6 (e).Nevertheless, even for the highest resolution, the artificial thermalization persists and likewise an approximate k −4 scaling observable in all simulations that sets in well above the kernel scale.The latter connects the resolved inertial ranges with the artificial thermal ranges and represents a SPH characteristic kinetic energy deficit also observed by other authors, e.g. in Ref. 13,29 .It clearly demonstrates that the WCSPH method produces coarse-grained LES solutions from the beginning and in accordance with the 1.Implication in Sec.III A. However, the explicit kernel fails to realize a sharp cutoff close to k kern , which could be due to the issue of Particle Duality resulting in a decreased effective resolution.This is especially evident in FIG. 7, in which the spectra of the DNS, the FVM LES and the WCSPH simulation with the highest resolution (TABLE I: Case   With these observations, the aim of the next section is to investigate whether the explicit consideration of a dissipative SFS model can reduce the thermalization in favor of scales larger than the kernel.This should not only lead to an improvement of the inertial range approximation but also in terms of the global dynamics.Consequently, we will seek for evidence for the 2. Implication in Sec.III A following the argumentation line of Sec.III B. We start with the SMAG model.

B. Effect of the Smagorinsky Model
Although the usage of a consistent LES model with an explicit SFS model should lead to an overall improvement of the solution, the investigation with the static SMAG model undermines this positive expectation (  6 (d), the consideration of the SMAG model only deteriorates the kinetic energy balance.This becomes evident from the temporal evolution of the averaged kinetic energy in FIG.6 (b).For the whole time range the level of the averaged kinetic energy is reduced compared to the reference WCSPH solution in FIG.6 (a).It holds true for all particle counts N analyzed.The deterioration of the solution must also be apparent in the corresponding dissipation rates in FIG.6 (d), which is indeed the case.A comparison of the WCSPH reference cases in FIG.6 (c) with the SMAG cases in FIG.6 (d) demonstrates in particular in the initial timeframe (t < 7.5 s), up to the point where the blue dissipation lines cross the black DNS line, that the excessive dissipation of the WCSPH solution is significantly enhanced by the SMAG model.For the lowest particle count of N = 128 3 it even results in a noticeable qualitative shift of the first dissipation peak from t ≈ 9 s to t ≈ 6 s.For the remaining particle counts the position of the first dissipation peak in FIG.6 (d) is quite robust, though a slight shift towards earlier times is perceptible.It is interesting to note that in this initial timeframe the vortex system is still in transition to HIT 1,68,69 , which prompts the eventuality that this deterioration of the solution might be linked to the drawbacks of the Smagorinsky model as explained in Refs. 58,60.In order to refute this eventuality, we will also present in Sec.V C an investigation employing the superior SIGMA model (Appendix A 4) .
From the discouraging global kinetic energy one might be tempted to that the described link between explicit LES and SPH in Sec.III might be flawed.However, the kinetic energy spectra at t = 14 s in FIG.6 (f) reveal that the anticipated causalities presented in Sec.III are correct.Most importantly the comparison of the spectra in FIG.6 (e) & FIG.6 (f) demonstrates the reduction of the artificial thermalization for a specific N, although not very effectively.The relative energy content below the kernel scales of the individual cases is mitigated by the SMAG model, which is consistent with a dissipative SFS model.This is also confirmed by an increase of the SNR metric (Eq.( 14)) for a given N, which are additionally listed in FIG.6 (e For the remainder of this section we will discuss whether scales larger than the kernel can profit from the mitigated artificial thermalization.The answer to this is ambivalent and we start with the positive aspects.As a matter of fact, the relative energy content of the first wavenumber shells (k ≤ 7 1/m) increases as the comparison of FIG. 6 (e) & 6 (f) demonstrates.Although the differences might seem extraneous, we want to emphasize that the plots are double logarithmic.Interestingly, this improved spectral signature can also be linked to physical space for particle counts N ≤ 256 3 by means of the backward finite-time Lyapunov exponent (FTLE) in the time range [11, 14]  x = π & y = π.This mirror symmetry is a characteristic of the vortex systems 68 and a vivid prove that a reduction of the artifical thermalization can positively influence the large scale coherent motion.As a sidenote, symmetry breaking of a similar kind was also observed in molecular approximations of the Taylor-Green system in the work of Gallis et al. 39 , albeit caused by physical thermal fluctuations.
Unfortunately, it must also be accentuated that the symmetry restoration is not discernible anymore for the highest resolution SPH runs.This is depicted in FIG. 10 (b) & (c).There, the vorticity field component ω z at z = π & t = 14 of the highest resolution SPH runs without and with SMAG model (TABLE I: Case 4 & 8) are illustrated.Obviously, the quantity behaves very similar to the FTLE field in FIG. 9.The lost symmetry restoration property is quite disadvantageous given that the N = 512 3 cases are the only runs, which reproduce a significant amount of inertial range scaling.Although the qualitative large scale pattern is matched in comparison with the symmetric FVM field in FIG. 10 (a), the SPH solutions evidently tend to smear out sharp features of the large scale vortices.This accordingly results in smaller vorticity values as the colorbars in FIG. 10 illustrate.
Despite the slight improvements described before, the overall effect of the SMAG model is rather malicious.The deterioration of the global dynamics in terms of the averaged kinetic energy in FIG.6 (b) and the corresponding dissipation rate in FIG.6 (d) strongly indicate that the negative aspects outweigh the positive ones.This is also reflected in the energy spectra in FIG.6 (f

C. Other Dissipative SFS Models and Model Calibration
In this part a comparison between the WCSPH solution without SFS model and with SMAG, SIGMA & SMAG-MCG model will be presented, as well as a sensitivity study of the Smagorinsky constant C S (TABLE I: Case 3 & 9-12).As the observations seem to be independent from the particle count N, only the N = 256 3 runs will compared.The results are depicted in FIG.11.All in all, the observations are very similar to those presented in Sec.V B. We will start with the influence of the explicit SFS models (TABLE I: Case 3, 7 & 9-10).
the superior SIGMA model nor the discrete angular momentum conserving SMAG-MCG model lead to an improvements in terms of the kinetic energy characteristics.In fact, both models are even worse.As depicted in FIG.11 (a), the averaged kinetic energy levels of both variants are slightly below the SMAG solution.The global dissipation rates in FIG.11  (c) confirm these results.Accordingly, this is also reflected by the spectra in FIG.11 (e).This is surprising for different reasons.For the SIGMA model a vanishing dissipation in the initial laminar phase by the nature of the SFS model 60,62 would be expected.However, as depicted in FIG.11 (c), the overall dissipation is increased from the beginning.One might be tempted to concluded that this is related to the missing zero order consistency of the SPH-LES model in Eq. ( 6), prohibiting the flow discrimination required for the SIGMA model.Nonetheless, similar observation were made using high-order Eulerian grid based schemes utilizing the same SFS model for Taylor-Green flows 76,77 .Consequently, this indicates that the Taylor-Green flow is a challenging problem for the SIGMA model independent of the numerical discretization scheme.The kinetic energy spectrum in FIG.11 (e So far our results can be criticized in terms of the fact that the utilized dissipative SFS models might be wrongly calibrated for the considered problem and discretization method 58 .For instance in the Smagorinsky model in Eq. (A24), one has the freedom to choose the constant C S and the filter width ∆.Since we set ∆ = D K /2 in Sec.IV, subsequently the influence of C S is analyzed (TABLE I: Case  11-12).Compared to our reference, the value of C S is either halved or doubled.We want to note that this is identical to the situation in which C S = 0.15 is fixed and ∆ changed accordingly.The results are depicted in the second column of FIG.11  in FIG.11 (b), which is also logically reflected by the dissipation rates in FIG.11 (d).For the highest value C S = 0.03, even an artificial dissipation rate plateau within the range t ∈ [6; 9] s is created.The dynamics increasingly departs from the case without SFS model for larger C S (or ∆).Same holds for the energy spectra in FIG.11 (f).It is interesting to note that the artificial thermalization is only slightly reduced compared to the range where the deficient energy scaling of k −4 prevails.This shows again the non-local character induced by the issue of Particle Duality.It becomes more pronounced with larger C S and incrementally represses the inertial range, showing that the best choice corresponds to C S → 0 (or ∆ → 0).
Clearly, from all these observations, it can be concluded that WCSPH without explicit SFS generates coarse-grained solution from the beginning, however at significant cost.This is likely due to the Particle Duality issue, which phenomenologically indicates that the method suffers from an decreased effective resolution.The intrinsic property also seems to prohibit improvements by means of standard dissipative SFS models.Although the models mitigate the artificial thermalization, which is believed to be the origin of numerical dissipation and the largest resolved scales can slightly profit from it, the SFS model dominantly remove kinetic energy from scales larger than the kernel.These scales are already badly resolved in terms of the energy cascade, which let us confirm the statement of Rennehen that dissipative SFS models overall degrade the SPH solution 29 .
Having verified that the best SPH-LES model relinquishes the explicit usage of standard dissipative SFS models, we finally proceed with an investigation of numerical dispersion effects for the pure WCSPH scheme.

D. Numerical Dispersion and Galilean Invariance
In this short paragraph, we will discuss whether numerical dispersion errors in the WCSPH scheme (Eq.( 6)) might influence the turbulent dynamics.Therefore, based on the idea of Yalla et al. 66 , we investigate whether our overall model is Galilean invariant.If the energy cascade is significantly altered, this could indicate that dispersion errors inhibit the transfer to the smallest scales.
However, before we compare the convected Case 13 with v z V = π and the quiescent Case 3, we will first study how Case 3 is influenced by the increase of p re f , which results from the requirement to keep initial Mach similarity (see Sec IV).The case is termed as "p-Boost".Evidently, as illustrated in FIG. 12, the global dissipation is heavily increased.Even a short initial phase is created in which kinetic energy is injected into the flow.This could be the consequence of a weak compressibilty effect.The overall increase of dissipation is comprehensible, considering that the value p re f is adjusting the magnitude of a force, which is the consequence of zero order errors only 6,70,78 .Its corresponding acceleration reads The dissipation of the convected case with v 0,z = π is strongly increased, however, this is only the consequence of the increase of p re f as the case p-Boost with v 0,z = 0 demonstrates.The lines for v 0,z = π and p-Boost perfectly overlap each other.
and physically causes an artificial momentum transport of particles that avoids the creation of void spaces.It ensures the stability of the method due to an inherent particle regularization, which reacts to irregular distributions, where the ideal condition ∑ N ngb j=1 ∇W h,i j V j = 0 is violated.It can be anticipated that this artificial transport will also take place perpendicular to the local main flow directions, causing numerical dissipation according to the ideas presented in Sec.III B. However, from the generic SPH dispersion study of Dehnen & Aly 65 , it is certain that not only numerical dissipation with p re f is increased but also the influence of dispersion errors.This is reflected in the spectra in FIG. 13.Beyond the increase of the artificial thermalization, which likely reflects the increase of numerical dissipation, the inertial range is even stronger suppressed in the p-Boost case.The deficient SPH scaling k −4 profits from increased p re f value.This could be an indication that dispersion errors inhibit the transfer of energy to the smallest scales, spoiling the development of an appropriate inertial range similiar to the study of Yalla et al. 66 .
Ensuing from the p-Boost case, the results of the convected Case 13 with v z V = π are comprehensible.As depicted in FIG. 12 & 13, the lines for v 0,z = π and p-Boost perfectly overlap each other.This shows that the WCSPH scheme is truly Galilean invariant, though only if the mean convective velocity is changed.Practically, the superposition of a mean flow also requires an adjustment of p re f to ensure stability, which evidently will spoil the Galilean invariance.The zero order term in Eq. (15) can drastically increase numerical dissipation and also inhibit the development of an inertial range due to the increase of dispersion errors.Consequently, the results obtained in our work might be rather positive compared 13. Kinetic energy spectra at t = 14 s obtained from the analysis of Galilean invariance.The spectrum of the convected case with v 0,z = π is strongly altered, however, this is only the consequence of the increase of p re f as the case p-Boost with v 0,z = 0 demonstrates.The lines for v 0,z = π and p-Boost perfectly overlap each other.
to cases in which a mean convection is superimposed onto the turbulent flow.It is certain that this interplay of numerical dissipation and dispersion caused by p re f requires a more thorough study, but from the selected tests analyzed in this section the first observations seem reasonable.

VI. CONCLUSION
Summarizing, the contributions of this work are numerous and especially important for WCSPH simulations trying to capture subsonic turbulent flows.Our main goal was to argue in Sec.II and Sec.III that, based on the Hardy theory from NEMD, SPH should be viewed as a non-local Lagrangian quadrature procedure intrinsically related to explicit LES.This implies on the one hand that subsonic turbulence captured by SPH will be correctly represented, in the best case, up to the kernel scale but also at significant cost, taking the convergence characteristics into account 9,53 .This is consistent with observations in the literature 13,53 , empirically supporting our LES perspective on SPH.On the other hand, it paves a potential way to mitigate SPH characteristic shortcomings by explicit consideration of a dissipative SFS model as explained in Sec.III B. The main idea followed herein was to replace the numerical dissipation, emerging from the artificial thermalization, by an explicit, dissipative SFS model.In order to test the hypothesis that a reduction of the artificial thermal range results in an improvement of the kinetic energy content of scales larger than the kernel, several simulations of freely decaying HIT at Re = 10 4 in accordance to Dairay et al. 1 were conducted and analyzed.However, it must be stated from the results presented in Sec.V that the explicit SFS model only leads to an marginal improvement of the largest coherent structures.This was vividly reflected by the symmetry restoration of the vortex system in the computed FTLE fields.Concerning the inertial range dynamics, the dissipative SFS models merely remove kinetic energy where SPH is already characterized by a spectral energy deficit.Eventually, it deteriorates the overall solution outweighing the positive effects.This is rooted in the non-local character of the Lagrangian quadrature, which can be explained by the concept of Particle Duality.The latter states that the SPH particles must simultaneously represent super fluid element approximants and fluid element surrogates at the same time, causing an non-physical increase of the effective particle interaction distance.Finally, our work allows to confirm Rennehen's expectation that explicit SFS models in a SPH framework only degrade the quality of the approximation for subsonic turbulent flows 29 and from our current understanding they should be disregarded.It seems to be the case that the excessive numerical dissipation mechanisms of SPH 15,16,30 outperform the explicit dissipative SFS models considered in this work.However, only adapted SFS models, originally developed for gridbased Eulerian methods, were tested.This possibly indicates that SPH native SFS models are needed that consider SPH specific characteristics.Furthermore, our work obviously shows that the understanding of numerical dissipation and dispersion in kernelbased methods like SPH and their influence on turbulent dynamics is very heuristic compared to conventional Eulerian grid-based methods.This lacking rigor must be definitely addressed in future works.As next step it would be interesting to study whether the issue of Particle Duality in terms of explicit dissipative SFS model can be either circumvented by higher-order schemes or by the use of the SPH-LES scheme developed by Di Mascio et al. 22 and Antuono et al. 23 .Although we currently believe that the issue of Particel Duality is a conceptual problem of SPH, there may be the chance that cross-effects between different modelling terms can restore the actual goal of standard explicit dissipative SFS models.a resolution of 2048 3 elements was required.Contrary, our results demonstrate that the SPH inertial ranges are controlled by numerical dissipation and dispersion effects representing an implicit SFS model.Therefore, we do not expect significant change with different Re as long as the kernel scale lies within the inertial range.This can be mathematically proved using the Chauchy-Schwartz inequality 32 and is confirmed in FIG. 14 for Case 3 in TABLE I with halved and doubled Re.For the given resolution the inertial range scaling prevails in the same wavenumber range.However, it is interesting to note that higher viscosities are beneficial in terms of the energy levels of the artificial thermalization.As a consequence, scales larger than the kernel relatively contain slightly more energy at lower Re, although this effect is evidently not very pronounced.

Appendix C: Influence of Kernel Size
In the main part of our work we only discuss cases with a kernel diameter of D K = 8∆l or radius of R K = 4∆l.This is motivated by the observation that adequate numerical convergence can hardly be observed for the given problem with smaller values, highlighting the paramount importance of N ngb .To give an idea, the temporal evolution of the density weighted averaged kinetic energy e v is depicted in FIG. 15 for smaller R K but otherwise same conditions as in Case 3 in TABLE I. Apparently, for a given N, the dissipation rate is strongly reduced with increasing N ngb .Likewise are the artificial compressibilty effects at the beginning of the simulation, which lead to a production of kinetic energy.In spectral space the effect is similarly pronounced, as shown in FIG.16.Only the R K = 4∆l case shows a pronounced inertial range scaling for N = 256 3 particles.These observations are compliant with the numerical convergence behaviour of SPH presented by Zhu et al. 53

FIG. 3 .
FIG. 3. Metrics of a SPH solution for a two dimensional turbulent Kolmogorov flow according to Rivera et al. 55,56 .(a) Snapshot of the nondimensional, noisy vorticity field, (b) Kinetic energy spectrum, (c) Snapshot of the nondimensional Frobenius norm of the SFS tensor and (d) Nondimensional bivariate probability density of the Frobenius norm of the SFS tensor (abscissa) and vorticity (ordinate).

re f v 2 0 16 ( 2 +
FIG. 5. Initialization of the Taylor-Green flow at Re = 10 4 with N = 192 3 particles.(a) Velocity magnitude field with highlighted shear flow planes and rotational direction for z = 2π.(b) Mass distribution corresponding to the initial pressure field.

FIG. 6 .
FIG. 6.Comparison of quantitative metrics for different particle counts N. The first column (a), (c) & (e) represents the WCSPH solutions without SFS model and the second column (b), (d) & (f) represents the solution with SMAG SFS model.The color coding is explained in (a).(a) & (b) Temporal evolution of the density weighted averaged kinetic energy e v .(c) & (d) Temporal evolution of the averaged dissipation rate ε t .(e) & (f) Kinetic energy spectra at t = 14 s.
4) are compared.Whereas the inertial range characteristics of the DNS (grey) and the FVM LES (red) are in perfect agreement and deviations are only evident beyond k DNS ≈ 50 1/m, the WCSPH (blue) solution fails to reproduce the inertial range up to k DNS .The spectrum is already damped at k SPH = 32 1/m = k kern /2, although k kern = 64 1/m exceeds the maximum inertial range wavenumber k DNS of the reference DNS.Compared to the FVM LES with only 384 3 cells, versus N = 512 3 particles and N ngb ≈ 250 neighbors, this result is rather disappointing.

FIG. 7 .
FIG.7.Comparison of kinetic energy spectra at t = 14 s.DNS of Dairay et al.1 is depicted in grey, the FVM LES in red and the WCSPH with the highest resolution in blue.

FIG. 8 .
FIG. 8. (a) Frobenius norm of the estimated SFS tensor without explicit SFS model and (b) with SMAG model.(c) Turbulent dynamic eddy viscosity resulting from the SMAG model.All distributions refer to the x = π plane at t = 14 s for N = 256 3 .
) & FIG.6 (f).Moreover, this reduction of the thermalization correlates with a reduction of the Frobenius norm of ||τ SFS || F , as anticipated in Sec.III B and illustrated in FIG. 8. There, the estimated ||τ SFS || F distribution is depicted for N = 256 3 at the plane x = π for t = 14 s, corresponding to the time instance of the energy spectra.Without explicit SFS (FIG.8 (a)) and with SMAG model (FIG.8 (b)) high values of ||τ SFS || F prevail in the shear zones of the flow (cp.FIG. 5).However, as a consequence of the SMAG model, the absolute SFS values are finally reduced.This damping is associated to the turbulent dynamic eddy viscosity η t distribution (FIG.8 (c)), which evidently correlates with the ||τ SFS || F distribution (FIG.8 (a) & FIG. 8 (b)).By construction of Eq. (A24), this is plausible as the turbulent eddy viscosity reacts to high values of spatial velocity gradients prevailing in the shear flow planes.
) further demonstrates that the qualitative effect of the SIGMA model and the SMAG model are similar.Compared to the SMAG model, the artificial thermalization below the kernel scale is slightly reduced, from which the first wavenumber shells (k ≤ 7 1/m) again slightly profit.This is also reflected by the FTLE field in FIG.9 (c), which compared to the WCSPH solution in FIG. 9 (a) is less tattered and, moreover, approximately mirror symmetric at the midplanes.However and most importantly, the dissipative SIGMA model still suffers from the issue of Particle Duality.A removal of kinetic energy from scales in the range k > 7 1/m and k < k kern is still present and even leads to an intensification of the observed deficient SPH energy scaling with E(k) ∼ k −n , n > 4. The SMAG-MCG model results tell a nearly identical story.Compared to the SMAG model the overall dissipation is enhanced.This becomes evident in FIG.11 (a) & FIG.11 (c).It is surprising as one would intuitively expect a general improvement related to the restoration of the angular momentum conservation property.Instead, the considered problem demonstrates that this comes at a certain cost.The kinetic energy spectrum in FIG.11 (e) once more confirms the already noted observations.The SMAG-MCG model is characterized by the strongest reduction of the artificial thermalization, which repeatedly has a positive effect on the first wavenumber shells (k ≤ 7 1/m) in the spectra, as well as the FTLE field in FIG. 9 (d).Nevertheless, the issue of Particle Duality for the SMAG-MCG model in the range k > 7 1/m and k < k kern is yet evident.Compared to the SMAG model, the angular momentum conservation property of the MCG form intensifies the observed deficient SPH scaling E(k) ∼ k −n , n > 4.
FIG. 11.Comparison of quantitative metrics for different N = 256 3 runs.The first column (a), (c) & (e) represents the influence of different explicit SFS models and the second column (b), (d) & (f) represents a sensitivity study of the Smagorinsky constant C S for the SMAG model.color coding is explained either in (a) or (b).(a) & (b) Temporal evolution of density weighted averaged kinetic energy e v .(c) & (d) Temporal evolution of the averaged dissipation rate ε t .(e) & (f) Kinetic energy spectra at t = 14 s.

)FIG. 12 .
FIG.12.Temporal evolution of the density weighted averaged kinetic energy e v obtained from the analysis of Galilean invariance.The dissipation of the convected case with v 0,z = π is strongly increased, however, this is only the consequence of the increase of p re f as the case p-Boost with v 0,z = 0 demonstrates.The lines for v 0,z = π and p-Boost perfectly overlap each other.

FIG. 15 .
FIG. 15.Temporal evolution of the density weighted averaged kinetic energy e v for different R K with N = 256 3 . .
Schematic of consecutive application of Hardy theory to different particle groups.Using a spherical kernel with size ∼ h, smaller particles can be transferred into coarse-grained particles by appropriate averaging.As a consequence, axiomatic conservation properties, e.g.mass, momentum, energy, are scale-independently satisfied.
Positions of specific Lagrangian particles are represented by y ∈ R 3 .Further, we assume a constant kinematic viscosity ν = const and Mach number Ma < 0.3.
The latter implies that we can simplify the viscous stress term to div[τ visc ] = div[2νρD] as ∇ • v 0 46,47 , with D denoting the symmetric strain rate tensor.Finally, the considered LES model in its Lagrangian form reads

TABLE I .
SPH cases for the considered Taylor-Green flow.

TABLE I :
Case 5-8).Especially in terms of the global dynamics, represented by FIG. 6 (b) & FIG.