Keywords: physically-based animation, smoothed particle hydrodynamics, deformable porous solids, fluid simulation, two-way coupling

1 Introduction↩︎

Most natural solids, like wood, soil, or wool, contain holes or small pores, which allow fluid to enter and flow through them. Even though this porous flow is very common, many simulations that couple solids and fluids ignore the porous aspects. This greatly reduces realism, since porous flow mechanics can have a huge influence on the macroscopic behavior of both phases.

The Smoothed Particle Hydrodynamics (SPH) method in particular is often used for multiphysics scenarios, since the purely particle-based simulation of all materials readily allows for unified discretization of interaction forces [1]. Still, we found the current state of the art to be lacking in regards to coupling solid and fluid phases through porous flow, which would further play into this strength of the SPH method. We want to simulate complex scenarios, driven by momentum conserving interaction forces between the fluid and solid phase. Additionally, we need to consider the available pore space, which is vital to determine how much fluid can be absorbed.

Previous SPH porous flow methods in computer graphics model absorption by scaling fluid particles based on how much of their mass is currently outside of the solid domain [2], [3]. Therefore, fluid particles can have vastly different sizes at the solid-fluid interface, that correspond to large mass differences between interacting particles. This leads to badly conditioned systems in implicit dual solvers [4] commonly used for SPH pressure computations [5]. The fluid movement inside the porous solid is then solved explicitly and simplified through the use of Darcy’s law, which describes stationary porous flow based on linear drag forces. These methods do not ensure balance of forces between the two phases and porous flow effects are modeled using heuristics for the pressure fields.

In this work, we propose a new SPH porous flow method, that does not require scaling of fluid particles and permits the usage of larger time steps by solving an implicit system. Our method allows fluid and solid particles to overlap by introducing a new density estimation that takes solid porosity into account. Further, we use an SPH pressure solver to guarantee fluid incompressibility even inside the porous domain. This allows us to model other porous flow effects with momentum conserving coupling forces, that are motivated by the physical forces acting between the fluid and solid, instead of using heuristics. Here, we take inspiration from granular flow models using the material point method (MPM) [6], [7], but derive new implicit interaction forces using state-of-the-art SPH solvers for effects like viscosity and surface tension. This allows us to present the, to the best of our knowledge, first implicit porous flow solver using SPH in computer graphics. For coupling with deformable solids, which we simulate using corotated linear elasticity, we propose modifications to include softening and bloating induced by absorbed fluid content and derive a strongly coupled implicit solver, which can handle large drag forces between solid and fluid phase. Finally, we propose a new model for capillary action using adhesion forces and include buoyancy based on fluid pressure, such that our method is able to simulate a wide range of porous flow effects (see Fig. 1).

In summary, our contributions are:

• the derivation of a physical model for porous flow, that includes drag, capillary action and buoyancy,

• an SPH framework that includes a pressure solver for overlapping particles in porous flow regions,

• modifications to linear stress computation to include porous solid effects, like softening and bloating,

• and a strongly coupled implicit linear formulation of the non-pressure forces.

2 Related Work↩︎

Our work investigates the macroscopic effects of the microscopic phenomenon of fluid moving within porous solids. We model this behavior using a continuum model, which is discretized using SPH. Based on the specific material properties, different approaches to model porous flow have already been developed. In this section, we will discuss related literature regarding general continuum porous-flow simulation using SPH-based methods in particular.

2.0.0.1 Smoothed Particle Hydrodynamics

The SPH method, being a truly mesh-free particle-based Lagrangian method, is well suited to model phenomena with complex and changing interfaces, such as free-surface or porous flow. In recent years, this method has seen many improvements regarding stability, efficiency, and the ability to model a wide range of physical effects. These include, among many others, fluid viscosity [8], [9], elasticity [10], [11] and surface tension [12]–[14]. These effects greatly influence porous phenomena, such as viscous drag or capillary action due to adhesion forces between the fluid and the solid phase. We therefore see great potential in using SPH to simulate porous flow and build on these models to derive strongly coupled interaction forces between the two phases.

To ensure that the fluid can only fill the available pore space, we use an SPH pressure solver to limit local fluid density. Modern pressure solvers ensure incompressibility by solving an implicit pressure Poisson equation [15]. [16] additionally show how density estimation can be adapted for the computation of pressure forces between non-dilutable fluid mixtures. We instead propose a new density estimation for overlapping particles based on porosity, which ensures fluid incompressibility while also prohibiting oversaturation of pores.

2.0.0.2 Porous Flow Simulation

Inspired by [2], many porous flow models in computer graphics use the SPH method for the non-absorbed part of the fluid phase, while the absorbed fluid mass is stored as a virtual saturation or wetness term on the solid discretization elements. In these methods, the number and size of fluid particles grows and shrinks to reflect how much fluid is currently outside of the porous object. Implementing this method is therefore challenging, as inserting particles in areas already filled with fluid is not a straightforward process and can lead to discontinuities [17]. The absorbed fluid on the other hand is transported between the solid elements based on diffusion equations. This approach can therefore be used with various solid simulation methods, like SPH [2], [18], the discrete element method (DEM) [19], discrete elastic rods (DER) [20] or mesh-based methods [21], [22].

Other approaches use multiphase methods, where each particle represents both solid and fluid at once [23], [24]. The exact local mixture is then tracked using volume fractions stored at the particles. These approaches are well suited to simulate dissolution processes, but the vastly different behavior of the two phases can make particle movement not well defined.

For other porous flow applications, an overlapping domain approach has therefore become more popular. Here, both phases are simulated based on their own constitutive models and respective discretization elements, which are allowed to occupy the same space as the other phase. This approach allows the use of specialized models for each phase and can make use of the inherent mass preserving property of particle-based simulation methods, as no mass transfer between discretization elements is required. [3] recently proposed a new SPH porous flow model that uses this kind of approach. While their method allows the coupling of porous materials with multiphase fluids, fluid particles entering the porous region is again made possible by reducing their size during the absorption process. [3] employ harmonic means for volume terms in the pressure solver to increase stability, but we still encountered issues in scenarios with fast porous flow where particles can rapidly change in size. In contrast to our approach, their method does not correct densities of absorbed fluid particles and only propagates them using Darcy’s law, using pore pressure heuristics in a way that does not consider momentum conservation between fluid and solid.

