Improved δ -SPH Scheme with Automatic and Adaptive Numerical Dissipation

: In this work we present a δ -Smoothed Particle Hydrodynamics (SPH) scheme for weakly compressible ﬂows with automatic adaptive numerical dissipation. The resulting scheme is a meshless self-adaptive method, in which the introduced artiﬁcial dissipation is designed to increase the dissipation in zones where the ﬂow is under-resolved by the numerical scheme, and to decrease it where dissipation is not required. The accuracy and robustness of the proposed methodology is tested by solving several numerical examples. Using the proposed scheme, we are able to recover the theoretical decay of kinetic energy, even where the ﬂow is under-resolved in very coarse particle discretizations. Moreover, compared with the original δ -SPH scheme, the proposed method reduces the number of problem-dependent parameters. consistency and accuracy of the presented approach has been studied by means of several two-dimensional benchmarks for different Reynolds numbers: Taylor–Green, lid-driven cavity and dam-break ﬂows. The results show that the proposed numerical method is able to simulate complex ﬂows. It alleviates the parameter dependency of δ -SPH methods and it is generally less dissipative than the δ -LES-SPH model. The obtained results suggest than the proposed approach could be a promising technique for LES computations of transitional ﬂows, where usual SGS models introduce excessive dissipation. This will be the subject of future research.


Introduction
Lagrangian methods have great interest in problems with large interface deformations or with strong dynamics. The Smoothed Particle Hydrodynamics (SPH) method is one of the most widely used Lagrangian methods for the simulation of fluid flows. It was first introduced for astrophysical applications [1], and nowadays it is applied to a wide variety of problems [2][3][4]. The SPH method has some interesting properties. For example, it easily handles free-surface flows, it has interesting conservation properties and presents very small numerical dissipation connected with convection [4].
Among the SPH methods, there are different ways of including numerical dissipation [3,5]. Thus, many SPH formulations use a term of artificial numerical dissipation to stabilize the schemes which may introduce excessive dissipation. A different possibility is using Riemann-solvers, which leads to a different family of methods [6][7][8][9][10][11][12][13].
In particular, this work focuses on the δ-SPH family models [14][15][16][17] which are a variant of SPH for weakly-compressible flows. Compared with standard SPH formulations, the main advantage of this scheme is the stabilization of the scheme avoiding the occurrence of a non-physical energy flux to high-frequency modes [15].

Governing Equations
In this work, we assume weakly compressible, Newtonian fluid flows in isothermal conditions. Under this hypothesis, the Navier-Stokes and displacement equations expressed in Lagrangian form read as where d(.) dt represents the material derivative following an infinitesimal fluid element (particle). ∇ is the nabla operator, ρ, p, u = (u, v), r, ν and g represent density, pressure, velocity vector, position vector, kinematic viscosity and the gravity acceleration, respectively.
Since the fluid is considered as weakly-compressible, the density fluctuations are supposed to be very small and a linear state equation is assumed [25] where c 0 is the speed of sound, and ρ 0 is the density along the free surface. In the case of confined fluids ρ 0 may represent the density in rest conditions. As a practical note about the numerical implementation, the fulfillment of the weakly-compressible hypothesis is attained through an adequate choice of the reference sound speed c 0 . Instead of using the physical sound velocity, it is imposed the following constraint on the speed of sound, in order to limit the variations of the density below the 1% [15] c 0 ≥ 10 max(U max , p max ρ 0 ) (2) where U max and p max are reference values of the velocity and pressure, which are case dependent. In practice, it is usual to estimate only the maximum velocity and not taken into account p max in Equation (2). This is the procedure that we have followed here.

