Dual-Time Smoothed Particle Hydrodynamics for Incompressible Fluid Simulation

TL;DR

This paper introduces a dual-time stepping scheme combined with EDAC for SPH, significantly improving accuracy and performance in incompressible fluid simulations, with open-source implementation and benchmark validation.

Contribution

It presents a novel dual-time SPH scheme with EDAC, enhancing accuracy and efficiency over traditional methods for incompressible fluid simulation.

Findings

Up to seven times better performance than standard SPH

Accurate and robust simulation results

Open-source implementation available

Abstract

In this paper we propose a dual-time stepping scheme for the Smoothed Particle Hydrodynamics (SPH) method. Dual-time stepping has been used in the context of other numerical methods for the simulation of incompressible fluid flows. Here we provide a scheme that combines the entropically damped artificial compressibility (EDAC) along with dual-time stepping. The method is accurate, robust, and demonstrates up to seven times better performance than the standard weakly-compressible formulation. We demonstrate several benchmarks showing the applicability of the scheme. In addition, we provide a completely open source implementation and a reproducible manuscript.

Figures

Tables

Table 1. Table 1 : CPU time taken for 2.5 2.5 2.5 secs of Taylor-Green simulation with 100 × 100 100 100 100\times 100 particles, with advection and without advection in pseudo time.

Scheme	CPU time (secs)
DTSPH frozen	326.87
DTSPH advect, no update neighbors	638.62
DTSPH advect, update neighbors	1934.21

Table 2. Table 2 : CPU time time taken for a simulation time of 2.5 2.5 2.5 secs with 100 × 100 100 100 100\times 100 particles with varying tolerance.

$ϵ$	CPU time (secs)
0.01	120.22
0.0001	133.87
1e-05	305.88

Table 3. Table 3 : CPU time time taken for a simulation time of 2.5 2.5 2.5 secs with 100 × 100 100 100 100\times 100 particles for various schemes.

Scheme	CPU time (secs)
DTSPH	124.49
WCSPH	239.70
EDAC	275.43
$\delta$-SPH	347.98

Table 4. Table 4 : CPU time time taken for different simulations of the 3D dam break problem.

Scheme	$ϵ$	CPU time (secs)
DTSPH	0.001	563.19
DTSPH	0.0001	2404.97
EDAC	NA	4035.25

Table 5. Table 5 : Speed-up obtained for different simulations when using the DTSPH scheme.

Problem	Scheme	Speed up
Dam break 3D	EDAC vs. DTSPH ($ϵ={10}^{-3}$)	7.16
Dam break 3D	EDAC vs. DTSPH ($ϵ={10}^{-4}$)	1.68
Dam break 2D	EDAC vs. DTSPH	3.57
Dam break 2D	EDAC vs. $\delta$-SPH	3.66
Cavity	TVF vs DTSPH	6.87

Equations

$$\frac{dp}{d\tau }=\frac{\partial p}{\partial \tau }+\tilde{\mathbf{V}}\cdot \text{grad}(p).$$

$$\frac{1}{{\beta }^2}\frac{\partial p^{\ast }}{\partial {\tau }^{\ast }}+\text{div}({\mathbf{V}}^{\ast })=0.$$

$$\frac{\partial p}{\partial \tau }=-\rho c^2\text{div}(\mathbf{V}).$$

$$\frac{dp}{d\tau }=\tilde{\mathbf{V}}\cdot \text{grad}(p)-\rho c^2\text{div}(\mathbf{V}).$$

$$\frac{dp}{d\tau }=\tilde{\mathbf{V}}\cdot \text{grad}(p)-\rho c^2\text{div}(\mathbf{V})+{\nu }_e{\nabla }^{2}p,$$

$${\nu }_e=\frac{hc}{16}.$$

$$\frac{\partial p}{\partial \tau }=-\rho c^2\text{div}(\mathbf{V})+{\nu }_e{\nabla }^{2}p.$$

$$\frac{\partial \mathbf{V}}{\partial \tau }+\frac{d\mathbf{V}}{dt}=-\frac{1}{\rho }\text{grad}(p)+\nu {\nabla }^2\mathbf{V}+\mathbf{f},$$

$$\frac{d\mathbf{V}}{d\tau }+\frac{d\mathbf{V}}{dt}=\tilde{\mathbf{V}}\cdot \text{grad}(\mathbf{V})-\frac{1}{\rho }\text{grad}(p)+\nu {\nabla }^2\mathbf{V}+\mathbf{f}.$$

$${\mathbf{V}}^{k=0}={\mathbf{V}}^n+\Delta t{\left(\frac{d\mathbf{V}}{dt}\right)}^n,$$

$${\mathbf{r}}^{k=0}={\mathbf{r}}^n+\Delta t{\mathbf{V}}^n+\frac{\Delta t^2}{2}{\left(\frac{d\mathbf{V}}{dt}\right)}^n,$$

$${\mathbf{r}}^{k+1/2}={\mathbf{r}}^k+\frac{1}{2}\Delta \tau {\tilde{\mathbf{V}}}^k,$$

$$p^{k+1/2}=p^k+\frac{1}{2}\Delta \tau {\left(\frac{dp}{d\tau }\right)}^k,$$

$${\mathbf{V}}^{k+1}={\mathbf{V}}^k+\Delta \tau {\left(\frac{d\mathbf{V}}{d\tau }\right)}^{k+1/2},$$

$${\mathbf{r}}^{k+1}={\mathbf{r}}^{k+1/2}+\frac{1}{2}\Delta \tau {\tilde{\mathbf{V}}}^{k+1},$$

$$p^{k+1}=p^{k+1/2}+\frac{1}{2}\Delta \tau {\left(\frac{dp}{d\tau }\right)}^{k+1/2}.$$

$${\tilde{\mathbf{V}}}^{k+1}={\left(\frac{d\mathbf{V}}{d\tau }\right)}^{k+1/2}\Delta t.$$

$${\left(\frac{d\mathbf{V}}{dt}\right)}^{n+1}=\frac{3{\mathbf{V}}^{n+1}-4{\mathbf{V}}^n+{\mathbf{V}}^{n-1}}{2\Delta t}+{O(\Delta t)}^2.$$

