Smoothed particle hydrodynamics method from a large eddy simulation perspective. Generalization to a quasi-Lagrangian model

The present work deals with some recent developments regarding the inclusion of the Large-Eddy Simulation (LES) in the weakly compress-ible Smoothed Particle Hydrodynamics (SPH) framework. Previously {see the work of Di Mascio et al. [Phys. Fluids 29 , 4 (2017)]}, this goal was achieved by applying a Lagrangian ﬁlter to the Navier–Stokes equations for compressible ﬂuids and, then, approximating the differential operators in a SPH fashion. Since the Lagrangian nature of the derived scheme turned out to be an obstacle for accurate simulations of high Reynolds number problems, the above approach is here modiﬁed to obtain a quasi-Lagrangian LES-SPH model. This relies on the addition of a small velocity deviation to the actual Lagrangian velocity based on the particle shifting technique and on the inclusion of the tensile instability control technique for eliminating the onset of the tensile instability in the ﬂuid regions characterized by large vorticity and negative pressure. The proposed model is successfully tested in both two-dimensional and three-dimensional frameworks by simulating the evolution of freely decaying turbulence problems and comparing the outputs with the available theoretical results and solutions from other numerical models.


I. INTRODUCTION
In the past years, an increasing number of studies have been dedicated to the extension of the Smoothed Particle Hydrodynamics (SPH) to the modeling of turbulent flows (see the works of Lo and Shao, 1 Dalrymple and Rogers, 2 Price, 3 and Mayrhofer et al. 4 ).This is a fundamental step for the development of a numerical method that aims at simulating engineering applications but, at the same time, it represents a great challenge in terms of both theoretical and numerical issues.
About the former aspect, the most straightforward approach to turbulence in the SPH framework seems to be the one at the basis of the Large-Eddy Simulations (LESs).In fact, this method relies on a filtering of the Navier-Stokes equations that resemble the one adopted for the derivation of the smoothed differential operators of the SPH.This resemblance was first noted by Bicknell 5 where the development of a LES model in the SPH framework was highlighted as a fundamental step forward for the SPH scheme.A preliminary effort in such a direction was made in the work of Di Mascio et al. 6 where a consistent Lagrangian LES-SPH scheme was proposed.This was obtained by applying a Lagrangian filter to the Navier-Stokes equations for compressible fluids and, then, approximating the differential operators in a SPH fashion.The model structure was similar to a general weakly compressible SPH scheme, apart from the presence of additional terms in both the continuity and momentum equations coming from the LES filtering procedure.
The model proposed by Di Mascio et al., 6 named δLES-SPH, performs well for moderate Reynolds number flows but, for larger values, it exhibits some drawbacks in those regions characterized by strong vorticity and negative values of the pressure field.In fact, the negative pressure leads to the onset of the so-called tensile instability (see the work of Swegle et al. 7 ), i.e., to the generation of large spurious voids in the fluid domain.In weakly compressible SPH schemes, this problem is usually avoided by adding a background pressure (BP) in the state equation to guarantee that the pressure field is positive everywhere.Unfortunately, in large Reynolds number flows, the background pressure induces a non-physical particle resettlement and a consequent generation of spurious high-frequency noise.The latter heavily affects the kinetic energy spectrum, causing a substantial deviation from the theoretical rate of dissipation in freely decaying turbulence simulations.
The aim of the present work is, therefore, to propose an improved δLES-SPH model that is able to simulate high Reynolds number flows without the onset of tensile instability and the use of the background pressure.This is achieved in two steps.The first step relies on the definition of a quasi-Lagrangian large-eddysimulation model (which extends the approach of Di Mascio et al. 6 ) where a small velocity deviation is added to the actual fluid velocity.This strategy allows us to cast the LES in the framework of the most advanced SPH schemes, which rely on a quasi-Lagrangian motion of the fluid particles (see the works of Xu et al. 8 and Nestor et al. 9 ).The second step is achieved by rearranging the LES equations in the SPH formalism.In particular, the velocity deviation is modeled through the Particle Shifting Technique (PST), similar to the δplus-SPH scheme derived by Sun et al. 10 This, along with the use of Tensile Instability Control (TIC) described by Sun et al., 11 allows for the definition of a stable and robust numerical scheme, which is suited to the simulations of high Reynolds number flows typical of the LES approach.
In the proposed δLES-SPH scheme, the presence of the velocity deviation with respect to the actual fluid velocity leads to the appearance of additional terms in the continuity and momentum equations that need proper turbulence closures.Specifically, the term in the momentum equation is represented through a classical LES closure, while that in the continuity equation is here modeled through the diffusive term of the δ-SPH scheme (see the work of Antuono et al. 12 ).This paper is organized as follows: Sec.I introduces the theoretical framework of the quasi-Lagrangian LES-SPH model and the closures adopted for the additional terms, Sec.III describes the derivation of the proposed δLES-SPH numerical scheme from the theoretical model, and, finally, Sec.IV displays the applications.Specifically, the proposed model is tested in both 2D and 3D frameworks by considering the evolution of a fluid at rest, which is initially put into motion by an external body force and, then, left free to evolve, generating freely decaying turbulence.The forcing term used in the above problem mimics the formation of the vortex patches of the Taylor solution (see the work of Taylor 13 ).The results obtained with the proposed model are compared with the solutions obtained through a mesh-based finite volume method.