Numerical Methods
SPH method discretizes the physical space into many discrete elements, usually called particles, without any connectivity among them. This method is based on the approximation of any physical scalar (or vector) field using a convolution. Numerically, it is performed by replacing the Dirac delta function with a regular smooth function, which is called kernel (W). This function must satisfy some conditions such as symmetry (even function), normalization, compactness of its support, among others. We refer the interested reader to [26] for more details. The kernel function used in this work is the C2-Wendland [27]. The kernel function depends on a parameter h, called the smoothing length, which defines the domain of influence of the kernel function. In this work, the smoothing length h is a constant which is chosen relative to the initial inter-particle distance δx 0 (h = 2δx 0 ) which allows the interpolation with about 50 neighboring particles when the kernel support is not truncated at boundaries. The initial particle volume is taken as V 0 = ∆x 0 d , where d is the space dimension number. The mass of each particle i is chosen to be constant and equal to m = ρ 0 V 0 during all the simulation time.

δ + -SPH and δ-SPH Methods
For the discretization of the governing equations system in the framework of weakly compressible SPH method, the δ + -SPH model [15,17] has been adopted. This model is characterized by the addition of a diffusive term to the continuity equation and particle shift δr i to the position calculated from governing displacement equation. The diffusive term is added to reduce the high-frequency noise to stabilize the density/pressure fields [14], while particle shift is added to overcome the so-called tensile instability and regularize the distribution of particles in order to achieve suitable numerical properties (such as accuracy, convergence and stability). The δ + -SPH model can be summarized as follows where the symbol ∇ i denotes the gradient operator applied to the particle i. The quantities ρ i , u i , p i , r i , r * i and V i represent the density, velocity vector, pressure, the advected position vector, the particle position vector after shifting and volume of the particle i, respectively. The quantities with sub-index j, represent the quantities related to the neighboring particle j of the particle i. Note that the sums in j are the sum in all the neighbors of particle i. The term Π µ i denotes the viscous acceleration vector. The notation W ij = W |r i − r j |, h denotes the kernel function.
The term D ij of the additional diffusive term is expressed as Here, ∇ρ L i and ∇ρ L j re-normalized density gradients [14] with respect to the particle i and its neighbor j, respectively.
where the re-normalization matrix L i is expressed as The dimensionless parameter δ in Equation (3) defines the magnitude of the density diffusion. In δ + -SPH and δ-SPH methods it is usually taken constant as δ = 0.1 [18].
In the momentum equation, the viscous part is defined as where the term π ij reads as and β is with K = 2(d + 2) and d is the number of spatial dimensions. Note that there is a term in β related to the molecular viscosity (µ) and another term related with numerical dissipation added for stability reasons. This artificial viscosity is controlled by the non-dimensional parameter α ij . In δ + -SPH and δ-SPH methods it is usually taken constant in the range [0.01-0.05]. On the other hand, the value of δ is kept constant with a value that depends on the problem. Taking δ = 0.1 is an usual choice. The main difference between δ + -SPH and δ-SPH methods is in the updating of the position of the particles. Thus, the δ + -SPH method uses the particle-shifting procedure as proposed in [15], where the position of the particles is corrected each time step in order to alleviate the tensile-instability phenomena. Thus, the position of the particles is updated as follows where CFL = c o ∆t h is the Courant-Friedrichs-Lewy number, Ma = U max c 0 is a reference Mach number and χ ij denotes the artificial pressure-like function described by Monaghan [28] and it is expressed as follows R and n are two constants set to be equal to 0.2 and 4, respectively. Note that particle shifting in the updating of the position of the particles is neglected (δr i = 0) in δ-SPH method [25].
In this work, the particle shifting technique is used only for confined flows cases, since for the free surface flow simulations, its application can accumulate the errors in time and causes non-physical behavior [29]. Thus, for free-surface simulations δr i = 0 and then the (δ + -SPH) scheme is reduced to the δ-SPH method [25].