$${\left(\frac{d\mathbf{V}}{dt}\right)}^{n+1}=\frac{3\Delta \tau }{2\Delta t}{\left(\frac{d\mathbf{V}}{d\tau }\right)}^{k+1/2}+\frac{3{\mathbf{V}}^k-4{\mathbf{V}}^n+{\mathbf{V}}^{n-1}}{2\Delta t}+{O(\Delta t)}^2.$$

$$\begin{array}{rl}{\left(\frac{d\mathbf{V}}{d\tau }\right)}^{k+1/2}\approx & \{\tilde{\mathbf{V}}\cdot \text{grad}(\mathbf{V})-\frac{1}{\rho }\text{grad}(p)+\nu {\nabla }^2\mathbf{V}+\mathbf{f}\\ & -\frac{(3{\mathbf{V}}^k-4{\mathbf{V}}^n+{\mathbf{V}}^{n-1})}{2\Delta t}\}\left(\frac{2\Delta t}{2\Delta t+3\Delta \tau }\right).\end{array}$$

$$p^{k+1/2}=p^k+\frac{1}{2}\Delta \tau {\left(\frac{\partial p}{\partial \tau }\right)}^k,$$

$${\mathbf{V}}^{k+1}={\mathbf{V}}^k+\Delta \tau {\left(\frac{\partial \mathbf{V}}{\partial \tau }\right)}^{k+1/2},$$

$$p^{k+1}=p^{k+1/2}+\frac{1}{2}\Delta \tau {\left(\frac{\partial p}{\partial \tau }\right)}^{k+1/2}.$$

$${\left(\frac{d\mathbf{V}}{dt}\right)}^{n+1}=\frac{3\Delta \tau }{2\Delta t}{\left(\frac{\partial \mathbf{V}}{\partial \tau }\right)}^{k+1/2}+\frac{3{\mathbf{V}}^k-4{\mathbf{V}}^n+{\mathbf{V}}^{n-1}}{2\Delta t},$$

$$\begin{array}{rl}{\left(\frac{\partial \mathbf{V}}{\partial \tau }\right)}^{k+1/2}\approx & \{-\frac{1}{\rho }\text{grad}(p)+\nu {\nabla }^2\mathbf{V}+\mathbf{f}\\ & -\frac{(3{\mathbf{V}}^k-4{\mathbf{V}}^n+{\mathbf{V}}^{n-1})}{2\Delta t}\}\left(\frac{2\Delta t}{2\Delta t+3\Delta \tau }\right).\end{array}$$

$${\mathbf{r}}^{n+1}={\mathbf{r}}^n+\frac{\Delta t}{2}\left({\mathbf{V}}^n+{\mathbf{V}}^{n+1}\right).$$

$$\frac{\partial p}{\partial \tau }+\frac{\partial p}{\partial t}+\mathbf{V}\cdot \text{grad}(\mathbf{V})=-\rho c^2\text{div}(\mathbf{V}).$$

$$\frac{\partial p}{\partial \tau }=-\mathbf{V}\cdot \text{grad}(\mathbf{V})-\rho c^2\text{div}(\mathbf{V}).$$

$$\frac{\partial \mathbf{V}}{\partial \tau }=-\mathbf{V}\cdot \text{grad}(\mathbf{V})-\frac{1}{\rho }\text{grad}(p)+\nu {\nabla }^2\mathbf{V}+\mathbf{f}.$$

$$\begin{array}{rl}\frac{d{p}_i}{d\tau }=& {\tilde{\mathbf{V}}}_i\cdot \sum _{j\in N(i)}{m}_j\left(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}\right)\nabla W_{ij}\\ +& \sum _{j\in N(i)}\frac{m_j{\rho }_i}{{\rho }_j}{c}^2{\mathbf{V}}_{ij}\cdot \nabla W_{ij}\\ +& \sum _{j\in N(i)}\frac{4m_j}{({\rho }_i+{\rho }_j){\rho }_j}{\nu }_e\frac{p_{ij}}{({\mathbf{r}}_{ij}^2+{\mathbf{h}}_{ij}^2)}\nabla W_{ij}\cdot {\mathbf{x}}_{ij},\end{array}$$

Full text

1 Introduction

The Smoothed Particle Hydrodynamics (SPH) method originated with the work of Gingold and Monaghan [1], and Lucy [2] as a method to simulate astrophysical problems. The method is grid-free and Lagrangian in nature. It has since become a very general purpose technique and applied to a variety of problems including incompressible fluid flow[3, 4, 5], solid mechanics[6], and fluid-structure interaction [7, 8].

There are several SPH schemes for simulating incompressible and weakly compressible fluid flows. The original weakly compressible SPH scheme (WCSPH) was proposed by Monaghan [3]. This scheme treats the fluid as weakly compressible with an artificial sound speed and a stiff equation of state. This allows the scheme to utilize a hyperbolic system of equations and integrate them in time. There are many significant variants of this scheme including a Transport Velocity Formulation (TVF) [9] which introduces a transport velocity to ensure particle homogeneity. The original WCSPH and their derivatives generally suffer from a large amount of pressure oscillations and the $\delta$-SPH scheme [10, 11] reduces these oscillations by introducing a dissipation into the continuity equation. Similarly, an Entropically Damped Artificial Compressibility SPH scheme (EDAC-SPH)[12] has been proposed which introduces entropy by diffusing the pressure. This approach is quite similar to the $\delta$-SPH scheme and both schemes produce superior pressure distributions. All of these schemes employ an artificial sound speed and this places severe time step limitations due to stability considerations.

Cummins and Rudman [4] proposed a family of projection based schemes for incompressible fluids. Shao and Lo [5] and Hu and Adams [13] proposed incompressible SPH (ISPH) schemes which satisfy incompressibility by solving a pressure-Poisson equation. These approaches eliminate the need for evolving the pressure at the sound speed and this significantly increases the allowed time steps. The difficulty with the projection and incompressible schemes is the requirement to solve a linear system of equations which can be time consuming and involved. Recently, a Predictive-Corrective ISPH (PCISPH) [14] has been proposed for use in the graphics community for rapid simulation of incompressible fluids. A more accurate and efficient scheme has been proposed called the Implicit-Incompressible SPH (IISPH) [15]. The IISPH is matrix-free, and very efficient. It has been shown to be close to an order of magnitude faster than traditional schemes. However, the IISPH can be more involved to implement than many of the traditional WCSPH-based schemes. Recently, another fast and matrix-free implementation of the ISPH method has been proposed [16].

