Improved δ-SPH Scheme with Automatic and Adaptive Numerical Dissipation

Abstract

In this work we present a $\delta$-Smoothed Particle Hydrodynamics (SPH) scheme for weakly compressible flows with automatic adaptive numerical dissipation. 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 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 flow is under-resolved in very coarse particle discretizations. Moreover, compared with the original $\delta$-SPH scheme, the proposed method reduces the number of problem-dependent parameters.

1. 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 $\delta$-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].

In $\delta$-SPH formulation, two extra diffusive terms are added, one in the continuity equation and the other in the momentum equation as a viscosity-like term. These terms are tuned through two parameters: $\delta$ for the continuity equation and $\alpha$ for the momentum equation. The additional dissipation introduced in the continuity equation leads to a reduction of the numerical high frequencies on the pressure field, leading to more regular distributions of particles. Thus, $\delta$-SPH family of schemes alleviate one of the main problems of Lagrangian schemes. Thus, in SPH schemes, the Lagrangian nature of the schemes leads to non-uniform distributions of particles, specially in flow regions with high velocity gradients, which causes problems of stability and accuracy of the numerical schemes [15]. However, in some conditions with certain irregularities of the particle spatial configurations, the $\delta$-SPH method does not prevent the tensile instability. Thus, a further improvement is presented in [15], where a particle shifting procedure is included in the $\delta$-SPH formulation, leading to the ${\delta }^{+}$-SPH method.

Unfortunately, there is not a single and well defined value for these parameters [16], as they are problem dependent. Moreover, it was also shown that different values of parameters could lead to similar values of total dissipated energy [18]. In [16], the $\delta$-SPH method is rewritten in terms of turbulence models, resulting in the $\delta$-LES-SPH scheme, in such a way that $\delta$ and $\alpha$ are computed each time step using the Smagorinsky subgrid scale model. This method represents a great improvement for $\delta$-SPH models since eliminates the need of ad-hoc definition of parameters for each problem. However, it is well known from grid-based methods that the use of the Smagorinsky closure model leads to excessive dissipation.

In this work, we propose a new method to define locally and automatically the value of the $\alpha$ parameter, to avoid the excessive dissipation caused by the classical Smagorinsky model. The $\delta$ parameter is kept fixed as suggested by [18]. Our proposal is based on the Automatic Dissipation Adjustment method (ADA), presented in [19] for low Mach computations using finite volume schemes and extended in [20] for supersonic and transonic computations. In the ADA method, the energy ratio (ER) criteria [21,22,23] are used to automatically adjust the amount of viscosity introduced by the numerical flux, increasing the numerical viscosity in those zones where the flow is under-resolved, and decreasing it in the other zones.

Our proposal for modifications to the $\delta$-SPH scheme is based on the ADA method in a weakly compressible Lagrangian approach. Using the proposed scheme, we are able to recover the theoretical decay of kinetic energy, even where the flow is under-resolved in very coarse particle discretizations. The fundamental idea of the proposed method is somewhat similar to that presented in [24], where an additional stress term depending on the content of high-frequencies of the solution.

The paper is structured as follows. In Section 2, the governing equations are presented, and in Section 3, different $\delta$-SPH formulations are presented, and then the new formulation is exposed. Section 4 is devoted to present some numerical tests to show the accuracy and robustness of the proposed methodology. Finally, the conclusions are drawn.

2. 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

$$\left\{\begin{array}{c}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla ·\mathbf{u}\hfill \\ \hfill \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t}=-\frac{\nabla p}{\rho }+\nu {\nabla }^2\mathbf{u}+\mathbf{g}\hfill \\ \hfill \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}=\mathbf{u}\hfill \end{array}\right.$$

(1)

where $\frac{\mathrm{d}(.)}{\mathrm{d}t}$ represents the material derivative following an infinitesimal fluid element (particle). ∇ is the nabla operator, $\rho$, p, $\mathbf{u}=(u,v)$, $\mathbf{r}$, $\nu$ and $\mathbf{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]

$$p=c_0^2(\rho -{\rho }_0)$$

