1. Introduction
The smoothed particle hydrodynamics (SPH) method was proposed in 1977 as a Lagrangian discretisation method for fluid dynamics (Gingold & Monaghan Reference Gingold and Monaghan1977; Lucy Reference Lucy1977) and matured significantly since then, as detailed in several reviews (Springel Reference Springel2010; Monaghan Reference Monaghan2012; Price Reference Price2012; Shadloo, Oger & Touze Reference Shadloo, Oger and Le Touzé2016; Ye et al. Reference Ye, Pan, Huang and Liu2019; Lind, Rogers & Stansby Reference Lind, Rogers and Stansby2020; Sigalotti, Klapp & Gesteira Reference Sigalotti, Klapp and Gesteira2021). Originally, the SPH method featured serious numerical convergence problems due to the fact that the consistency of spatial derivative operators is strongly affected by the local particle arrangement, which only can be compensated by a drastic increase in the number of neighbour particles
$N_{\textit{ngb}}$
(Zhu, Hernquist & Li Reference Zhu, Hernquist and Li2015). This is especially problematic in strong shear flows and subsonic turbulence, resulting in zeroth-order errors related to excessive numerical dissipation for small
$N_{\textit{ngb}}$
(Ellero, Español & Adams Reference Ellero, Español and Adams2010; Bauer & Springel Reference Bauer and Springel2012; Colagrossi et al. Reference Colagrossi, Souto-Iglesias, Antuono and Marrone2013; Hopkins Reference Hopkins2015).
Pioneered by the work of Vila (Reference Vila1999), current SPH methods, as recently compared by Eirís et al. (Reference Eirís, Ramírez, Couceiro, Fernández-Fidalgo, París and Nogueira2023), mostly eliminated this convergence issue by the use of at least one of the following two strategies:
-
(i) a consistent approximation of spatial derivatives by either a reproducing kernel (RK) or moving least squares (MLS) approach (Hopkins Reference Hopkins2015; Frontiere, Raskin & Owen Reference Frontiere, Raskin and Owen2017)
-
(ii) an arbitrary Lagrangian–Eulerian (ALE) formulation with transport velocity as noise mitigation technique (Oger et al. Reference Oger, Marrone, Le Touzé and de Leffe2016; Antuono et al. Reference Antuono, Sun, Marrone and Colagrossi2021b ), also rigorously incorporating particle shifting (Xu, Stansby & Laurence Reference Xu, Stansby and Laurence2009; Lind et al. Reference Lind, Xu, Stansby and Rogers2012).
In line with these positive developments, the confidence in the ability of current SPH methods to capture incompressible turbulence increases. Especially in complex multiphase flow situations, where SPH can play to its strengths, these methods are nowadays optimistically combined with a large eddy simulation (LES) perspective (Colagrossi et al. Reference Colagrossi, Marrone, Colagrossi and Touzé2021; Lai et al. Reference Lai, Li, Yan, Liu and Zhang2022; King et al. Reference King, Lind, Rogers, Stansby and Vacondio2023; Meringolo et al. Reference Meringolo, Lauria, Aristodemo and Filianoti2023). In the following, we will refer to this combination of SPH methods with LES simply as SPH-LES for brevity. As already mentioned by Bicknell (Reference Bicknell1991), such a confluence is intuitive since the SPH kernel and the LES low-pass filter, building the foundation of both methods, are mathematically congruent. Certainly, the development of rigorous combined SPH-LES theories is a very recent topic; e.g. Di Mascio et al. (Reference Di Mascio, Antuono, Colagrossi and Marrone2017), Antuono et al. (Reference Antuono, Marrone, Mascio and Colagrossi2021a ), Okraschevski et al. (Reference Okraschevski, Buerkle, Koch and Bauer2022).
One central remaining issue is that all these current SPH-LES studies intuitively model one of the central LES objects, namely the subfilter stress (SFS) tensor
$\boldsymbol{\tau }_{\textit{SFS}}$
, by classical, functional eddy viscosity approaches, employing the Boussinesq hypothesis (Schmitt Reference Schmitt2007). In Okraschevski et al. (Reference Okraschevski, Buerkle, Koch and Bauer2022), we could argue that for classical SPH, physically reinterpreting this method from a spatial coarse-graining perspective, eddy viscosity modelling must fail due to the non-locality of the SPH method (Du & Tian Reference Du and Tian2020; Vignjevic, DeVuyst & Campbell Reference Vignjevic, DeVuyst and Campbell2021; Yao, Zhou & Qian Reference Yao, Zhou and Qian2022). The resulting incompatibility is illustrated in figure 1 and manifests in a spectral mismatch where the classical SFS model is introduced. Yet one might object that this is just a consequence of the classical SPH approach considered in our former study. This is where the following work comes into play showing that our coarse-graining perspective generally applies to current SPH methods. (The aforementioned convergence issues with SPH, the countermeasures listed above and the subsequent rationale strictly apply only to Lagrangian SPH methods, inherently including particle disorder. For regular particle distributions in an Eulerian frame, the convergence issues to begin with can be fully eliminated (Hopkins Reference Hopkins2015; Lind & Stansby Reference Lind and Stansby2016) with proper numerical schemes (Appendix A), and accordingly, the spectral peculiarities in figure 1. However, since we believe that the Lagrangian character is still a key argument for current SPH methods, we will subsequently focus on this specific reference frame.) Hence classical eddy viscosity modelling, as already indicated by Rennehan (Reference Rennehan2021), is compromising the most accurate prediction of turbulence possible, and necessitates the development of completely new and specific SFS models in the SPH-LES context. This perspective could accentuate even more the advantages of current SPH methods over traditional grid-based approaches in application areas such as multiphase flows encompassing turbulence.
Figure 1. Typical distribution of spectral energy density obtained with SPH methods for incompressible turbulence. The properly resolved range with large eddies passes into an energy deficit range that is non-locally caused and followed by a Lagrangian noise range. From an optimal SFS model, we would expect a reduction of the Lagrangian noise in favour of the deficit range. However, with incompatible classical SFS models, the noise is barely reduced and the deficit range is exacerbated due to non-locality.
However, should such an improved SFS model not be available soon, we definitively advise operating current SPH methods for aforementioned flows using no explicit SFS model. Thus we will broadly argue in favour of the notion that current Lagrangian SPH methods operate intrinsically as implicit LES. The latter is an established concept in the grid-based community, and relies on properly designed discretisation errors to provide an implicit SFS contribution; e.g. Grinstein, Margolin & Rider (Reference Grinstein, Margolin and Rider2007), Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017), Moura et al. (Reference Moura, Mengaldo, Peiró and Sherwin2017), Fehn et al. (Reference Fehn, Kronbichler, Munch and Wall2022), Volpiani (Reference Volpiani2024). Heuristically, it seems that these implicit SFS in SPH methods emerge from statistical physics principles (Posch, Hoover & Kum Reference Posch, Hoover and Kum1995; Ellero et al. Reference Ellero, Español and Adams2010; Borreguero et al. Reference Borreguero, Bezgin, Adami and Adams2019; Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2021a ), which we will confirm in more detail below using Hardy’s theory (Hardy Reference Hardy1982).
2. Novelty and implications
In our former works (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2021a
,
Reference Okraschevski, Hoffmann, Stichling, Koch and Bauerb
) we laid the theoretical foundation to be summarised in § 3 with Hardy’s theory (Hardy Reference Hardy1982) at its core. Based on this, we could demonstrate that SPH-LES with classical SPH seems fundamentally incompatible with explicit viscosity models due to the non-local characteristic of the discretisation (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022; Okraschevski Reference Okraschevski2024). Here, we report evidence for the first time that this incompatibility also holds for SPH-LES with current SPH methods, which are not plagued by the classical SPH problems as described in § 1. As a consequence, we believe that classical eddy viscosity models for the SFS tensor
$\boldsymbol{\tau }_{\textit{SFS}}$
are detrimental in SPH-LES of incompressible turbulence, and that novel models matching the discretisation characteristics must be developed. This insight will particularly improve the predictive power of current SPH methods in complex multiphase flows, when incompressible turbulence is an inevitable aspect of the considered flow.
3. Coarse-graining perspective on SPH methods
Although current SPH methods might not suffer from the same convergence issue as their original ancestor, at their heart they still employ quasi-Lagrangian particles. (This term is usually used within ALE frameworks (Vogelsberger et al. Reference Vogelsberger, Sijacki, Kereš, Springel and Hernquist2012; Oger et al. Reference Oger, Marrone, Le Touzé and de Leffe2016; Antuono et al. Reference Antuono, Marrone, Mascio and Colagrossi2021a
). Here, we also include pseudo-Lagrangian particles (Vogelsberger et al. Reference Vogelsberger, Sijacki, Kereš, Springel and Hernquist2012) – also called purely Lagrangian particles in other works (Oger et al. Reference Oger, Marrone, Le Touzé and de Leffe2016; Antuono et al. Reference Antuono, Marrone, Mascio and Colagrossi2021a
) – in the definition.) These are connected by a spherical, positive and monotonously decaying kernel
$W: \mathbb{R}^3 \to \mathbb{R}$
with compact support
$V_{\boldsymbol{x}}$
, being centred at the quasi-Lagrangian particles
$\boldsymbol{x} \in \mathbb{R}^3$
. Hence even current SPH methods intrinsically contain two resolution scales, namely the mean particle distance
$\Delta l$
, and the larger kernel diameter
$D_K$
. We termed this property particle duality in our former work (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022). This is a characteristic peculiarity compared to grid-based discretisation techniques, and results in strongly overlapping discretisation elements, i.e. a non-local discretisation. Despite the fact that numerical convergence in current SPH methods can be reached using a constant ratio
$D_K/\Delta l=\textit {O}(1)$
(Vila Reference Vila1999; Hopkins Reference Hopkins2015) resulting in a fixed number of neighbours
$N_{\textit{ngb}}$
inside the kernel, one still might wonder which flow scales can be effectively resolved. Taking a conservative point of view, it must be expected that flow scales can be captured maximally up to the kernel diameter
$D_K$
. By means of such a rationale, one implicitly interprets SPH methods from a spatial coarse-graining perspective at the effective scale
$D_K$
. Such a coarse-graining perspective is not only a convenient footing in the following, but also the physical foundation of the LES community (Eyink & Drivas Reference Eyink and Drivas2018), and even more so a general perspective employed by fluid dynamicists (Irving & Kirkwood Reference Irving and Kirkwood1950; Okraschevski et al. Reference Okraschevski, Hoffmann, Stichling, Koch and Bauer2021b
; Eyink Reference Eyink2024). In the LES community, the coarse-graining is also synonymously called the filtering approach (Germano Reference Germano1992). Although the different terminologies describe equivalent mathematical operations, we believe that the term coarse-graining raises the awareness for a geometric interpretation in terms of a hierarchical clustering of Lagrangian particles (figure 2). The aforementioned can be vividly unravelled by generalising the theory of Hardy from statistical physics (Hardy Reference Hardy1982; Okraschevski et al. Reference Okraschevski, Hoffmann, Stichling, Koch and Bauer2021b
), and highlights the conceptual similarity to the Lagrangian discretisation techniques of interest. Hence we anticipate that SPH methods aim to solve an effective field equation and intrinsically operate as Lagrangian LES. By defining the Lagrangian derivative as
\begin{equation} \frac {\mathrm{d}}{\mathrm{d}t}:= \partial _t + \tilde {\boldsymbol{v}} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \end{equation}
and the spatial coarse-graining of a scalar field
$f: \mathbb{R}^3 \times \mathbb{R}^+_0 \to \mathbb{R}$
over
$V_{\boldsymbol{x}}$
as
\begin{equation} \overline {f} (\boldsymbol{x}, t) = \int _{V_{\boldsymbol{x}}} f(\boldsymbol{y}, t)\, W(\boldsymbol{x} - \boldsymbol{y})\,\mathrm{d}\boldsymbol{y}, \end{equation}
we will subsequently consider barotropic flows in a Lagrangian reference frame
\begin{equation} \frac {\mathrm{d} \overline {\rho }} {\mathrm{d}t}(\boldsymbol{x}, t) = -\overline {\rho } (\boldsymbol{x}, t)\, \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\tilde {\boldsymbol{v}}(\boldsymbol{x}, t), \end{equation}
\begin{align} \overline {\rho } (\boldsymbol{x}, t) \frac {\mathrm{d} \tilde {\boldsymbol{v}}} {\mathrm{d}t}(\boldsymbol{x}, t)& = - \boldsymbol{\nabla} _{\boldsymbol{x}} \overline {p}(\boldsymbol{x}, t) + \text{div}_{\boldsymbol{x}} \left [\tilde {\boldsymbol{\tau }}_{v\textit{isc}} + \boldsymbol{\tau }_{\textit{SFS}} \right ](\boldsymbol{x}, t), \end{align}
\begin{align} \overline {p}(\boldsymbol{x}, t) & = \overline {p}_{\textit{ref}} + c_s^2 ( \overline {\rho }(\boldsymbol{x}, t) - \overline {\rho }_{\textit{ref}} ) .\end{align}
Figure 2. Illustration of spatial coarse-graining emerging from the generalisation of Hardy’s theory (Hardy Reference Hardy1982; Okraschevski et al. Reference Okraschevski, Hoffmann, Stichling, Koch and Bauer2021b ). Adapted from Okraschevski et al. (Reference Okraschevski, Buerkle, Koch and Bauer2022).
Moreover, we will assume for (3.3), (3.4) and (3.5) a weakly compressible, low Mach number flow (
$ \textit{Ma} \lt 0.3$
). Hence bulk viscous stresses are neglected. The fields
$\overline {\rho }$
and
$\overline {p}$
denote the coarse-grained density and pressure, whereas
$\tilde {\boldsymbol{v}} = \overline {\rho \boldsymbol{v}}/\overline {\rho }$
is the density-weighted coarse-grained velocity as proposed by Reynolds (Reference Reynolds1895), nowadays called the Favre averaged velocity (Bilger Reference Bilger1975). For a Newtonian fluid, then, the viscous stress tensor reads
\begin{equation} \tilde {\boldsymbol{\tau }}_{v\textit{isc}}=\eta \left(\unicode{x1D645}_{\tilde {\boldsymbol{v}}} +\unicode{x1D645}_{\tilde {\boldsymbol{v}}}^{{\rm T}} - \frac {2}{3}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\tilde {\boldsymbol{v}}\right), \end{equation}
with
$\unicode{x1D645}_{\tilde {\boldsymbol{v}}}$
the Jacobian of
$\tilde {\boldsymbol{v}}$
and
$\unicode{x1D645}_{\tilde {\boldsymbol{v}}}^{\rm T}$
its transpose. The dynamic viscosity
$\eta$
, the reference density
$\overline {\rho }_{\textit{ref}}$
and pressure
$\overline {p}_{\textit{ref}}$
, and the speed of sound
$c_s$
are dealt with as constant parameters to be specified.
The most important object emerging from the spatial coarse-graining at the arbitrary scale
$D_K$
is the SFS tensor
$\boldsymbol{\tau }_{\textit{SFS}}$
. It can be formally defined as (Vreman, Geurts & Kuerten Reference Vreman, Geurts and Kuerten1994; Okraschevski et al. Reference Okraschevski, Hoffmann, Stichling, Koch and Bauer2021b
)
\begin{equation} \boldsymbol{\tau }_{\textit{SFS}} (\boldsymbol{x}, t) :=- \int _{V_{\boldsymbol{x}}} \rho (\boldsymbol{y}, t)\, \boldsymbol{w}(\boldsymbol{x}, \boldsymbol{y}, t)\,\boldsymbol{w}^{{\rm T}}(\boldsymbol{x}, \boldsymbol{y}, t)\, W(\boldsymbol{x} - \boldsymbol{y}) \, \mathrm{d}\boldsymbol{y} =- \overline {\rho \boldsymbol{w}\boldsymbol{w}^{{\rm T}}}(\boldsymbol{x}, t), \end{equation}
with
$\boldsymbol{w}$
denoting the peculiar velocity. The peculiar velocity is a relative velocity connecting coarse-grained velocities
$\tilde {\boldsymbol{v}}(\boldsymbol{x}, t)$
with associated continuum fluid element velocities
$\boldsymbol{v}(\boldsymbol{y}, t)$
(figure 3). By convention, we label and distinguish the spatial coordinates of fluid elements by
$\boldsymbol{y} \in \mathbb{R}^3$
, and of coarse-grained elements by the spatial coordinate
$\boldsymbol{x} \in \mathbb{R}^3$
. Consequently, the peculiar velocity is a non-local quantity by definition, and emerges from
\begin{equation} \boldsymbol{w}(\boldsymbol{x}, \boldsymbol{y}, t) = \boldsymbol{v}(\boldsymbol{y}, t) - \tilde {\boldsymbol{v}}(\boldsymbol{x}, t) . \end{equation}
Figure 3. Visualisation of the velocity decomposition in (3.8). Adapted from Okraschevski (Reference Okraschevski2024).
The velocity decomposition is illustrated in figure 3. We argue that
$\boldsymbol{w}$
is the physically appropriate fluctuating velocity field in the coarse-graining framework. It satisfies
$\overline {\rho \boldsymbol{w}}=0$
by construction (Okraschevski et al. Reference Okraschevski, Hoffmann, Stichling, Koch and Bauer2021b
), and does not require the introduction of generalised central moments to identify the averaging invariance of the turbulent equations (Germano Reference Germano1992). Since the velocities
$\boldsymbol{v}(\boldsymbol{y}, t)$
of the fluid elements are unknown at the coarse-grained level, we practically face the well-known closure problem for the SFS tensor
$\boldsymbol{\tau }_{\textit{SFS}}$
in (3.7).
For turbulent flows at finite resolution, the closure problem is often resolved by explicit modelling of
$\boldsymbol{\tau }_{\textit{SFS}}$
with functional eddy viscosity approaches (Silvis, Remmerswaal & Verstappen Reference Silvis, Remmerswaal and Verstappen2017; Moser, Haering & Yalla Reference Moser, Haering and Yalla2021) assuming dominant physical SFS error over discretisation error. However, as long as the filter width and resolution width are similar, it is well known in the grid-based LES community that these errors are at least of the same order of magnitude, or in fact that the latter exceeds the former (Ghosal Reference Ghosal1996; Dairay et al. Reference Dairay, Lamballais, Laizet and Vassilicos2017). This insight motivated the development of the nowadays well-established implicit LES approaches, in which the discretisation error is designed to provide the SFS contribution; e.g. Grinstein et al. (Reference Grinstein, Margolin and Rider2007), Moura et al. (Reference Moura, Mengaldo, Peiró and Sherwin2017), Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017), Fehn et al. (Reference Fehn, Kronbichler, Munch and Wall2022), Volpiani (Reference Volpiani2024). In this light, it is natural to ask how the explicit and implicit SFS contributions in current SPH methods interact. Is it possible to reduce the significant implicit SFS in current SPH methods, emerging from statistical physics principles (Posch et al. Reference Posch, Hoover and Kum1995; Ellero et al. Reference Ellero, Español and Adams2010; Borreguero et al. Reference Borreguero, Bezgin, Adami and Adams2019; Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2021a
), using explicit models for
$\boldsymbol{\tau }_{\textit{SFS}}$
? This is the leading theme of this work, and in the spirit of similar works in the grid-based LES community (Ghosal Reference Ghosal1996; Dairay et al. Reference Dairay, Lamballais, Laizet and Vassilicos2017).
4. Methods
There is a large variety of SPH methods today, as contrasted by Eirís et al. (Reference Eirís, Ramírez, Couceiro, Fernández-Fidalgo, París and Nogueira2023). Hence care must be taken in the choice of the SPH method for the verification of our incompatibility hypothesis concerning SPH-LES with eddy viscosity models. We decided to use the locally conservative and second-order accurate meshless finite-mass method (MFM), developed and made publicly available in the open source code GIZMO by Hopkins (Reference Hopkins2015). Since the MFM belongs to the large class of current SPH methods termed MLS-SPH-ALE (Eirís et al. Reference Eirís, Ramírez, Couceiro, Fernández-Fidalgo, París and Nogueira2023), it inherently incorporates both strategies mentioned in § 1 to eliminate the convergence issues of classical SPH. The MFM is based on an ALE formulation without particle shifting, which can be operated in either fully Eulerian or quasi-Lagrangian mode, although we will focus on the latter. (In Appendix A, we perform tests on a Cartesian grid in Eulerian mode to investigate whether the incompatibility is generally related to the SPH method or also influenced by the chosen reference frame. Eventually, we realised that by removing the Lagrangian noise and the induced implicit SFS (figure 1), numerical stability becomes an issue hindering the drawing of a final conclusion. Yet we hypothesise that classical eddy viscosity models will withdraw turbulent kinetic energy mostly from scales larger than the kernel even in an Eulerian frame. This would be the natural consequence of the non-locality that we criticise.) We are convinced that the Lagrangian character is still a key argument for the discretisation of fluid flows with current SPH methods, giving natural access to Lagrangian flow properties such as Lagrangian coherent structures (Haller Reference Haller2015; Dauch et al. Reference Dauch, Rapp, Chaussonnet, Braun, Keller, Kaden, Koch, Dachsbacher and Bauer2018). Applying the MFM method to (3.3), (3.4) and (3.5) with a second-order-accurate MLS approximation of spatial derivatives, one arrives, for
$i \in \{1,\ldots ,N \}$
particles with mass
$M_i$
, at
\begin{align} \frac {\mathrm{d} M_i} {\mathrm{d}t}& = 0 \ \implies \overline {\rho }_i = M_i \sum _{j=1}^{N_{\textit{ngb}}}W(\boldsymbol{x}_{i}-\boldsymbol{x}_{j}),\end{align}
\begin{align} M_i \frac {\mathrm{d} \tilde {\boldsymbol{v}}_i} {\mathrm{d}t}& = \sum _{j=1}^{N_{\textit{ngb}}} - \overline {p}_{\textit{ij}}^*\boldsymbol{A}_{\textit{ij}}^{\mathit{eff}} + \left [\tilde {\boldsymbol{\tau }}_{v\textit{isc},\textit{ij}}^* + \boldsymbol{\tau }_{\textit{SFS},\textit{ij}}^* \right ]\boldsymbol{A}_{\textit{ij}}^{\mathit{eff}},\end{align}
\begin{align} \overline {p}_i &= \overline {p}_{\textit{ref}} + c_s^2 ( \overline {\rho }_i - \overline {\rho }_{\textit{ref}} ) ,\end{align}
where we use standard particle notation. Thus the single index
$i$
indicates the numerical proxy of the corresponding field at
$\boldsymbol{x}_{i}$
. The discretised momentum equation (4.2) can be interpreted as a Lagrangian finite-volume formulation with fluxes to be approximated at effective interface areas
$\boldsymbol{A}_{\textit{ij}}^{\mathit{eff}} \in \mathbb{R}^3$
between particles
$i$
and
$j$
. These interface areas depend on the local particle configuration and the chosen kernel, subsequently the pairing-stable Wendland C4 with
$N_{\textit{ngb}}=128$
as our default (Dehnen & Aly Reference Dehnen and Aly2012). The interface fluxes, namely
$\overline {p}_{\textit{ij}}^* \in \mathbb{R}$
and
$\tilde {\boldsymbol{\tau }}_{v\textit{isc},\textit{ij}}^*, \boldsymbol{\tau }_{\textit{SFS},\textit{ij}}^* \in \mathbb{R}^{3 \times 3}$
, are computed with approximate Riemann solvers that are slope- and flux-limited, and embedded into an explicit single-stage second-order-accurate time integration scheme. More details on the computation of the effective interface areas
$\boldsymbol{A}_{\textit{ij}}^{\mathit{eff}}$
, the Harten–Lax–van Leer contact solver for
$\overline {p}_{\textit{ij}}^*$
, the Harten–Lax–van Leer solver for
$\tilde {\boldsymbol{\tau }}_{v\textit{isc},\textit{ij}}^*, \boldsymbol{\tau }_{\textit{SFS},\textit{ij}}^*$
, and aspects beyond that can be found in the works of Hopkins (Reference Hopkins2015, Reference Hopkins2016). We note that the resulting system (4.1), (4.2) and (4.3) is very similar to the well-known classical weakly compressible SPH (WCSPH). Thus we will term it accordingly as weakly compressible MFM (WCMFM). As a consequence of the explicit weakly compressible approach, we expect more Lagrangian noise with stronger implicit SFS due to the allowed acoustic waves, comparing with a stable projection-based incompressible approach at the same resolution (Xu et al. Reference Xu, Stansby and Laurence2009; Lind et al. Reference Lind, Xu, Stansby and Rogers2012). This expectation is not SPH agnostic but seems to be a generally anticipated effect, as demonstrated by Vittoz et al. (Reference Vittoz, Oger, de Leffe and Le Touzé2019) in a grid-based high-order finite-volume context. After all, the resulting WCMFM discretisation is obviously non-local, realising that the differential operators in the momentum balance of (3.4) are transferred to
$N_{\textit{ngb}}$
flux evaluations per particle on the kernel scale
$D_K$
. This is in the spirit of classical SPH (Du & Tian Reference Du and Tian2020; Vignjevic et al. Reference Vignjevic, DeVuyst and Campbell2021; Yao et al. Reference Yao, Zhou and Qian2022). As will be demonstrated in § 5, it seems to be exactly this non-locality that leads to the incompatibility of SPH-LES with eddy viscosity models for incompressible turbulence (figure 1).
As a canonical benchmark we will consider a Taylor–Green flow (Taylor & Green Reference Taylor and Green1937) on the periodic domain
$\varOmega =[0,2\pi ]^3$
for three different particle counts
$N\in \{128^3, 256^3, 512^3\}$
. This is a precious test to evaluate the dissipation characteristics of a numerical solver and its ability to resolve incompressible turbulence; e.g. (Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983; Drikakis et al. Reference Drikakis, Fureby, Grinstein and Youngs2007; Dairay et al. Reference Dairay, Lamballais, Laizet and Vassilicos2017; Moura et al. Reference Moura, Mengaldo, Peiró and Sherwin2017; Pereira et al. Reference Pereira, Grinstein, Israel, Rauenzahn and Girimaji2021; Fehn et al. Reference Fehn, Kronbichler, Munch and Wall2022). As in our former work with classical SPH (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022), we will use the direct numerical simulation (DNS) solution of Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017) at
$Re=10^4$
as a reference. It was computed with the sixth-order finite difference code Incompact3d (Laizet & Li Reference Laizet and Li2011). We will initialise and evaluate the WCMFM simulations exactly as we did in Okraschevski et al. (Reference Okraschevski, Buerkle, Koch and Bauer2022) for classical WCSPH. Thus we will specify the initial root mean square Mach number as
$ \textit{Ma}_{\textit{rms}}(t=0\,\text{s})=\sqrt {2e_v(t=0\,\text{s})}/c_s = 0.1$
such that
$c_s=5\,\text{m}\,\text{s}^{-1}$
,
$\overline {\rho }_{\textit{ref}}=1\,\text{kg}\,{\text{m}^{-3}}$
and
$\overline {p}_{\textit{ref}}=\overline {\rho }_{\textit{ref}} c_s^2/4 = 6.25\,\text{Pa}$
in (4.3). Consequently, the dynamic viscosity in (3.6) must be
$\eta =0.0001\,\rm {Pa\,s}$
to reach the targeted Reynolds number. We will study and compare the temporal evolution of the averaged kinetic energy
$e_v$
, its corresponding averaged dissipation rate
$\epsilon _t=-({\mathrm{d}e_v}/{\mathrm{d}t})$
, and the spectral energy density
$E(k)$
at
$t=14\,\text{s}$
. The time instance was selected in accordance with Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017), which ensures that the turbulence is developed and exhibits the expected inertial range scaling
$E(k)\sim k^{-5/3}$
(Kolmogorov Reference Kolmogorov1941; Obukhov Reference Obukhov1941; Onsager Reference Onsager1945; Heisenberg Reference Heisenberg1948). (Note that there is recent experimental doubt about the quantitative correctness of this scaling in the inertial range; see Küchler et al. (Reference Küchler, Bewley and Bodenschatz2023).) For the evaluation of the spectra, we use the validated methodology of Bauer & Springel (Reference Bauer and Springel2012) in combination with considerations by Durran, Weyn & Menchaca (Reference Durran, Weyn and Menchaca2017). This preserves physically important small-scale features of the flow, and guarantees the validity of the discrete Parseval relation.
Finally, to work out the interplay between explicit and implicit SFS, we need to specify an explicit SFS model and think about how the implicit SFS contribution can be estimated. For the former, we decided to opt for the
$\sigma$
model by Nicoud et al. (Reference Nicoud, Toda, Cabrit, Bose and Lee2011), which is one of the most sophisticated static eddy viscosity models. It eliminates artificial dissipation in two-dimensional flows, laminar shear zones and solid body rotation, but likewise shows proper asymptotic scaling near walls (Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011; Silvis et al. Reference Silvis, Remmerswaal and Verstappen2017; Moser et al. Reference Moser, Haering and Yalla2021). The explicit model reads in continuous representation
\begin{equation} \boldsymbol{\tau }_{\textit{SFS}}^{exp}=\eta _{\textit{SFS}}\left(\unicode{x1D645}_{\tilde {\boldsymbol{v}}} +\unicode{x1D645}_{\tilde {\boldsymbol{v}}}^{{\rm T}} - \frac {2}{3}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\tilde {\boldsymbol{v}}\right), \quad \eta _{\textit{SFS}}:=\overline {\rho }(C_\sigma \varDelta )^2 \frac {(\sigma _1-\sigma _2)(\sigma _2-\sigma _3)\sigma _3}{\sigma _1^2}, \end{equation}
where
$\sigma _k,\ k\in \{1,2,3\}$
, are the singular values of the tensor
$\unicode{x1D645}_{\tilde {\boldsymbol{v}}}^{{\rm T}}\unicode{x1D645}_{\tilde {\boldsymbol{v}}}$
,
$C_{\sigma }=1.35$
is a model constant, and the filter width is
$\varDelta =D_K$
. The latter is an unambiguous choice emerging from our coarse-graining perspective and a matter of debate in the SPH-LES context (Rennehan Reference Rennehan2021; King et al. Reference King, Lind, Rogers, Stansby and Vacondio2023). We will elaborate on it more closely in § 5. For the estimation of the implicit SFS tensor, we assume that density changes of the fluid elements over the kernel scale
$D_K$
are much weaker than the corresponding velocity changes. This is reasonable for a weakly compressible flow developing characteristics of incompressible turbulence, and implies that
$\tilde {\boldsymbol{v}}=\overline {\boldsymbol{v}}$
. Then with the spatial coarse-graining operation in (3.2) and linearisation of the continuum velocity field at a position
$\boldsymbol{z}\in \mathbb{R}^3$
, one obtains
$\tilde {\boldsymbol{v}}(\boldsymbol{x}=\boldsymbol{z},t) = \boldsymbol{v}(\boldsymbol{y}=\boldsymbol{z},t) + \textit {O}(D_K^2)$
, which gives a consistency preserving (with respect to the MFM discretisation) second-order proxy for the peculiar velocity in (3.8), namely
\begin{equation} \boldsymbol{w}(\boldsymbol{x}_1, \boldsymbol{x}_2, t) \approx \tilde {\boldsymbol{v}}(\boldsymbol{x}_2, t) - \tilde {\boldsymbol{v}}(\boldsymbol{x}_1, t), \end{equation}
with two different coarse-grained coordinates
$\boldsymbol{x}_1,\boldsymbol{x}_2$
. Inserting (4.5) into the SFS tensor in (3.7) gives, after discretisation of the integral into finite-mass elements, the following estimator for the implicit SFS in particle notation:
\begin{equation} \boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}} \approx - \sum _{j=1}^{N_{\textit{ngb}}} (\tilde {\boldsymbol{v}}_j-\tilde {\boldsymbol{v}}_i)(\tilde {\boldsymbol{v}}_j-\tilde {\boldsymbol{v}}_i)^{{\rm T}}\,W(\boldsymbol{x}_i-\boldsymbol{x}_j)\,M_j. \end{equation}
Since incompressible turbulence is a convection-driven phenomenon, we will evaluate the local importance of the implicit SFS on a particle
$i$
in comparison to the coarse-grained convective stress tensor
$\overline {\rho }_i\tilde {\boldsymbol{v}}_i\tilde {\boldsymbol{v}}_i^{{\rm T}}$
. Therefore, we define the
$R$
-index as
\begin{equation} R_i:=\frac {\| \boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}} \|_F}{\|\boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}}\|_F + \|\overline {\rho }_i\tilde {\boldsymbol{v}}_i\tilde {\boldsymbol{v}}_i^{{\rm T}}\|_F}, \end{equation}
with
$\|\cdot\|_F$
denoting the Frobenius norm.
5. Results and discussion
In this section, we will present and critically discuss the results of our study. We start with a qualitative verification of our implementation of the
$\sigma$
model (Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011) in (4.4) into the code GIZMO (Hopkins Reference Hopkins2015). Therefore, the uniquely coloured individual particle IDs are utilised as flow tracers and visualised for the purpose of structure identification (figure 4). Two snapshots before and after the well-known dissipation peak (Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983), namely at
$t_{1,2}=9\pm 3\ \text{s}$
, are shown in figures 4(a,c) and 4(b,d), respectively, for the highest resolution runs
$N=512^3$
. They vividly render the development of primary instabilities and turbulence in the flow field. Figure 4(a,b) display the case without explicit SFS model (WCMFM), and figure 4(c,d) with explicit SFS model (WCMFM + SIGMA). In figure 4(e,f), the central quantity of the
$\sigma$
model, namely the eddy viscosity field according to (4.4), is shown. It is scaled by the dynamic viscosity of the flow. The ratio, denoted as
$\eta ^*$
, visually correlates with the flow structures and is evidently non-negligible in the shear flow planes (figure 4
a) where incompressible turbulence develops. In this region, the eddy viscosity is dominant over the dynamic viscosity by up to an order of magnitude. This is an anticipated consequence as the coarse-grained viscous stress tensor is bounded from above by the Cauchy–Schwarz inequality, and it can be proven that the bound scales with
$1/D_K$
(Eyink & Drivas Reference Eyink and Drivas2018). Hence for regions of under-resolved turbulence, in which the kernel scale
$D_K$
is larger than viscous length scale, one would expect
$\eta ^* \gg 1$
. A comparison of the flow structures predicted by WCMFM and WCMFM + SIGMA reveals that the explicit SFS model does not inhibit the dynamics of the primary instabilities before the dissipation peak (figure 4
a versus figure 4
c). However, it apparently damps noisy small-scale features in the turbulent flow field after the dissipation peak (figure 4
b versus figure 4
d). This could be incorrectly interpreted as removal of numerically dissipative Lagrangian noise (artificial thermalisation). As extensively discussed by Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017) for grid-based methods, the removal of such an artificial thermalisation should be the ultimate goal of an explicit SFS model, eventually resulting in a quantitative improvement for the dissipation characteristics in physical and spectral space. We will show now that the explicit SFS model completely fails with respect to such a quantitative analysis, although visually it seems to perform well.
Figure 4. Qualitative verification of the implementation of the
$\sigma$
model (Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011) for
$N=512^3$
. Flow structures before and after the dissipation peak (a,b) for the case without explicit SFS model (WCMFM), and (c,d) for the case with explicit SFS model (WCMFM + SIGMA). (e,f) The scaled eddy viscosity field.
Therefore, we compare the temporal evolution of the averaged kinetic energy
$e_v$
, the averaged dissipation rate
$\epsilon _t$
and the spectral energy density
$E(k,t= 14\,\text{s})$
in figure 5. The latter is scaled with the corresponding kinetic energy value
$e_v(t= 14\,\text{s})$
and
$L_c=1\,\text{m}$
, such that integration over the wavenumber shells always results in unity. Cases without explicit SFS model (WCMFM) are displayed as solid blue lines, whereas cases with explicit SFS model (WCMFM + SIGMA) are displayed as red dashed lines. It is evident for the kinetic energy evolution in figure 5(a) that qualitative convergence towards the DNS (solid black line) for increasing particle count can be obtained. After the dissipation peak, as soon as incompressible turbulence develops, the theoretically predicted Saffman decay rate (Skrbek & Stalp Reference Skrbek and Stalp2000)
$e_v\sim (t-t_0)^{-1.2}$
can be matched (dashed black line). Here,
$t_0$
denotes a time shift parameter that accounts for the earlier transitions at lower resolution. These observations are independent of the explicit SFS model and also reflected by the dissipation rate profiles in figure 5(b). However, the results clearly demonstrate the detrimental effect of the eddy viscosity model on the dissipation characteristics in physical space. Coincidentally, for the chosen configuration, it seems that the explicit SFS model leads to a fallback of the accuracy by approximately a whole resolution step. The WCMFM cases for
$N=128^3$
and
$N=256^3$
behave very similar to the WCMFM + SIGMA cases for
$N=256^3$
and
$N=512^3$
. In other words, under ideal code scaling with CFL restriction, to achieve the same result with an explicit SFS model, at least a
$2^4=16$
times higher computational effort is required. This is a quite drastic finding. Note that for a quasi-Lagrangian particle method, the highest-resolution WCMFM run (
$N=512^3$
) without explicit SFS model leads to a comparably accurate prediction of the dissipation characteristics in physical space. In particular, the sharper prediction of the global dissipation peak and the local peak nearby in figure 5(b) compared to WCSPH with twice as many neighbours (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022) is prominent.
Figure 5. Quantitative effect of the
$\sigma$
model (Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011) in physical and spectral space for different resolutions: (a) averaged kinetic energy; (b) averaged dissipation rate; (c) scaled spectral energy density at
$t=14\,\text{s}$
for DNS and WCMFM run (
$N=512^3$
) without explicit SFS model; (d) scaled spectral energy density at
$t=14\,\text{s}$
. For orientation, the kernel scale for
$N=512^3$
is included.
Before we proceed with the comparison on the spectral dissipation characteristics for all resolutions, we will first elaborate on the highest-resolution WCMFM run in figure 5(c). Provably, WCMFM as current SPH method is able to predict the inertial range scaling of incompressible turbulence by approximately an order of magnitude in wavenumber space. However, in comparison to the DNS case, the inertial range already terminates at
$k_{\textit{MFM}}\lt k_{\textit{DNS}}$
(by
$k_{\textit{MFM}}$
we denote the effective wavenumber up to which a qualitatively correct spectral behaviour can be observed and not the wavenumber corresponding to the mean particle diameter
$\Delta l$
), then passes into the energy deficit range
$k\in [k_{\textit{MFM}},k_{kern}]$
and is followed by the Lagrangian noise range (artificial thermalisation) resulting from kernel scale errors. Qualitatively, this is similar to classical SPH, and happens although the kernel wavelength satisfies
$k_{kern}\gt k_{\textit{DNS}}$
(Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022). It is indicative for the non-local character of the method, the emerging implicit SFS according to (4.6) and a reaction to peculiar velocities on the kernel scale, which manifest as artificial thermalisation. In the next paragraph this will be detailed, but prior to that we want to show in figure 5(d) that the explicit SFS model also deteriorates the situation in spectral space due to the non-locality of the method. Instead of removing a significant part of kinetic energy from the artificial thermalisation, it dominantly withdraws kinetic energy from the energy deficit range and the partially resolved inertial range. Hence it erroneously attacks scales that are already badly resolved, and not the numerically dissipative ones. While these observations should be sufficient to advise against the usage of classical eddy viscosity models in the SPH-LES context, we want to unravel the interplay of the implicit SFS and explicit SFS models in the following paragraph, and underpin these spectral observations in physical space.
Ergo, we will visualise the implicit SFS relative to the convective stress using the
$R$
-index defined in (4.7), and investigate how the field is affected by resolution and the explicit SFS model. We will focus on the plane
$x=\pi$
(figure 4
a) for the time
$t=14\,\text{s}$
, exactly corresponding to the spectra in figures 5(c) and 5(d). The resulting fields are depicted in figure 6 and contain a kernel element at the given resolution in the upper right corner to assess the extent of emerging coherent structures. We will start with the influence of the spatial resolution for WCMFM without explicit SFS model shown in figure 6(a,c,e). Note that due to the definition of the
$R$
-index, the exact ratio of implicit SFS to the convective stress is given by
$\| \boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}} \|_F/\|\overline {\rho }_i\tilde {\boldsymbol{v}}_i\tilde {\boldsymbol{v}}_i^{{\rm T}}\|_F=R_i/(1-R_i)$
, and the upper limit of the colour bar indicates that the implicit SFS is a third of the convective stress. Interestingly, for all resolutions, the implicit SFS is relevant only around the shear flow planes where incompressible turbulence develops and forms a coherent network that surpasses the kernel scale; although the coherent network becomes more delicate for higher resolution (cf. figure 6
a,e), structures with significant contribution remain larger in extent than the given kernel element. Considering that the implicit SFS is a consequence of the peculiar velocities evaluated at the kernel scale but unfolds its effect well beyond the kernel scale, we believe we can see the deterministic reason for how the artificial thermalisation causes the energy deficit range in spectral space (figure 5
c,d). Hence the
$\| \boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}} \|_F$
field unravels how the kernel scale effects propagate to larger scales due to the non-local character of the method. Focusing now on figure 6(b,d,f), in which the WCMFM + SIGMA cases are shown, a reflection of the spectral behaviour in figure 5(d) can be clearly seen as well. Even though we would expect a working explicit SFS model to significantly reduce
$\| \boldsymbol{\tau }_{\textit{SFS},i}^{\textit{imp}} \|_F$
with respect to
$R$
on the kernel scale, we see a strong non-local damping. The network itself is destroyed, instead of diminishing its amplitude. Recalling the similarity of the WCMFM cases for
$N=128^3$
and
$N=256^3$
, and WCMFM + SIGMA cases for
$N=256^3$
and
$N=512^3$
in physical space highlighted above, our requirements on a working explicit SFS model can be specified more precisely. It manifests in the comparison of figures 6(a) versus 6(d), and figures 6(c) versus 6(f). There, the network structure is mainly preserved, and the amplitude of
$R$
diminished. Together with the spectral statistics in figure 5(d), it is evident that this observation correlates with a diminished artificial thermalisation and gain in the energy deficit range. This is what we would expect from an explicit SFS model in the SPH-LES context in a much stronger form for a fixed resolution, but unfortunately see that even sophisticated eddy viscosity models, such as the
$\sigma$
model by Nicoud et al. (Reference Nicoud, Toda, Cabrit, Bose and Lee2011), fail.
Figure 6. Implicit SFS measured by the
$R$
-index in (4.7) at the plane
$x=\pi$
for the time
$t=14\,\text{s}$
. Different resolutions are shown without explicit SFS model (WCMFM) and with explicit SFS model (WCMFM + SIGMA).
Taking this new evidence for a current SPH method into account with our former, congruent results for classical SPH (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022), we are convinced that SPH-LES with classical eddy viscosity models is highly detrimental for the prediction of incompressible turbulence. In our opinion, this is due to the mismatch of discretisation characteristics, namely quasi-Lagrangian particles and non-locality, with the explicit SFS model therefore being spectrally introduced in the already problematic energy deficit range. We observe the well-known issue from grid-based LES that an a priori correct model can perform badly in simulations a posteriori (Park, Yoo & Choi Reference Park, Yoo and Choi2004; Dairay et al. Reference Dairay, Lamballais, Laizet and Vassilicos2017). From our coarse-graining perspective, current SPH methods operate intrinsically as Lagrangian LES with implicit SFS for incompressible turbulence. We want to stress that the chosen terminology is not indicative of the quality of this implicit LES approach, which manifests in the spectral energy density in figure 5(c). Obviously, it is plagued by an energy deficit range due to artificial thermalisation. We have shown already in our former work that the inertial range can be reproduced with grid-based finite-volume Smagorinsky LES at much lower resolution
$N=384^3$
and lower computational cost (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022).
Before we finally conclude our work, we want to address two more aspects that are no less important: first, the effect of the filter width
$\varDelta$
when using an explicit SFS model; and second, the effect of the neighbour particles
$N_{\textit{ngb}}$
for simulations without explicit SFS model. The results are displayed in figure 7 for all cases considering
$N=256^3$
.
Figure 7. (a,c) Influence of the filter width
$\varDelta$
for WCMFM + SIGMA, and (b,d) the neighbour particles
$N_{\textit{ngb}}$
for WCMFM without explicit SFS model. All simulations were performed with
$N=256^3$
.
The first aspect is an ongoing matter of debate in the SPH-LES context (Rennehan Reference Rennehan2021; King et al. Reference King, Lind, Rogers, Stansby and Vacondio2023). Although from our coarse-graining perspective the choice
$\varDelta =D_K$
is unambiguous, we have considered the two further cases with
$\varDelta =D_K/2$
(Rennehan Reference Rennehan2021) and
$\varDelta =2D_K$
. The case
$\varDelta =0$
corresponds to the run without an explicit SFS model. Obviously from figures 7(a) and 7(c), we observe a monotonic trend with increasing filter width. The effect is detrimental in physical space as well as in spectral space, but starts to get irrelevant for
$\varDelta \leqslant D_K/2$
. This is congruent with the observation of Rennehan (Reference Rennehan2021), who studied forced subsonic turbulence with MFM and different explicit SFS models using a dynamic procedure. It prompts the eventuality that SPH-LES studies with
$\Delta$
equal to or smaller than the kernel radius – e.g. Antuono et al. (Reference Antuono, Sun, Marrone and Colagrossi2021b
), Colagrossi et al. (Reference Colagrossi, Marrone, Colagrossi and Touzé2021), Lai et al. (Reference Lai, Li, Yan, Liu and Zhang2022), King et al. (Reference King, Lind, Rogers, Stansby and Vacondio2023) – just run (slightly) more expensive simulations in which the influence of the explicit SFS model is negligible.
The second aspect concerns the influence of the number of neighbours
$N_{\textit{ngb}}$
for WCMFM without explicit SFS model. In classical SPH, better convergence towards the DNS can be obtained by increasing
$N_{\textit{ngb}}$
(Zhu et al. Reference Zhu, Hernquist and Li2015; Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022), which is in the spirit of explicitly filtered LES (Lund Reference Lund2003; Bose, Moin & You Reference Bose, Moin and You2010). However, for current SPH methods as MFM one would expect minor influence of
$N_{\textit{ngb}}$
by construction (Vila Reference Vila1999; Hopkins Reference Hopkins2015). Indeed, a negligible influence is found according to figures 7(b) and 7(d) for the well-resolved scales. Merely in spectral space, we see a loss of kinetic energy from the artificial thermalisation and the energy deficit range that seems irrelevant for the physical space dynamics. This is contrary to classical SPH, practically eliminates
$N_{\textit{ngb}}$
as a crucial calibration parameter, and is beneficial for the computational cost at finite resolution. Coincidentally, we observe for
$N_{\textit{ngb}}\neq 128$
that the associated kernel scale in figure 7(d) is not perfectly separating the energy deficit range and the artificial thermalisation as in the former cases.
6. Concluding remarks
In this study, we presented evidence for the incompatibility of current SPH methods and classical eddy viscosity models for scale-resolved incompressible turbulence. This result is of particular importance for current SPH methods, which are often advantageously applied in complex multiphase flows potentially encompassing incompressible turbulence (Colagrossi et al. Reference Colagrossi, Marrone, Colagrossi and Touzé2021; Lai et al. Reference Lai, Li, Yan, Liu and Zhang2022; King et al. Reference King, Lind, Rogers, Stansby and Vacondio2023; Meringolo et al. Reference Meringolo, Lauria, Aristodemo and Filianoti2023). With our coarse-graining perspective, we could argue and show for MFM, as a representative from the class of MLS-SPH-ALE approaches (Eirís et al. Reference Eirís, Ramírez, Couceiro, Fernández-Fidalgo, París and Nogueira2023), that it intrinsically operates as Lagrangian LES with significant implicit SFS. Even a sophisticated eddy viscosity model such as the
$\sigma$
model by Nicoud et al. (Reference Nicoud, Toda, Cabrit, Bose and Lee2011) is not able to reduce this implicit SFS due to the non-locality of the discretisation method. For the unambiguous filter width
$\varDelta =D_k$
, the explicit SFS model attacks dominant scales larger than the kernel that are already underresolved. Hence the explicit SFS model is either highly detrimental physically, or for choices
$\varDelta \leqslant D_K/2$
, is likely irrelevant and introduces computational overhead. The latter complies with results of Rennehan (Reference Rennehan2021). To our knowledge, such a study of the interplay between implicit SFS and explicit SFS models for a current SPH method is presented for the first time, revealing the familiar a priori versus a posteriori dilemma in grid-based LES; e.g. Park et al. (Reference Park, Yoo and Choi2004), Dairay et al. (Reference Dairay, Lamballais, Laizet and Vassilicos2017). It elucidates that SPH-LES approaches, even with current SPH methods, must focus on the mitigation of the implicit SFS to improve the predictive power of the Lagrangian LES. In particular, one core application area of SPH, namely multiphase flows, could benefit from it if turbulence matters. This progress can be realised either directly by completely new explicit SFS models, which match the discretisation characteristics of current SPH methods, or indirectly by numerical noise mitigation techniques.
Hypothetically, in this work there exist two disregarded mainstream aspects of current SPH methods compliant with the latter strategy. Particle shifting (Xu et al. Reference Xu, Stansby and Laurence2009; Lind et al. Reference Lind, Xu, Stansby and Rogers2012) and density diffusion (Antuono et al. Reference Antuono, Colagross, Marrone and Molteni2010; Marrone et al. Reference Marrone, Antuono, Colagrossi, Colicchio, Le Touzé and Graziani2011) are established numerical noise mitigation techniques, also in the SPH-LES context (Antuono et al. Reference Antuono, Marrone, Mascio and Colagrossi2021a , Reference Antuono, Sun, Marrone and Colagrossib ; Colagrossi et al. Reference Colagrossi, Marrone, Colagrossi and Touzé2021; Meringolo et al. Reference Meringolo, Lauria, Aristodemo and Filianoti2023). It is likely that both reduce the peculiar velocities on the kernel scale according to (4.5), and indirectly have positive feedback on the implicit SFS according to (4.6). However, we believe that the density diffusion, which can be heuristically rationalised by a coarse-graining without the density-weighted Favre velocity (Di Mascio et al. Reference Di Mascio, Antuono, Colagrossi and Marrone2017), will also be plagued by non-local effects of the discretisation scheme. Very likely, it will mitigate not only noisy but also physically under-resolved scales. Nevertheless, this needs to be tested in follow-up studies.
Another concluding point that needs to be stressed is that our coarse-graining theory is valid only for
$D_K\approx \text{const.}$
and
$\rho \approx \text{const.}$
, hence when compressibility effects are negligible. This is what we ensured by our initial root mean square Mach number choice
$ \textit{Ma}_{\textit{rms}}(t=0\,\text{s}) = \sqrt {2e_v(t=0\,\text{s})}/c_s = 0.1\lt 0.3$
(Jakobsen Reference Jakobsen2014) for the decaying Taylor–Green flow at
$Re=10^4$
. Considering this aspect, an essential difference between our results and the results of Rennehan (Reference Rennehan2021) should be highlighted. The energy deficit range observed in this study and our former one for classical SPH (Okraschevski et al. Reference Okraschevski, Buerkle, Koch and Bauer2022) is there replaced by an energy pile-up known from highly accurate discontinuous Galerkin methods (Moura et al. Reference Moura, Mengaldo, Peiró and Sherwin2017; Fehn et al. Reference Fehn, Kronbichler, Munch and Wall2022). We believe that this difference is rooted in the Mach number at which the spectral statistics are computed. Whereas
$ \textit{Ma}_{\textit{rms}}\approx 0.3$
is controlled by forcing in Rennehan (Reference Rennehan2021), the Mach number in our decaying case can be estimated to be
$ \textit{Ma}_{\textit{rms}}(t=14\,\text{s}) = \sqrt {2e_v(t=14\,\text{s})}/c_s\approx 0.06$
. It was already shown by Hopkins (Reference Hopkins2015) that the accuracy of the MFM method deteriorates with lower Mach number, possibly explaining the switch from the energy pile-up to the energy deficit also known from classical SPH. Such low Mach numbers are not unusual in engineering applications, and favourably comply with our coarse-graining theory.
For completeness, we would like to bring up that alternative Lagrangian LES approaches exist that are much better suited to accurately predict incompressible turbulence, if the latter is the exclusive goal. These are the so-called vortex particle methods (Alvarez & Ning Reference Alvarez and Ning2024), which solve the Navier–Stokes equations in their velocity–vorticity form by a decomposition of the fields into self-adaptive Lagrangian vortex elements. Interestingly, these methods suffer from the exact opposite picture developed in this work for current SPH methods, namely insufficient implicit SFS to obtain stability. In order to restore the latter, explicit SFS models seem inevitable (Mansfield et al. Reference Mansfield, Knio and Meneveau1998, Reference Mansfield, Knio and Meneveau1999; Alvarez & Ning Reference Alvarez and Ning2024). This said, one might wonder whether there is a chance to develop a Lagrangian LES method to retain the high accuracy of the vortex particle methods for incompressible turbulence and flexibility and robustness of current SPH methods required for complex multiphase flows. A promising machine-learning based route towards such an ultimate goal could be the one presented by Woodward et al. (Reference Woodward, Tian, Hyett, Fryer, Stepanov, Livescu and Chertkov2023) and Tian et al. (Reference Tian, Woodward, Stepanov, Fryer, Hyett, Livescu and Chertkov2023).
Acknowledgements
The authors acknowledge support by the state of Baden-Württemberg through bwHPC. The authors would like to thank S. Joshi for critical feedback on the manuscript as well as for valuable discussions.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Incompatibility in the Eulerian reference frame?
In order to elaborate on the influence of the chosen reference frame, we repeated the numerical experiments presented in § 5 in a fully Eulerian mode for
$N=256^3$
. The GIZMO code permits such a course of action as it is methodologically based on an MLS-SPH-ALE framework. The underlying idea is to remove the Lagrangian noise (Hopkins Reference Hopkins2015; Lind & Stansby Reference Lind and Stansby2016) and the corresponding implicit SFS by construction to clearly isolate the effect of the explicit SFS model. This allows us to analyse whether the incompatibility of classical eddy viscosity models with current SPH methods, rooted in the non-local discretisation, persists in an Eulerian reference frame.
The results for the same parameters and initial conditions are shown in figure 8(a). Obviously, the numerical dissipation solely introduced by the approximate Riemann solvers for the interface fluxes (§ 4) is not sufficient to obtain a stable solution, even leading to a blow-up of kinetic energy beyond
$t\gt 10\,\text{s}$
in the Eulerian case. Although the
$\sigma$
model by Nicoud et al. (Reference Nicoud, Toda, Cabrit, Bose and Lee2011) can prevent the blow-up (Eulerian + SFS), it seems not appropriate to eradicate the actual cause of the instability as the kinetic energy starts to plateau beyond
$t\gt 10\,\text{s}$
, also developing oscillatory behaviour. The instability is rooted in constructively interfering acoustic waves leading to density changes of 30 % deviation from
$\overline {\rho }_{\textit{ref}} = 1\,\text{kg}\,\text{m}^{-3}$
(not shown herein), which is in contrast to the allowed magnitude of our explicit weakly compressible approach for the targeted
$ \textit{Ma}_{\textit{rms}}\leqslant 0.1$
(§ 4). We will now show that these interfering acoustic waves cannot be eliminated in the Eulerian reference frame with numerical adaptions in GIZMO, unfortunately hindering the drawing of a final conclusion in terms of generalisation of the incompatibility that we observe in the Lagrangian mode. This interference remains even if we switch to the most dissipative numerical scheme in GIZMO (figure 8
f).
Figure 8. Comparison of Lagrangian and Eulerian results for
$N=256^3$
. (a) Averaged kinetic energy with same parameters. (b) Averaged kinetic energy with conservative slope limiter. (c) Averaged kinetic energy with Kurganov–Tadmor (KT) scheme and conservative slope limiter. (d) Velocity magnitude at
$t=14\,\text{s}$
from the Lagrangian WCMFM run without SFS. (e) The Eulerian run with Kurganov–Tadmor scheme and conservative slope limiter and SFS. (f) The corresponding density field of the latter.
As noted by Hopkins (Reference Hopkins2015), the slope-limiting procedure specifically designed and calibrated for his MLS-SPH-ALE methods in the Lagrangian reference frame must be likely adapted to more conservative parameters in the context of Cartesian grids in an Eulerian mode. Using the most conservative choice in GIZMO, available as a pre-processor directive (SLOPE_LIMITER_TOLERANCE = 0), results in the averaged kinetic energy evolution in figure 8(b). Although this gives a stable result, we obtain seemingly sub-dissipative behaviour in the kinetic energy after
$t\gt 6\,\text{s}$
, and adding the explicit SFS model does not lead to qualitative improvement. However, realising that the Eulerian run without SFS can perfectly match the kinetic energy characteristics up to
$t\leqslant 6\,\text{s}$
, we additionally exchanged the Harten–Lax–van Leer contact flux approximation by the central Kurganov–Tadmor scheme (Kurganov & Tadmor Reference Kurganov and Tadmor2000; Panuelos, Wadsley & Kevlahan Reference Panuelos, Wadsley and Kevlahan2020) without the novel dissipation switch developed by Panuelos et al. (Reference Panuelos, Wadsley and Kevlahan2020). This should give the most dissipative Eulerian run on a Cartesian grid in GIZMO, and hopefully eliminate the seemingly sub-dissipative behaviour observed before. Evidently from figure 8(c), the overall dissipation increases from the beginning, but the sub-dissipative behaviour persists while it is shifted to
$t\gt 9\,\text{s}$
. Again, the explicit SFS model does not qualitatively improve this situation. To rationalise these effects, the magnitude of the velocity fields from the least dissipative Lagrangian run (figure 8
d) and the most dissipative Eulerian run (figure 8
e) can be compared at
$t=14\,\text{s}$
, where incompressible turbulence should be present. Apparently, taking all dissipation measures into account, the Eulerian velocity field is much smoother than the Lagrangian reference, even showing a laminarisation of the dynamics in the shear flow planes. This is surprising, given that the averaged kinetic energy levels exceed the Lagrangian ones after
$t\gt 9\,\text{s}$
(figure 8
c). Actually, a decrease of the (effective) Reynolds number should result in a monotonic decrease in the averaged kinetic energy levels (Fehn et al. Reference Fehn, Kronbichler, Munch and Wall2022). Hence kinetic energy must be introduced by another mechanism. These are the constructively interfering acoustic waves mentioned above, which manifest in the corresponding density field in figure 8(f). Although their magnitude is reasonable in terms of the weakly compressible approach taking all dissipation measures into account, they are the dominant, unphysical feature in the density field. Unfortunately, they also introduce oscillatory behaviour in the velocity field (figure 8
e) and are the root of instability in figure 8(a) and the seemingly sub-dissipative behaviour in figure 8(b,c).
From these observations, we extract three main conclusions.
-
(i) To perform implicit LES of incompressible turbulence with the MLS-SPH-ALE methods of Hopkins (Reference Hopkins2015), the Lagrangian reference frame with its Lagrangian noise and the related implicit SFS seem inevitable. In the Eulerian frame, the numerical schemes lack a direct dissipation mechanism for the detrimental acoustic waves.
-
(ii) If an implicit LES in an Eulerian reference frame were possible, then we would expect that the incompatibility of classical eddy viscosity models and current SPH methods persists as it is rooted in the non-local discretisation.
-
(iii) The
$\sigma$
model by Nicoud et al. (Reference Nicoud, Toda, Cabrit, Bose and Lee2011) does not resolve the acoustic wave issue due to a lack of awareness. This could also be interpreted as a sort of incompatibility, but is different from the non-local incompatibility to which we refer in our work.
Open access