In this paper we propose a new scheme for weakly-compressible fluid flows. Our paper takes inspiration from the Artificial Compressiblility-based Incompressible SPH (ACISPH) scheme proposed by Rouzbahani and Hejranfar [17]. Our scheme uses a different formulation that is also very efficient. The original scheme was not noted in particular for its efficiency. We propose an original derivation and suggest many improvements that make the proposed scheme efficient. The performance is significantly better than that of the traditional WCSPH schemes and comparable to that of the ISPH schemes without sacrificing any accuracy or being unduly hard to implement. Our approach employs the classic artificial compressibility of Chorin [18] in a dual-time stepped framework. We call the resulting scheme, DTSPH for Dual-Time stepped SPH.

Fatehi et al. [19] propose density-based dual-timestepping schemes for the SPH method. They propose two different formulations in order to update the pressure in pseudo-time. They perform an accurate discretization of the derivatives. While their method is in principle similar to the present work, it does not demonstrate any significant performance advantage over the traditional WCSPH schemes. Moreover, they do not demonstrate their method with any free-surface problems. Our proposed implementation and method is however significantly faster than the WCSPH scheme and has been demonstrated to work well for free-surface and three-dimensional problems. Zhang et al. [20] on the other hand propose a dual-criteria timestepping strategy that is quite different from the method proposed here or in [19, 17]. They use an “acoustic timestep” to relax the pressure and an “advective timestep” to update the verlet lists used for neighbor computation. This optimization results in a performance improvement without significant complexity albeit at the cost of increased memory. The present method demonstrates much improved performance and does not involve a significant memory penalty.

The new scheme is designed to be implemented as an extension of the classic weakly-compressible schemes. The advantage of this is that these are relatively easy to implement, boundary conditions may be easily enforced, and there are several well-established schemes that may be used. Although we have not done so in this manuscript, it is important to note that this approach may also be employed for solid mechanics problems where the artificial speed of sound is usually very large.

The new scheme can be adapted to any WCSPH formulation which uses a density or pressure evolution equation based on a continuity equation. We demonstrate the scheme with the EDAC scheme [12]. We note that it can be easily applied to other schemes like the $\delta$-SPH. We show how our scheme can be used to obtain steady state solutions, although this is not particular to the new scheme and can be easily performed for a variety of other schemes. Obtaining steady state solutions in the context of SPH simulations is useful in different contexts. For example, in the case of the impulsively started flow past a complex geometry, an initial potential flow solution is useful and this may be easily obtained using this approach.

In this manuscript we provide a new formulation as compared to the work of [17], explore several important details for the implementation of the scheme, and more importantly provide a high-performance, open source implementation of the scheme. Our implementation uses the open source PySPH framework [21, 22] and all the code related to the manuscript is available at https://gitlab.com/pypr/dtsph. In order to facilitate reproducible research, this entire manuscript is completely reproducible and every figure in this paper is automatically generated using automan [23].

In the next section we discuss the proposed DTSPH scheme in a general setting and then discuss the SPH discretization. We briefly discuss the stability requirements of the scheme. We show how the resulting scheme is efficient and then proceed to simulate various standard benchmark problems. We perform comparisons with the TVF [9], the $\delta$-SPH scheme [10], and the EDAC scheme [12] where relevant to demonstrate the accuracy and efficiency of the proposed scheme.

2 The dual-time SPH method

In dual-time stepping schemes, a new time dimension called the “dual-time”, denoted by the variable $\tau$, is introduced. We have the following important considerations to keep in mind. If $\mathbf{r}$ is the position vector of a particle, then the real velocity of the particle is defined as, $\mathbf{V}=\frac{d\mathbf{r}}{dt}$. On the other hand, if the particle were to move in pseudo-time, we define the velocity in pseudo-time as $\mathbf{\tilde{V}}=\frac{d\mathbf{r}}{d\tau}$.

If we consider a property of the particle, $p$, then we can define a material derivative in pseudo-time as,

[TABLE]

The ACISPH formulation [17] uses a non-dimensionalized form of the equation,

[TABLE]

Here we use a star for the non-dimensional terms. By assuming that $p=(\rho-\rho_{0})c^{2}$ where $c$ is an artificial sound speed and $\rho_{0}$ is a reference density, the above can be written in a dimensional form as,

[TABLE]

By rewriting the non-dimensional form in equation (2) with suitable dimensional quantities, it is easy to see that $\beta=1/M$ where $M=V_{ref}/c$ is the Mach number of the reference speed of the flow, $V_{ref}$.

If we move the particles in pseudo-time, we introduce a material derivative to get,

[TABLE]

With the EDAC formulation [12], the pressure evolution equation also includes a diffusive term and we get,

[TABLE]

where $\nu_{e}$ is an artificial viscosity parameter which is typically chosen as,

[TABLE]

Here $h$ is the SPH kernel radius. In the above, we have assumed that the particles move in pseudo-time. In later sections, we show that it is computationally efficient to not move the particles in pseudo-time. This results in a simpler equation for the pressure evolution without a time derivative as,

[TABLE]

The momentum equation can be written similarly by adding a time derivative of velocity in pseudo-time,

[TABLE]

where $\nu$ is the kinematic viscosity of the fluid and $\mathbf{f}$ is the body force. Again, if we choose to move the particles in pseudo-time, we can write this in terms of a material derivative in pseudo-time as,

[TABLE]

Note the key difference here from what is proposed in [17] is that they have used $\mathbf{V}$ where they should have only used $\mathbf{\tilde{V}}$.

There are two possible approaches we can choose for implementation:

Move particles in real time and pseudo-time and use the equations (4) and (9). Note that when using the EDAC formulation we would add the pressure diffusion term as shown in the right hand side of (5).

• 2.

Move the particles only in real time and use equations (3) (or (7) in the case of the EDAC scheme) and (8).