δ-LES-SPH Method
In [18], the δ-LES-SPH method [16] is written in a similar fashion to the δ-SPH method. Using the same expressions that were presented in the previous section, we can define the δ-LES-SPH method as follows. First, instead of having a constant α ij , this parameter can vary locally for each pair of interacting particles as [16] where the parameter α i is usually defined using a subgrid-scale model [30,31]. For the δ-LES-SPH model presented in [18], the Smagorinsky SGS model is used, and thus α i can read as follows where ν T i is the kinematic eddy viscosity that can be expressed as C s = 0.12 is the constant of the Smagorinsky model [32] and I LES is a reference length which has been set to be equal to the kernel radius (I LES = 2h = 4∆x) [16]. Moreover D i = √ 2D i : D i is a rescaled Frobenius norm for the particle i, and the strain-rate tensor is computed as where the symbol ⊗ presents the outer product between two vectors. Similarly, in the δ-LES-SPH method, the value of δ ij is not constant, and it varies locally for each pair of interacting particles as well as α ij Here, the constant of the Smagorinsky SGS model is taken as C δ = 1.5. The variation of α i as well as δ i is limited to 0.2 in order to avoid stability issues [18].

A New δ-ADA-SPH Scheme
The Adaptive Dissipation Adjustment (ADA) method was recently developed as an implicit SGS model. It was applied for the computation of low mach flows [19] in a finite volume framework. The ADA method is based on the local Energy Ratio (ER) introduced by Tantikul and Domaradzki in the context of the Truncated Navier-Stokes (TNS) procedure [22]. In this work, we propose to extend its application to Lagrangian methods to eliminate the dependency of δ and α parameters. Thus, we propose to fix the value of δ and substitute the parameter (α ij ) in Equation (7) by the a new parameter which depends on the high-frequency content of the solution. The key idea of the proposed method is to link the magnitude of the dissipation part with the Energy Ratio (ER i ) of the particle i [21] In Equation (14), u i is the velocity field obtained as a result of the SPH computation. Moreover, u i and u i are two filtered velocity fields, obtained through filtering of u i using a low-pass filter with different widths. Thus, in Equation (14) u i and u i are computed using two different values of the filter width. Differently from what is performed in [19,22] where a top-hat filter is used, in this work we use kernel filtering. In particular, we use the Wendland-C2 kernel [27] with a cut-off radius set at 2h, where h is the smoothing length. Thus, the filtered velocity fields u i and u i will respectively read as and where h and h denote the smoothing lengths that define the Wendland-C2 kernel radius and they are evaluated respectively as h = h = 2∆x 0 and h = 2h = 4∆x 0 . As stated in Equation (14), ER i can be seen as a ratio of the spatial high-frequency components of the velocity field for two different filters. In this work, we use the ER i value to determine if the energy content of the solution is excessive or not. If it is, the numerical dissipation is reduced. This is achieved through the introduction of a multiplicative parameter in the viscous part of the momentum Equation (6). The value of this multiplicative parameter is automatically updated every time step depending on the value of ER i .
In order to automatically adjust the parameter, we define the value i associated to a particle i following the rule proposed in [19] Here, a value of φ = 0.001 is used to adjust the value of i continuously and gradually. The maximum value of i is determined following [33].
It is important to remark that we have chosen different values than those presented in [19] to define the range of ER i . The reason is that we have used different filters than those presented in previous works. However, it is important to note that it is possible to get similar results using another configuration of the filters if an adequate ER i interval is found, since the range of validity is completely dependent on the filter chosen [23]. Note that the values of these parameters are fixed for all the examples presented in this work, and have been calibrated to reproduce the decay kinetic turbulent energy in the two dimensional Taylor-Green vortex. Thus, since in our formulation α is substituted by (which is automatically adapted), we reduce the number of problem-dependent parameters compared with the base scheme.
The multiplicative parameter is defined as the mean value of We introduce this approach in the SPH scheme (3), by substituting the dimensionless parameter α ij in (8) by defined as in Equation (18).
Moreover, following [17], the value of δ is kept constant as δ = 0.1 for all simulation cases. This choice is justified since the range of variability of δ is narrow [33], and it has been shown that with this choice, the errors in energy conservation are negligible [18]. Moreover, it has also been shown in [18] that keeping this value of δ ensures that the dissipation varies with α. Thus, in the framework of the δ-ADA-SPH method, we ensure that the dissipation is totally controlled by .
To ensure the stability of the method, the time step (δt) must be chosen to fulfill the kinetic, body force [34] In this work, low storage 4th-order Runge-Kutta scheme [35] is used to integrate the governing equation system. The Courant-Friedrichs-Lewy constant CFL is set equal to 1.5.