II. QUASI-LAGRANGIAN LES-SPH SCHEME
In this section, we describe the derivation of the quasi-Lagrangian LES-SPH model.Let us consider the following Navier-Stokes equations for a barotropic weakly compressible Newtonian fluid: where u is the flow velocity, p and ρ denote the pressure and density fields, respectively, and F represents the state equation.The hypothesis that the fluid is weakly compressible corresponds to the following requirement: where δp indicates the variation of the pressure field and c = c(ρ) is the sound speed (see the works of Monaghan 14 and Marrone et al. 15 ).The viscosity coefficients ν and λ ′ indicate the ratios between the Lamé constants μ and λ and the density ρ.Since the fluid is weakly compressible, ν and λ ′ are assumed constant.Now, let us define a generic filter in R 3 × R as follows: The above filter is supposed to have a compact support to depend only on ∥xp(t) − y∥ and |t − τ| and to be an even function with respect to both arguments.Here, xp(t) indicates the position of a quasi-Lagrangian point that moves in the fluid domain according to the following equation: where δũ is a (small) arbitrary velocity deviation (to be specified later), while ũ is given by the following definition: Hereinafter, we refer to xp and ũ as the filtered position and velocity, respectively.Accordingly, the main idea is to rewrite system (2) in terms of the filtered quantities and obtain a quasi-Lagrangian LES scheme.With respect to this point, we observe that, since the state equation is generally nonlinear, F(ρ) is different from F(ρ) and, consequently, the filtering procedure cannot be applied to both pressure and density.To avoid inconsistency, when we refer to filtered pressure, we mean Under this hypothesis, we apply the filter in (3) to the Navier-Stokes equations for weakly compressible flows and, integrating over where the total time derivatives d/dt are obtained with respect to the velocity ũ + δũ.The tensor T ℓ is given by and is equivalent to the sub-grid stress tensor.Following the derivation shown by Di Mascio et al., 6 we now rearrange system (7) in the framework of the Smoothed Particle Hydrodynamics (SPH) scheme.In particular, we highlight that the latter scheme is based on a smoothing procedure that somehow recalls the one used to obtain the Lagrangian LES.The main differences between these filtering procedures are that SPH approaches adopt a filter ("kernel" in the SPH terminology) that depends only on the spatial variables (xp − y); furthermore, the smoothing procedure is not applied to the Navier-Stokes equations, but to the differential operators, that are replaced by their smoothed counterpart (see the work of Colagrossi et al. 16,17 ).It is possible, nevertheless, to reinterpret the Lagrangian LES through the SPH approach.To this purpose, we split the filter ϕ into where W indicates the SPH kernel.In SPH notation, the smoothing procedure of a generic scalar field f is indicated as Here, a time filtering alone is introduced by the symbol Using the above definitions, it is easy to prove that f ≙ ⟨ f ⟩ and that, generally, the time and space filters do not commute, i.e., ⟨ f ⟩ ≠ ⟨ f ⟩ [see the work of Di Mascio et al. 6 ].Since the time filter is the inner one, the overall LES-SPH scheme may be regarded as a spatial Lagrangian filter applied to a set of time-averaged variables.In this sense, the time filter may be thought as an implicit filter whose presence is accounted for through the modeling of the additional terms.We stress here that the dependence on the time of the filter ϕ is fundamental for the derivation of the model (see Appendix A for details).Furthermore, it is conceptually well-posed, since turbulence is characterized by both time and space fluctuations.Now, suppose that we want to model a high Reynolds number flow for which LES filtering is required.Then, we need the filtered variables ũ, p, ρ for each fluid particle at positions xp; at the same time, we want to approximate the operators in Eq. (7) in the SPH fashion.Using the above definitions and the fact that the spatial filter and the gradient commute, we can rearrange the gradient of f as follows: where the first term in the right-hand side is the SPH operator, while the latter term accounts for small scale "fluctuations" in space, hereinafter denoted through f ′ ≙ f − f .For confined flows, the noncommutability of filtering and differentiation must be taken into account for a rigorous extension of the filtering close to the boundaries (i.e., in those points whose distance from the boundaries is smaller than the kernel radius).The above procedure can be applied to all the remaining operators.By doing so, in SPH formalism, system (7) reads where Here, C 1 and M 1 come from the SPH approximation procedure and require a SPH closure, whereas C 2 and M 2 include all terms from the Lagrangian LES and require a LES closure.Incidentally, we note that the term C 2 is not present in LES models for compressible flows where Favre-filtered variables are adopted (see the work of Moin et al. 18 ).Finally, C 3 and M 3 come from the use of the generic deviation velocity δũ and are new terms in comparison to the model defined by Di Mascio et al. 6 Incidentally, we highlight that system (12) shows close similarities with the ALE (Arbitrary Lagrangian-Eulerian) models (see the work of Oger et al. 19 ).As described by Antuono et al., 20 the same peculiarities are shared by a class of SPH models based on the use of the shifting algorithm.