We discuss these two approaches next.

2.1 Moving particles in pseudo-time

In this section we show how the above equations are integrated in real time and pseudo-time. The following equations apply to each particle, $i$. We suppress the subscript $i$ in the following to simplify the notation. We use the index $k$ to denote pseudo-time iterations and $n$ for the real time. Before iterating in pseudo time the particles are updated to a guessed new state ($k=0$) for the next real time ($n+1$) using,

[TABLE]

[TABLE]

where ${(\frac{d\mathbf{V}}{dt})}^{n}$ is given by the momentum equation (8) without the partial derivative of velocity in pseudo time (i.e. considering $\frac{\partial\mathbf{V}}{\partial\tau}=0$). Then the integration in pseudo time proceeds in the following fashion, with $\mathbf{\tilde{V}}^{k=0}=0$ as the starting value,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In addition to these we have,

[TABLE]

This is to ensure consistency of the motion. A detailed proof for equation (17) is provided in the appendix at the end of this manuscript.

We need an expression for the term $\frac{d\mathbf{V}}{dt}$ in the momentum equation (9). We obtain this using an implicit three-point backward difference scheme to discretize the real-time derivative,

[TABLE]

Substitute $\mathbf{V}^{k+1}$ instead of $\mathbf{V}^{n+1}$ in (18) and add and subtract the term $3\mathbf{V}^{k}$ in the numerator and rearrange terms to get,

[TABLE]

If we use equation (9), we can rewrite the above as,

[TABLE]

We may now discretize the right hand side using SPH and find the acceleration to the velocity in pseudo-time. This is done in section 3.

2.2 Fixed particles in pseudo-time

In this section we consider the case where the particles are not moved in pseudo-time. In this case, we can perform the integration as follows,

[TABLE]

[TABLE]

[TABLE]

Further, starting from equation (22), the equation corresponding to (19) becomes,

[TABLE]

where, the material derivative in the first term on the right hand side is replaced with a partial pseudo-time derivative. The equation (25) also changes appropriately to no longer include the pseudo-time advective term, resulting in

[TABLE]

We note that usually the velocity of a particle in pseudo-time, $\mathbf{\tilde{V}}$, is very small and this makes the changes to the position even smaller. This means that it is computationally efficient to fix the position of the particles in pseudo-time. In addition, even if we were to move the particles, we do not need to recompute the neighbor information.

Once the pseudo-time iterations are completed, we update the particle positions using,

[TABLE]

2.3 Steady state solutions

We can use the dual-time to seek a solution to steady state problems. To do this we set the partial derivative in time to zero and retain only the pseudo time derivative. Further, we do not move the particles at all. We start with the following form of the continuity equation,

[TABLE]

We set the partial derivative in time to zero and this results in using the following form of equation for the evolution of the pressure,

[TABLE]

We note that we have not used the EDAC scheme in this case. Similarly, the equation (8) for momentum reduces to,

[TABLE]

Here, we have moved the convection term to the right side and removed any partial time derivatives. This can be easily solved purely in pseudo-time while keeping the particles fixed in space. Technically, we could replace $\tau$ with $t$. However, the dual-time offers a convenient perspective for seeking a steady state solution iteratively.

The above approach to obtain steady state solutions is simple and not tied to any particular SPH scheme. Any scheme that uses a density or pressure evolution equation that is dependent on the divergence of the velocity will work. The approach is numerically efficient as it does not require any re-computation of neighbors. While simulating the steady-state problem, we iterate until the changes in pseudo-time become small enough.

3 SPH discretization

The basic scheme discussed in the previous section should work for any particular SPH discretization of the momentum and pressure equations. In the following, we use a WCSPH formulation for the SPH discretization. We keep density fixed as per the original problem. We consider two different cases one where we move the particles in pseudo-time and the other where the particles are fixed in pseudo-time.

3.1 Moving particles in pseudo-time

When the particles move in pseudo-time, the pressure evolution is computed using equation (5). This is discretized as,

[TABLE]

where $i$ denotes the $i$’th particle, $\mathbf{V}_{ij}=\mathbf{V}_{i}-\mathbf{V}_{j}$, $p_{ij}=p_{i}-p_{j}$, $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $\mathbf{x}_{ij}=\mathbf{x}_{i}-\mathbf{x}_{j}$, $W_{ij}=W(|\mathbf{r_{i}}-\mathbf{r_{j}}|,h)$, is the SPH smoothing kernel with $h$ as the smoothing length of the kernel and $\nu_{e}$ is the artificial viscosity coefficient used for the EDAC scheme [12] as shown in equation (6). The kernel is compact so the summation is over all the nearest neighbor particles that influence the particle $i$, $N(i)$. We note that the second term in the right hand side of equation (5) corresponds to a density diffusion term employed in the $\delta$-SPH scheme[10].

For the momentum equation where the particles move in pseudo-time given in (9), we can write,

[TABLE]

where $\Pi_{ij}$ is the artificial viscosity term [24] added to the momentum equation is given by,

[TABLE]

[TABLE]

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$, $h_{ij}=(h_{i}+h_{j})/2$, $\bar{\rho}_{ij}=(\rho_{i}+\rho_{j})/2$, $\bar{c}_{ij}=(c_{i}+c_{j})/2$, $\eta=0.01$, and $\alpha$ is the artificial viscosity parameter. Note that in the above equations the viscosity is modeled using the discretization proposed by Morris et al. [25].

We now write out the final form of the rate of change of the velocity in pseudo-time for the case where the particles are moving which we get by substituting equation (19) in the momentum equation (9), and discretizing the equations using the SPH formulation to get,

[TABLE]

It is important to note here that while we have used standard WCSPH discretizations, the DTSPH formulation would work just as well with any other SPH discretization that is based on a weakly-compressible formulation. While in principle this may be extended to Riemann-solver based SPH formulations [26], this is not explored in the present work.

We observe that in equation (34), we require the velocity at the current time and the previous time. When starting the simulations, if we do not have an exact solution, we assume that $\mathbf{V}^{-1}=\mathbf{V}^{0}$.

3.2 Fixed particles in pseudo-time