Viscous Two-Dimensional Taylor-Green Flow
The first test case is the 2D Taylor-Green flow with Reynolds numbers Re = 100 and Re = 1000. This example is aimed to show that the proposed methodology is able to compute low Mach flows without adding excessive dissipation. The analytic solution of this test case is [36] u(x, y, t) = −Ue bt cos(2πx) sin(2πy), v(x, y, t) = Ue bt sin(2πx) cos(2πy) (21) where the reference velocity U is set to U = 1 and b = −8π 2 /Re is the decay rate of velocity field. The density is assumed to be the unity, and the initial pressure is given by In this test case, the velocity is initiated using t = 0 in Equation (21). This configuration leads to a periodic array of 2 vortices. We use a discretization with 2500 particles. The initial position of the particles is obtained by using the particle packing algorithm described in [37]. We use the value U max = 1 for computing Equation (2).
We compare the results using the proposed scheme with those obtained with the original δ + -SPH scheme, the δ-LES-SPH method, and a modified δ-LES-SPH method keeping a constant value of δ = 0.1. The value of β is set to match the molecular viscosity in the δ-SPH original scheme, and the term related with the added numerical dissipation is neglected (α ij = 0). In δ-LES-SPH and δ-ADA-SPH models, the term (α ij ) is added to the molecular viscosity, as in Equation (8). Figures 1 and 2 show the results for the value of the viscous parameters ( and α) and the vorticity field at t = 0.5 a for the Re = 100 case and t = 3 for the Re = 1000 case. The vorticity field is in good agreement with the one obtained in [24,38] in the two cases. Concerning the numerical viscosity, it is shown that the proposed methodology introduces a reduced amount of dissipation, compared with the δ-LES-SPH.
The decay of the maximum velocity for the two test cases is plotted in Figure 3 and compared with the theoretical solution [39]. Note that, for the case Re = 1000, the theoretical solution remains valid until t = 10 [16,29]. The results for the three schemes are practically coincident for the Re = 100 case, and are in agreement with the exact solution. The reason is that the different numerical viscosity introduced by the numerical schemes is not relevant considering the value of the molecular viscosity. The situation changes for Re = 1000. In this case, the dissipation introduced by the δ-LES-SPH is excessive. It is also seen that for this test case, using a δ parameter constant or variable, marginally affects the results. However, the result obtained of the δ-ADA-SPH method is practically coincident with the theoretical solution [39], since = 0 in most of the domain, as shown in Figure 2. For this test case, the original δ-SPH scheme lies slightly over the theoretical solution.

