Position Based Fluids (PBF)
|
|||||||||||||||||||||||||||
PBF algorithm:This example shows the PBF method introduced by Macklin and Müller [MM13]:
1. Time integrationThe particle positions are advected using a symplectic Euler method to obtain predicted positions $\mathbf x^*$: $$\begin{align*} \mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a_\text{ext}(t) \\ \mathbf x^* &= \mathbf x(t) + \Delta t \mathbf v^*, \end{align*}$$ where $\mathbf a_\text{ext}$ are accelerations due to external forces.2. 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).3. Solver loopPosition corrections are computed and applied in a loop until the density constraint is satisfied (within a given tolerance) or a maximum number of iterations is reached.4. Compute Lagrange multipliersThe general form to compute the Lagrange multipliers in position-based dynamics is $$\mathbf J \mathbf M^{-1} \mathbf J^T \boldsymbol \lambda = -\mathbf C(\mathbf x^*),$$ where $\mathbf J = (\partial C / \partial \mathbf x)^T$ is the Jacobian of the constraint function and $\mathbf M$ denotes the mass matrix. Since the density constraint is a scalar function the Lagrange multiplier is determined as $$\lambda_i = - \frac{C_i}{\sum_j \frac{1}{m_j} \| \partial C_i / \partial \mathbf x_j\|^2 + \varepsilon},$$ where $\varepsilon$ is a small constant which is used to avoid a division by zero. To compute the Lagrange multipliers the SPH discretization of the density constraint is required: SPH discretization of constraint:$$C_i(\mathbf x) = \frac{\rho_i}{\rho_0} - 1 = \frac{1}{\rho_0} \sum_j m_j W_{ij} - 1$$SPH constraint gradients:$$\begin{eqnarray*} \frac{\partial C_i}{\partial \mathbf x_i} &= \frac{1}{\rho_0} \sum_j m_j \nabla W_{ij} \\ \frac{\partial C_i}{\partial \mathbf x_j} &= - \frac{1}{\rho_0} m_j \nabla W_{ij} \\ \end{eqnarray*}$$5. Position correction:A symmetric position correction is determined using the Lagrange multipliers of particle $i$ and its neighbors $j$: $$\Delta \mathbf x_i = \frac{1}{m_i \rho_0} \sum_j m_j (\lambda_i + \lambda_j) \nabla W_{ij}, $$ where the Lagrange multipliers are computed for the predicted positions $\mathbf x^*$. Then the correction is applied to update $\mathbf x^*$ as $$\mathbf x_i^* := \mathbf x_i^* + \Delta \mathbf x_i.$$6. Velocity update:The final positions and velocities are computed as: $$\begin{align*} \mathbf x_i(t+\Delta t) &= \mathbf x_i^* \\ \mathbf v_i(t+\Delta t) &= \frac{1}{\Delta t} (\mathbf x_i(t+\Delta t) - \mathbf x_i(t)) \end{align*} $$7. 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}.$$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
|