In the case the particles are fixed in pseudo-time. The right hand side of the equation (3) is discretized as,

[TABLE]

In the case of the EDAC scheme, the right hand side of the equation (7) is discretized as,

[TABLE]

For the momentum equation given in (8), we can write,

[TABLE]

By substituting equation (37) in equation (24) or equation (31) in equation (19), we get an equation for either $\frac{\partial\mathbf{V}}{\partial\tau}$ or $\frac{d\mathbf{V}}{d\tau}$. We can then use this to integrate the set of equations (12) – (17) or (21) – (23). Note that we first update the velocity and position of the particle as per equation (10) and (11).

The above equations govern the velocity and pressure of the fluids in the simulation. In order to satisfy the boundary conditions when solids are present we use an implementation of the boundary conditions presented in [27] where the pressure and the ghost velocity of the solid walls are set. Furthermore, following the work of [28], we ensure that the pressure is always positive on the solid walls to prevent particles from sticking to them in our free surface simulations. We do not impose any specific free-surface boundary conditions.

3.3 Steady state solutions

For the case of the steady state simulations, we do not move the particles in time and use the original distribution of particles. This is reasonable as steady solutions are usually sought where the geometry and boundaries are fixed. As discussed earlier, this leads to a very efficient solution procedure. For the pressure evolution we simply use equation (35). For the momentum equation we use the following discretization for equation (29),

[TABLE]

By solving these until the pseudo-time derivatives are small, we can obtain steady state solutions. This can be implemented very efficiently. The neighbors can be computed once and never need to be updated. The time step restrictions though continue to be as per the original weakly-compressible scheme.

3.4 Stability and convergence

It is important to choose the real and pseudo-timesteps carefully. We choose $\Delta t$ such that $\Delta t\mathbf{V}_{max}=Ch$, and $C$ is around 0.25. This is similar to the timesteps used in the ISPH schemes.

We choose $1<\beta<20$, recall that $\beta=1/M=c/V_{ref}$. We choose $\Delta\tau$ as $\Delta t/\beta$. The choice of these parameters is due to the following observations:

The real time step is limited by the amount of permitted motion of the particles in one time step.

• 2.

The pseudo-time step can be seen to be essentially similar to the original weakly-compressible scheme and is therefore limited by the speed of sound. The pressure waves travel at the speed of sound and therefore the pseudo-timesteps should be limited to around $\Delta t/\beta$.

We use the following approach to decide when to stop iterating in pseudo-time. The user specifies a particular tolerance, $\epsilon$, which determines the termination of the pseudo-time iterations. During every pseudo-time iteration we compute the mean rate of change in the pressure, let us call this quantity $\delta p/\delta\tau$. We also compute the mean value of $|\mathbf{\tilde{V}}|$ for all particles. When checking for convergence we ensure the following,

[TABLE]

Note that we multiply the rate of change of pressure by $\Delta t$ in order to ensure that the change over several pseudo-iterations would be accounted for.

Due to the inaccuracies of the SPH approximations and particle disorder, it is likely that the divergence does not become less than the tolerance and that the derivatives do not reduce. In order to prevent needless iterations we keep track of the changes in each pseudo-time iteration and stop iterations if the peak-to-peak relative changes in the last 3 or 4 iterations is less than 5%. This ensures that if the pressure and velocity do not change with increasing iterations we stop the iterations. This works very well in practice. We also stop iterating if there are more than 1000 iterations. In practice for reasonable tolerance values (larger than $10^{-5}$) we typically have far less than 50 iterations per real time step. Our default tolerance is $\epsilon=10^{-3}$.

3.5 Particle shifting

We use the shifting algorithm proposed by Lind et al. [29] which is based on Fickian diffusion. All the the test cases use shifting unless stated explicitly. The positions are shifted to their new position,

[TABLE]

where,

[TABLE]

where the present implementation uses $\Delta x=h$, the particle smoothing length. After shifting, all the fluid properties are corrected using Taylor series approximation as follows,

[TABLE]

where $\varphi_{i}$ and $\varphi_{i^{\prime}}$ denotes the fluid property before and after shifting respectively for $i^{\text{th}}$ fluid particle.

We apply the free-surface identification as suggested by Lind et al. [29] where a particle is considered as a free-surface particle if $\nabla\cdot\mathbf{r}_{i}<d-1/2$ where $d$ is the spatial dimension of the problem. Such free-surface particles are not shifted.

3.6 Complete algorithm

We summarize the proposed scheme using pseudo-code for clarity. The case of particles moving in pseudo-time is shown in Algorithm 1. The case of particles fixed in pseudo-time is shown in Algorithm 2.

4 Results and discussion

In this section we perform various numerical experiments using standard benchmark problems. We explore the following specific questions using the Taylor-Green problem which has an exact solution:

Is it worth moving the particles in pseudo-time or can we freeze the particles? This has significant performance implications.

• 2.

What possible values of $\beta$ can be used?

• 3.

What suitable values of the tolerance $\epsilon$ can be chosen and what does this imply for accuracy?

Once these are explored, a suite of test problems are simulated with the DTSPH scheme and compared with other schemes like the standard WCSPH [3], transport velocity formulation (TVF) [9], Entropically damped artificial compressibility (EDAC)  [12], and the $\delta$-SPH scheme [10]. Except where noted, the DTSPH scheme uses the EDAC formulation in all the simulations below i.e. the pressure evolution uses equation (36).

The TVF, EDAC, $\delta$-SPH, and WCSPH schemes are part of the PySPH [21, 22] framework. All the results presented below are automated and the code for the benchmarks is available at https://gitlab.com/pypr/dtsph. The tools used to automate the results are described in detail in [23]. This allows us to automatically reproduce every figure and table in this manuscript.

All the simulations are performed on a four core Intel (R) Core (TM) i$5-7400$ CPU with a clock speed of 3.00GHz. The problems are executed on four cores using OpenMP. We use the DTSPH scheme along with the EDAC in all results below unless explicitly mentioned otherwise.

4.1 Taylor-Green problem

The Taylor-Green problem is a classical problem which is periodic along both $x$ and $y$ axes and has an exact solution. This is a particularly challenging problem for SPH [9, 17, 19] since the particles move along the streamlines towards a stagnation point leading to particle disorder.