A. Closures
Here, we propose some possible closures for the terms C 1 , C 2 , and C 3 in the continuity equation and for the terms M 1 , M 2 , and M 3 in the momentum equation.
First, assuming that the filtering procedure is performed far from the boundary of the fluid domain and using the symmetry ARTICLE scitation.org/journal/phfproperties of the kernels W and θ, we derive the following relations (see the work of Di Mascio et al. 6 ): where the length and time scales σ and t only depend on the kernels W and θ and are defined as follows: Combining the above relations, it immediately follows We highlight that the dependence of f on the time scale is inspected here for the first time, while the work of Di Mascio et al., 6 only the dependence on the spatial scales was studied.Now, choosing a reference scale f 0 for the scalar function f and a reference length, x 0 , and time, t 0 , we can rewrite the above expression as where the starred variables indicate dimensionless values (e.g., σ * = σ/x 0 and t * = t/t 0 ).Physically, (x 0 /t 0 ) = U 0 , where U 0 indicates the reference fluid velocity.Then, choosing the kernel θ such that σ/t = c (we recall that c is the sound velocity), we obtain where Ma = U 0 /c is the Mach number.Under the hypothesis that the flow is weakly compressible [see Eq. ( 2)], we can assume Ma ≪ 1 and simplify Eq. ( 19) as follows: Note that the condition Ma ≪ 1 leads to t ≪ σU 0 , which implies that the kernel θ is very narrow.This means that we can generally neglect the action of the time filter in the terms {Ci, Mi}.
More in depth, if we require that the term O(M 2 a σ 2 * ) is consistent with the order of approximation of Eq. ( 21) (which is fourth order in σ * ), we obtain a specific requirement on the Mach number, namely, Ma ≙ O(σ * ).From a different point of view (also going back to dimensional variables), this corresponds to σ ≙ O(Max 0 ) or, equivalently, to where x 0 can be regarded as the smallest spatial scale that we expect to solve in the present model, while shorter lengths are modeled through the LES closures.Taking this into consideration, we assume Ma ≙ O(σ * ) hereinafter and simplify the expressions for the additional terms in the continuity and momentum equation accordingly.
Under the hypotheses described above, the weakly compressibility assumption implies that the gradients of the density field are of order O(M 2 a ) ≙ O(σ 2 * ).If we additionally require δu 0 /U 0 ≙ O(σ * ), where δu 0 denotes the order of magnitude of the deviation δũ, it is possible to prove that C * 3 , M * 3 ≙ O(σ 5 * ) (see Appendix B), and therefore, they can be neglected.
About C 1 and M 1 , these may be regarded as a sort of defiltering of the leading order SPH operators.In fact, using the expansions in (18) and the assumptions for weakly compressible flows, it immediately follows In the work of Di Mascio et al., 6 where q 2 represents the turbulent kinetic energy, ν T is the turbulent kinetic viscosity, and D is the strain-rate tensor, that is, D ≙ (∇ũ + ∇ũ T )/2.Similar to the work of Yoshizawa, 21 we assume where where ν δ is assumed to be a function of D. In agreement with the usual approaches adopted in the LES framework, this is equivalent to model C 2 as a diffusive term.In the present work, we adopt a closure that is expected to be more effective on smaller spatial scales.In particular, this is defined through the difference between the actual (i.e., Lagrangian) density field and its filtered value as follows: Using Eq. ( 22), we obtain and, consequently, the following closure for C 2 : In the SPH framework, a simple closure for C 2 is obtained by using the diffusive term proposed by Antuono et al. 12 for the δ-SPH scheme, since this contains fourth-order spatial derivatives of the density field.In particular, far from the boundary of the fluid domain, such a term reduces to a bi-Laplacian of the density field (see Appendix C).In Sec.III, we introduce the numerical model of the proposed LES-SPH model and describe the modeling of each additional term in detail.