Other methods that allow particle overlapping instead simulate fluid transport through the porous domain using momentum-conserving interaction forces. This approach is mostly used for the mixture of granular materials and fluids [25], [26]. Our method follows this idea for general porous flow scenarios, while both phases are simulated using SPH.

We adapt the IISPH pressure solver [27] to allow overlapping particles, while still taking the presence of the other material into account. To the best of our knowledge, this is a novel contribution, as other porous flow approaches that use SPH do not account for the available pore space in the overlapping regions [3] or only use explicit pressure forces [26], that require small time step sizes.

Aside from SPH, particle-grid hybrid methods like MPM [28] are a popular alternative to simulate porous effects. MPM methods have been used to great effect to simulate porous materials, such as fluid interacting with hair [29], cloth [30], or granular materials [6], [7], just to mention a few. However, our goal is to enable physically consistent porous flow simulation within SPH frameworks, where MPM approaches cannot be directly applied due to the different discretization techniques.

3 Method↩︎

If the internal structure of a solid contains a connected network of pores, fluid from the outside can freely enter into and move through the resulting open space, giving rise to the porous flow phenomenon. Here, the very small scale of the pores results in macroscopic behavior not usually seen in other fluid simulation contexts. One well-known example is capillary rise, caused by surface tension forces, where fluid is pulled upwards against gravity. While porous flow properties often depend on the specific internal solid structure or fluid properties, we did not tune our model to specific materials. We instead provide a general porous flow framework that includes the most common associated effects and results in visually intriguing simulations with great artistic control.

In our work, we simulate porous flow from a macroscopic perspective, since it would be vastly inefficient to resolve all the microscopic pore details. We therefore adopt a continuum approach and account for the pore volume inside of a solid just by using the porosity value \(\phi \in [0,1)\). This value is a measure for the free space inside the porous object, such that a non-porous material has a porosity of \(\phi = 0\), while a porosity of \(\phi = 1\) corresponds to pure void space. Conceptually, this can be seen in Fig. 2, where a solid particle \({s}\) with mass \({{m}_{{s}}}\) has a volume of \({{V}_{{s}}^{0}}\), which is only partially occupied by the solid phase itself.

We use SPH to simulate both fluid and solid, where each particle represents a material parcel from one of the two phases. These two sets of particles behave according to their own equations of motion, but interact with particles of the other set through momentum conserving interaction forces. In this section, we first derive a new way to estimate particle densities for porous mixtures, which accounts for the pore space in regions where solid and fluid particles overlap. Based on these density estimations, we propose formulations for pressure accelerations that counteract compression and ensure that fluid absorption is limited by the porosity.

Next, we introduce implicit formulations for the different forces required to solve the equations of motion, to which we add solid-fluid coupling forces for drag, buoyancy and capillary action. These forces are assembled in a linear system for the particle velocities, that strongly couples the interaction forces with other non-pressure forces, like elastic deformation response and fluid viscosity.

3.1 Porous Densities↩︎

In the SPH method, a continuous material is sampled with particles, where each particle corresponds to a parcel with constant mass. Given the sampling volume \(V_i^0\) and the constant rest density \(\rho_i^0\) of the material, the mass \(m_i\) of particle \(i\) is \(m_i = \rho_i^0 V_i^0.\) The current density \(\rho_i\) is approximated using the mass of all nearby particles \(j\) of the same phase: \[\label{eqn:SPH95density95gen} \rho_i = \sum_j m_j W_{ij},\tag{1}\] where \(W_{ij} = W(\boldsymbol{\mathbf{x}}_i - \boldsymbol{\mathbf{x}}_j, h)\) is a normalized kernel with smoothing length \(h\), that assigns a weight based on the distance between the particle positions \(\boldsymbol{\mathbf{x}}_i\) and \(\boldsymbol{\mathbf{x}}_j\). This equation directly links particle spacing to density and \(\rho_i \approx \rho_i^0\) if each particle covers a volume that is equal to its sampling volume \(V_i^0\), in which case \(\sum_j V_j^0 W_{ij} \approx 1\) holds. We refer to the tutorial by [15] for the derivation and more information on SPH.

The density estimation plays a critical role in pressure solvers, which ensure that the distances between particles correspond to their physically accurate volume. While Eq. 1 gives a decent estimate if only one phase occupies the region around particle \(i\), in the case of porous flow we also need to consider the presence of the other phase. We therefore introduce a new density estimation to account for porosity, based on our solid sampling. To differentiate particles, we will use the following indices based on their phase: fluid particles \({i}\), fluid neighbors \({j}\), solid particles \({s}\), and solid neighbors \({t}\).

3.1.1 Porous Solid Density↩︎

We assume that the fluid phase and the solid phase itself (but not necessarily the pores) are incompressible. Both phases then have to fulfill the continuity equation, which links the material derivative of the density \(\frac{D \rho}{D t}\) to the velocity field \(\boldsymbol{\mathbf{v}}\): \[\label{eqn:conti} \frac{D \rho}{D t} = -\rho \nabla \cdot \boldsymbol{\mathbf{v}} \equiv 0.\tag{2}\]

Since we want to model porous behavior from a macroscopic perspective, we consider porous solids as continuous materials that include both the solid phase and the so-called void regions that can be filled with a fluid phase. We further differentiate between two types: those that allow compression of pores and those that do not. The first group includes materials like sponges, which can be compressed until all pores have collapsed, while the solid phase itself is incompressible. The second group includes materials like wood or compact soil, which are resistant to any compression due to the structural integrity of the solid skeleton. In this section, we derive an estimation for the density of the porous object, while the different behavior of the two groups is achieved using SPH pressure solver properties, as detailed in Sec. 3.2.