The exact solution for the Taylor-Green problem is given by

[TABLE]

where $U$ is chosen as $1m/s$, $b=-8\pi^{2}/Re$, $Re=UL/\nu$, and $L=1m$.

The Reynolds number $Re$ is set to 100 and various cases are tested to better understand the scheme. For all the simulations, the quintic spline kernel is used with $h/\Delta x=1.0$, no artificial viscosity is used. The following cases are considered,

a comparison of results when particles are either advected or frozen in pseudo-time;

• 2.

the effect of changing the artificial speed of sound, $\beta$;

• 3.

the effect of changing the convergence tolerance, $\epsilon$;

• 4.

comparison of results with different schemes; and

• 5.

comparison of the results with different number of particles and with different Reynolds numbers.

The results are compared against the exact solution. Particle plots are also shown wherever necessary as the error plots do not always reveal any particle disorder, particle clumping, or voiding occurring in the flow. In addition, we compare the performance of the schemes where the difference is noticeable.

We first discuss a simple method of initializing the particles that we use to compare all the schemes in a consistent and fair manner.

4.1.1 Perturbation in initial positions

As discussed earlier, it is not always easy to obtain good results for the Taylor-Green problem. When the particles are initially distributed in a uniform manner, they tend to move towards the stagnation points and this often leads to severe particle disorder. In [12], it was found that this problem can be reduced by introducing a small amount of noise in the initial particle distribution. To this end, a small random displacement is given to the particles with a maximum displacement of $\Delta x/10$. The random numbers are drawn from a uniform distribution and the random seed is kept fixed leading to the same distribution of particles for all cases with the same number of particles. This allows us to perform a fair comparison. Initially the particles are arranged in a $100\times 100$ grid, with smoothing length of $h/\Delta x=1$, $\Delta t=0.00125$, $\epsilon=10^{-4}$, and simulated for $2.5s$ with $\beta=5$.

The decay rate is computed by computing the magnitude of maximum velocity $|\mathbf{V}_{\max}|$ at each time step, the $L_{1}$ error is computed as the average value of the difference between the exact velocity magnitude and the computed velocity magnitude, given as

[TABLE]

where $\mathbf{V}_{i}$ is computed at the particle positions for each particle $i$ in the flow.

Fig. 1(a) shows the decay of the velocity compared with the exact solution for the case where there is no initial perturbation and with a small amount of perturbation (of at most $\Delta x/10$). As can be clearly seen, the introduction of the perturbation significantly improves the results. This is clearly seen in the $L_{1}$ norm of the error in the velocity magnitude in Fig. 1(b). While we have not shown this, the results are similar for simulations made using most other schemes including the new scheme, TVF, and the EDAC. Given this, we henceforth use a small initial perturbation for the results for the Taylor-Green vortex problem.

4.1.2 Advection of particles in pseudo-time

As discussed earlier, it is important to study the effect of moving the particles in pseudo-time as against keeping them frozen in pseudo-time. If we find that there is no significant advantage gained by moving the particles in pseudo-time, we can simplify the implementation of the scheme as well as improve its performance considerably.

We consider the Taylor Green problem and compare the results of simulations where we advect the particles and where we hold them frozen. The rest of the parameters are held fixed. The same small perturbation is added to the initial particle position. Parameters used for this simulation are, initial particle spacing $\Delta x=\Delta y=0.01$, Reynolds number $Re=100$, the value of $\beta=5$, time step $\Delta t=0.00125$, and tolerance of $\epsilon=10^{-6}$.

We recall that when we advect the particles in pseudo-time, we need to update the neighbors, however the displacements are very small and this is not necessary. We perform simulations to see if the differences are significant. Fig. 2(a) shows the decay rate for the case with and without advection in pseudo-time while updating the neighbours. There are no noticeable differences in the results and the plots for each case lie on each other. This is also seen in Fig. 2(b) which shows the $L_{1}$ error. Fig. 3(a) and Fig. 3(b) show the decay rate and the $L_{1}$ error in the velocity magnitude while not updating the neighbours resulting in a similar conclusion that movement of particles in pseudo-time is too small to significantly influence the results.

While the accuracy is unaffected, the performance is significantly different as can be seen from Table 1. This shows that advection of particles reduces performance by close to a factor of two. This increase in performance is largely due to the fact that we re-calculate the neighbour particles when we advect them. There is also some increase due to the additional computations required for the advection. Fig. 4 shows the particle plots with color representing velocity magnitude for the case where the particles are advected and frozen. The results look identical. Based on these results, we do not advect the particles in pseudo-time for any of the other simulations.

4.1.3 The influence of $\beta$

The parameter $\beta$ is the ratio of $c/V_{ref}$, as discussed earlier. The pseudo-timestep is also determined such that $\Delta t=\beta\Delta\tau$. In this section we consider the Taylor-Green problem simulated at $Re=100$ using $100\times 100$ particles, using the new scheme with different values of $\beta$ chosen between 2 and 20 for a tolerance $\epsilon=10^{-4}$, and run for a simulation time of $t=1$ sec.

Figs. 5(a) and 5(b) show the decay rate and the $L_{1}$ error in the velocity for the different cases. From these it appears that $\beta$ between 5 – 20 works well. We use a value of $\beta=10$ in all simulations unless explicitly mentioned.

4.1.4 Changing the convergence tolerance parameter $\epsilon$

We next choose $\beta=10$ and vary the tolerance from $10^{-2}$ to $10^{-5}$. Fig. 6(a) shows the decay rates as the tolerance is changed and Fig. 6(b) shows the $L_{1}$ error in the velocity. These results clearly show that reducing the tolerance improves the accuracy. However, we do see that the solutions are by-and-large robust to changes in $\epsilon$ over a very large range. As expected, increase in the tolerance leads to increase in the simulation time taken as seen in Table. 2.

4.1.5 Varying Reynolds number

We simulate the problem at $Re=1000$ using the most appropriate parameters based on the previous results. The DTSPH scheme is used with a quintic spline kernel, maximum initial random particle displacement of $\Delta x/10$, with a tolerance of $\epsilon=10^{-3}$, and $\beta=5$ for different initial particle arrangement of $25\times 25$ to $200\times 200$ and simulated for $2.5$ secs.