III. THE δLES-SPH SCHEME
To write system (7) in the discrete formalism, we rely on the work of Sun et al. 10 where the additional δũ-terms are included in the SPH framework in a consistent way.Here, the parameter σ is replaced with ℓ, where the latter represents a reference length for the spatial filter (to be specified later) and ℓ ≙ O(σ) (to be consistent with the theoretical framework).In particular, following the work of Sun et al. 10 for the discretization of the differential operators, we obtain where mi is the ith particle mass (assumed to be constant) and Vi is its volume, while Wij ≙ W(xi − xj).Here, Pij represents the argument of the pressure gradient.Following the work of Sun et al., 11 the Tensile Instability Control (TIC) is added to the present model to avoid the onset the tensile instability when large negative pressure regions arise in the fluid domain (for example, inside the cores of strong vortical structures).This corresponds to a switch from the "plus" formulation (namely, Pij ≙ pj + pi) to the "minus" form (that is, Pij ≙ pj − pi) when the pressure pi is negative, namely, Pij ≙ pj +|pi|.More details can be found in the work of Sun et al. 11 The coefficient of the viscous term is K = 2(n + 2), where n is the number of spatial dimensions, while its arguments are In 2D, CS is set equal to 0.12, as suggested by Lo and Shao, 1 whereas for 3D problems, it is set equal to 0.18. 22The first contribution in the expression of αij represents the actual fluid viscosity (see the works of Monaghan 23 and Colagrossi et al. 17 ), while the latter one is the LES closure for turbulence.Note that all contributions related to the fluid compressibility have been neglected, since they are negligible in comparison to the leading order stress tensor.The symbol ψij is the argument of the diffusive term in the work of Antuono et al., 12 namely, Here, the superscript L indicates that the gradient is evaluated through the renormalized gradient formula 24 as follows: Incidentally, we observe that the right-hand side of the expression (32) is of order O(σ 2 ) (see the work of Antuono et al. 12 for more details) and, therefore, it is consistent with the closures made in Sec.II A. As done by Di Mascio et al. 6 and by Meringolo et al., 25 the diffusive term is not multiplied by an external parameter, but the latter is included directly inside the summation and modeled as a viscous-like coefficient following a standard LES approach.In particular, we choose in which C δ is a dimensionless constant.This has been set equal to 6 after a tuning analysis in order that the term C 2 has approximately the same order of magnitude of the diffusive term in the δ-SPH model.
As usually done for weakly compressible fluids, the state equation is linearized around the reference density ρ 0 , leading to pi ≙ c 2 0 (ρi where c 0 = c(ρ 0 ) is a numerical sound speed that satisfies the following requirement: Here, Umax and δpmax indicate the maximum fluid velocity and the maximum pressure variation expected during the simulation.The above inequality allows the density variations to be below 1% (see the work of Monaghan 23 ).
For what concerns the velocity deviation δûi, this is defined following the work of Sun et al. 10 as follows: where and Δx is the initial mean particle distance.Here, the constants χ and n are set equal to 0.2 and 4, respectively.The expression in (37) has to be further corrected for particles close to the domain boundaries to avoid a non-physical particle motion.For further details, refer to the work of Sun et al. 10 31) is integrated in time by using a fourth-order Runge-Kutta scheme with the Courant-Friedrichs-Lewy number equal to 1.5 (see the work of Meringolo et al. 25 for more details).