For the SPH discretization, we sample porous objects in their rest configuration using particles \(s\) with volumes \({{V}_{{s}}^{0}}\) (see Fig. 2). Given the constant rest density of the solid phase \({{\rho}_{{s}}^{0}}\), the particle \({s}\) has an effective rest density of \((1-\phi) {{\rho}_{{s}}^{0}} \leq {{\rho}_{{s}}^{0}}\) due to the constant rest porosity \(\phi\). We therefore define the solid particle mass as \({{m}_{{s}}} = (1 - \phi) {{\rho}_{{s}}^{0}} {{V}_{{s}}^{0}}.\) Inserting this particle mass definition into Eq. 1 and summing only over solid neighbors \(t \in \mathcal{N}^{S}\) results in: \[\label{eqn:SPH95density95solid} \begin{align} {{\rho}_{{s}}} &= (1 - \phi) {{\rho}_{{s}}^{0}} \sum_{{t}\in \mathcal{N}_{{{s}}}^S} {{V}_{{t}}^{0}} {W_{{s}{t}}}, \end{align}\tag{3}\] where we exploit that the solid has a constant rest density (\({{\rho}_{{s}}^{0}} = {{\rho}_{{t}}^{0}}\)). Due to our sampling choice, \({{\rho}_{{s}}}\) approximates the density \((1 - \phi) {{\rho}_{{s}}^{0}}\) of the porous object and not that of the solid phase itself. Therefore, \({{\rho}_{{s}}}\) should be constant for incompressible porous objects. In the case of compressible porous objects, the pores can be squeezed out such that densities greater than \((1 - \phi) {{\rho}_{{s}}^{0}}\) are allowed. This corresponds to solid particles moving closer together, as Eq. 3 increases for denser particle configurations. In this case, the particle’s estimated density \({{\rho}_{{s}}}\) should never exceed the solid phase density \({{\rho}_{{s}}^{0}}\), at which point the porous object does not contain void space anymore. Further volume loss corresponds to the solid phase itself being compressed, which we prohibit in our method. These density constraints can be simulated using SPH pressure solvers, as will be discussed later.

3.1.2 Fluid Density↩︎

We sample the continuous fluid phase with particles \({i}\), where we assume that the sampling volume \({{V}_{{i}}^{0}}\) is completely filled by the fluid phase. This is different from our definition of the porous solid phase, as shown in Fig. 2. Outside of the porous domain, we can use Eq. 1 to estimate the fluid particle densities \({{\rho}_{{i}}}\). However, when the fluid particle overlaps with the solid particles we need to account for the space already occupied by the solid phase.

We take inspiration from the number density approach [16], where all neighbors of particle \(i\) are included in the density estimation while all particles are treated as having the same phase as \(i\). This approach allows for consistent pressure coupling between different phases but prohibits particle overlapping. For porous flow, we instead want fluid particles to be able to occupy the same space as the solid particles, based on their porosity. We achieve this by only including the part that corresponds to the solid phase itself (without pores) in the density estimation of fluid particles. Since the porosity \(\phi\) measures pore space inside the porous solid, we can exclude this space from the solid particle sampling volumes \({{V}_{{s}}^{0}}\) by multiplying with the factor \(1-\phi\). We then approximate the density of fluid particle \({i}\) as: \[\label{eqn:SPH95density95covered} {{\hat{\rho}}_{{i}}} = {{\rho}_{{i}}^{0}} \sum_{{j}\in \mathcal{N}_{{{i}}}^F} {{V}_{{j}}^{0}} {W_{{i}{j}}}+ {{\rho}_{{i}}^{0}} \sum_{{t}\in \mathcal{N}_{{{i}}}^S} (1 - \phi) {{V}_{{t}}^{0}} {W_{{i}{t}}}.\tag{4}\] This results in a term for the solid particles that is similar to our proposed solid density estimation (Eq. 3 ), except that the sum approximating the local solid volume fraction is multiplied with the fluid rest density instead of the solid one. If no solid particles are close to the fluid particle, the term for the solid particles vanishes. In that case, Eq. 4 only contains terms from the fluid neighborhood \(\mathcal{N}^F\), and we recover Eq. 1 . Therefore, \({{\hat{\rho}}_{{i}}}\) is consistent both outside and inside the porous domain. The porosity \(\phi\) now limits how much fluid is allowed to overlap with the solid. This also influences fluid particle spacing, as shown in Fig. 3.

3.2 Pressure Solver↩︎

We enforce incompressibility by solving a pressure Poisson equation: \[\label{eqn:PPE} \Delta t\nabla^2 p = \frac{\rho^0 - \rho^*}{\Delta t},\tag{5}\] where \(p\) is the pressure, \(\rho^*\) is the predicted density based on the non-pressure forces and \(\rho^0\) is the rest density of the material.

As shown by [5], most implicit SPH pressure solvers use an equation like Eq. 5 to counteract compression. These solvers require the density change based on velocities, to be able to predict \(\rho^*\), and a formulation for the accelerations resulting from pressure gradients. The former we derive by taking the temporal derivative of Eq. 3 and Eq. 4 , where \(\boldsymbol{\mathbf{v}}_{{i}{j}} = \boldsymbol{\mathbf{v}}_{{i}} -\boldsymbol{\mathbf{v}}_{{j}}\): \[\label{eqn:SPH95density95change} \begin{align} \frac{D {{\rho}_{{s}}}}{D t} = &-(1 - \phi){{\rho}_{{s}}^{0}} \sum_{{t}\in \mathcal{N}_{{{s}}}^S} {{V}_{{t}}^{0}} \boldsymbol{\mathbf{v}}_{{t}{s}} \cdot \nabla {W_{{s}{t}}}, \\ \frac{D {{\hat{\rho}}_{{i}}}}{D t} = &-{{\rho}_{{i}}^{0}} \sum_{{j}\in \mathcal{N}_{{{i}}}^F} {{V}_{{i}}^{0}} \boldsymbol{\mathbf{v}}_{{j}{i}} \cdot \nabla {W_{{i}{j}}}\\ &- {{\rho}_{{i}}^{0}} \sum_{{t}\in \mathcal{N}_{{{i}}}^S} (1 - \phi) {{V}_{{t}}^{0}} \boldsymbol{\mathbf{v}}_{{t}{i}} \cdot \nabla {W_{{i}{t}}}. \end{align}\tag{6}\] Note that these equations are SPH approximations of the velocity divergence, which is consistent with Eq. 2 .