Fig. 7(a) shows the maximum velocity decay as well as the $L_{1}$ error of the velocity magnitude for the case of $Re=100$. Fig. 7(b) shows the same for $Re=1000$. It can be seen that in the case of $Re=1000$ that there is a clear reduction in the errors as the resolution is increased. In the case of $Re=100$, it appears that the errors are lower for the $50\times 50$ resolution and increase by a small amount as the resolution is increased. We believe that this occurs because of possible issues with the convergence of the discretization of the diffusive terms in the governing equations. We point out that at $Re=500$ we obtain similar convergence as in the case of $Re=1000$ showing that issue is not with the convergence of the DTSPH scheme per-se. These results show that the new scheme performs very well.

4.1.6 Comparision with other schemes

Here we simulate the problem for $t=2.5s$ for $Re=100$ using the new DTSPH scheme comparing it with WCSPH, $\delta$-SPH, and EDAC. The quintic spline kernel is used for all the schemes with $h=\Delta x$. For all the cases the particles are perturbed by atmost $\Delta x/10$. For DTSPH, we use $\beta=10$ with a tolerance of $\epsilon=10^{-4}$. We use an initial configuration of $100\times 100$ particles.

As can be seen from the results shown in Fig. 8, the new scheme is more accurate than the standard WCSPH scheme. The scheme is more accurate than the standard $\delta$-SPH scheme. We note that we do not employ any form of shifting for the $\delta$-SPH and WCSPH scheme cases. The scheme is not more accurate than the EDAC scheme as in the EDAC scheme, the particles are also regularized using the transport velocity formulation which significantly improves the results. As can be seen from the Table 3, the new scheme is anywhere from 1.9 to 2.8 times faster than the other schemes.

In Fig. 9, we show the streamlines as well as the particle plot with the color indicating pressure at $t=2.5$secs for the DTSPH case for $Re=100$. As can be seen the particle distribution is smooth.

4.2 The lid-driven-cavity problem

We next consider the classic lid-driven-cavity problem. This is a fairly challenging problem to simulate with SPH [9, 30, 31]. The fluid is placed in a unit square with a lid moving with a unit speed to the right. The bottom and side walls are treated as no-slip walls. The Reynolds number of the problem is given by $Re=\frac{V}{\nu}$, where $V$ is the lid velocity. We use a quintic spline kernel with $h=\Delta x$. The problem is simulated at $Re=100$ using a $50\times 50$, $100\times 100$, and $150\times 150$ grid for a simulation time of $t=10s$ until there is no change in the kinetic energy of the system. For the DTSPH scheme we use a $\beta=10$ with a tolerance $\epsilon=10^{-4}$. The results are compared with those of the TVF scheme[9] and the established results of Ghia et al. [32].

Fig. 10, shows the centerline velocity profiles for $u$ vs. $y$ and $v$ vs. $x$ for different resolutions of particles. It is seen that the TVF scheme produces better results as expected. However, the results of the new scheme are in good agreement. In Fig. 11, we show the particle distribution with the velocity magnitude on the left and streamlines on the right.

4.2.1 Steady Lid-driven cavity

In order to show that we are able to obtain steady state results, we employ the steady state equations discussed in Section 2.3 to solve the lid-driven-cavity problem. We solve the problem until there is no change in the kinetic energy of the system. We simulate the problem using a quintic spline kernel for $Re=100$ and $Re=1000$ using a $50\times 50$, $100\times 100$ and $150\times 150$ particle grid. For $Re=100$ we simulate the problem up to $\tau=10$ and for $Re=1000$ we simulate up to $\tau=50$.

Fig. 12 shows the velocity profiles for the $Re=100$ case and Fig. 13 shows velocity profiles for the $Re=1000$ case.

These results show that we are able to simulate internal flows very well using the new DTSPH scheme. We have also demonstrated that the steady-state equations also work very well. We next consider problems that involve a free-surface.

4.3 Square patch

The square patch problem [33, 34, 35] is a free surface problem where a square patch of fluid of side $L$ is subjected to the following initial conditions,

[TABLE]

[TABLE]

where $X^{*}=x+L/2$ and $y^{*}=y+L/2$.

We simulate this problem for $t=3s$ using the DTSPH, and EDAC schemes for comparison. In this case, the EDAC simulations also employ particle shifting [29]. The quintic spline kernel with $h/\Delta x=1.3$ is used for all the schemes, artificial viscosity $\alpha=0.1$ is used for all the schemes. For the DTSPH scheme, $\beta=10$ with a tolerance $\epsilon=10^{-3}$ is used. Two different initial configurations of $50\times 50$ and $200\times 200$ particles are used.

The particle distribution for each scheme at the end of $t=3s$ is shown in Fig. 14. The plots of DTSPH and EDAC are in good agreement with each other showing that the new scheme is as good as the EDAC scheme.

4.4 Elliptical drop

The elliptical drop problem was first solved in the context of the SPH by Monaghan [3]. This problem is also solved in the context of truly incompressible SPH [34, 36, 37]. In this problem an initially circular drop of inviscid fluid having unit radius is subjected to the initial velocity field given by $-100x\hat{i}+100y\hat{j}$. The outer surface is treated as a free surface. Due to the incompressibility constraint on the fluid there is an evolution equation for the semi-major axis of the ellipse.

This problem is simulated using the DTSPH, $\delta$-SPH (without shifting), and EDAC (with shifting) respectively. An artificial viscosity parameter of $\alpha=0.15$ is used for all the schemes. An error tolerance of $\epsilon=10^{-4}$ is used for the DTSPH and scheme. $\beta=10$, $\Delta x=0.02$, $h=1.3\Delta x$ and a quintic spline kernel is used for all the schemes. The simulation is run for $t=0.0076s$.

Fig. 15 shows the distribution of particles for different schemes. The colors indicate the pressure. As can be seen, the DTSPH and EDAC results are similar. It is important to note that all the pressure values are in a similar range with none of the schemes exhibiting severe noise in the pressure.

