An efficient, open source, iterative ISPH scheme

TL;DR

This paper introduces a simple, robust, and efficient open-source iterative ISPH scheme suitable for CPU and GPU implementation, featuring a matrix-free approach and novel boundary condition handling, demonstrated through various benchmarks.

Contribution

The paper presents a new iterative, matrix-free ISPH method that improves boundary conditions and is easily implementable on GPUs, with open-source code and reproducible results.

Findings

The method is fast and GPU-compatible.

It maintains accuracy comparable to existing approaches.

It demonstrates robustness across multiple benchmark problems.

Abstract

In this paper a simple, robust, and general purpose approach to implement the Incompressible Smoothed Particle Hydrodynamics (ISPH) method is proposed. This approach is well suited for implementation on CPUs and GPUs. The method is matrix-free and uses an iterative formulation to setup and solve the pressure-Poisson equation. A novel approach is used to ensure homogeneous particle distributions and improved boundary conditions. This formulation enables the use of solid wall boundary conditions from the weakly-compressible SPH schemes. The method is fast and runs on GPUs without the need for complex integration with sparse linear solvers. We show that this approach is sufficiently accurate and yet efficient compared to other approaches. Several benchmark problems that illustrate the robustness, performance, and wide range of applicability of the new scheme are demonstrated. An open source implementation is provided and the manuscript is fully reproducible.

Tables

Table 1. Table 1 : Average number of iteration taken for Jacobi for various tolerances and problems.

Problem	Tolerance	Jacobi (avg no of. iterations)
Cavity	0.001	3.0
Dam-break 2D	0.001	5.0
Cavity	0.01	2.0
Dam-break 2D	0.01	2.0

Equations

$$\nabla \cdot \mathbf{u}=0$$

$$\frac{d\mathbf{u}}{dt}=-\frac{1}{\rho }\nabla p+\nu \nabla \cdot (\nabla \mathbf{u})+\mathbf{f}$$

$${\mathbf{u}}^{\ast }=\mathbf{u}+\frac{\Delta t}{\rho }\nabla p$$

$$\nabla \cdot (\frac{1}{\rho }\nabla p)=\frac{\nabla \cdot {\mathbf{u}}^{\ast }}{\Delta t}$$

$$\frac{\tilde{d}\mathbf{u}}{dt}=-\frac{1}{\rho }\nabla p+\nu \nabla \cdot (\nabla \mathbf{u})+\frac{1}{\rho }\nabla \cdot \mathbf{A}+\mathbf{f}$$

$${\rho }_i=\sum _{j\in N(i)}{m}_j{W}_{ij}$$

$$\frac{\tilde{d}{\mathbf{u}}_i}{dt}={\mathbf{f}}_{\text{p}}+{\mathbf{f}}_{\text{visc}}+{\mathbf{f}}_{\text{b}}+{\mathbf{f}}_{\text{avisc}}+{\mathbf{f}}_{\text{astress}}$$

$${\mathbf{f}}_{\text{p, symm}}=-\sum _{j\in N(i)}{m}_j\left(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}\right)\nabla W_{ij}$$

$${\mathbf{f}}_{\text{p, asymm}}=-\sum _{j\in N(i)}\frac{m_j}{{\rho }_i{\rho }_j}(p_j-p_i)\nabla W_{ij}$$

$${\mathbf{f}}_{\text{visc}}=\sum _{j\in N(i)}{m}_j\frac{4\nu }{({\rho }_i+{\rho }_j)}\frac{{\mathbf{r}}_{ij}\cdot \nabla W_{ij}}{(|{\mathbf{r}}_{ij}{|}^2+\eta h_{ij}^2)}{\mathbf{u}}_{ij}$$

$${\mathbf{f}}_{\text{avisc}}=\sum _{j\in N(i)}{\Pi }_{ij}\nabla W_{ij}$$

