|
| 1 | +# [Implicit Incompressible SPH](@id iisph) |
| 2 | + |
| 3 | +Implicit Incompressible SPH (IISPH) as introduced by [Ihmsen et al. (2013)](@cite Ihmsen2013) |
| 4 | +is a method that achieves incompressibility by solving the pressure Poisson equation. |
| 5 | +The resulting linear system is iteratively solved with the relaxed Jacobi method. |
| 6 | +Unlike the [weakly compressible SPH method](@ref wcsph), incompressible methods determine |
| 7 | +pressure by enforcing the incompressibility constraint rather than using an equation of |
| 8 | +state. |
| 9 | + |
| 10 | +```@autodocs |
| 11 | +Modules = [TrixiParticles] |
| 12 | +Pages = [joinpath("schemes", "fluid", "implicit_incompressible_sph", "system.jl")] |
| 13 | +``` |
| 14 | +## Derivation |
| 15 | +To derive the linear system of the pressure Poisson equation, we start by discretizing the |
| 16 | +continuity equation |
| 17 | +```math |
| 18 | +\frac{D\rho}{Dt} = - \rho \nabla \cdot \bm{v}. |
| 19 | +``` |
| 20 | + |
| 21 | +For a particle ``i``, discretizing the left-hand side of the equation with a forward |
| 22 | +difference yields |
| 23 | + |
| 24 | +```math |
| 25 | +\frac{\rho_i(t+ \Delta t) - \rho_i(t)}{\Delta t}. |
| 26 | +``` |
| 27 | + |
| 28 | +The divergence in the right-hand side is discretized with the SPH discretization for |
| 29 | +particle ``i`` as |
| 30 | +```math |
| 31 | +-\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \nabla W_{ij}, |
| 32 | +``` |
| 33 | +where ``\bm{v}_{ij} = \bm{v}_i - \bm{v}_j``. |
| 34 | + |
| 35 | +Together, the following discretized version of the continuity equation for a particle ``i`` |
| 36 | +is achieved: |
| 37 | +```math |
| 38 | +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij}. |
| 39 | +``` |
| 40 | + |
| 41 | +The right-hand side contains the unknown velocities ``\bm{v}_{ij}(t + \Delta t)`` at the |
| 42 | +next time step, making it an implicit formulation. |
| 43 | + |
| 44 | +Using the semi-implicit Euler method, we can obtain the velocity in the next time step as |
| 45 | +```math |
| 46 | +\bm{v}_i(t + \Delta t) = \bm{v}_i(t) + \Delta t \frac{\bm{F}_i^\text{adv}(t) + \bm{F}_i^p(t)}{m_i}, |
| 47 | +``` |
| 48 | + |
| 49 | +where ``\bm{F}_i^{\text{adv}}`` denotes all non-pressure forces such as gravity, viscosity, surface |
| 50 | +tension and more, while ``\bm{F}_i^p``denotes the unknown pressure forces, which we |
| 51 | +want to solve for. |
| 52 | + |
| 53 | +Note that the IISPH is an incompressible method, which means that the density of the |
| 54 | +fluid remain constant over time. By assuming a fixed reference density ``\rho_0`` for all |
| 55 | +fluid particle over the whole time of the simulation, the density value at the next time |
| 56 | +step ``\rho_i(t + \Delta t)`` also has to be this rest density. So ``\rho_0`` can be plugged in |
| 57 | +for ``\rho_i(t + \Delta t)`` in the formula above. |
| 58 | + |
| 59 | +The goal is to compute the pressure values required to obtain the pressure acceleration that |
| 60 | +ensures each particle reaches the rest density in the next time step. At the moment these |
| 61 | +pressure values are unknown in ``t``, but all the non-pressure forces are already known. |
| 62 | +Therefore a predicted velocity can be calculated, which depends only on the non-pressure |
| 63 | +forces ``\bm{F}_i^{\text{adv}}``: |
| 64 | + |
| 65 | +```math |
| 66 | +\bm{v}_i^{\text{adv}}(t+\Delta t)= \bm{v}_i(t) + \Delta t \frac{\bm{F}_i^{\text{adv}}(t)}{m_i}, |
| 67 | +``` |
| 68 | +Using this predicted velocity and the continuity equation, a predicted density can be defined |
| 69 | +in a similar way as |
| 70 | + |
| 71 | +```math |
| 72 | +\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t). |
| 73 | +``` |
| 74 | + |
| 75 | +To achieve the rest density, the unknown pressure forces must counteract the compression |
| 76 | +caused by the non-pressure forces. In other words, they must resolve the deviation between |
| 77 | +the predicted density and the reference density. |
| 78 | + |
| 79 | +Therefore, the following equation needs to be fulfilled: |
| 80 | +```math |
| 81 | +\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. |
| 82 | +``` |
| 83 | + |
| 84 | +This expression is derived by substituting the reference density ``\rho_0`` for |
| 85 | +``\rho_i(t+\Delta t)`` in the discretized continuity equation and inserting the definitions of |
| 86 | +the predicted velocity ``\bm{v}_{ij}(t+\Delta t)`` and predicted density ``\rho_i^{\text{adv}}``. |
| 87 | + |
| 88 | +The pressure force is defined as: |
| 89 | + |
| 90 | +```math |
| 91 | +\bm{F}_i^p(t) = -\frac{m_i}{\rho_i} \nabla p_i = -m_i \sum_j m_j \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_j(t)}{\rho_j^2(t)} \right) \nabla W_{ij}. |
| 92 | +``` |
| 93 | + |
| 94 | +Substituting this definition into the equation yields a linear system ``\bm{A}(t) \bm{p}(t) = \bm{b}(t)`` |
| 95 | +with one equation and one unknown pressure value for each particle. |
| 96 | +For ``n`` particles, this is a linear system with ``n`` equations and ``n`` unknowns: |
| 97 | + |
| 98 | +```math |
| 99 | +\sum_j a_{ij} p_j = b_i = \rho_0 - \rho_i^{\text{adv}}. |
| 100 | +``` |
| 101 | + |
| 102 | +To solve for the pressure values, a relaxed Jacobi scheme is used. |
| 103 | +This is an iterative numerical method to solve a linear system ``Ap=b``. In each iteration, the |
| 104 | +new values of ``p`` are obtained as computed as |
| 105 | + |
| 106 | +```math |
| 107 | +p_i^{(k+1)} = (1-\omega) p_i^{(k)} + \omega \left( \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} p_j^{(k)} \right)\right). |
| 108 | +``` |
| 109 | + |
| 110 | +Substituting the right-hand side, we obtain |
| 111 | + |
| 112 | +```math |
| 113 | +p_i^{l+1} = (1-\omega) p_i^l + \omega \frac{\rho_0 - \rho_i^{\text{adv}} - \sum_{j \neq i} a_{ij}p_j^l}{a_{ii}}. |
| 114 | +``` |
| 115 | + |
| 116 | +Therefore the diagonal elements ``a_{ii}`` and the sum ``\sum_{j \neq i} a_{ij}p_j^l`` need to |
| 117 | +be determined. |
| 118 | +This can be done efficiently by separating the pressure force acceleration into two |
| 119 | +components: One that describes the influence of particle ``i``'s own pressure on its |
| 120 | +displacement, and one that describes the displacement due to the pressure values of the |
| 121 | +neighboring particles ``j``. |
| 122 | + |
| 123 | +The pressure acceleration is given by: |
| 124 | +```math |
| 125 | +\Delta t^2 \frac{\bm{F}_i^p}{m_i} = -\Delta t^2 \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right)\nabla W_{ij} = \underbrace{\left(- \Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij} \right) }_{d_{ii}} p_i + \sum_j \underbrace{-\Delta t^2 \frac{m_j}{\rho_j^2} \nabla W_{ij}}_{d_{ij}}p_j. |
| 126 | +``` |
| 127 | + |
| 128 | +The ``d_{ii}p_i`` value describes the displacement of particle ``i`` because of the particle ``i`` |
| 129 | +and ``d_{ij}p_j`` describes the influence from the neighboring particles ``j``. |
| 130 | +Using this new values the linear system can be rewritten as |
| 131 | + |
| 132 | +```math |
| 133 | +\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \nabla W_{ij}, |
| 134 | +``` |
| 135 | + |
| 136 | +where the first sum over ``k`` loops over all neighbor particles of ``i`` and |
| 137 | +the second sum over ``k`` loops over the neighbor particles of ``j`` |
| 138 | +(which is a neighbor of ``i``). |
| 139 | +So the sum over the neighboring pressure values ``p_j`` also includes the pressure values |
| 140 | +``p_i``, since ``i`` is a neighbor of ``j``. |
| 141 | +To separate this sum, it can be written as |
| 142 | + |
| 143 | +```math |
| 144 | +\sum_k d_{jk} p_k = \sum_{k \neq i} d_{jk} p_k + d_{ji} p_i. |
| 145 | +``` |
| 146 | + |
| 147 | +With this separation, the equation for the linear system can again be rewritten as |
| 148 | + |
| 149 | +```math |
| 150 | +\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \nabla W_{ij}. |
| 151 | +``` |
| 152 | + |
| 153 | +In this formulation all coefficients that are getting multiplied with the pressure value ``p_i`` |
| 154 | +are separated from the other. The diagonal elements ``a_{ii}`` can therefore be defined as: |
| 155 | + |
| 156 | +```math |
| 157 | +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij}. |
| 158 | +``` |
| 159 | + |
| 160 | +The remaining part of the equation represents the influence of the other pressure values ``p_j``. |
| 161 | +Hence, the final relaxed Jacobi iteration takes the form: |
| 162 | + |
| 163 | +```math |
| 164 | +p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right). |
| 165 | +``` |
| 166 | + |
| 167 | +Because interactions are local, limited to particles within the kernel support defined by |
| 168 | +the smoothing length, each particle's pressure ``p_i`` depends only on its own value and |
| 169 | +nearby particles. |
| 170 | +Consequently, the matrix ``A`` is sparse (most entries are zero). |
| 171 | + |
| 172 | +The diagonal elements ``a_{ii}`` get computed and stored at the beginning of the simulation step |
| 173 | +and remain unchanged throughout the relaxed Jacobi iterations. The same applies for the |
| 174 | +``d_{ii}``. The coefficients ``d_{ij}`` are computed, but not stored, as part of the calculation |
| 175 | +of ``a_{ii}``. |
| 176 | +During the Jacobi iterations, two loops over all particles are performed: |
| 177 | +The first updates the values of ``\sum_j d_{ij}`` and the second computes the updated |
| 178 | +pressure values ``p_i^{l+1}``. |
| 179 | + |
| 180 | +The final pressure update follows the equation of the relaxed Jacobi scheme as stated above, with two |
| 181 | +exceptions: |
| 182 | +#### 1. Pressure clamping |
| 183 | +To avoid negative pressure values, one can enable pressure clamping. In this case, the |
| 184 | +pressure update is given by |
| 185 | + |
| 186 | +```math |
| 187 | +p_i^{l+1} = \max \left(0, (1-\omega) p_i^l + \omega \frac{\rho_0 - \rho_i^{\text{adv}} - \sum_{j \neq i} a_{ij}p_j^l}{a_{ii}}\right). |
| 188 | +``` |
| 189 | + |
| 190 | +#### 2. Small diagonal elements |
| 191 | +If the diagonal element ``a_{ii}`` becomes too small or even zero, the simulation may become |
| 192 | +unstable. This can occur, for example, if a particle is isolated and receives little or no |
| 193 | +influence from neighboring particles, or due to numerical cancellations. In such cases, the |
| 194 | +updated pressure value is set to zero. |
| 195 | + |
| 196 | +There are also other options, like setting ``a_{ii}`` to the threshold value if it is beneath |
| 197 | +and then update with the usual equation, or just don't update the pressure value at all, and |
| 198 | +continue with the old value. The numerical error introduced by this technique remains small, |
| 199 | +as only isolated or almost isolated particles are affected. |
| 200 | + |
| 201 | + |
| 202 | +## Boundary Handling |
| 203 | + |
| 204 | +The previously introduced formulation only considered interactions between fluid particles |
| 205 | +and neglected interactions between fluid and boundary particles. To account for boundary |
| 206 | +interactions, a few modifications to the previous equations are required. |
| 207 | + |
| 208 | +First, the discretized form of the continuity equation must be adapted for the case in which |
| 209 | +a neighboring particle is a boundary particle. From now on, we distinguish between |
| 210 | +neighboring fluid particles (indexed by ``j``) and neighboring boundary particles (indexed by |
| 211 | +``b``). |
| 212 | + |
| 213 | +The updated discretized continuity equation becomes: |
| 214 | + |
| 215 | +```math |
| 216 | +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \nabla W_{ib}. |
| 217 | +``` |
| 218 | + |
| 219 | +Since boundary particles have zero velocity, the difference between the fluid |
| 220 | +particle's velocity and the boundary particle's velocity simplifies to just the fluid |
| 221 | +particle's velocity ``\bm{v}_{ib}(t+\Delta t) = \bm{v}_{i}(t+\Delta t)``. |
| 222 | +Accordingly, the predicted density ``\rho^{\text{adv}}`` becomes: |
| 223 | + |
| 224 | +```math |
| 225 | +\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \nabla W_{ib}(t). |
| 226 | +``` |
| 227 | + |
| 228 | +This leads to the following updated formulation of the linear system: |
| 229 | + |
| 230 | +```math |
| 231 | +\Delta t^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. |
| 232 | +``` |
| 233 | + |
| 234 | +Note that, since boundary particles are fixed, the force ``F_b^p`` is zero and does not appear in this equation. |
| 235 | + |
| 236 | +The pressure force acting on a fluid particle is computed as: |
| 237 | + |
| 238 | +```math |
| 239 | +\bm{F}_i^p(t) = -\sum_j m_j \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_j(t)}{\rho_j^2(t)} \right) \nabla W_{ij}(t) - \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). |
| 240 | +``` |
| 241 | + |
| 242 | +This also leads to an updated version of the equation for the diagonal elements: |
| 243 | +```math |
| 244 | +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_b m_b (-d_{bi}) \nabla W_{ib}. |
| 245 | +``` |
| 246 | + |
| 247 | +From this point forward, the computation of the coefficients required for the Jacobi scheme |
| 248 | +(such as ``d_{ii}``, ``d_{ij}`` etc.) depends on the specific boundary density evaluation method |
| 249 | +used in the chosen boundary model. |
| 250 | + |
| 251 | + |
| 252 | + |
| 253 | +### Pressure Mirroring |
| 254 | +When using pressure mirroring, the pressure value ``p_b`` of a boundary particle in the equation |
| 255 | +above is defined to be equal to the pressure of the corresponding fluid particle ``p_i``. |
| 256 | +In other words, the boundary particle "mirrors" the pressure of the fluid particle interacting |
| 257 | +with it. As a result, the coefficient that describes the influence of a particle's own |
| 258 | +pressure value ``p_i`` must also include contributions from boundary particles. Therefore, |
| 259 | +the equation for calculating the coefficient ``d_{ii}`` must be adjusted as follows: |
| 260 | + |
| 261 | +```math |
| 262 | +d_{ii} = -\Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij} - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. |
| 263 | +``` |
| 264 | + |
| 265 | +The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: |
| 266 | + |
| 267 | +```math |
| 268 | +\begin{align*} |
| 269 | +p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} |
| 270 | + - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj}p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right. \\ |
| 271 | +& \quad - \left. \sum_b m_b \sum_j d_{ij} p_j^l \nabla W_{ib} \right). |
| 272 | +\end{align*} |
| 273 | +``` |
| 274 | + |
| 275 | +### Pressure Zeroing |
| 276 | +If pressure zeroing is used instead, the pressure value of a boundary particle ``p_b`` |
| 277 | +is assumed to be zero. Consequently, boundary particles do not contribute to the pressure |
| 278 | +forces acting on fluid particles. |
| 279 | +In this case, the computation of the coefficient ``d_{ii}`` remains unchanged and is given by: |
| 280 | + |
| 281 | +```math |
| 282 | +d_{ii} = -\Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij}. |
| 283 | +``` |
| 284 | + |
| 285 | +The equation for the relaxed Jacobi iteration remains the same as in the pressure mirroring |
| 286 | +approach. However, the contribution from boundary particles vanishes due to their zero |
| 287 | +pressure: |
| 288 | + |
| 289 | +```math |
| 290 | +\begin{align*} |
| 291 | +p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} |
| 292 | + - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj}p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right. \\ |
| 293 | +& \quad - \left. \sum_b m_b \sum_j d_{ij} p_j^l \nabla W_{ib} \right). |
| 294 | +\end{align*} |
| 295 | +``` |
0 commit comments