WCSPH algorithm:
This example shows the WCSPH method introduced by Becker and Teschner [BT07].
- perform neighborhood search
- compute particle densities $\rho_i$
- compute particle pressure values $p_i$ using an equation of state
- compute pressure forces
- compute viscosity (XSPH)
- time integration
1. Neighborhood search
The neighborhood search is performed using spatial hashing. In each step all particles are added to a spatial grid using a cell size that equals the support radius of the SPH kernel function. Hence, the neighbor particles of a particle in cell (i,j) lie in one of the 9 neighboring cells: (i±1, j±1).
2. SPH density computation
Using the SPH formulation the density is determined as:
$$\rho_i = \sum_j \frac{m_j}{\rho_j} \rho_j W_{ij} = \sum_j m_j W_{ij}.$$
Hence, the density only depends on the masses and positions of the particles.
3. Equation of State
Pressure is determined locally by the density:
$$p_i = \frac{\kappa \rho_0}{\gamma} \left ( \left ( \frac{\rho_i}{\rho_0} \right )^\gamma - 1 \right ),$$
where he constants $\kappa, \gamma$ define the stiffness.
4. Pressure force
The pressure force of a particle is determined by a symmetric SPH formulation to preserves linear and angular momentum:
$$\mathbf f_i = - m_i \sum_j m_j \left ( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right ) \nabla W_{ij}.$$
5. Viscosity (XSPH)
Artificial viscosity is introduced by smoothing the velocity field as:
$$\hat {\mathbf v}_i = \mathbf v_i + \alpha \sum_j \frac{m_j}{\rho_j} (\mathbf v_j - \mathbf v_i) W_{ij}.$$
6. Time integration
Finally, the particles are advected using a symplectic Euler method:
$$\begin{align*}
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t) \\
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t),
\end{align*}$$
where $\mathbf a$ are the accelerations corresponding to the sum of non-pressure and pressure forces.
Boundary handling
In order to implement a rigid-fluid coupling the SPH equation for the density computation is extended by a second sum over the contributing neighboring boundary particles $k$ [AIA+12]:
$$\rho_i = \sum_j m_j W_{ij} + \sum_k \Psi_k W_{ik}.$$
Each boundary particle $k$ has a pseudo-mass $\Psi$ which considers the density of the boundary sampling points:
$$\Psi_k = \rho_0 V_k = \rho_0 \frac{m_k}{\rho_k} = \rho_0 \frac{m_k}{\sum_l m_k W_{kl}} = \rho_0 \frac{1}{\sum_l W_{kl}}$$
where $l$ denotes the boundary particle neighbors of particle $k$.
References
- [BT07] Markus Becker and Matthias Teschner. Weakly compressible SPH for free surface flows. In Proceedings of ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007. Eurographics Association.
- [AIA+12] Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler, and Matthias Teschner, "Versatile rigid-fluid coupling for incompressible SPH", ACM Transactions on Graphics 31(4), 2012
|