IV. APPLICATIONS
In this section, we apply the proposed scheme to the study of freely decaying turbulence in two and three dimensions.In all the simulations, we adopt a C 2 Wendland kernel, as described by Wendland. 26Consistent with the constraint in Eq. ( 36), we choose Ma = 0.1.Under this hypothesis, the length x 0 described in Eq. ( 23) is of the same order of the diameter of the spatial kernel.Specifically, x 0 /D = 1.31 and x 0 /D = 1.58 in two-and threedimensions, respectively, where D is the diameter of the Wendland kernel.
According to the space dimension, we consider bi-and threeperiodic squared domains containing a fluid, which is initially at rest.The fluid is then put to motion by applying an external forcing term during a finite time window.This is realized by multiplying an external body force (which changes according to the space dimension) for the following time ramp: for t * ∈ ∥0.1, 0.9), 10(1 − t * ) for t * ∈ ∥0.9, 1), where t * = tU/L is the dimensionless time.Here, U is the maximum velocity attained during the forcing stage and L is the length of the squared domain.After the forcing term turns off (namely, for t * > 1), the flow rapidly becomes unstable and leads to the onset of freely decaying turbulence.Hereinafter, the Reynolds number is indicated by Re = UL/ν, where ν is the kinematic viscosity.
Before proceeding, we clarify some fundamental aspects of the above approach.Basically, the use of external forces to put into motion a still fluid is motivated by the necessity of a careful inspection of the modeled and resolved kinetic energies in the numerical scheme.Indeed, the assignment of an initial fluid field would immediately correspond to an assignment on both the modeled and resolved components of the kinetic energies and, consequently, this would make the analysis of the energy balance of the model rather difficult.Incidentally, we highlight that in both two and three dimensions, the adopted external body forces are divergencefree.

A. Two dimensions
In two dimensions, we consider Re = 10 4 , 10 5 , and 10 6 .The spatial filter length ℓ is set equal to the initial particle distance, hereinafter denoted by Δx, and the kernel radius is R = 4Δx.The forcing term is defined by using the two-dimensional solution described by Taylor, 13 namely, F 2d ≙ r(t * )A(sin(8πx * ) cos(8πỹ * ), − cos(8πx * ) sin(8πỹ * )), (40) where x * ≙ x/L and A = 1.3U 2 /L.The time histories of the kinetic energy at the different Reynolds numbers are reported in the left panel of Fig. 1, where it is possible to observe the action of the forcing term during the early stages (i.e., for t * ≤ 1) and subsequent free decay for t * > 1.The right panel of the same figure displays the horizontal velocity field at t * = 1 for Re = 10 6 , showing that the forcing term defined in (40) induces a pattern of eight vortex cores in the squared domain.In this latter case, the spatial resolution is L/Δx = 1200.
Before showing the results of the proposed model, we briefly describe the reasons that led us to define the present quasi-Lagrangian LES-SPH scheme.When we use a Lagrangian LES-SPH model as that described by Di Mascio et al., 6 the simulation of large Reynolds number problems represents a great challenge.Specifically, the occurrence of regions with large negative pressure values leads to the onset of the so-called tensile instability (see Swegle et al. 7 ) and to a consequent generation of voids inside the fluid domain.In the SPH framework (also for confined flows), this issue is usually avoided by introducing a background pressure so that, during the whole evolution, the pressure field always maintains positive.
Unfortunately, this strategy induces short-length oscillations (due to particle resettlement) that affect both the pressure and velocity fields and that lead the energy spectrum to deviate from the theoretical decay.This behavior is briefly summarized in Fig. 2 where the left panel displays a snapshot of the kinetic energy field at tU/L = 3 and the right panel shows the corresponding energy spectrum, hereinafter denoted by E k , scaled by E 0 = ρ 0 U 2 L 3 (details about the computation of the SPH energy spectrum are given in Appendix E).The deviation from the theoretical rate of decay occurs at wavelengths that are comparable with the kernel radius, heuristically confirming the particle resettlement (in Fig. 2, the wave number associated with the kernel radius is indicated by kR, while k f denotes the wave number of the forcing term).In this case, the use of the Tensile Instability Control (TIC) instead of the background pressure leads to even worse results because the particle distribution becomes disordered and this strongly affects the pressure field (see the left panel of Fig. 3).
The results shown above convinced us about the necessity of including the TIC along with the Particle Shifting Technique (PST) described by Sun et al. 10 in the δ-LES-SPH scheme in order to get rid of the tensile instability and, at the same time, to obtain a regular particle distribution and a consequent reduction of the spurious noise.This led to the derivation of the quasi-Lagrangian scheme described in Secs.II and III.Incidentally, we observe that the use of the PST alone (namely, the implementation of the present model without the addition of TIC) is not enough to avoid the onset tensile instability.About this point, the right panel of Fig. 3 shows the generation of void regions inside the fluid domain in correspondence with the negative pressure regions.
The behavior of the proposed δ-LES-SPH model is summarized in Fig. 4. In the left panel, the energy spectrum for Re = 10 6 , L/Δx = 1200, and tU/L = 3 is compared with the results obtained through a finite-volume scheme (details are reported in Appendix D), showing a good agreement between the numerical outputs.In this panel, we also drew the wave numbers related to the forcing term (i.e., k f ) to the kernel radius (that is, kR) and the positions of the wave numbers k λ and kD related to the Taylor and Kolmogorov scales, that is, where ϵ is the rate of energy dissipation and E is the kinetic energy (the flow is assumed to be fully turbulent).Generally, the length scales that are larger than the Taylor microscale are not strongly affected by viscosity while, on the contrary, length scales smaller than the Kolmogorov microscale are completely dissipated by viscous effects.In particular, these scales are used to identify the inertial range of turbulence.The fact that the Kolmogorov scale is smaller than the minimum simulated length and the Taylor scale is close to the length of the forcing term confirms that the LES simulation is correctly applied.The right panel of the same figure displays the kinetic energy field at the same instant (namely, tU/L = 3).Finally, the pressure and vorticity fields are shown in Fig. 5, proving that both are free from spurious short-length noise.The behavior of the energy spectra for Re = 10 4 and 10 5 is reported in Fig. 6.In the latter case, the numerical output is in good agreement with the direct and inverse cascades for turbulence in 2D, while a slight deviation from the theoretical profile is observed for the case at Re = 10 4 close to kR.