Fig. 16, shows the evolution of the kinetic energy, the results are similar for DTSPH and EDAC schemes. This is to be expected. However it is interesting to note that the kinetic energy of the $\delta$-SPH scheme decays faster than the other schemes. Fig. 17 shows the error in the semi-major axis as compared to the exact solution. The DTSPH and EDAC both perform slightly better than the $\delta$-SPH scheme. The results indicate that the new scheme performs well in comparison with state of the art weakly-compressible schemes.

4.5 Dam-break in 2 dimensions

A two dimensional dam-break over a dry bed [31, 11, 36] is considered next. The DTSPH and the standard EDAC schemes are compared. The simulation is performed for $1s$. The quintic spline kernel is used with $h/\Delta x=1.0$, and an artificial viscosity of $\alpha=0.1$ is used for all the schemes. A tolerance of $\epsilon=10^{-4}$ is used for DTSPH.

The problem considered is described in [38] with a block of fluid column of height $h=2m$, width $w=1m$. The block is released under gravity which is assumed to be $-9.81m/s^{2}$.

For the DTSPH and EDAC schemes, the particle distribution is shown in Fig. 18 at various times with color indicating pressure. Fig. 19 shows the particle distribution with color indicating velocity magnitude at various times. The results of the new scheme seem largely comparable with that of the EDAC scheme. The results also show the improvements obtained by the addition of diffusive term in the pressure evolution equation. This significantly reduces the noise. The standard EDAC results are very similar to those of the $\delta$-SPH and are hence not shown. Fig. 20 plots the position of the toe of the dam versus time as compared with the results of the Moving Point Semi-implicit scheme of [39]. The results of the DTSPH scheme are in good agreement with those of the EDAC scheme.

4.6 Dam-break in three dimensions

A three dimensional case is shown to demonstrate the performance of the new scheme as compared to the EDAC scheme. This is an important case as the previous problems only require a smaller number of particles. We consider a three-dimensional dam break over a dry bed with an obstacle. A cubic spline kernel is used for both the new scheme and the EDAC with $h/\Delta x=1.3$ and artificial viscosity $\alpha=0.1$. The problem is simulated for a total time of 1 second. We do not use any particle shifting in this case.

The problem considered is described in [38] with a block of fluid column of height $h=0.55m$, width $w=1.0m$ and length $l=1.228m$. The container is $3.22m$ long. The block is released under gravity with an acceleration of $-9.81m/s^{2}$. Both schemes are simulated with a fixed time step. For both schemes we use a CFL of 0.25. The speed of sound for the EDAC case is set to $10\sqrt{2gh}$, where $h$ is the height of the water column. For the DTSPH, we use a time step of $\frac{0.25h}{\sqrt{2gh}}$, and use $\beta=10,\epsilon=10^{-3}$. The particle spacing, $\Delta x=0.02$, leading to around 240000 particles in the simulation.

Fig. 21 shows the particle distribution at various times for both the schemes. This indicate that the new scheme produces good results. Table 4 shows the time taken for the different 3D dam break simulations. Depending on the tolerance chosen, we are able to obtain between a 1.7 to 7.16 fold improvement in performance as compared to the EDAC scheme. We note that as the tolerance $\epsilon$ is reduced, the DTSPH requires more iterations in pseudo-time in order to attain convergence and this reduces the performance. However even with lower tolerance values we are able to obtain very good results. In the next section we demonstrate the performance achievable with the new scheme for different problems.

4.7 Performance

The performance of DTSPH is compared with that of other schemes for various problems. In Table 5 we list the different problems and the speedup obtained by using the new scheme.

As we have seen before, for the dam-break problem in three dimensions it can be seen that DTSPH (with a tolerance of $10^{-3}$) can be up to 7.16 times faster than the standard EDAC scheme. When a very low tolerance is used, the scheme is about 1.7 times faster this is only to be expected as lower tolerances require many more iterations to obtain a converged pressure. For the two-dimensional dam break cases, the new scheme is about 3.5 times faster than the EDAC or the $\delta$-SPH. For the unsteady cavity problem at $Re=100$ with a $150\times 150$ grid of particles, we get up to a 7 times performance improvement. This is quite significant since in these cases we have compared these cases where the effective Mach number is 0.1. Given this, the primary advantage with the DTSPH is that it can take a time step that is 10 times smaller. Hence in this case a speed-up of close to 7 suggests that it is an efficient scheme.

The suite of benchmarks considered shows that the new scheme is robust, simulates a variety of problems, and is as accurate as the EDAC scheme. In addition it is very efficient and can be as much as seven times faster than the EDAC scheme. Indeed, it is possible to improve the performance even more by caching the values of the kernel and kernel gradients during the pseudo-time iterations but we have not done this in the present work.

5 Conclusions

In this paper we propose a scheme called Dual-Time stepped SPH (DTSPH) that employs a dual-time stepping approach for incompressible fluid flow simulations. The method has been demonstrated with the EDAC formulation [12]. There are a few recent developments in the literature that employ a similar approach [17, 19] but these implementations are not efficient despite the ability to use much smaller timesteps than the WCSPH scheme. We find that by not moving the particles in pseudo-time we are able to improve performance significantly without loss of accuracy as attested by our simulations. We show that the scheme is robust and accurate. Through several benchmarks in two and three dimensions we show that the scheme produces results that are as accurate as the $\delta$-SPH scheme [10] as well as the EDAC scheme [12] while being up to seven times faster. The method is matrix-free and may be implemented in the context of any explicit SPH scheme. The DTSPH method does introduce a few new parameters in the form of the term $\beta=\Delta t/\Delta\tau$ and the tolerance $\epsilon$. We discuss how these parameters can be set based on rational considerations. For a reasonable choice of parameters, the performance of the method is comparable to that of incompressible SPH schemes. An open source implementation of the new scheme is provided and the manuscript is fully reproducible. Given the results presented with this method, a promising line of future work would be to apply this to the case of solid mechanics problems where the timesteps are very small. A method of this kind would provide significant computational advantages. This would also be of considerable importance in the area of fluid-structure-interaction which is an important emerging area of investigation with the SPH method.

Acknowledgments

We are grateful to the anonymous reviewers for their comments that have improved the quality of this manuscript.