The pressure accelerations on the other hand result from the pressure gradient. For fluid particles, we treat the solid particles similar to boundary particles and adapt the method presented by [31] to our porous density estimation, with corrected sampling volumes as in Eq. 4 : \[\label{eqn:pressure95accel95fluid} \begin{align} {{\boldsymbol{\mathbf{a}}}_{{i}}^{\text{press}}} = &- {{\rho}_{{i}}^{0}} \sum_{{j}\in \mathcal{N}_{{{i}}}^F} {{V}_{{j}}^{0}} \left( \frac{{{p}_{{i}}}}{{{\rho}_{{i}}^{2}}} + \frac{{{p}_{{j}}}}{{{\rho}_{{j}}^{2}}} \right) \nabla {W_{{i}{j}}}\\ &- {{\rho}_{{i}}^{0}} \sum_{{t}\in \mathcal{N}_{{{i}}}^S} (1 - \phi) {{V}_{{t}}^{0}} \frac{{{p}_{{i}}}}{{{\rho}_{{i}}^{2}}} \nabla {W_{{i}{t}}}. \end{align}\tag{7}\] Here, we chose the symmetric formula for the gradient estimation, which is commonly used in SPH pressure solvers. This formulation produces pressure forces that encourage ordered particle patterns, that in turn enhance the accuracy of all other SPH approximations [32]. For solid particles, we do not consider any pressure forces acting from nearby fluid particles at this stage and instead use: \[\label{eqn:pressure95accel95solid} {{\boldsymbol{\mathbf{a}}}_{{s}}^{\text{press}}} = - (1 - \phi){{\rho}_{{s}}^{0}} \sum_{{t}\in \mathcal{N}_{{{s}}}^S} {{V}_{{t}}^{0}} \left( \frac{{{p}_{{s}}}}{{{\rho}_{{s}}^{2}}} + \frac{{{p}_{{t}}}}{{{\rho}_{{t}}^{2}}} \right) \nabla {W_{{s}{t}}}.\tag{8}\] The interaction with fluid particles is handled separately as an additional coupling force to simulate buoyancy. Since elastic forces also counteract volume change, which is not considered in the pressure solver, we achieve more plausible results by coupling this interaction force directly with the elasticity solver.

In SPH, it is common to use pressure clamping [15], such that only compression is considered as a violation of the continuity equation. Hence, the right-hand side of Eq. 5 is required to be negative for pressure forces to be active. For solid particles, we can use this property to switch between compressible and incompressible porous objects. In the first case, all densities \({{\rho}_{{s}}}\) smaller than \({{\rho}_{{s}}^{0}}\) are acceptable. For incompressible porous materials, we instead only allow densities under \((1 - \phi) {{\rho}_{{s}}^{0}}\), to reflect that we do not allow reduction of the pore space. This is achieved by replacing the rest density term in Eq. 5 . For fluid particles, we always require that their density is smaller or equal to the rest density \({{\rho}_{{i}}^{0}}\). More details on this and the derivation of the density constraints are given in the supplementary document.

We used a modified IISPH pressure solver [27] to solve Eq. 5 , using the updated density change and pressure acceleration terms. This pressure solver uses the predicted density changes to compute pressure values for each particle, such that the resulting accelerations remove any density deviations. The resulting implicit linear system is solved using the relaxed Jacobi method, where in each iteration the pressure values are clamped to be positive, such that they only counteract compression. For boundary handling we use volume maps [33] or the approach by [34], but our method does not have strict requirements in this regard.

3.3 Porous Interaction Forces↩︎

In the previous section, we presented how we can allow fluid and solid particles to overlap, while still adhering to the continuity equation. The pressure forces do not depend on any specific pore structure or fluid properties, since they directly follow from incompressibility and the available pore space given by the porosity value. Next, we need to cover all other forces that act on the fluid and solid particles according to their equations of motion. We start with the porous interaction forces, which are less well defined and greatly depend on the material properties. Here, we chose force formulations that not only cover most of the common porous flow effects but can also be combined into a unified implicit solver.

For a solid particle \({s}\), we use a modified Cauchy momentum equation, which includes external forces \({{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{ext}}}\) and the porous coupling forces \({{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{pore}}}\): \[\label{eqn:solid95momentum} \frac{D {{\boldsymbol{\mathbf{v}}}_{{s}}}}{D t} = \frac{1}{ {{\rho}_{{s}}}} \nabla \cdot \boldsymbol{\sigma} + \frac{1}{{{m}_{{s}}}} {{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{pore}}} + \frac{1}{{{m}_{{s}}}} {{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{ext}}},\tag{9}\] where \(\boldsymbol{\sigma}\) is the local stress tensor. For a fluid particle \({i}\), we instead add the porous forces to the Navier-Stokes momentum equation: \[\label{eqn:navier95stokes} \frac{D {{\boldsymbol{\mathbf{v}}}_{{i}}}}{D t} = -\frac{1}{{{\hat{\rho}}_{{i}}}} \nabla p + \frac{{{\mu}_{{i}}}}{{{\hat{\rho}}_{{i}}}} \nabla^2 \boldsymbol{\mathbf{v}} + \frac{1}{{{m}_{{i}}}} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{pore}}} + \frac{1}{{{m}_{{i}}}} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{ext}}}.\tag{10}\] For our method, we formulate \(\boldsymbol{\mathbf{f}}^{\text{pore}}\) as the sum of the most common effects observed in porous flow, which are capillary action \(\boldsymbol{\mathbf{f}}^{\text{cap}}\), drag \(\boldsymbol{\mathbf{f}}^{\text{drag}}\), and buoyancy \(\boldsymbol{\mathbf{f}}^{\text{buo}}\). Often, the drag forces inside the porous solid are very strong, which makes explicit solvers highly unstable. Instead, we propose implicit formulations for all non-pressure forces, that depend linearly on the velocities \(\boldsymbol{\mathbf{v}}^{n+1}\) of the next time step \(n+1\) and can easily be combined into one strongly coupled system.

3.3.1 Capillary Action↩︎

Capillary forces allow fluid to rise inside the porous material, due to adhesive forces acting between the solid and the fluid. To model this effect, we simplify the adhesion model by [14] and formulate attraction force between fluid and solid using a capillary coefficient \(C^{\text{cap}}\) as: \[\label{eqn:adhesion} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{cap}}} = -\sum_{{t}\in \mathcal{N}_{{{i}}}^S} C^{\text{cap}}({{S}_{{t}}^{n}}) \frac{\bar{m}_{{i}{t}}}{\bar{\rho}_{{i}{t}}^{n}} \boldsymbol{\mathbf{x}}_{{i}{t}}^{n+1} W_{{i}{t}}^n,\tag{11}\] where \(\bar{m}_{ij} = \frac{m_i + m_j}{2}\) and \(\bar{\rho}_{it} = \frac{\rho_i + (1 - \phi) \rho_t^0}{2}\) are averaged to ensure momentum balance, \(S\) is the saturation of the solid particle and \(\boldsymbol{\mathbf{x}}_{i t} = \boldsymbol{\mathbf{x}}_{i} - \boldsymbol{\mathbf{x}}_{t}\). For a constant \(C^{\text{cap}}\), this formulation yields a vector which points into the direction with the most solid neighbors, such that the resulting acceleration will pull fluid particles inside the porous body. To include increased attraction forces towards unsaturated regions we further weight them based on solid particle saturation.

