Divergence-Free SPH (DFSPH)
|
|||||||||||||||||||||||||||||||||
DFSPH algorithm:This example shows the DFSPH method introduced by Bender and Koschier [BK17]:
In each simulation step DFSPH first solves the pressure Poisson equation (PPE) for the constant density source term: $$\Delta t \nabla^2 p = \frac{\rho_0 - \rho^*}{\Delta t},$$ where $\rho^*$ is the predicted density after applying all non-pressure forces. Discretizing the PPE using the SPH formulation yields a linear system $$\mathbf{A} \mathbf{p} = \mathbf{s}$$ for the unknown pressure values $\mathbf{p}$ and the source term $\mathbf{s}$. The first solve enforces that the density stays constant over time. After the time integration with the resulting forces the divergence of the velocity field should be zero. However, due to the numerical integration this is not guaranteed. Therefore, DFSPH solves a second PPE for the divergence source term: $$\Delta t \nabla^2 p = \rho \nabla \cdot \mathbf v.$$ Finally, by the combination of both solvers we get a constant density and a divergence-free velocity field. 1. Neighborhood searchThe 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 computationUsing 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. Compute diagonal matrix elements $a_{ii}$The diagonal elements of the system matrix $\mathbf{A}$ are determined as $$a_{ii} = -\frac{\Delta t}{\rho_i^2} \left (\sum_j \left \| m_j \nabla W_{ij}\right \|^2 + \left \| \sum_j m_j \nabla W_{ij}\right \|^2 \right ).$$4. Viscosity (XSPH)Artificial viscosity is introduced by smoothing the velocity field as: $$\mathbf {a}^{\text{visco}}_i = \frac{\alpha}{\Delta t} \sum_j \frac{m_j}{\rho_j} (\mathbf v_j - \mathbf v_i) W_{ij}.$$5. Integrate velocities only considering non-pressure forcesThe particles velocities are updated using an Euler method: $$\begin{align*} \mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a^\text{np}(t) \\ \end{align*}$$ where $\mathbf a^\text{np}$ are the accelerations corresponding to the sum of non-pressure forces (e.g. viscosity).6. Compute constant density source termDetermine the source term on the right hand side of the linear system: $$s_i = \frac{\rho_0 - \rho_i^*}{\Delta t},$$ where $$\rho_i^* = \rho_i + \Delta t \frac{D \rho_i}{D t} = \rho_i + \Delta t \sum_j m_j (\mathbf{v}^*_i - \mathbf{v}^*_j) \cdot \nabla W_{ij}$$7. Solver loopThe PPE is solved using relaxed Jacobi iterations with a relaxation factor of $\omega = 0.5$ until the density error is smaller than a given tolerance or a maximum number of iterations is reached.8. Compute pressure accelerationsThe pressure acceleration of a particle is determined by a symmetric SPH formulation to preserves linear and angular momentum: $$\mathbf a^\text{p}_i = - \sum_j m_j \left ( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right ) \nabla W_{ij}.$$9. Determine $\mathbf{A} \mathbf{p}$Compute the product of the system matrix $\mathbf{A}$ and the current pressure vector: $$(\mathbf{A} \mathbf{p})_i = \Delta t \sum_j m_j \left (\mathbf{a}^\text{p}_i - \mathbf{a}^\text{p}_j \right ) \cdot \nabla W_{ij}.$$10. Update pressure valuesIn each solver iteration we update the pressure value for each particle $i$ as $$p_i := p_i + \frac{\omega}{a_{ii}} \left (s_i - (\mathbf{A} \mathbf{p})_i \right ).$$11. Time integration with pressure forcesFinally, the particles are advected using a symplectic Euler method: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a^\text{p}(t) \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t), \end{align*}$$ where $\mathbf a^\text{p}$ are the accelerations corresponding to the sum of pressure forces.16. Compute divergence source termDetermine the source term on the right hand side of the linear system: $$s_i = \rho_i \nabla \cdot \mathbf v_i,$$ where $$\rho_i \nabla \cdot \mathbf v_i = -\frac{D \rho_i}{D t} = -\sum_j m_j (\mathbf{v}_i - \mathbf{v}_j) \cdot \nabla W_{ij}.$$ Note that in contrast to the constant density source term, we do not need predicted velocities here since we want to make the current velocity field divergence-free.Boundary handlingIn 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
|