Inviscid Two-Dimensional Taylor-Green Flow
In this case, we solve the inviscid 2D Taylor-Green vortex. This is equivalent to consider Re = ∞. Thus, the part of β corresponding with the molecular viscosity is set to zero, and the only viscosity introduced in the problem is related to the α parameter, which is adaptive in δ + -ADA-SPH and δ + -LES-SPH methods and constant for the original δ + -SPH. For the latter case, we use a value of α = 0.01. The use of a constant value of the explicit dissipation is an approach often used for the simulation of inviscid flows with the δ-SPH and δ + -SPH methods [15,25]. Again, we use the value U max = 1 for computing Equation (2).
Thus, this test case is intended to determine if the proposed approach is able to recover the physical dynamics in a severely under-resolved simulation. The base case is computed using a resolution of 10,000 particles on a periodic [0, 1] × [0, 1] square domain. The initial position of the particles is obtained by using the particle packing algorithm described in [37].  In Figure 4, the value of given by Equation (17) is plotted, together with the resulting vorticity field. It is observed that the scheme is able to capture certain fine structures of the flow. Very interestingly, it is observed that the artificial dissipation is only introduced in regions where the smallest structures of the flow are present. This is specially clear for t = 25. The core of the final two vortices is well captured by the spatial discretization, and the artificial dissipation is null there. However, out of this region, the resolution is not enough to capture all the fine structures all the flow and the dissipation reaches the maximum value defined. It is also observed that the method is also able to recover the final equilibrium state of two-dimensional decay turbulence, in which a pair of big vortices with opposite sign and with a size comparable to the computational domain appear as a result of merging and pairing of the smaller vortices [24,40,41].   (17). On the right, vorticity field. Figures (a,c) show the results at t = 2 and figures (b,d) show the results obtained at t = 25, using a resolution of 10,000 particles.
In Figure 5 it is clearly observed that the results obtained for the δ + -LES-SPH and δ + -SPH are clearly over-dissipative. Thus, for this spatial resolution, these methods are not able to recover the fine features of the flow. This test case is representative of two-dimensional turbulence. Thus, it tests the ability of a new method for simulating this kind of problems. Two-dimensional turbulence has different energy-transfer characteristics than three-dimensional turbulence. Two-dimensional turbulence is characterized by a typical direct entropy cascade with −3 scaling of the energy spectrum [24]. In Figure 6 (left), it is shown that the proposed scheme is able to reproduce the right slope in most of the represented scales, independently of the resolution of the distribution of particles. The slope is not accurate for the highest scales, but this error is quickly reduced as the number of particles is increased. Moreover, in Figure 6 (right), it is shown that the proposed δ + -ADA-SPH scheme, greatly improves the results obtained by other δ + -SPH schemes. a b Figure 6. Inviscid two-dimensional Taylor-Green flow. Figure (a): Energy spectrum obtained using the proposed δ + -ADA-SPH method with different particle resolutions . Figure (b): comparison of the results using different δ + -SPH schemes using a resolution of 10,000 particles (right). Figure 7 shows the evolution of the kinetic energy. Ideally, it should remain constant and equal to the initial value, since at infinite Reynolds number there is no viscous dissipation. Thus, all the dissipation observed is introduced by the numerical method. When the proposed δ + -ADA-SPH scheme is used, the results are less dissipative than those obtained by the other δ + -SPH schemes. As the grid is refined, the proposed methodology also tends to reduce the dissipation, and, for the resolution of 40,000 particles only a slight decrease of the kinetic energy is observed. These results can be explained by examination of and α values. Figure 8 shows the values of these parameters at several times. It is seen that the values of are consistently lower than those of α, which indicates that the numerical scheme is introducing more dissipation. c d e f Figure 8. Inviscid two-dimensional Taylor-Green flow. Results for and α fields. On the left, the results obtained using the δ + -ADA-SPH for t = 0.5 (a), t = 2 (c) and t = 20 (e). On the right, the results obtained using the δ + -LES-SPH [16] approach for t = 0.5 (b), t = 2 (d) and t = 20 (f). Results computed using a resolution of 10,000 particles.