The saturation \({{S}_{{s}}}\) measures how much the local pores are already filled with fluid, which we estimate by comparing the fraction of the local fluid and solid volume with the porosity: \[\label{eqn:saturation} {{S}_{{s}}} = \frac{1}{\phi} \frac{\sum_{{j}\in \mathcal{N}_{{{s}}}^F} {{V}_{{j}}^{0}} {W_{{s}{j}}}}{\sum_{{t}\in \mathcal{N}_{{{s}}}^S} \frac{{{m}_{{t}}}}{{{\rho}_{{t}}}} {W_{{s}{t}}}}.\tag{12}\] To avoid issues arising from interpolation errors at the solid-fluid interface, we clamp saturation to ensure \({{S}_{{s}}} \in [0,1]\). We then decrease the capillary potential for increasing saturation, based on a falloff factor \(\eta_{\text{cap}} \in [0,1]\): \[\label{eqn:adhesion95mod} C^{\text{cap}}({{S}_{{s}}}) = C^{\text{cap},0} (1 - {{S}_{{s}}} \eta_{\text{cap}}),\tag{13}\] where the adhesion \(C^{\text{cap},0}\) determines the maximum height fluid can reach based on capillary action.

As capillary pressure depends on the channel radius, smaller pores can pull fluid much higher than the larger pores can, such that only a portion of the solid will be fully saturated [35]. Our model allows direct control over the final height of this fully saturated region, based on the falloff factor \(\eta_{\text{cap}}\). As shown in Fig. 4, this front will move upwards and closer to the maximum height when decreasing \(\eta_{\text{cap}}\). [2] also propose to decrease capillary potential based on saturation, but use a power function which always reduces the adhesion force to zero for fully saturated regions. For \(\eta_{\text{cap}} < 1\), we instead include adhesion forces for saturated solid particles, which then still attract fluid particles that are close to the porous solid boundary. This absorption force in turn results in pressure forces acting on already absorbed fluid, which is pushed further inside towards regions with available pore space, if there are no opposing forces stronger than the adhesion. This allows gravitational, pressure and capillary forces to balance out in each time step, resulting in stable configurations even when the fully saturated region lies above the outer fluid surface.

For the averaged density \(\bar{\rho}_{it}\) we use \({{\rho}_{{i}}}\) instead of \({{\hat{\rho}}_{{i}}}\) for the fluid particle density. Since the sum in Eq. 11 approximates an integral over the whole volume around the fluid particle, not just the integral over a single phase, we require the approximated particle influence volumes \(\frac{m}{\rho}\) to reflect this. We therefore use the densities averaged over the whole surrounding volume, not just the volume occupied by the phase.

While the solid particle densities \({{\rho}_{{s}}}\) already include the pore volumes, we could directly use them in the force estimation. But for large adhesion forces, we observed simulation artifacts when using the approximated solid density. Since solid particles at the surface of the porous object have a smaller estimated density (due to missing neighbors), their volume is greatly overestimated and fluid particles close to but still outside the porous domain experience unnaturally strong capillary forces. We instead use the rest density \((1 - \phi) \rho_s^0\) here to solve this problem.

The forces acting on the porous solid follow from momentum balance, since \(\boldsymbol{\mathbf{f}}^{\text{cap}}_{{i}\leftarrow {s}} = - \boldsymbol{\mathbf{f}}^{\text{cap}}_{{s}\leftarrow {i}}\) must hold: \[\label{eqn:adhesion95mirror} {{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{cap}}} = -\sum_{{j}\in \mathcal{N}_{{{s}}}^F} C^{\text{cap}} ({{S}_{{s}}^{n}}) \frac{\bar{m}_{{s}{j}}}{\bar{\rho}_{{s}{j}}^{n}} \boldsymbol{\mathbf{x}}_{{s}{j}}^{n+1} W_{{s}{j}}^n.\tag{14}\] This formulation of capillary adhesion does not directly include the velocities, but the particle positions. Since we aim for formulations linear in \(\boldsymbol{\mathbf{v}}^{n+1}\), we approximate \(\boldsymbol{\mathbf{x}}^{n+1} = \boldsymbol{\mathbf{x}}^{n} + \Delta t \boldsymbol{\mathbf{v}}^{n+1}\) and split the force into two terms \(\mathbf{f} = \mathbf{f}(\mathbf{v}^{n+1}) + \mathbf{b}^n\), e.g. for fluid particles: \[\label{eqn:adhesion95semi95implicit} \begin{align} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{cap}}} \left(\boldsymbol{\mathbf{v}}^{n+1} \right) &= - \Delta t \sum_{{t}\in \mathcal{N}_{{{i}}}^S} C^{\text{cap}}({{S}_{{t}}^{n}}) \frac{\bar{m}_{{i}{t}}}{\bar{\rho}_{{i}{t}}^n} \boldsymbol{\mathbf{v}}_{{i}{t}}^{n+1} W_{{i}{t}}^n, \\ {{\boldsymbol{\mathbf{b}}}_{{i}}^{\text{cap},n}} &= - \sum_{{t}\in \mathcal{N}_{{{i}}}^S} C^{\text{cap}}({{S}_{{t}}^{n}}) \frac{\bar{m}_{{i}{t}}}{\bar{\rho}_{{i}{t}}^n} \boldsymbol{\mathbf{x}}_{{i}{t}}^{n} W_{{i}{t}}^n. \end{align}\tag{15}\] The term \(\mathbf{b}^{\text{cap},n}\) only depends on known values and will contribute to the right-hand side of the implicit system.

3.3.2 Drag↩︎

As proposed by [26], the drag forces between the fluid and solid particles can be simulated as viscous friction. In our work, we express these drag forces as: \[\label{eqn:SPH95viscosity95porous} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{drag}}} \left(\boldsymbol{\mathbf{v}}^{n+1} \right) = \tilde{d} \sum_{{t}\in \mathcal{N}_{{{i}}}^S} \frac{\mu^{\text{por}}}{1 - \phi} \frac{{{m}_{{i}}} {{m}_{{t}}}}{{{\rho}_{{i}}^{n}} {{\rho}_{{t}}^{0}}} \frac{\boldsymbol{\mathbf{v}}_{{i}{t}}^{n+1} \cdot \boldsymbol{\mathbf{x}}_{{i}{t}}^{n}}{\|\boldsymbol{\mathbf{x}}_{{i}{t}}^{n} \|^2 + 0.01 h^2} \nabla W_{{i}{t}}^{n},\tag{16}\] where \(\mu^{\text{por}}\) is the porous viscosity coefficient and \(\tilde{d} = 2(d+2)\) where \(d\) is the number of space dimensions. Eq. 16 is an approximation of the relative velocity Laplacian, leading to a force similar to the fluid viscosity term in the Navier-Stokes equation (Eq. 10 ). The drag forces acting on the solid particles then follow from force balance. The implicit formulation is inspired by the viscosity method by [9], which we modified to achieve symmetric forces between particles of different mass.