B. Three dimensions
In three dimensions, we consider the problems with Re = 10 6 and ℓ is set equal to the kernel radius, which, in this case, is set equal to 3Δx.The value of ℓ changes in 3D as a consequence of the change in the ratio R/Δx between 2D and 3D simulations (here, R indicates the kernel radius).In fact, in 2D, we usually adopt a value of R/Δx that guarantees about 50 neighboring particles in the kernel support.In 3D, we are bound to reduce R/Δx in order to maintain the neighbor particles to about 100 (otherwise, using the same value as in 2D, we would obtain about 260 neighboring particles in 3D, leading to a huge increase in the computational costs).In turn, the use of a smaller value of R/Δx in 3D makes the simulations slightly more noisy than 2D ones.Then, to recover similar qualitative behavior, we use a larger value of ℓ.
The forcing term is defined through the three-dimensional solution described by Antuono. 27In particular, we write where ua,t is the time derivative of the analytical velocity field (evaluated at t = 0) and A * = 1.25 ⋅ 10 −3 (here, A * is dimensionless).For the sake of brevity, we do not include the explicit formula, and refer to the work of Antuono. 27The above-mentioned solution (which can be regarded as an extension in three dimensions of that by Taylor 13 ) is tri-periodic and fully three-dimensional.This implies that shorter times are required to reach a fully developed homogeneous and isotropic turbulent flow in comparison to those works that adopt an initial three-dimensional condition characterized by a null velocity component in one direction (see the work of Goldstein, 28 Orszag, 29 and Sharma and Sengupta 30 ).The forcing term in Eq. ( 42) generates a tri-periodic motion with kx = ky = kz = 4π/L, where L indicates the length of the squared domain and κ = (kx, ky, kz) is the wave number.This implies that two vortex modules are considered in each coordinate direction.In the left plot of Fig. 7, we display a snapshot of the L 2 -norm of the vorticity at t * = 1 (that is, when the ramp is switched off).As described by Antuono, 27 the Y-shaped structures (i.e., the isocontours at ||ω|| = 3L/U) enclose the stagnation points.The spatial resolution is L/Δx = 128.
As shown in the left plot of Fig. 8, the finite volume code and the present δLES-SPH model show a different transition to homogeneous turbulence because of their different structures.Indeed, the resolution L/Δx = 128 here adopted is rather low with respect to those attained in the 2D test cases.As a consequence, the flow evolution is significantly affected also by the inherent numerical diffusion of the two models.In this sense, the ways in which numerical diffusion acts are driven by the substantial difference between the two codes: to mention just the most relevant, the SPH is a quasi-Lagrangian solver formally using a second-order spatial differential operator (which can locally degrade down to the zeroth order), whereas the FV is a pure Eulerian scheme formally employing a fifthorder spatial operator.It is, therefore, reasonable that these aspects play an important role in the specific development of the initial flow instability and on the following evolution.Thus, we compare their behavior when the amount of the kinetic energy is approximately the same in both the simulations.In particular, we compare the output of the finite volume code at t * = 4.5 and that of the δLES-SPH model at t * = 6 (see the left plot of Fig. 8).As displayed in the right panel of the same figure, notwithstanding the different time histories of the simulations, the energy spectra of both solutions have a quite similar behavior and, furthermore, are in good agreement with the theoretical rate of decay.In this figure, we also show the wave numbers of the external forcing (namely, k f ) and the kernel radius (i.e., kR).
Finally, in Fig. 9, we plot the L 2 -norm of the vorticity as predicted by the δLES-SPH model (left panel) and by the finite volume solver (right panel).Both plots show a similar distribution and intensity of ||ω||.This is further confirmed by the isocontour at ||ω|| = 15L/U of the vorticity displayed in Fig. 10 for both the numerical schemes.Remarkably, the SPH vorticity field appears free of high-frequency noise and coherent vortex structures are clearly visible, similar to the two-dimensional case.