where $c_0$ is the speed of sound, and ${\rho }_0$ is the density along the free surface. In the case of confined fluids ${\rho }_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\ge 10max(U_{max},\sqrt{\frac{p_{max}}{{\rho }_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.

3. 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 $\delta x_0$ ($h=2\delta 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={\mathsf{\Delta }{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={\rho }_0{V}_0$ during all the simulation time.

3.1. ${\delta }^{+}$-SPH and $\delta$-SPH Methods

For the discretization of the governing equations system in the framework of weakly compressible SPH method, the ${\delta }^{+}$-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 $\delta {\mathit{r}}_{\mathit{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 ${\delta }^{+}$-SPH model can be summarized as follows

$$\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}=-{\rho }_i\sum _j({\mathit{u}}_j-{\mathit{u}}_i)·{\nabla }_i{W}_{ij}{V}_j+h{c}_0\delta \sum _j{\mathcal{D}}_{ij}·{\nabla }_i{W}_{ij}{V}_j}\\ {\displaystyle \frac{\mathrm{d}{\mathit{u}}_i}{\mathrm{d}t}=-\frac{1}{{\rho }_i}\sum _j(p_i+p_j){\nabla }_i{W}_{ij}{V}_j+{\mathsf{\Pi }}_{\mathit{i}}^{\mu }+{\mathbf{g}}_i}\\ \frac{\mathrm{d}{\mathit{r}}_i}{\mathrm{d}t}={\mathit{u}}_i;p_i=c_0^2({\rho }_i-{\rho }_0);{\mathit{r}}_i={\mathit{r}}_{\mathit{i}}^{*}+\delta {\mathit{r}}_{\mathit{i}};V_i=\frac{m}{{\rho }_i};\end{array}\right.$$

(3)

where the symbol ${\nabla }_i$ denotes the gradient operator applied to the particle i. The quantities ${\rho }_i$, ${\mathit{u}}_i$, $p_i$, ${\mathit{r}}_i$, ${\mathit{r}}_{\mathit{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 ${\mathsf{\Pi }}_{\mathit{i}}^{\mu }$ denotes the viscous acceleration vector. The notation $W_{ij}=W\left(|{\mathit{r}}_i-{\mathit{r}}_j|,h\right)$ denotes the kernel function.

The term ${\mathcal{D}}_{ij}$ of the additional diffusive term is expressed as

$${\mathcal{D}}_{ij}=2\left[\left({\rho }_j-{\rho }_i\right)-\frac{1}{2}\left({\langle \nabla \rho \rangle }_j^L+{\langle \nabla \rho \rangle }_i^L\right)·\left({\mathbf{r}}_j-{\mathit{r}}_i\right)\right]\frac{\left({\mathit{r}}_j-{\mathit{r}}_i\right)}{{\left|{\mathit{r}}_j-{\mathit{r}}_i\right|}^2}$$

(4)

Here, ${\langle \nabla \rho \rangle }_i^L$ and ${\langle \nabla \rho \rangle }_j^L$ re-normalized density gradients [14] with respect to the particle i and its neighbor j, respectively.

$${\langle \nabla \rho \rangle }_i^L=\sum _j\left({\rho }_j-{\rho }_i\right){\mathbb{L}}_i{\nabla }_i{W}_{ij}{V}_j$$

(5)

where the re-normalization matrix ${\mathbb{L}}_i$ is expressed as

$${\mathbb{L}}_i={\left[\sum _j\left({\mathit{r}}_j-{\mathit{r}}_i\right)\otimes {\nabla }_i{W}_{ij}{V}_j\right]}^{-1}$$

The dimensionless parameter $\delta$ in Equation (3) defines the magnitude of the density diffusion. In ${\delta }^{+}$-SPH and $\delta$-SPH methods it is usually taken constant as $\delta =0.1$ [18].

In the momentum equation, the viscous part is defined as

$${\displaystyle {\mathsf{\Pi }}_{\mathit{i}}^{\mu }=\frac{1}{{\rho }_i}\sum _j{\beta \pi }_{ij}{\nabla }_i{W}_{ij}{V}_j}$$

(6)

where the term ${\pi }_{ij}$ reads as

$${\pi }_{ij}=\frac{\left({\mathit{u}}_j-{\mathit{u}}_i\right)·\left({\mathit{r}}_j-{\mathit{r}}_i\right)}{{\left({\mathit{r}}_j-{\mathit{r}}_i\right)}^2}$$

(7)

and $\beta$ is

$$\beta =\left({\rho }_{0}h{c}_0{\alpha }_{ij}+K{\mu }_i\right)$$

(8)

with $K=2(d+2)$ and d is the number of spatial dimensions. Note that there is a term in $\beta$ related to the molecular viscosity ($\mu$) and another term related with numerical dissipation added for stability reasons. This artificial viscosity is controlled by the non-dimensional parameter ${\alpha }_{ij}$. In ${\delta }^{+}$-SPH and $\delta$-SPH methods it is usually taken constant in the range [0.01–0.05]. On the other hand, the value of $\delta$ is kept constant with a value that depends on the problem. Taking $\delta =0.1$ is an usual choice.

The main difference between ${\delta }^{+}$-SPH and $\delta$-SPH methods is in the updating of the position of the particles. Thus, the ${\delta }^{+}$-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

$$\delta {\mathit{r}}_{\mathit{i}}=-\mathrm{CFL}·\mathrm{Ma}·2·h^2·\sum _j\left[1+{\chi }_{ij}\right]{\nabla }_i{W}_{ij}{V}_j$$

(9)

where $\mathrm{CFL}=c_o\frac{\mathsf{\Delta }t}{h}$ is the Courant–Friedrichs–Lewy number, $\mathrm{Ma}=\frac{U_{max}}{c_0}$ is a reference Mach number and ${\chi }_{ij}$ denotes the artificial pressure-like function described by Monaghan [28] and it is expressed as follows

$${\chi }_{ij}=R{\left(\frac{W_{ij}}{W\left(\mathsf{\Delta }{x}_0\right)}\right)}^n$$

(10)

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 ($\mathbf{\delta }{\mathit{r}}_{\mathit{i}}=0$) in $\delta$-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 $\mathbf{\delta }{\mathit{r}}_{\mathit{i}}=0$ and then the (${\delta }^{+}$-SPH) scheme is reduced to the $\delta$-SPH method [25].

3.2. $\delta$-LES-SPH Method

In [18], the $\delta$-LES-SPH method [16] is written in a similar fashion to the $\delta$-SPH method. Using the same expressions that were presented in the previous section, we can define the $\delta$-LES-SPH method as follows. First, instead of having a constant ${\alpha }_{ij}$, this parameter can vary locally for each pair of interacting particles as [16]

$${\alpha }_{ij}=2\frac{{\alpha }_i{\alpha }_j}{{\alpha }_i+{\alpha }_j}$$

(11)

where the parameter ${\alpha }_i$ is usually defined using a subgrid-scale model [30,31]. For the $\delta$-LES-SPH model presented in [18], the Smagorinsky SGS model is used, and thus ${\alpha }_i$ can read as follows

$${\alpha }_i=\frac{K{\nu }_i^T}{c_{0}h}$$

(12)

where ${\nu }_i^T$ is the kinematic eddy viscosity that can be expressed as

$${\nu }_i^T={\left(C_s{I}_{LES}\right)}^2\parallel {\mathbb{D}}_i\parallel$$

(13)

$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\mathsf{\Delta }x$) [16]. Moreover $\parallel {\mathbb{D}}_i\parallel =\sqrt{2{\mathbb{D}}_i:{\mathbb{D}}_i}$ is a rescaled Frobenius norm for the particle i, and the strain-rate tensor is computed as

$${\mathbb{D}}_i=\frac{1}{2}\sum _j\left[\left({\mathit{u}}_j-{\mathit{u}}_i\right)\otimes \left({\mathbb{L}}_i·\nabla W_{ij}\right)+\left({\mathbb{L}}_i·\nabla W_{ij}\right)\otimes \left({\mathit{u}}_j-{\mathit{u}}_i\right)\right]V_j$$

where the symbol ⊗ presents the outer product between two vectors.

Similarly, in the $\delta$-LES-SPH method, the value of ${\delta }_{ij}$ is not constant, and it varies locally for each pair of interacting particles as well as ${\alpha }_{ij}$

$${\delta }_{ij}=2\frac{{\delta }_i{\delta }_j}{{\delta }_i+{\delta }_j}$$

where

$${\delta }_i=\frac{{\nu }_i^{\delta }}{c_{0}h};{\nu }_i^{\delta }={\left(C_{\delta }{I}_{LES}\right)}^2\parallel {\mathbb{D}}_i\parallel$$

Here, the constant of the Smagorinsky SGS model is taken as $C_{\delta }=1.5$. The variation of ${\alpha }_i$ as well as ${\delta }_i$ is limited to $0.2$ in order to avoid stability issues [18].

3.3. A New $\delta$-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 $\delta$ and $\alpha$ parameters. Thus, we propose to fix the value of $\delta$ and substitute the parameter (${\alpha }_{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 ($E{R}_i$) of the particle i [21]

$$E{R}_i=\frac{{({\mathit{u}}_{\mathit{i}}-\widetilde{{\mathit{u}}_{\mathit{i}}})}^2}{{({\mathit{u}}_{\mathit{i}}-\widehat{{\mathit{u}}_{\mathit{i}}})}^2}$$

(14)

In Equation (14), ${\mathit{u}}_{\mathit{i}}$ is the velocity field obtained as a result of the SPH computation. Moreover, $\widetilde{{\mathit{u}}_{\mathit{i}}}$ and $\widehat{{\mathit{u}}_{\mathit{i}}}$ are two filtered velocity fields, obtained through filtering of ${\mathit{u}}_{\mathit{i}}$ using a low-pass filter with different widths. Thus, in Equation (14) $\widetilde{{\mathit{u}}_{\mathit{i}}}$ and $\widehat{{\mathit{u}}_{\mathit{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 $\widetilde{{\mathit{u}}_{\mathit{i}}}$ and $\widehat{{\mathit{u}}_{\mathit{i}}}$ will respectively read as

$$\widetilde{{\mathit{u}}_{\mathit{i}}}=\frac{{\sum }_j{\mathit{u}}_{\mathit{j}}W(r_{ij},\tilde{h})V_j}{{\sum }_{j}W(r_{ij},\tilde{h})V_j}$$

(15)

and

$$\widehat{{\mathit{u}}_{\mathit{i}}}=\frac{{\sum }_j{\mathit{u}}_{\mathit{j}}W(r_{ij},\widehat{h})V_j}{{\sum }_{j}W(r_{ij},\widehat{h})V_j}$$

(16)

where $\tilde{h}$ and $\widehat{h}$ denote the smoothing lengths that define the Wendland-C2 kernel radius and they are evaluated respectively as $\tilde{h}=h=2\mathsf{\Delta }{x}_0$ and $\widehat{h}=2h=4\mathsf{\Delta }{x}_0$.

As stated in Equation (14), $E{R}_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 $E{R}_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 $E{R}_i$.

In order to automatically adjust the parameter, we define the value ${ϵ}_i$ associated to a particle i following the rule proposed in [19]

$$\left\{\begin{array}{cc}E{R}_i<0.5,\hfill & {ϵ}_i=\max [({ϵ}_i-\varphi ),0]\hfill \\ E{R}_i>0.55,\hfill & {ϵ}_i=\min [({ϵ}_i+\varphi ),0.01]\hfill \\ 0.5\le E{R}_i\le 0.55,\hfill & {ϵ}_i\mathrm{does}\mathrm{not}\mathrm{change}\hfill \end{array}\right.$$

(17)

Here, a value of $\varphi =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 $E{R}_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 $E{R}_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 $\alpha$ 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

$$ϵ=\frac{{ϵ}_i+{ϵ}_j}{2}$$

(18)

We introduce this approach in the SPH scheme (3), by substituting the dimensionless parameter ${\alpha }_{ij}$ in (8) by $ϵ$ defined as in Equation (18).

$${\displaystyle {\mathsf{\Pi }}_{\mathit{i}}^{\mu }=\frac{1}{{\rho }_i}\sum _j{\left({\rho }_{0}h{c}_0ϵ+K{\mu }_i\right)\pi }_{ij}{\nabla }_i{W}_{ij}{V}_j}$$

(19)

Moreover, following [17], the value of $\delta$ is kept constant as $\delta =0.1$ for all simulation cases. This choice is justified since the range of variability of $\delta$ 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 $\delta$ ensures that the dissipation varies with $\alpha$. Thus, in the framework of the $\delta$-ADA-SPH method, we ensure that the dissipation is totally controlled by $ϵ$.

To ensure the stability of the method, the time step ($\delta t$) must be chosen to fulfill the kinetic, body force [34]

$$\delta t=min\left\{CFL\frac{h}{c_0},0.25{\left(\frac{h}{\parallel \mathbf{g}\parallel }\right)}^{1/2}\right\}$$

(20)

In this work, low storage 4th-order Runge–Kutta scheme [35] is used to integrate the governing equation system. The Courant–Friedrichs–Lewy constant $\mathrm{CFL}$ is set equal to $1.5$.

4. Numerical Examples

4.1. 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)=-U{e}^{bt}cos\left(2\pi x\right)sin\left(2\pi y\right),v(x,y,t)=U{e}^{bt}sin\left(2\pi x\right)cos\left(2\pi y\right)$$

(21)

where the reference velocity U is set to $U=1$ and $b=-8{\pi }^2/Re$ is the decay rate of velocity field. The density is assumed to be the unity, and the initial pressure is given by

$$P(x,y,0)=-\frac{{\rho }_0}{4}(cos\left(4\pi x\right)+cos\left(4\pi y\right))$$

(22)

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 ${\delta }^{+}$-SPH scheme, the $\delta$-LES-SPH method, and a modified $\delta$-LES-SPH method keeping a constant value of $\delta =0.1$. The value of $\beta$ is set to match the molecular viscosity in the $\delta$-SPH original scheme, and the term related with the added numerical dissipation is neglected (${\alpha }_{ij}=0$). In $\delta$-LES-SPH and $\delta$-ADA-SPH models, the term (${\alpha }_{ij}$) is added to the molecular viscosity, as in Equation (8).

Figure 1 and Figure 2 show the results for the value of the viscous parameters ($ϵ$ and $\alpha$) 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 $\delta$-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 $\delta$-LES-SPH is excessive. It is also seen that for this test case, using a $\delta$ parameter constant or variable, marginally affects the results. However, the result obtained of the $\delta$-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 $\delta$-SPH scheme lies slightly over the theoretical solution.

4.2. Inviscid Two-Dimensional Taylor–Green Flow

In this case, we solve the inviscid 2D Taylor–Green vortex. This is equivalent to consider $Re=\infty$. Thus, the part of $\beta$ corresponding with the molecular viscosity is set to zero, and the only viscosity introduced in the problem is related to the $\alpha$ parameter, which is adaptive in ${\delta }^{+}$-ADA-SPH and ${\delta }^{+}$-LES-SPH methods and constant for the original ${\delta }^{+}$-SPH. For the latter case, we use a value of $\alpha =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 $\delta$-SPH and ${\delta }^{+}$-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]\times [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].

In Figure 5 it is clearly observed that the results obtained for the ${\delta }^{+}$-LES-SPH and ${\delta }^{+}$-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 ${\delta }^{+}$-ADA-SPH scheme, greatly improves the results obtained by other ${\delta }^{+}$-SPH schemes.

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 ${\delta }^{+}$-ADA-SPH scheme is used, the results are less dissipative than those obtained by the other ${\delta }^{+}$-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 $\alpha$ 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 $\alpha$, which indicates that the numerical scheme is introducing more dissipation.

5. 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\times 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].

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 $\alpha$ for the $Re=3200$ case. Note that, to facilitate the comparison, we have used the same scale for $ϵ$ and $\alpha$, but the maximum value achieved by $\alpha$ 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 $\delta$-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 ($\mathit{g}={\left(0,-9.81\right)}^T$). Thus, the part of $\beta$ related to the molecular viscosity is set to zero. A pressure probe $\mathit{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/\mathsf{\Delta }{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\sqrt{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 $\delta$-SPH scheme of Marrone et al. [25]. The results obtained are compared with those of the $\delta$-LES-SPH (variable $\alpha$ and $\delta$) and those of the original $\delta$-SPH scheme [25].

Figure 13 and Figure 14 show snapshots of the evolution of pressure, $ϵ$ and vorticity fields for $H/\mathsf{\Delta }{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.

In order to analyze the water wave front position, we define the adimensional position $X_f=\frac{x_f-2H}{H}$, 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 $\delta$-LES-SPH for $H/\mathsf{\Delta }{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.

In order to study the intrinsic dissipation of the proposed numerical methodology, we analyze the evolution of the mechanical energy ($\mathsf{\Delta }E$). Here, we use the expression proposed by Marrone et al. [25] to compute $\mathsf{\Delta }E$

$$\mathsf{\Delta }E=\frac{E_{kin}\left(t\right)+E_p\left(t\right)-E_p^0}{E_p^0-E_p^{\infty }}$$

(23)

where $E_{kin}\left(t\right)$ is the kinetic energy, $E_p\left(t\right)$ the potential energy, $E_p^0$ the initial potential energy and $E_p^{\infty }$ 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\sqrt{g/H}=20$ which corresponds well with the experimental results of Buchner [49].

Figure 17 compares the evolution of $\mathsf{\Delta }E$ before the first splash-up, obtained using the $\delta$-ADA-SPH scheme with the results obtained with the $\delta$-LES-SPH and the results obtained with the original $\delta$-SPH scheme [25], for $H/\mathsf{\Delta }{x}_0=80$. It is clearly observed that for this case, the proposed $\delta$-ADA-SPH scheme outperforms the results of the original $\delta$-SPH scheme and it is slightly less dissipative than the $\delta$-LES-SPH.

6. 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 $\delta$-SPH and ${\delta }^{+}$-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 $\delta$-SPH methods and it is generally less dissipative than the $\delta$-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.