3.3.3 Buoyancy↩︎

Some porous materials float due to the pressure acting on them by the surrounding fluid. We achieve this effect using pressure forces acting from fluid on solid particles. Including this force in the implicit pressure solver would allow fluid particles to push away solid particles, without regard for the elastic forces keeping them in their relative configurations. We therefore consider the solid particles to be fixed when computing the pressure forces for the fluid particles, which results in a more accurate estimate for interactions with relatively stiff solid objects with open pores, where pressure forces are rather resolved by fluid displacement than solid deformation.

Based on the fluid pressure acceleration in Eq. 7 , the sum of the pressure forces acting from solid particle \({s}\) on fluid particles \({j}\) results in the mirrored force: \[\label{eqn:buoyancy95interaction95forces} \boldsymbol{\mathbf{b}}_{{s}}^{\text{buo},n} = - (1 - \phi) {{V}_{{s}}^{0}} \sum_{{j}\in \mathcal{N}_{{{s}}}^F} {{m}_{{j}}} {{\rho}_{{j}}^{0}} \left( \frac{{{p}_{{j}}}}{{{\hat{\rho}}_{{j}}^{2}}} \right) \nabla W_{{s}{j}}.\tag{17}\] Here, we use the pressure forces as determined in the most recent pressure solve step to approximate the fluid pressure gradient, such that this force does not depend on particle velocities. As this interaction force is already included in the pressure solver for the fluid particles, those do not experience a buoyancy force (\(\boldsymbol{\mathbf{f}}_{{i}}^{\text{buo}} = \boldsymbol{\mathbf{0}}\)).

3.4 Non-Coupling Forces↩︎

For the other terms in the momentum balance equations (Eq. 9 and Eq. 10 ), there already exist many specialized SPH solvers. In this section we will show how to adapt the implicit fluid viscosity solver by [9] and the corotated linear elasticity solver by [10] to our porous flow model. We chose those solvers, since they both use force formulations which are linear in velocity and can therefore be easily combined into a strongly coupled implicit system with our porous interaction forces.

3.4.1 Fluid Viscosity↩︎

For the fluid viscosity, we use a variation of Eq. 16 to achieve symmetric forces between the particles: \[\label{eqn:SPH95viscosity} {{\boldsymbol{\mathbf{f}}}_{{i}}^{\text{vis}}} \left(\boldsymbol{\mathbf{v}}^{n+1} \right) = \tilde{d} \sum_{{j}\in \mathcal{N}_{{{i}}}^F} \mu^{\text{vis}} \frac{{{m}_{{i}}} {{m}_{{j}}}}{{{\hat{\rho}}_{{i}}^{n}} {{\hat{\rho}}_{{j}}^{n}}} \frac{\boldsymbol{\mathbf{v}}_{{i}{j}}^{n+1} \cdot \boldsymbol{\mathbf{x}}_{{i}{j}}^{n}}{\|\boldsymbol{\mathbf{x}}_{{i}{j}}^{n} \|^2 + 0.01 h^2} \nabla W_{{i}{j}}^{n}.\tag{18}\] Here we use the fluid phase densities \(\hat{\rho}\) instead of \(\rho\), as this is an internal fluid force and we therefore only want to integrate over the domain occupied by the fluid.

3.4.2 Solid Elasticity↩︎

For the porous elastic solid, we assume that the elasticity coefficients are defined to model the behavior of the whole porous body, including pores. Then we can apply the elastic solid solver by [10] directly, using \({{V}_{{s}}^{0}}\) for the particle rest volumes and \({{\rho}_{{s}}}\) for the density, to model elastic forces: \[\label{eqn:SPH95elasticity} \frac{1}{{{m}_{{s}}}}{{\boldsymbol{\mathbf{f}}}_{{s}}^{\text{elast}}} = \frac{1}{ {{\rho}_{{s}}}} \nabla \cdot \boldsymbol{\sigma}_{{s}} \left( \boldsymbol{\mathbf{x}}^n, \boldsymbol{\mathbf{v}}^{n+1} \right).\tag{19}\] We refer to the work by [10] for details on how to compute the stress tensor divergence, which can be split into a term that depends linearly on particle velocities \(\boldsymbol{\mathbf{v}}^{n+1}\) and a term \(\mathbf{b}^{\text{elast},n}\) that only depends on the current particle positions \(\boldsymbol{\mathbf{x}}^{n}\).

3.4.3 Porous Solid Effects↩︎

We further modify the stress tensor computation to include porous material effects like bloating and saturation dependent elastic behavior. This results in the modified computation of the stress tensor \({{\boldsymbol{\sigma}}_{{s}}}\): \[\label{eqn:stress} {{\boldsymbol{\sigma}}_{{s}}} = 2 \mu({{S}_{{s}}}) {{\boldsymbol{\epsilon}}_{{s}}} + \lambda({{S}_{{s}}}) \text{tr} ( {{\boldsymbol{\epsilon}}_{{s}}} ) \mathbb{1}+ {{\boldsymbol{\sigma}}_{{s}}^{\text{bloat}}}.\tag{20}\]

One possible cause for bloating is the increase of volumetric pressure inside filled pores, as the absorbed fluid has a higher pressure than air. This in turn can cause the pores to expand for soft solids. For biological materials, like wood or natural sponges, expansion is instead more prevalently caused by hygroscopic swelling [35]. Here, water not only flows through the pores, but is also absorbed into the cells of the solid skeleton. The solid itself then increases in volume, observable as a volumetric strain \(\epsilon_h = \eta_{h} S\) based on the hygroscopic expansion coefficient \(\eta_{h}\). Following [2], instead of modifying the strain tensor, we add a bloating term to the stress tensor: \[\label{eqn:bloating} {{\boldsymbol{\sigma}}_{{s}}^{\text{bloat}}} = - \eta_{\text{bloat}} {{S}_{{s}}} \mathbb{1},\tag{21}\] controlled by the bloating factor \(\eta_{\text{bloat}}\), which can be understood as a hygroscopic expansion coefficient that is premultiplied with the elasticity parameters. This term modifies the rest configuration of the solid particles and is therefore independent of the current deformation measured by the strain tensor \({{\boldsymbol{\epsilon}}_{{s}}}\). In our implicit solver for the particle velocities, given in Sec. 3.5, it is only included on the right-hand side of the system.