Lid-Driven Cavity Flow
The lid-driven cavity flow is a classical benchmark widely-studied to test new numerical methods [42][43][44][45][46]. Here, we use this test case to assess the accuracy of the proposed method for wall-bounded flows at different flow regimes. This problem consists of a two-dimensional square cavity filled with an incompressible fluid and covered by four rigid walls of length L. The top wall (lid) moves laterally at a constant velocity U that subsequently creates a fluid motion within the cavity under viscosity effect. The U max velocity for Equation (2) is estimated equal to one.
Here, 50 × 50 particles are used to solve the Re = 400, Re = 1000 and Re = 3200 test cases. The velocity profiles along vertical and horizontal centerlines obtained by the three different schemes are shown in Figure 9 compared with the data of [42].   The solutions are compared with those of Ghia [42]. Figures (a,b) show the results for Re = 400, figures (c,d) plot the results for Re = 1000 and figures (e,f) show the results obtained for Re = 3200.
In general, a good agreement is observed for the three schemes at low Reynolds number. However, as the Reynolds number is increased, the proposed method obtains more accurate results. We also note that, the present 2D simulations obtains a solution that remains steady. Some authors [43] have raised the question about the steadiness and two-dimensionality of this flow at Reynolds numbers higher than 1000.
The analysis of the distribution of the viscous parameters helps to understand these results. Figure 10 shows the values of and α for the Re = 3200 case. Note that, to facilitate the comparison, we have used the same scale for and α, but the maximum value achieved by α is more than three times larger than that of . We notice the large differences in the top boundary and in the top-right corner of the cavity, and the larger dissipation introduced by the δ-LES-SPH around the main vortex.  Figure 11 shows the results obtained using the proposed scheme as the grid is refined. As expected, the numerical solution converges to the reference one and no spurious artifacts are observed.