$${\Pi }_{ij}=\{\begin{array}{ll}\frac{-\alpha h_{ij}{\bar{c}}_{ij}{\varphi }_{ij}}{{\bar{\rho }}_{ij}},& {\mathbf{u}}_{ij}\cdot {\mathbf{r}}_{ij}<0\\ 0,& {\mathbf{u}}_{ij}\cdot {\mathbf{r}}_{ij}\ge 0.\end{array}$$

$${\mathbf{f}}_{\text{astress}}=\sum _{j\in N(i)}\left(\frac{{\mathbf{A}}_i}{{\rho }_i^2}+\frac{{\mathbf{A}}_j}{{\rho }_j^2}\right)$$

$${\tilde{\mathbf{u}}}_i^{n+1}={\mathbf{u}}_i+\Delta t\left(\frac{\tilde{d}{\mathbf{u}}_i}{dt}+{\mathbf{f}}_{i,\text{GTVF}}\right)$$

$${\mathbf{f}}_{i,\text{GTVF}}=-p_i^0\sum _{j\in N(i)}\frac{m_j}{{\rho }_j^2}\nabla W({\mathbf{r}}_{ij},{\tilde{h}}_{ij}).$$

$$\nabla \cdot {\left(\frac{1}{\rho }\nabla p\right)}_i=\frac{\nabla \cdot {{\mathbf{u}}^{\ast }}_i}{\Delta t}$$

$$\sum _j\left(\frac{4m_j}{{\rho }_i({\rho }_i+{\rho }_j)}\right)\frac{p_{ij}{\mathbf{r}}_{ij}\cdot {\nabla }_i{W}_{ij}}{|r_{ij}^2|+\eta h_{ij}^2}$$

$${\mathbf{r}}_i^{\ast }={\mathbf{r}}_i^n+\Delta t{\tilde{\mathbf{u}}}_i^n$$

$${\mathbf{u}}_i^{\ast }={\mathbf{u}}_i^n+\Delta t({\mathbf{f}}_{i,\text{visc}}^n+{\mathbf{f}}_{i,\text{avisc}}^n+{\mathbf{f}}_{i,\text{astress}}^n+{\mathbf{f}}_{i,\text{body}}^n).$$

$${\mathbf{u}}_i^{n+1}={\mathbf{u}}_i^{\ast }+\Delta t{\mathbf{f}}_{i,\text{p}}^n.$$

$${\mathbf{r}}^k={\mathbf{r}}^{(k-1)}+\Delta \tau {\tilde{\mathbf{u}}}^{k-1}+\frac{1}{2}(\Delta \tau {)}^2\mathbf{f}({\mathbf{r}}^{k-1}{)}_{\text{GTVF}}$$

$${\tilde{\mathbf{u}}}^k={\tilde{\mathbf{u}}}^{(k-1)}+\Delta \tau \mathbf{f}({\mathbf{r}}^{(k-1)}{)}_{\text{GTVF}}$$

$${\tilde{\mathbf{u}}}_i^{n+1}={\mathbf{u}}_i^{n+1}+\frac{{\mathbf{r}}^K-{\mathbf{r}}^0}{\Delta t}$$

$${\mathbf{r}}_i^{n+1}={\mathbf{r}}_i^n+\Delta t\left(\frac{{\tilde{\mathbf{u}}}_i^{n+1}+{\tilde{\mathbf{u}}}_i^n}{2}\right)$$

$$D_{ii}{p}_i^{k+1}={\text{RHS}}_i-\sum _j{\text{OD}}_{ij}{p}_j^k,$$

$$D_{ii}$$

$${\text{OD}}_{ij}$$

$$p_i^{k+1}=\omega \frac{({\text{RHS}}_i-{\sum }_j{\text{OD}}_{ij}{p}_j^k)}{D_{ii}}+(1-\omega )p_i^k$$

$$\sum _i\frac{|p_i^{k+1}-p_i^k|}{\max (\Lambda ,{\sum }_i|p_i^{k+1}|)}<ϵ,$$

$$\Delta t_{\text{cfl}}=0.25\frac{h_{ij}}{|\mathbf{U}|}$$

$$\Delta t_{\text{visc}}=0.25\frac{h_{ij}^2}{\nu }$$

Full text

Program summary

Program title: SISPH

Licensing provisions: BSD 3-Clause

Programming language: Python

External routines/libraries: PySPH (https://github.com/pypr/pysph), scipy (https://pypi.org/project/scipy/), matplotlib (https://pypi.org/project/matplotlib/), compyle (https://pypi.org/project/compyle/), automan (https://pypi.org/project/automan/).

Nature of problem: Incompressible SPH solvers solve a sparse system of equations to obtain a solution to the pressure-Poisson equation. This is often implemented by solving a sparse matrix, and can be complex to implement on GPUs. SPH accuracy deteriorates when particles become disordered. Ensuring homogeneous particle distributions requires particle shifting algorithms which requires tuning of parameters. We propose the use of matrix-free algorithm and an accurate and general purpose way to ensure particle homogeneity.

Solution method: We use an explicit iterative approach to solve the pressure-Poisson equations which makes it easy to implement in parallel on GPUs. This also makes it easy to implement the boundary conditions. We propose an accurate way to ensure homogeneous particle distribution by using a background pressure.

Additional comments: The source code for this repository can be found at https://gitlab.com/pypr/sisph.

1 Introduction

The Smoothed Particle Hydrodynamics method was originally proposed to simulate astrophysical hydrodynamic problems by Lucy [1] and Gingold and Monaghan [2]. It is a mesh-free, Lagrangian, particle-based method that has since been used to simulate incompressible fluids. The Weakly-Compressible SPH (WCSPH) [3] scheme was proposed to simulate incompressible flows by treating the fluid as weakly compressible and using a stiff equation of state. This method has been used to simulate fluid flows with a free-surface. One significant difficulty with the WCSPH method is the fact that it introduces an artificial sound speed which is typically around 10 times the maximum speed of the flow. This introduces severe timestep limitations in the scheme.

The Incompressible SPH (ISPH) schemes have their origin in the work of Cummins and Rudman [4] who proposed a projection method where a pressure-Poisson equation (PPE) is solved to ensure a divergence free velocity field. In this, the Laplacian of the pressure is related to the divergence of the velocity field. Shao and Lo [5] proposed a slightly different variant which relates the Laplacian of the pressure to a change in density. Later, Hu and Adams [6] proposed a method that combines both of these approaches to enforce incompressibility. The significant advantage with the class of ISPH schemes is that they do not involve the sound speed and hence can use timesteps that are an order of magnitude larger than the WCSPH schemes.

The ISPH schemes have the disadvantage that they require the solution to large (but sparse) linear systems. This makes implementing them in parallel and on GPUs fairly difficult. The WCSPH implementations are generally more popular as they are much easier to implement and parallelize. The recent work of Chow et al. [7] discusses many of the challenges in implementing traditional ISPH schemes for large scale computing on GPUs.

An explicit version of the ISPH (EISPH) was first proposed by Hosseini et al. [8] which focuses on simulating non-Newtonian flows. This approach removes the requirement to solve a linear system and instead sets up an explicit (non-iterative) equation to compute the pressure. A similar explicit approach was used by Rafiee and Thiagarajan [9] for solving fluid-structure interaction problems. Barcarolo [10, 11] validates and explores the issues with EISPH and compare the scheme with the WCSPH in the context of internal and free surface flows. Nomeritae et al. [12] assess the performance of EISPH with a comparison with WCSPH and $\delta$-SPH schemes [13]. Basser et al. [14] simulate multi-fluids with porous media using EISPH and compared their results with experimental data. While Hosseini et al. [8], Rafiee and Thiagarajan [9], and Barcarolo et al. [11] solve for the pressure by satisfying a constant density condition [5], Nomeritae et al. [12] and Basser et al. [14] use the divergence free condition [4] to obtain the pressure. It is important to note that in all these cases, the authors perform a single step to solve the PPE and do not perform any iterations thus making the method a truly explicit ISPH method.

In the graphics community, the Implicit ISPH (IISPH) scheme [15] has been developed that provides an iterative solution to the ISPH formulation. This formulation is tied to the way in which the pressure forces are approximated between pairs of particles. The method does not work well with negative pressures. Our own implementation of this scheme demonstrates some sensitivity to changes in the timestep or choice of smoothing kernels. The method is however, matrix free, and very fast.

In the present work, we consider the projection-based ISPH scheme of Cummins and Rudman [4]. We use an iterative approach to solve the PPE as constrasted with the EISPH approach. This makes the formulation completely matrix-free, and easy to derive and implement. This approach is not new and is similar to the approach used in the EISPH [8, 10, 12, 14] with the key difference being that we iteratively solve the same PPE. This makes a significant difference when we simulate problems at higher Reynolds numbers. We use a successive-over-relaxation (SOR) procedure to accelerate the convergence of our iterations. Our implementation is efficient, and works comfortably on a GPU architecture. This reformulation also allows us to implement boundary relatively easily as is done in WCSPH schemes and this is harder to do with the traditional ISPH. We note that despite the existence of the similar EISPH for many years, the traditional approach of solving the system of equations using a separate linear solver continues to be popular [7].

Both the WCSPH and ISPH schemes suffer from inaccuracies when the particles become disordered. In flows with significant shear, the particles can become significantly disordered leading to poor accuracy and particle clumping in extreme cases. The ISPH community has developed a few particle shifting techniques [16, 17, 18] that add a small amount of particle motion to ensure uniform distributions. These improve the particle distribution at the cost of a negligible amount of diffusion. It does requires some tuning to handle free surfaces.

For the WCSPH schemes, the Transport Velocity Formulation (TVF) [19] and the Generalized TVF [20] provide a slightly different and more accurate way of generating homogeneous particle distributions. The method relies on the fact that the naive SPH gradient of a constant pressure field will generate forces that tend to push particles into a uniform configuration. Hence, the particles are moved with a “transport velocity” which is due to the fluid forces along with this constant pressure force. The background pressure is subtracted from the main momentum equation through an artificial stress term. The Generalized TVF [20], extends the TVF to free-surface flows. In this paper we apply the GTVF to work with ISPH schemes. Since the timesteps are larger, we use a small modification to move the particles during the regularization.

Although Barcarolo et al. [11] explore internal flows and flow past a circular cylinder, they do so for low Reynolds numbers where the chances of particle disorder are less. In fact, in [10], it is found that there is some particle voiding present for the case of a lid driven cavity at a Reynolds number of 3200. Moreover, in EISPH, the PPE is solved using a single step. This approach does not work well for high Reynolds number cases and we demonstrate these issues in the current work. Previous works employing the EISPH do not explore the performance on GPUs where the scheme has a high potential for rapid, accurate simulation. Furthermore, they do not simulate the Taylor-Green problem which typically fails to work unless a careful implementation of particle shifting is included. Our implementation, on the other hand works well for this problem and also demonstrates good performance on GPUs.

In order to handle solid boundaries accurately, we modify the method proposed by Adami et al. [21] to work with the ISPH schemes. The original method does not work too well with the ISPH due to the larger timesteps that the method allows. We provide a simple way to compute the normals on a solid body with just the point information, and use this to improve the solid boundary condition. The resulting method works effectively for a variety of test cases.

In summary, the proposed method is easy to implement, accurate, matrix-free, fast, and works on GPUs. We apply the transport velocity formulation to ensure that the particle distribution is uniform and suggest an improved solid wall boundary condition. We demonstrate the method with several standard benchmarks. These include internal flows, external free-surface problems, wind-tunnel type problems requiring an inlet and an outlet and also demonstrate good performance on a GPU. We demonstrate the performance of the new scheme and compare this with the performance of the traditional matrix-based ISPH. We also demonstrate the performance achieved on different GPUs as well as multi-core CPUs.

Our entire implementation is open source. We employ the open source PySPH [22, 23] framework for the implementation of the scheme. In the interest of reproducible research, our entire manuscript is reproducible. Every figure presented in the results section of this manuscript is automated [24] and the code for the computations is made available at https://gitlab.com/pypr/sisph.

2 The ISPH method

The basic formulation is that of Cummins and Rudman [4]. For an incompressible fluid, the continuity equation is,

[TABLE]

where $\mathbf{u}$ is the velocity field, and the momentum equation is given by,

[TABLE]

where $\frac{d\mathbf{u}}{dt}$ is the material derivative of the velocity field, $\rho$ is density, $p$ is pressure, $\nu$ is kinematic viscosity and $\mathbf{f}$ is the external force.

Pressure is obtained by the projection method, which uses Hodge decomposition to project any velocity field, $\mathbf{u}^{*}$ into a divergence-free component, and a curl-free component,

[TABLE]

where $\mathbf{u}$ is the divergence free velocity field, $\nabla p$ is the curl-free component. the pressure is then obtained by taking divergence of equation (3), which reads,

[TABLE]

where $\mathbf{u}^{*}$ is the intermediate velocity obtained by integrating the momentum equation without considering the pressure gradient terms [4]. After solving for the pressure, $p$ (curl-free component), the divergence free velocity $\mathbf{u}$ is obtained by subtracting the gradient of pressure from the intermediate velocity $\mathbf{u}^{*}$.

The ISPH like most SPH methods suffers from particle disorder. We consider the Generalized Transport Velocity Formulation (GTVF) by Zhang et al. [20] to remedy the particle disorder. GTVF extends the Transport Velocity Formulation (TVF) by Adami et al. [19] to free surface flows and solid dynamics by using variable background pressure instead of a constant global background pressure. In the GTVF formulation, positions and velocity are advected using the transport velocity, which introduces an additional artificial stress term to the momentum equation. The transport velocity is regulated by variable background pressure, which keeps the particles in a uniform configuration. The momentum equation for GTVF is,

[TABLE]

where $\frac{\tilde{d}\mathbf{u}}{dt}=\frac{\partial\mathbf{u}}{\partial t}+\tilde{\mathbf{u}}\cdot\nabla\mathbf{u}$ is the material derivative of velocity field which is advected by the transport velocity, $\tilde{\mathbf{u}}$, and $\mathbf{A}=\rho\mathbf{u}\otimes(\tilde{\mathbf{u}}-\mathbf{u})$ is the artificial stress term.

2.1 SPH discretization

We apply SPH discretization to the scheme discussed above. In all the simulations we do not keep the density constant, density is computed by the summation density,

[TABLE]

where the subscript $i$ denotes the $i^{\text{th}}$ particle, $W_{ij}=W(|\mathbf{r}_{ij}|,h_{ij})$ is the SPH smoothing kernel, where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $h_{ij}=(h_{i}+h_{j})/2$, $\mathbf{r}$ is the position of the particle and $h$ the smoothing length of the kernel. The summation is carried over $N(i)$ neighbors, of the $i^{\text{th}}$ particle. We use the density in order to provide a better estimate of the particle volume given by $m/\rho$.

The SPH discretization of the momentum equation (5) is given by,

[TABLE]

where $\mathbf{f}_{\text{p}}$ is the force due to pressure gradient, $\mathbf{f}_{\text{visc}}$ is the force due to viscous forces, $\mathbf{f}_{\text{b}}$ is the body force, $\mathbf{f}_{\text{avisc}}$ is the force due to artificial viscosity and $\mathbf{f}_{\text{astress}}$ is the force due to artificial stress which is a consequence of the GTVF formulation.

In the literature two ways of computing pressure gradient are encountered. One is a symmetric, momentum-preserving form [6, 4] and the other is asymmetric and does not preserve momentum [16, 18, 17]. In the present work the symmetric form is labelled as “symm” and the asymmetric form as “asymm”. We have used both the forms in this work and report the differences between them in our results. The SPH discretization of the symmetric form [25] is,

[TABLE]

and the asymmetric form [25] is,

[TABLE]

where $\nabla W_{ij}$ is the gradient of the smoothing kernel. Viscous forces are discretized as [4],

[TABLE]

where $\mathbf{u}_{ij}=\mathbf{u}_{i}-\mathbf{u}_{j}$, $\nu$ is the kinematic viscosity, $\eta=0.01$. Artificial viscosity [26] is added to the momentum equation wherever necessary given by,

[TABLE]

where, $\Pi_{ij}$ is computed by,

[TABLE]

The artificial stress due to the GTVF formulation [20] is evaluated using,

[TABLE]

where, ${A}_{i}=\rho_{i}\mathbf{u}_{i}\otimes(\tilde{\mathbf{u}}_{i}-\mathbf{u}_{i})$. The transport velocity, $\tilde{\mathbf{u}}$ due to GTVF [20] is updated by adding an additional background pressure to the simulation,

[TABLE]

where, $\frac{\tilde{d}\mathbf{u}_{i}}{dt}$ is obtained from equation (7) and

[TABLE]

On numerical investigation, we suggest $\tilde{h}=0.5h_{ij}$ for external and free surface flows, and $\tilde{h}_{ij}=h_{ij}$ for internal flows, $p^{0}_{i}=\min(10|p_{i}|,p_{\text{ref}})$ for external flows, and $p^{0}_{i}=p_{\text{ref}}$ is kept constant for internal flows. The choice of reference pressure, $p_{\text{ref}}$, is typically taken to be $\rho c^{2}$ for external flows and $2\max(p_{i})$, i.e., twice the maximum pressure in the entire flow after the first iteration for internal flows. Here, we assume that $c=10|\mathbf{u}|_{\text{max}}$ and this is simply a reference pressure.

The pressure-Poisson equation that is solved is,

[TABLE]

where $\Delta t$ is the timestep. The approximate projection form for the PPE operator in equation (16) given by Cummins and Rudman [4] is used, resulting in the discretized form as,

[TABLE]

where $p_{ij}=p_{i}-p_{j}$.

2.2 Time Integration

Due to its simplicity we use the time integration as given by Cummins and Rudman [4]. However, the higher order method proposed in Nair and Tomar [27] could also be used. Time integration is performed by first calculating the intermediate position using the transport velocity at the current time,

[TABLE]

The intermediate velocity is obtained by integrating the momentum equation neglecting the force due to the pressure gradient,

[TABLE]

The pressure-Poisson equation given by (17) is solved to obtain the pressure. The divergence free velocity is then found by correcting the intermediate velocity as,

[TABLE]

We use modified GTVF to redistribute particles. The GTVF formulation is first used in the context of weakly compressible flows, where the time step is much lower than ISPH schemes. The timesteps are large in the ISPH and applying the GTVF correction in one timestep introduces stability issues. We modify the GTVF formulation by applying the GTVF force ($\mathbf{f}_{\text{GTVF}}$, equation (15)) to shift the particles iteratively in one time step. We then compute the transport velocity using the initial and final positions after the iteration.

The GTVF acceleration is a function of positions $\mathbf{r}$, background pressure $p_{\text{ref}}$ and smoothing length $h_{ij}$. In all our simulations $h_{ij}$ is fixed and for a given time the background pressure remains unchanged. We use $K$ sub-timesteps to integrate the GTVF accelerations without recomputing any neighbors such that $\Delta\tau=\Delta t/K$. If we denote $\mathbf{r}^{k},\mathbf{\tilde{u}}^{k}$ as the position and GTVF velocity of a particle at the $k$’th sub-iterate, we first start with $\mathbf{r}^{0}$ set to the current particle position, and set $\mathbf{\tilde{u}}^{0}=0$.

[TABLE]

[TABLE]

After the GTVF iterative step is done, we compute the transport velocity by looking at the inital and final positions

[TABLE]

particles are then advected to the next time step using transport velocity of current and previous time steps given by,

[TABLE]

3 Iterative formulation

Traditionally, equation (17) is solved by setting up a sparse matrix corresponding to the coefficients of $p_{i}$ and $p_{j}$ in the equation and computing the right-hand side. This system is then solved using any fast, sparse matrix solver. In the EISPH [8, 10, 12, 14] a simple explicit formula is provided to update the pressure.

We use this approach to write this system in terms of a Jacobi iteration we introduce a time index, $k$, and write $p^{k}$ as the pressure at the $k$th iteration. We rewrite equation (17) in the form,

[TABLE]

where $\text{OD}_{ij}$ represents the off-diagonal coefficients which are multiplied by pressure at the earlier iteration $p^{k}$ and $\text{RHS}_{i}$ is the right-hand-side of equation (17). We can easily write the following expressions for the coefficients in equation (25) as,

[TABLE]

To begin with, we set $p^{k=0}(t+\Delta t)=p(t)$, that is our first iteration starts with the existing value of pressure of the particle. We find the pressure for the next iteration using Successive-Over-Relaxation (SOR) as,

[TABLE]

The default value of $\omega=0.5$. We note that the entire term $\sum_{j}\text{OD}_{ij}p_{j}^{k}$, needs to be accumulated into a single term and hence only requires a single number per particle for storage. Thus, per particle we require storing three terms, viz. $D_{ii},\sum_{j}\text{OD}_{ij}p_{j}^{k}$, and $\text{RHS}_{i}$. We note that if a particle has no neighbors, then $D_{ii}$ can be zero, in which case we set the pressure of the particle to zero. Similarly, to handle free surfaces we compute the summation density of the particles and if the ratio $\rho_{i}/\rho_{0}<0.8$, where $\rho_{0}$ is the reference or rest density, we treat the particle as a free-surface particle and set its pressure to zero. Note that during these pressure iterations we do not move the particles. We continue to iterate until the relative change in the pressure is less than some user-specified amount,

[TABLE]

where $\Lambda=\max(RHS_{i}/D_{ii})$. When the mean pressure in the flow is less than one, we use the largest value of $RHS_{i}/D_{ii}$ as a measure of the mean pressure and this conveniently avoids zeros in the denominator while also using a relative scale for the error. The minimum number of iterations performed is set to two.

We vary the tolerance parameter, $\epsilon$, in our numerical studies. We consider the suitable choice of this tolerance parameter next.

3.1 Choice of the tolerance parameter

While iterating for convergence, we note that using very small values of $\epsilon$ is not necessary. We first note that the discretized terms in equation (26) involve errors of at best $O(h^{2})$ so there is no major advantage to solving the linear system to high accuracy.

Furthermore, for many timesteps, the pressure difference at a point between two timesteps is a small quantity and since we start with an accurate pressure, we expect that obtaining the new pressure will take very few iterations. By setting a tolerance of say $0.001$, we ensure that if there are any changes to the pressure during the Jacobi iterations that are less than 0.1% that we stop iterations. In practice, we find that most often only a few iterations (typically 2–10) are necessary. It is important to note that Chorin [28], in his original fractional step work, makes a similar suggestion. In addition as mentioned in the introduction, the EISPH [8, 10, 12, 14] only uses a single iteration and obtains acceptable results for a wide variety of problems. However, in the present work, we find that a single iteration is insufficient especially in the case of high-Reynolds number problems and we explore this in section 4.2. As such, we recommend choosing the $\epsilon$ to be of $O(h)$ or $O(\Delta t)$.

In Table 1, we show the average number of Jacobi iterations for different problems and tolerances. These results seem to show that using Jacobi iterations is very effective.

Although we do not discuss this in detail, our repository contains an implementation of a BiCGStab [29] iterative algorithm to solve the PPE. A BiCGStab is needed since when we apply the boundary conditions, the resulting matrix may become non-symmetric. We find that implementing the BiCGStab correctly is much more complex as compared to the Jacobi iterations. Further, even though it converges faster, we find that the iterations do not seem to produce significantly better results. We find that in practice, it is much faster to use the Jacobi approach since most of the time we only need a few iterations to converge. This suggests that using more accurate iterative sparse linear solvers may be unnecessary.

We therefore recommend the use of the Jacobi iterations approach for its efficiency and simplicity.

3.2 The Algorithm

For clarity, the Algorithm (1) shows the steps involved in the proposed scheme. This simple formulation works very well in practice as borne out by our numerous benchmarks. The formulation allows us to satisfy boundary conditions that are traditionally not easy to do with a matrix-based formulation.

The time steps are chosen based on the following criteria, the CFL criterion is,

[TABLE]

where $|\mathbf{U}|$ is the magnitude of maximum velocity in the simulation, the viscous criterion is,

[TABLE]

the condition due to external force (if applicable) is,

[TABLE]

time steps are then chosen to be the minimum of above conditions,

[TABLE]

Adaptive time steps can be used in the following way. At each time step the maximum velocity in the simulation is calculated, and used in the CFL condition,

[TABLE]

where $\max{|\mathbf{u}^{n}_{i}|}$ is the maximum magnitude of velocity in the domain at the current time step. Similarly, the time step restriction due to force is calculated as,

[TABLE]

where $\mathbf{f}^{n}_{i}=\mathbf{f}^{n}_{\text{p}}+\mathbf{f}^{n}_{\text{visc}}+\mathbf{f}^{n}_{\text{body}}+\mathbf{f}^{n}_{\text{avisc}}+\mathbf{f}^{n}_{\text{astress}}$. The restriction due to viscosity has no effect as viscosity is held constant in all our simulations. Hence the time step for the next iteration is given by,

[TABLE]

3.3 Implementation details

Here we outline the procedure in solving the iterative formulation using the PySPH framework [22, 23]. PySPH is a Python framework which implements various SPH formulations. The user code is written in Python from which high performance serial (OpenMP) or parallel (OpenCL or CUDA C) code is automatically generated. PySPH also supports multi-CPU execution using MPI.

Listing 1 shows the Python code used to define the inter-particle interactions for a Taylor-Green problem. This is done in two separate stages. Each stage is defined as a list of equations. Each equation has a dest and sources keyword argument. The dest refers to the particles on which the equation is to be solved and sources is the list particles that influence the dest particles. stage1 is a list of equations that are to be solved before the first integration step similarly stage2 are solved before the second integration step. A Group is a PySPH construct which updates all the particles with a given set of equations. A Group is solved using the current particle properties. For example, in the stage1, the summation density is evaluated and updated density is found for all the particles before moving on to the next group of equations. In the Listing 1, we see that the first stage performs the steps 2-4 in Algorithm 1. The second stage performs the remaining computations. Step 5 in the algorithm is performed by the VelocityDivergence equation, next the equations of steps 7 and 8 are performed in an iterated group. The iterated group iterates until the convergence criterion is satisfied or if predefined maximum iterations are completed (in practice this limit is never reached). Finally, the last group performs steps 9 and 10. A suitable integrator updates the velocities and positions. More details on PySPH are available in [22].

An example equation is shown in Listing 2. This shows how the equation for the pressure coefficients, i.e. equations (26) and (27), are written. The initialize method of the class initializes the variables to zero, d_idx is the index of dest particles. The loop method is called for every pairwise interaction with the neighbors, where s_idx is the index of the source (neighbor) particles. The source particle properties are prefixed with s_ and the destination with d_, for example s_pk is an array of the pressure of the source particles at the $k$th iteration, $p_{k}$, and d_diag is the diagonal term, $D_{ii}$ of the destination particle. Various quantities like DWIJ which is the gradient vector for the current particle, and XIJ is the distance between destination particle and source particle are available to each method. The listings demonstrate the ease with which the new scheme can be implemented in the PySPH framework.

3.4 Solid wall boundary conditions

In order to satisfy solid wall boundary conditions, the method proposed by Adami et al. [21] is slightly modified to work with the ISPH scheme. This method uses 3–4 layers of dummy (or ghost) solid particles and then performs a Shepard interpolation of the pressure and also accounts for any fluid acceleration to set the solid wall pressure. As discussed in the original paper the pressure of the dummy wall particles are set using,

[TABLE]

where the subscript $w$ denotes a wall particle and $f$ a fluid particle, $\mathbf{a}_{w}$ denotes the acceleration of the wall and $\mathbf{g}$ the acceleration due to gravity if relevant. In the case of free-surface flows, we ensure that the wall does not have any negative pressures as this often causes particles to stick to and sometimes penetrate the body. Hence if the pressure on the wall is negative, we set it to zero.

In our iterative scheme to solve for the pressure, during each iteration, we perform a Shepard interpolation of the fluid pressure using equation (37) and use those values in evaluating the diagonal and off-diagonal terms. Thus the solid particles are treated just as fluids but with their pressure evaluated based on the extrapolated fluid pressure from the previous iteration. The dummy particle velocities are normally set by first calculating a Shepard-extrapolated fluid velocity on the ghost particles using,

[TABLE]

where $h_{ij}=(h_{i}+h_{j})/2$. Note that we use a smaller kernel radius than the original scheme to ensure that the closest fluid particles have the most influence on the wall velocity. The velocity of the dummy particles is then set by,

[TABLE]

where $\mathbf{u}_{i}$ is the physical velocity of the wall itself. If a no-slip boundary condition is satisfied, the above wall velocity is used for the wall particles when calculating the viscous forces. For the ISPH method, due to the larger allowed timesteps, this basic approach generally leads to some amount of leakage. We modify this boundary condition and ensure that the resulting ghost velocity $\mathbf{u}_{w}$ has no component into the solid itself. This is done by using the normal of the solids. We compute the normal of the solid bodies and if the component of the ghost velocity in the direction of the normal is into the solid, we subtract that component,

[TABLE]

Note that this is not done when the direction of the dummy particle velocity is away from the solid body.

When solving the PPE, the boundary condition on $\mathbf{u}^{*}$ is typically set such that the walls do not have any slip. However, we find that when simulating high-Reynolds number flows this introduces unnecessary divergence on the fluid particles leading to noisy results. This can be mitigated by using a slip velocity for the $\mathbf{u}^{*}$. This is illustrated for the case of the lid-driven cavity problem in section 4.3.

3.4.1 Computation of normals

The normal is itself computed purely using the particle positions in a novel way by first computing the following quantity for the solid particles,

[TABLE]

If the magnitude of the resulting vector is less than $\frac{1}{4h_{i}}$, then the $\mathbf{n}^{*}$ is set to zero otherwise we normalize the vector by its magnitude. Then, we smooth these normals using an SPH approximation,

[TABLE]

These normal vectors are then made into unit normals. This produces smooth normals with the positions of the particle alone.

3.5 Inlet and outlet boundary conditions

It is relatively straightforward to simulate wind-tunnel type problems with our scheme as well. This requires that the inlet and outlet boundary conditions be set appropriately.

For the inlet, we set the prescribed velocity, this is often a constant or a particular velocity profile. The pressure of the inlet region is also solved for, based on the iterative pressure-Poisson equation.

For the outlets we use a modified “do-nothing” type condition as elaborated in Negi et al. [30]. As a fluid particle enters the outlet region, its properties $\mathbf{u}$ and $p$ are frozen. The particle is advected in the outlet region using a Shepard extrapolated velocity from the current fluid particles along the direction of the outlet. For example if the outlet is oriented perpendicular to the $x$-axis, then we take the fluid’s $x$-component of velocity $\mathbf{u}$ and perform a Shepard interpolation of that on the fluid at each timestep and move the outlet particles. This simple approach works very well as seen from the results presented later.

4 Results and discussion

We simulate several standard benchmarks illustrating the proposed method. These benchmarks include internal flows with exact solutions, external flows with a free-surface, two and three dimensional cases and finally the flow past a circular cylinder. We also demonstrate the performance of the new scheme on a GPU. Where possible we compare our results with those produced using the traditional ISPH implementations, and with weakly-compressible SPH using the TVF [19] or EDAC-SPH [31] formulations. Our implementation of the scheme is made using the PySPH framework [22, 23]. In the interest of reproducible research, all the code for producing the results in this manuscript are made available, furthermore every figure is automatically generated using a simple automation framework [24]. The code is available at https://gitlab.com/pypr/sisph.

4.1 Taylor Green Vortex

The Taylor-Green Vortex problem in two-dimensions is periodic and has an exact solution given by,

[TABLE]

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

We simulate this problem by setting the initial condition at $t=0$ for a given $Re$. We choose $Re=100$ for our initial tests and compare the results with the exact solution. We use the quintic spline with $h/\Delta x=1.0$.

The decay rate of the velocity is computed by plotting the maximum velocity $|\mathbf{u}_{\max}|$ at each time. We compute the $L_{1}$ error in the velocity magnitude as,

[TABLE]

where $\mathbf{u}_{i,exact}$ is found at the position of the $i$’th particle. We compute the $L_{1}$ error in the pressure using,

[TABLE]

Here $p_{i,avg}$ is the average pressure due to the particle and its neighbors. This is done to avoid any constant pressure bias in any of the schemes.

Since the Taylor-Green problem has an exact solution, we also take the opportunity to test various parameters in the new scheme. In particular we explore variations in the following,

Comparison with the original matrix-based formulation of ISPH, with and without the use of shifting, vs. GTVF and perturbation.

• 2.

Different tolerance value, $\epsilon$.

• 3.

The use of the symmetric form of the pressure gradient.

• 4.

Two different resolutions with different Reynolds numbers.

We then compare the results with the WCSPH scheme, the IISPH scheme and the EDAC scheme. We add a small amount of perturbation (of at most $\Delta x/5$), in the initial position of the particles for these schemes to ensure uniform particle distribution. But, the mass and density remain unchanged.

4.1.1 Comparison with ISPH

We compare the proposed scheme hereinafter called as SISPH, to the original matrix form of ISPH [4], we use a particle shifting technique [16] to maintain uniform particle distribution, as without particle shifting the ISPH method is unstable for higher Reynolds number due to particle disorder [16]. Here the tolerance is set to, $\epsilon=10^{-2}$. The “asymm” form of pressure gradient is used. For the pressure coefficient matrix in ISPH we construct a sparse matrix, and use biconjugate gradient stabilized method from SciPy [32] sparse matrix library to solve the pressure-Poisson matrix.

In Fig. 1, the decay of the velocity and the error in the pressure are plotted. As can be seen, the decay in velocity of the new scheme is close to the exact solution. The pressure plots also reveal that the new scheme performs very well. Fig. 2 shows the particle plots for ISPH, ISPH with Fickian diffusion based shifting [17, 18], and SISPH schemes. As said before ISPH is unstable without shifting [16]. ISPH with Fickian based shifting and SISPH both result in a uniform distribution of particles. This shows that the new scheme performs better than the original ISPH and is comparable with the case of ISPH along with shifting.

4.1.2 Change in form of pressure gradient.

Here we observe the effect of the two ways of computing the pressure gradient. The tolerance used for iteration is $\epsilon=10^{-2}$. We also note that the “symm” form of pressure gradient works well even without the use of the GTVF formulation, but the “asymm” pressure gradient does not work without the addition of the GTVF. However, the “asymm” pressure gradient with the GTVF produces very good results, as can be seen in Fig. 3. These result are better than the “symm” form for the pressure and decay rates with and without an initial perturbation. Hence for the Taylor-Green problem we use the “asymm” form to compute pressure gradient.

4.1.3 Change in tolerance

We next study the variation in the iteration tolerance for the SISPH scheme. We vary the tolerance in the range $\epsilon=[0.1,0.05,10^{-2},5\times 10^{-3},10^{-3},5\times 10^{-4}]$. The “asymm” pressure gradient form with a $100\times 100$ grid of particles is used. As can be seen in Fig.4, we see that the results are the same even as we vary the tolerance showing that it is not necessary to use very low tolerance when the spatial resolution does not demand it.

4.1.4 Change in Re with varying resolution.

We next compare the Taylor-Green problem with different Reynolds numbers, $Re=100,1000$, with different resolutions of $50\times 50$ and $100\times 100$. The tolerance is, $\epsilon=10^{-2}$. “asymm” pressure gradient is used here. Fig. 5 shows the results. These indicate that the scheme does work well for higher Reynolds numbers.

4.1.5 Comparison with other schemes

Here we compare with other established schemes such as WCSPH, IISPH and EDAC. All the schemes have a perturbed initial condition where the particles are displaced by a maximum of $\Delta x/5$ randomly about their original position with a uniform distribution. All the schemes use the same initial particle distribution. We use a tolerance of $\epsilon=10^{-2}$ for the SISPH scheme, “asymm” form of pressure gradient with a $100\times 100$ grid of initial particles.

In Fig. 6(a) we note that EDAC and SISPH formulation matches the exact decay rate, and in the Fig. 6(b) the errors in pressure, $p_{L_{1}}$, are much better compared to WCSPH and IISPH.

In conclusion SISPH with a tolerance of $\epsilon=10^{-2}$ works well with the Taylor-Green problem and hence we choose this for simulating the subsequent problems.

4.2 Comparison with EISPH

Here we compare the Explicit ISPH (EISPH) [12, 8, 9, 11] with SIPSH. In EISPH, the PPE is solved with an explicit equation which is equivalent to one single iteration of the Jacobi iterations in each time step. In SISPH the PPE is solved iteratively to reach a desired tolerance per every time step. This is important in the case of high Reynolds number simulations. We compare the results using both the methods for the Taylor-Green problem at $Re=10000$. We perturb the initial particle positions by a maximum of $\Delta x/10$ without changing the mass or density.

Figure 7 shows the particle plots of EISPH vs SISPH where the color indicates velocity. This shows that the using EISPH approach leads to an inaccurate velocity. This shows the importance of multiple iteration per time step for high Reynolds number flows.

4.2.1 Performance of SISPH.

Here we compare the performance of SISPH scheme with ISPH scheme on a CPU for a single-core and show the performance of SISPH on multi-cores. Fig. 8 shows the time taken versus the number of particles, $N$, for the different cases. When $N$ is large the SISPH method is an order of magnitude faster than the matrix based ISPH method. Fig. 9 compares the scale up of SISPH with ISPH scheme on single core, the single core performance of SISPH is increase to 4 times, and the multi-core (with 4 cores) is nearly 12 times as fast as single core ISPH for large $N$. The matrix ISPH cannot be run on GPUs without a significant amount of programming effort whereas SISPH can be executed on the GPUs very easily as demonstrated for a 3D problem later.

4.3 Lid Driven Cavity

The next problem considered is the classical lid-driven cavity problem. This problem allows us to test the solid wall boundary condition and the ability of the new scheme to handle the solid wall boundaries.

The problem involves a unit square vessel filled with incompressible fluid. The upper boundary over the fluid is made to move rightward with a unit speed, $U$. The other walls are treated as no-slip walls. For this problem the Reynolds number is defined as $Re=UL/\nu$, where $L$ is the length of a side. We simulate the problem at $Re=100$ with two grid sizes, $50\times 50$ and $100\times 100$. The results are compared with those of Ghia et al. [33] where the velocity profile along a line passing horizontally and vertically through the center is shown after the simulation has stabilized. We also compare the results with those from the TVF [19] scheme. We use the parameter $h/dx=1.0$ and employ the quintic spline kernel, the tolerance for SISPH is, $\epsilon=10^{-2}$, and since this is an internal flow the background pressure in GTVF formulation is fixed to a constant reference pressure which is $2\max(p_{i})$ after the first iteration. The “symm” form of pressure gradient is used.

Fig. 10, plots the distribution of $u$ and $v$ along the centerline along the horizontal and vertical. The results are in very good agreement. In this simulation it is observed that, either use of “symm” form of pressure gradient or the “asymm” form doesn’t affect the solution significantly, but the “symm” form works well even when no GTVF is used whereas the “asymm” form fails to work without the GTVF.

Due to the larger time steps the ISPH method allows, the new scheme is about 5 times faster than the TVF simulation. The $Re=100$ case with $100\times 100$ fluid particles takes only around 87 seconds to execute for 10 seconds of simulation time (on a quad core Intel i5-7400) with the new scheme. With the TVF scheme the same problem takes about 433 seconds. This shows that the new method works well and is much more efficient.

We next compare the results for $Re=10000$. The velocity profiles show a good agreement with Ghia et al. [33]. It should however be noted that the simulation of Barcarolo [10] shows voiding with $Re=3200$. In our case the GTVF and an accurate pressure leads to very good results.

Fig. 12, plots the velocity magnitude of the cavity problem at $Re=10000$. It highlights the difference between the boundary condition that is applied to the intermediate velocity ($\mathbf{u}^{*}$) when solving the PPE. In the plot on the left, the boundary condition on $\mathbf{u}^{*}$ is set to ensure a no-slip wall, whereas, on the right, the boundary condition is set to ensure a slip wall. The slip wall boundary condition on the intermediate velocity reduces the noise introduced in the velocity field due to the unnecessary divergence created at the wall.

4.4 Evolution of a Square Patch

The simulation of an initially square patch is a free-surface benchmark introduced by Colagrossi [34]. Here, an initially square patch of fluid having length $L$ is given the following initial conditions,

[TABLE]

[TABLE]

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

This problem is simulated for 3 seconds using the new scheme, and the WCSPH scheme, for the discussion on the results of IISPH scheme see Ramachandran et al. [35]. The quintic spline kernel with $h/\Delta x=1.3$ is used for all schemes. An artificial viscosity coefficient of $\alpha=0.15$ is used for the WCSPH and the new scheme. The tolerance chosen for the new scheme is $\epsilon=10^{-2}$. Initial particle distribution of $100\times 100$ is used for all the schemes.

A noticeable difference is observed when using the different forms of pressure gradient, the “asymm” gradient of pressure simulates the problem significantly better than the other, where even though the core structure remains, the particles along the free surface are deviated further. We also observed that computing the density after each time step using summation density is very crucial for this problem, without it i.e., using a constant density through out makes the simulation unstable. A comparison with WCSPH scheme is made, and the effect of pressure gradient is shown in the Fig. 13. It is clear that the new scheme with the “asymm” form of pressure gradient performs very well.

4.5 Dam break in two-dimensions

The dam-break problem in two-dimensions is a classic problem involving free-surfaces and solid walls. We simulate the standard problem as described in [36]. This involves a vessel of width 4m with a block of fluid of width 1m and height 2m on the left side. The block of water is released from rest and falls under the influence of gravity ($g=9.81m/s$).

The problem is simulated with the proposed scheme using a quintic spline kernel. Artificial viscosity of $\alpha=0.05$ is used for all the schemes. The fluid is not treated as viscous. The $h/\Delta x=1.3$ used for all the schemes, the default particle spacing, $\Delta x=0.01$, is used which results in 37k particles. The time step is fixed and chosen as $\Delta t=0.125h/\sqrt{2gh_{w}}$, where $g$ is the acceleration due to gravity and $h_{w}$ is the initial height of the water which is 2m. The tolerance value for the new scheme is, $\epsilon=10^{-2}$. The “symm” form of pressure gradient is used, although no significant difference is observed with either of the forms.

The particle distribution of the flow at different times is plotted in Fig. 14. The results show that the scheme works very well. Fig. 15 plots the toe of the breaking water versus time. The results are compared with the WCSPH scheme, EDAC scheme and the numerical simulations in Koshizuka and Oka [37]. As can be seen, the results are in very good agreement.

4.6 Dam break in three-dimensions

A three-dimensional dam break is considered in order to show that the scheme works in three dimensions and demonstrates the performance of the scheme on a GPU. The problem considered is a small modification of the problem described in [36], wherein an obstacle is placed in the path of the flow. We do not use an obstacle in our example. A rectangular block of fluid having height $h_{w}=0.55m$, width of $1m$ and length of $1.228m$ is at rest at one end of a container. The container is a long rectangular container of length $3.22m$ and open at the top. Four layers of boundary particles are used to enforce the no-penetration boundary condition. A quintic spline kernel is used with $h/\Delta x=1.0$. The reference velocity for the flow is $v_{ref}=\sqrt{2gh_{w}}$, where $g=9.81m/s^{2}$, and the tolerance is, $\epsilon=0.01$. The particle spacing is chosen as $\Delta x=0.01m$, this results in around 664k fluid particles and around 627k solid wall particles. The time step is chosen to be $h/(8v_{ref})$. An artificial viscosity is used with the parameter $\alpha=0.25$. The results at different times are shown in Fig. 16. There is a large amount of splashing initially but the scheme produces good results and are similar to those of Chow et al. [7].

We execute this on a GPU using OpenCL. The code is executed on an NVIDIA 1070 Ti processor with single precision. It is also executed with double precision on a quad core Intel i5–7400 CPU running at 3.0 Ghz. The user-code is identical as PySPH [22, 23] automatically executes the code using OpenCL or OpenMP depending on the options passed. On the GPU, the execution for 3 seconds of simulation takes 4531 seconds and the same simulation on the CPU with OpenMP takes 40264 seconds suggesting that the GPU execution is around 8.9 times faster than the CPU (quad core) for this size of problem.

The GPU implementation is not the most highly optimized but the above clearly shows that our new algorithm can be executed effectively on a GPU without requiring any additional effort as normally needed for ISPH schemes [7].

4.6.1 Performance of SISPH

Here we show the performance of SISPH scheme on single-core, and multi-core CPU (Intel i5-7400) with double precision, and NVIDIA 1050Ti, 1070Ti GPUs with single precision. In Fig. 17 the increase in number of particles, $N$, vs time taken by SISPH on various platforms is show, CPUs outperform GPUs when $N$ is small as can be expected. For sufficiently large $N$, at around 100k particles, we can start to see the improved performance of the GPUs (albeit with single precision), and for 1M particles the time is reduced by an order of magnitude, also the “scale up” is linear. Fig. 18 shows speed up of SISPH on multi-core and GPUs compared to that of single core CPU, while the performance on multi-core (with 4 cores) expectedly peaks around 4 times that of single core, performance of GPUs increases with number of particles and reaches the peak of around 16 times that of single CPU core for 1050Ti and around 32 times for 1070Ti, showing the increase in performance with the proposed scheme.

4.7 Flow past a circular cylinder

As a final example, the flow past a circular cylinder is considered. A cylinder of diameter $D=2m$ is considered, the problem setup is as shown in Fig. 19. Note that the particles discretizing the cylinder are placed at a distance of $\Delta x/2$ from the outer circumference in order to effectively simulate the correct location of the boundary of the cylinder. The cylinder is placed a distance of $5D$ from the inlet and the outlet is placed $10D$ from the center of the cylinder. The inlet is set to a uniform velocity, $U$, along the x-axis of 1m/s. The sides of the wind-tunnel are set to slip walls. These walls are placed at a distance of $7.5D$ from the center of the cylinder. The viscosity is set to ensure a Reynolds number, $Re=UD/\nu=200$. No artificial viscosity is used. The flow is started impulsively and simulated for $200$ seconds. A quintic spline kernel is used with $h/\Delta x=1.2$. Boundary conditions are implemented as discussed in Sections 3.4 and 3.5. The convergence tolerance is set as $\epsilon=0.01$. The particle spacing is chosen as $\Delta x=D/30$, which results in 200k fluid particles. The “symm” form is used for the calculation of pressure gradient. Since the resolution is low, the geometry curvature is captured by placing particles along the circumference of the cylinder such that the volume occupied by the particles is approximately constant as done in [30]. For this case we encountered particle voiding when we use the TVF reference pressure of $2\max(p_{i})$. This was resolved by choosing $p_{\text{ref}}=\rho c^{2}$, where $c$ is $10|\mathbf{U}|_{\text{max}}$.

In Fig. 16, we plot the velocity and pressure contour by taking the snapshot at times 10s, 153s and 185s. The pressure and velocity contour shows an excellent match with the result presented in [38, 39]. The pressure and velocity variations have much less noise than the WCSPH schemes even with a low resolution of particles as compared to [30].

The drag and lift coefficients, $c_{d},c_{l}$ are plotted as a function of time in Fig. 21. It must be noted that unlike WCSPH schemes, the data presented has not gone through any filtering to remove high frequency oscillations. In order to calculate the forces on the cylinder, we compute acceleration on the solid due to the pressure and velocity gradients using the following equation,

[TABLE]

where $\mathbf{f}_{\text{solid}}$ is the force on the cylinder particles, using SPH discretization [19, 31], the above equation is written as,

[TABLE]

where,

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{u}_{w\,i}$ is the velocity of the solid wall particles as obtained from equation (40), $\eta_{i}=\rho_{i}\nu_{i}$.

We obtain an average coefficient of drag, $c_{d}=1.609$ (once the vortex shedding has been established) and the maximum coefficient of lift, $c_{l}=0.804$. In the Fig. 21, $c_{d}$ and $c_{l}$ show clear oscillations without any change in the amplitude, which is consistent with the results of others. With this we have demonstrated the suitability and robustness of the proposed scheme for a variety of problems.

5 Conclusions and Future work

In this paper we introduce a simple, iterative, incompressible SPH scheme that is based on the original projection formulation of Cummins and Rudman [4]. The method is matrix-free, fast, and suitable for execution on GPUs. The formulation is simple and easy to implement. We introduce a novel technique to ensure a homogeneous distribution of particles. In addition we introduce a modified solid wall boundary condition and show how we can implement simple inlet and outlet boundary conditions. We demonstrate the accuracy and efficiency of the new scheme with a suite of benchmark problems in two and three dimensions involving internal and external flows.