The elastic response of a solid then also depends on the Lamé coefficients \(\mu\) and \(\lambda\). Many porous materials change their elastic behavior when absorbing fluid. One example is wood, which often becomes more malleable with increasing moisture [36]. For our model, we include this property of porous solids by adjusting the Lamé coefficients based on saturation: \[\label{eqn:elasticity95change} \mu ({{S}_{{s}}}) = (1 + \eta_m {{S}_{{s}}}) \mu^0, \quad \lambda ({{S}_{{s}}}) = (1 + \eta_l {{S}_{{s}}}) \lambda^0,\tag{22}\] with change factors \(\eta_m, \eta_l\). Positive change factors increase and negative ones decrease resistance against deformation, while in the latter case the Lamé coefficients have to be clamped to ensure positivity. Even though this simplified linear model for elasticity change is motivated purely phenomenologically, it can achieve visually plausible results. Note that these adjustments can be replaced with more fitting ones for materials where the coefficient change is not monotone or very non-linear. This is the case for sand, where a very specific amount of added water can greatly increase pile stability, while further increasing saturation causes the grains to flow like a liquid [37].

3.5 Solver↩︎

At the beginning of each time step, we first compute particle densities (Eq. 1 , Eq. 3 and Eq. 4 ) and solid particle saturation (Eq. 12 ). Then we solve all non-pressure forces to predict the new particle positions. Solving the porous interaction forces from Section 3.3 in an explicit fashion becomes highly unstable for large drag forces. We therefore propose to use an implicit solver.

Given the mass matrix \(\boldsymbol{\mathbf{M}}\) and forces \(\boldsymbol{\mathbf{f}}\) which are linear in \(\boldsymbol{\mathbf{v}}^{n+1}\), we write the velocity update as: \[\label{eqn:velocity95update} \begin{align} \boldsymbol{\mathbf{M}} \boldsymbol{\mathbf{v}}^{n+1} - \Delta t \boldsymbol{\mathbf{f}}\left( \boldsymbol{\mathbf{v}}^{n+1}\right) \approx \boldsymbol{\mathbf{M}} \boldsymbol{\mathbf{v}}^{n} + \Delta t \boldsymbol{\mathbf{b}}^{n}. \end{align}\tag{23}\] At timestep \(n\) we solve this linear system for \(\boldsymbol{\mathbf{v}}^{n+1}\), which is the stacked vector of all particle velocities for the next time step. The total force \(\boldsymbol{\mathbf{f}}\) and right-hand side \(\boldsymbol{\mathbf{b}}\) are the sums of all the corresponding terms in Section 3.3 and Section 3.4. Due to our choice of discretization, these forces can readily be combined into one linear system and solved using a matrix-free conjugate gradient method. More details about this linear system are given in our supplementary document.

As proposed by [9], we add the velocity updates \(\boldsymbol{\mathbf{v}}^{n} - \boldsymbol{\mathbf{v}}^{n-1}\) from the last timestep to the current velocities \(\boldsymbol{\mathbf{v}}^{n}\) to serve as an initial guess when solving the system for \(\boldsymbol{\mathbf{v}}^{n+1}\). After the new velocities are obtained, particle densities are predicted using Eq. 6 , which are corrected by pressure accelerations, as described in Section 3.2. These accelerations in conjunction with the velocities \(\boldsymbol{\mathbf{v}}^{n+1}\) are then integrated using the symplectic Euler method, resulting in new particle positions for the next time step.

4 Results↩︎

The following results were simulated on an AMD Ryzen Threadripper PRO 5975WX CPU with 32 cores, 3.60 GHz base core frequency and 256 GB of RAM. We have implemented our method into SPlisHSPlasH [38], an SPH fluid simulation library written in C++, and plan to release our code as open source. To increase performance, we employ AVX vectorization and parallelization using OpenMP. Surface reconstructions were done using the splashsurf library [39] or mesh skinning as described by [11]. We assume that the porous object is incompressible if not stated otherwise. Additional simulation parameter choices as well as runtime measurements are given in the supplementary document.

4.1 Stability↩︎

Using the explicit formulation of Eq. 16 (i.e., replacing \(\boldsymbol{\mathbf{v}}^{n+1}\) with \(\boldsymbol{\mathbf{v}}^{n}\)) we found that for a given time step size of \(1 \si{\milli \second}\) and \(\mu^{\text{por}} = 10 \si{\pascal \second}\), the resulting force cannot sufficiently reduce the fluid particle speed upon contact with the porous solid and they instead immediately permeate far inside the block, as shown in Fig. 5. Increasing the viscosity coefficient does not solve this issue and instead leads to instabilities. Implicit forces on the other hand can notably slow down the fluid particles even with the lower coefficient and are still stable for the higher one.

In Fig. 6 we show simulations of a cube (\(\phi = 0.6\)) placed between two shallow pools of water, where only one half is shown. We chose a significant capillary action (\(C^{\text{cap},0} = 1000 \si{\newton \per \meter}\)) to compare our work with the model by [3] (\(p_0 = 10000 \si{\pascal}\)). While our method can stably simulate both slow (\(k = 1e^{-6}\) or \(\mu^{\text{por}} = 50 \si{\pascal \second}\)) and fast (\(k = 1e^{-5}\) or \(\mu^{\text{por}} = 3 \si{\pascal \second}\)) porous flow, the method by [3] suffers from pressure spikes originating in dense fluid particle clumps at the solid-fluid interface, where the particles are scaled according to the absorption model. This results in explosions when particles rapidly increase in size, which happens more frequently in scenarios with fast porous flow. While this issue could be addressed with stabilization techniques, e.g., by introducing an upper limit on how much particles are allowed to grow in a time step, our method does not have these concerns since we completely avoid severe fluid particle scaling.