Two-Dimensional Dam-Break Flow
Finally, we consider the two-dimensional dam break problem. The aim of this test case was to assess the ability of the proposed method to deal with complex free-surface flows. This test case has been experimentally investigated by several authors such as Martin et al. [47], Zhou et al. [48] or Buchner [49], among others. This test case addresses a very complex flow with large deformations, and it is one of the most used benchmark to validate SPH and other particle-based methods [25,[50][51][52][53]. The initial configuration of the problem is illustrated in Figure 12. A water column of height H and length of 2H is located in a rectangular tank of length 5.366H. Here, H is taken equal to 0. 6 [m] as in the experimental test of Buchner [49]. Following [13,25], the water flow is considered inviscid with a density value of 1000 [kg/m 3 ], and the effect of gravity is also considered (g = (0, −9.81) T ). Thus, the part of β related to the molecular viscosity is set to zero. A pressure probe P located at the downstream wall at y/H = 0.19, measured from the bottom of the tank (see Figure 12). This position of the probe is slightly different to the one used in the experiments, following the recommendations of [54] to reduce some uncertainties in the measurements.  The generalized wall boundary condition method [50,55] is used to simulate this test case. Three different particle resolutions (H/∆x 0 = 120, 80, 60) are investigated. The pressure field is initialized from the solution of a Poisson equation on the water column [56]. The water column is considered at rest at t = 0 and c 0 = 10 2gH. The U max velocity required for Equation (2) is estimated using Torricelli's law as in [55].
We note that for this example, the particle shifting technique is not used. As shown in [29], the application of this technique to violent flows causes non-physical behavior, due to accumulation of errors in time. Then the ADA methodology is used combined with the δ-SPH scheme of Marrone et al. [25]. The results obtained are compared with those of the δ-LES-SPH (variable α and δ) and those of the original δ-SPH scheme [25]. Figures 13 and 14 show snapshots of the evolution of pressure, and vorticity fields for H/∆x 0 = 120. We note that the pressure field is relatively smooth, and no numerical artifacts are detected. Concerning the value of , it is shown that there are regions with a very small value of this parameter. This indicates a reduction of the numerical viscosity introduced by the numerical method. Moreover, the vorticity field shows the presence of fine scales of the flow. . Two-dimensional dam-break flow: Snapshots of the computational results at different times t g/H for H/∆x 0 = 120. On the left column the pressure field is presented, field is placed on the right. Figures (a,b) shows the results for t g/H = 2.35, (c,d) shows the results at t g/H = 4.04, whereas (e,f) plots the results for t g/H = 7.12. Figures (g,h) shows the results obtained at t g/H = 7.12.
a b c d Figure 14. Cont.
e f g h Figure 14. Two-dimensional dam-break flow: Snapshots of the computational results for different times t g/H with H/∆x 0 = 120. On the left column the vorticity field is presented, field is placed on the right. Figures (a,b)  In order to analyze the water wave front position, we define the adimensional position X f = , where x f is the x-coordinate of the water front. Figure 15 (left) shows the results obtained by the present methodology for three different particle resolutions, compared with the analytical solution from shallow water theory [57], and the experimental data obtained by Buchner et al. [49] and Lobovskỳ et al. [58]. A good agreement with the theoretical line is obtained with the three discretizations used. A quantitative agreement with the experimental data is also observed, although it should be noted that a slower front wave is observed in the experiments. This could be caused by several factors such as uncertainties of the measurements and wall roughness, as exposed in [59]. Figure 15 (right) shows a comparison of the time evolution of the pressure signals at the probe between the proposed method and the δ-LES-SPH for H/∆x 0 = 80. Both numerical models are in quantitative agreement with the experimental data. Note that the different location of the peak of the numerical solutions and experiments, is explained in [51] as the influence of the air-cushioning cavity on the sensor measurement which is not taken into account in the numerical simulation. a b Figure 15. Two-dimensional dam-break flow. Figure (a) shows the time evolution of the water front. The green solid line shows the analytical solution from shallow water theory [57] and the black triangles shows the experimental data given by Buchner [49]. The blue diamond shows the experimental data given by Lobovskỳ et al. [58]. Magenta (solid), blue (dashed) and red (dash-dot) lines shows the results obtained by the proposed scheme with H/∆x 0 = 120, 80 and 60, respectively. Figure (b) shows the time evolution of the pressure signals recorded at P is presented. The blue dashed line presents the results obtained by using the presented δ-ADA-SPH model. The magenta dotted line presents the results obtained with δ-LES-SPH approach. The particle resolution is H/∆x 0 = 80.
In order to study the intrinsic dissipation of the proposed numerical methodology, we analyze the evolution of the mechanical energy (∆E). Here, we use the expression proposed by Marrone et al. [25] to compute ∆E where E kin (t) is the kinetic energy, E p (t) the potential energy, E 0 p the initial potential energy and E ∞ p the potential energy when the flow reaches a hydrostatic steady state.
Results are plotted in Figure 16. It is observed that the use of the ADA method reduces the differences in the dissipation obtained for different discretizations. This is different to the usual behavior of SPH methods [13,25], and the explanation is that the dissipation is only introduced where the flow is underesolved. Moreover, it is also observed that the flow tends to stabilize and reach the rest state at t g/H = 20 which corresponds well with the experimental results of Buchner [49]. a b Figure 16. Two-dimensional dam-break flow. Figure (a) shows the time evolution of mechanical energy. Figure (b) plots the results obtained by the proposed method with different particle resolution. A zoom of the region of the curve before the first splash-up is displayed on the right. Figure 17 compares the evolution of ∆E before the first splash-up, obtained using the δ-ADA-SPH scheme with the results obtained with the δ-LES-SPH and the results obtained with the original δ-SPH scheme [25], for H/∆x 0 = 80. It is clearly observed that for this case, the proposed δ-ADA-SPH scheme outperforms the results of the original δ-SPH scheme and it is slightly less dissipative than the δ-LES-SPH.  [25] using δ-SPH scheme. All the results are computed with H/∆x 0 = 80.

Conclusions
In this work, an SPH scheme for weakly compressible flows with automatic adaptive numerical dissipation (ADA) is developed. The resulting scheme is a meshless self-adaptive method, in which the introduced artificial dissipation is designed to increase the dissipation in zones where the flow is under-resolved by the numerical scheme, and to decrease it where dissipation is not required. The proposed methodology is applied in the framework of δ-SPH and δ + -SPH formulations. The validation in terms of consistency and accuracy of the presented approach has been studied by means of several two-dimensional benchmarks for different Reynolds numbers: Taylor-Green, lid-driven cavity and dam-break flows. The results show that the proposed numerical method is able to simulate complex flows. It alleviates the parameter dependency of δ-SPH methods and it is generally less dissipative than the δ-LES-SPH model. The obtained results suggest than the proposed approach could be a promising technique for LES computations of transitional flows, where usual SGS models introduce excessive dissipation. This will be the subject of future research.