V. CONCLUSIONS
An improved version of the δ-LES-SPH model described by Di Mascio et al. 6 is derived based on the following steps: (i) the definition of a quasi-Lagrangian LES-SPH model, (ii) a novel closure for the term C 2 based on the use of a fourth-order diffusive term, (iii) the modeling of the velocity deviations from the Lagrangian velocity through the PST proposed by Sun et al., 10 and (iv) the implementation of the tensile instability control term.
The proposed model is tested in two and three dimensions by simulating initially forced motions followed by freely decaying turbulence.These test cases resemble the generation of vortices similar to those predicted by the Taylor solution. 13Thanks to the straightforward inclusion of the particle shifting technique and the tensile instability control, the proposed model overcomes some of the wellknown issues that affect SPH simulations of high Reynolds number flows, such as the generation of high-frequency spurious noise and the onset of the tensile instability.The comparisons of the energy spectra with the theoretical rate of decay and with results by a finite volume solver confirm the accuracy and reliability of the proposed quasi-Lagrangian scheme.
In future works, further efforts will be dedicated to (i) the validation of the proposed δ-LES-SPH model against experimental or numerical data for turbulent flows characterized by Reynolds numbers larger than 10 6 , (ii) the definition of proper wall-functions to deal with solid walls, and (iii) the extension of the scheme to represent free-surface flows.
Incidentally, it is worth mentioning that the present model may also be extended to the Riemann-ALE SPH variant (see the work of Oger et al. 19 ).However, this extension could unlikely result in a significant improvement, as the diffusive term of the δ-SPH has been shown to be equivalent to a simplified Riemann solver (see the work of Green et al. 31 ) and has similar properties regarding numerical diffusion and accuracy of the differential operators (see also the recent work by Hammani et al. 32 where δ-SPH is compared to a Riemann solver).Those are the main aspects to be tackled to enhance the present scheme and, in this sense, a higher-order approach could improve the present SPH results.details of the derivation.In brief, substituting the expression of the total derivative of the velocity [see Eq. ( 1)] in the first integral and using integration by parts, we obtain For what concerns the second integral in the right-hand side of (A4), we find where the tensor T ℓ is given by Finally, the last integral in the right-hand side of Eq. ( A4) is new and is equivalent to Collecting all the results together, we obtain For what concerns the filtered density, we find Substituting and integrating by parts, we find which, after some algebra, can be recast in the following form: Here, we prove that C 3 is of order O(σ 5 ).The same outcome is obtained for M 3 by following the same proof with minor changes.For these reasons and for the sake of brevity, we omit the details of this latter term and only consider C 3 .First, using the dimensionless variables as in Sec.II A, we write Using the properties of the spatial filter, we expand the above equation as follows: (B2) Expanding Eq. (B2), we find out that the terms that do not contain derivatives of ρ cancel out.All the remaining terms are of order O(σ 4 * δu 0 /U 0 ), since the derivatives of ρ * are of order O(M 2 a ) ≙ O(σ 2 * ) because of the weak-compressibility assumption.Then, the expression (B2) simplifies as follows: Then, choosing δũ in such a way that δu 0 /U 0 ≙ O(σ * ), we finally obtain This confirms that the term C 3 can be straightforwardly neglected in the proposed model.
APPENDIX C: THE DIFFUSIVE TERM OF ANTUONO ET AL. 12 As described by Antuono et al., 12 if we neglect the presence of the boundary of the fluid domain, the diffusive term of the δ-SPH model can be expressed as follows: q k q l qmqn q ∂W ∂q dVy + O(σ 4 ), (C1) where q = (y − x)/σ and the subscripts k, l, m, and n indicate the vector components and q = ||q||.We observe that ∂W/∂q is a radial negative kernel with a compact support.Because of the isotropy of the integral tensor in the Eq.(C1), we can rearrange the expression as follows: where Note that α and β do not depend on the specific choice of the indices k and l because of the isotropy (symmetry) of the integral.As proved by Violeau and Fonty 33 (see Appendix A of that paper), α = 3β and, consequently, we can rearrange the expression (C2) as follows: where Δ 2 indicate the bi-Laplacian, namely, Δ 2 ρ ≙ Δ(Δρ).