On the other hand, their method is faster, which can be attributed to the missing pressure computations for absorbed fluid particles and the only one-way coupled explicit solver for the porous flow. For the example with slow flow, our method took on average \(194 \si{\milli \second}\) per time step, while the the method by [3] only took \(126 \si{\milli \second}\). Note that, while our reimplementation of their method closely follows code published by the authors, we simulated these examples on the CPU without AVX vectorization, such that their method could potentially perform even faster.

4.2 Validation↩︎

We further investigated the well-researched use case of fluid transport through a trapezoidal dam, shown in Fig. 7. This scenario does not exhibit obvious capillary action, which is difficult to simulate using the method by [3], since we could not control porous flow speed independently of capillary height due to their definition of the pore pressure. We therefore extend their approach with a version of Darcy’s law that includes gravitational forces: \(\boldsymbol{\mathbf{v}}_{fs} = - \frac{k}{\phi \mu}(\nabla p - \rho_f \boldsymbol{\mathbf{g}})\), for which the derivation is given in the supplementary document. While this allows control over the upwards pull on the absorbed fluid through the rest pore pressure \(p_0\), we still could not find parameters that produce the expected result, which we determine using the Casagrande solution [40] for the fluid surface. Additionally, the fluid does not leave the dam at the opposite side and is instead greatly compressed at the bottom, due to the lack of proper pressure forces acting on absorbed fluid particles. Our method on the other hand can easily handle the common case of porous flow without capillary action and produces the expected result.

Furthermore, we showcase the simulation of objects sinking due to absorbed fluid mass. Fig. 8 shows porous spheres being dropped in a water bath, where our method correctly captures their expected floating behavior based on the solid phase density. Only the first sphere floats, while the third sphere sinks faster than the second one due to the smaller buoyancy force.

4.3 Parameter Studies↩︎

To demonstrate the effect of our model parameters, we simulate fluid seeping down into a porous block. In particular, we are interested in the influence of the porosity \(\phi\), the viscous drag coefficient \(\mu^\text{por}\) and capillary action controlled by \(C^\text{cap,0}\). First, we analyze the porosity (see Fig. 9), which we vary between \(\phi \in \{ 0.4, 0.6, 0.8 \}\). The final state showcases the expected result that the amount of absorbed liquid is linearly proportional to the porosity value.

Next, we investigate the effects of viscous drag and capillary action, while keeping the porosity (\(\phi = 0.5\)) constant. The results in Fig. 10 show how we can control the uniform seeping speed using the viscous drag coefficient \(\mu^{\text{por}}\), which slows down fluid particles when they interact with solid particles. The capillary force on the other hand introduces a pull towards unsaturated regions, which increases fluid velocity and leads to a smoother transition between dry and wet regions. Of these simulations, which each use 110.6k fluid and 110.6k solid particles, the version with adhesion took the longest to compute with an average time of \(373.7 \si{\milli \second}\) per time step. Here, \(69 \%\) of the time was spent for the pressure solver and \(20 \%\) to compute the strongly coupled non-pressure forces.

To show the effect of the capillary force parameters, we place a cube (\(\phi = 0.6\)) between two shallow pools of water, shown in Fig. 11. By varying the capillary force parameter \(C^\text{cap,0}\) we can influence both the speed and final height of the absorption. We first show the simulation without capillary effect (\(C^{\text{cap},0} = 0 \si{\newton \per \meter}\)), where the fluid permeates only the bottom of the porous medium. Increasing the capillary force pulls the fluid upwards, such that it can reach the top of the object. We also show that the falloff \(\eta_\text{cap}\), which determines how much the capillary potential diminishes with saturation, controls how much fluid is pulled up into the lower regions of the cube (see also Fig. 4), while the total height of the wet region only depends on \(C^{\text{cap},0}\).

4.4 Visual Showcases↩︎

Lastly, we demonstrate that our method can produce complex simulation results with great imaginative liberty. Shown in Fig. 1 is the simulation of a compressible sponge, which is first dropped into a water bath and then subjected to external compression. We cause the sponge to increase its volume by setting the bloating factor to \(\eta_\text{bloat} = 0.01\). The compression causes solid particles to move closer together, as shown in Fig. 12, and the corresponding reduction in available pore space causes the absorbed water to be ejected.

Fig. 13 shows fusilli-shaped porous objects being submerged in water and then spilled into a container. The objects are initially very stiff, but setting the shear resistance change factor \(\eta_m = -2.0\), while ensuring \(\mu(S_s) \geq 0.001 \mu^0\), causes the material to become soft as it absorbs water.

For simulations that require more turbulent behavior, our solver can also be used in conjunction with SPH fluid vorticity methods. For the simulation in Fig. 14 we use the micropolar model by [41] to better capture the behavior of fluid accelerated by a rotating propeller. Three porous bunnies with a low density are dropped into the fluid bath, where they interact with the surrounding flow field, creating interesting dynamic movements.

We also show capillary action with a more artistic scene, using letters made of porous material (see Fig. 15). The solid is fixed in space, while the fluid enters from the left and is transported up into the letters through capillary action.

5 Conclusion↩︎

While previous SPH porous flow methods in computer graphics suffer from stability issues caused by severe runtime particle scaling, our method not only avoids these issues entirely, but further increases stability with the use of implicit formulations and strong coupling of non-pressure forces. We achieve this by considering porosity within the pressure solver, which allows solid and fluid particles to overlap while still taking incompressibility into account. We can simulate scenarios that feature the most common porous flow effects, employing viscous drag forces and capillary action due to adhesion. Overall, our SPH porous flow method not only improves on the current state of the art in terms of stability, but also in physical soundness, due to the use of momentum conserving interaction forces and consistent pressure computations.

Still, our solver has some limitations. So far, we only incorporated corotated linear elastic solids. Investigating the use of a different solver for the solid phase is an interesting course of further research, like the use of a non-linear elasticity solver (e.g., [42]) or one for granular materials (e.g., [43]). We believe our approach of coupling fluid and solid phase has great potential to include even more porous materials. Additionally, we have not yet investigated non-linear drag models, as we already managed to achieve plausible results with the linear one. While non-linear drag forces would further increase the physical soundness of our model, deriving an efficient non-linear solver is not trivial and therefore remains as future work.

In our experiments, we observed that the pressure solver is a bottleneck regarding simulation speed. Using a more efficient pressure solver like divergence-free SPH [44] could increase our performance, but requires further modifications to allow intended velocity divergences occurring in porous flow absorptions and ejections.

References↩︎