Focused on the Best Accuracy ISPH

Hello! This is E8IGHT!

Today, I want to introduce the ISPH, developed to supplement the weak portion of SPH by adapting the formula which is normally used in mesh based CFD.

1. ISPH

The original formulation of SPH, also known as WCSPH (Weakly Compressible SPH), uses an equation of state to link density with pressure. During this calculation, density is permitted to vary within about 1% of the reference value. But this leads to fluctuations in pressure and velocity field. To solve this problem, Chorin’s projection method (1), which conventional mesh based CFD is originally using, is adapted to SPH. This is called Incompressible SPH(ISPH). ISPH is highly accurate and enforce the incompressibility, but it has two major problems, slow computational time due to large matrix system and limited solving scale due to the lack of memory.

Fig 1. Wave impact on a circular cylinder

[Source = Skillen, A., et al. (2013). "Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised Fickian smoothing applied to body–water slam and efficient wave–body interaction." Computer Methods in Applied Mechanics and Engineering 265: 163-173.]

1.1 Helmholtz – Hodge decomposition

Helmholtz-Hodge decomposition states that any vector field on a simply connected domain can be resolved into the sum of an irrotational (curl-free) vector field and a solenoidal (divergence-free) vector field.

Fig 2. Closed surface and surfaces with boundary (left): closed surface (right): surfaces with boundary

For arbitrary function u in simply connected domain, this could resolved into solenoidal part and irrotational part.

Taking the divergence of equation yields

If the vector fireld u is known, the above equation can be solved for the scalar function ϕ and the divergence-free part of u can be extracted using the relation

1.2 Chorin’s projection method

The incompressible Navier-Stokes equation can be written as

From the above equation, one first computes an intermediate velocity u^(*) explicitly using the momentum equation by ignoring the pressure gradient term.

Where u^n is the velocity at n time step. In the second half of the algorithm, the projection step, we correct the intermediate velocity to obtain the final solution of the time step u^(n+1).

Computing the right-hand side of the second half step requires knowledge of the pressure at the n+1 time step. This is obtained by taking the divergence and requiring that ∇∙u^(n+1)=0 which is the divergence (continuity) condition, thereby deriving the following Poisson equation.

To satisfy the Helmholtz-Hodge decomposition’s condition, Neumann boundary condition should be adapted to boundary of fluid domain.

1.3 SPH projection method

To adapt Chorin’s projection method into SPH, viscous force and Laplacian of pressure can be calculated using following equations.

In the SPH, because of the symmetry of the kernel for pairwise particle interactions, the matrix is initially symmetric without boundary conditions. The inclusion of Dirichlet boundary condition, as well as Neumann boundary condition (eq. 8), creates a non-symmetric matrix which can be solved iteratively.

1.4 CSR format (Compressed Sparse Row)

The matrix which stores the pressure Poisson equation is sparse matrix since the interaction only calculates the neighborhood particles. To store this matrix, CSR format is used due to the sparsity of pressure Poisson equation matrix. This can reduce the memory requirements by storing only non-zero entries. CSR format comprises three arrays. One contains the non-zero entries of the matrix, another contains the corresponding column index of each non-zero entry, the other contains pointers to the memory location of entries that start a new row in the matrix.

Fig 3. Sparse matrix from ISPH pressure Poisson equation

[Source = Chow, A. D. (2016). Converting DualSPHysics to solve strictly incompressible SPH. 3rd DualSPHysics User Workshop.]

2. ISPH in NFLOW

As ISPH is a large-scale complicated calculation method, it is slow and causes memory matters to save the pressure Poisson formula of all steps of calculations. However, the GPU technology adopted NFLOW can solve the speed problem. Also, the strong rendering functions and simplified pre-processes of CFD helps users in many aspects.

So far, we learned the basic of ISPH.

It is amazing that the speed problem of ISPH can be solved when it is applied to NFLOW!

We will come back with more technical information next time.

Thank you.

References

(1) Chorin, A. J. (1967). "A numerical method for solving incompressible viscous flow problems." Journal of Computational Physics 2(1): 12-26.

(2) Skillen, A., et al. (2013). "Incompressible smoothed particle hydrodynamics (SPH) with

reduced temporal noise and generalised Fickian smoothing applied to body–water slam

and efficient wave–body interaction." Computer Methods in Applied Mechanics and

Engineering 265: 163-173.

(3) Cummins, S. J. and M. Rudman (1999). "An SPH Projection Method." Journal of

Computational Physics 152(2): 584-607.

(4) Chow, A. D., et al. (2018). "Incompressible SPH (ISPH) with fast Poisson solver on a

GPU." Computer Physics Communications 226: 81-103.

New Flow Leads Our World!

www.e8-global.com

Tel. 1811-2707

Email. info@e8-korea.com