APPENDIX D: FINITE VOLUME CODE
The finite-volume code used for the simulations is a general purpose in-house code that contains several physical models and numerical discretizations.For all the simulations reported in the present paper, the weakly compressible approximation of the Navier-Stokes equations (the same adopted for the SPH simulations) was used.Spatial discretization was performed by means of the classical fifth order WENO approach by Jiang and Shu 34 for the convection and pressure parts of the equations, while a standard central second order approximation is adopted for the viscous terms.Time integration was performed by means of the third order Total Variation Diminishing (TVD) Runge-Kutta scheme developed by Gottlieb and Shu. 35As to LES modeling, the standard model by Smagorinsky 22 was used for consistency with the analogous SPH simulations.Examples of application of the code to complex turbulent flow simulations can be found in the works of Muscari et al. 36 and Magionesi et al. 37

APPENDIX E: ENERGY SPECTRA OF SPH FIELDS
In order to evaluate the kinetic energy spectrum of the SPH results, the procedure described by Shi et al. 38 has been adopted.Specifically, the magnitude of the velocity fields has been interpolated on a uniform Cartesian mesh with spacing equal to the particle size Δx.For the interpolation of the SPH scattered data, a Moving Least Square (MLS) technique (see the work of Fries and Matthies 39 ) has been used, which in the work of Shi et al. 38 is shown, through numerical experiments, to be the most appropriate for extracting energy spectra.The MLS interpolation has been performed using about 50 particles in the kernel support for both two-and threedimensional data.This corresponds to a kernel radius equal to 4Δx and 2.3Δx, respectively, for 2D and 3D.Then, in order to obtain the kinetic energy spectra, the same FFT algorithm used for the finite volume solution has been applied to the SPH interpolated data.

FIG. 8 .FIG. 10 .
FIG. 8. Left panel: time histories of the kinetic energy as predicted by the finite volume scheme and the δLES-SPH model.Right panel: energy spectra.Symbols k f and k R indicate the wave numbers of the external forcing and the kernel radius, respectively.
6hese terms were directly included inside the SPH differential operators as de-filtering contributions of the main filtered field.In any case, their influence is generally negligible in comparison with C 2 and M 2 and, consequently, they are dropped hereinafter.Finally, C 2 and M 2 represent typical terms coming from the Lagrangian LES-filtering procedure and require specific LES closures.In particular, M 2 is modeled as done by Di Mascio et al.,6 The dimensionless parameters CY and CS are called the Yoshizawa and Smagorinsky constants, respectively (see Ref. 21).For what concerns C 2 , the closure proposed by Di Mascio et al.