Simulation of weak‐inertia single‐phase flow in porous materials using Smoothed Particle Hydrodynamics

Intrinsic permeability of a porous material is a crucial material parameter in various application fields like e.g. geosciences, materials science, biological sciences and mechanical engineering. In contrast to the permeability often used in purely linear models, however, it is also a (nonlinear) function of the Reynolds number (Re) with the onset of the influence of inertial forces. In coarse‐grained continuum models the Forchheimer equation describes this incipient influence of viscous momentum exchange but it is based on the adaptation of experimentally obtained data and its range of validity remains a subject for discussion. Performing three‐dimensional Direct Numerical Simulations for single‐phase flow through porous media this work aims to show a smooth transition of effective intrinsic permeability within the interval 0.1 < Re < 1000. We use the fully‐Lagrangian meshless Smoothed Particle Hydrodynamics method (SPH) to solve the weakly compressible Navier‐Stokes equations for pore‐scale resolved flow in the so‐called weak‐inertia regime. Our implementation utilizes HOOMD‐blue thus allows us to perform massively parallel computations which are crucial for these applications.


Introduction
For traditional grid-based flow simulations on the continuum scale, e.g. in groundwater research, scalar values for the permeabilities that can be determined from experiments are often sufficient. In addition, however, the understanding of processes at the pore scale it is becoming more and more important. For Direct Numerical Simulations (DNS) of flow through regular porous materials and XRCT-scanned sandstone samples, we employ the kernel based fully-Lagrangian particle method SPH [3].

Model
Mathematical Model: The identity for evaluating the vector-valued field function f (x) (see Eq. (1) 1 ) is the basis of SPH. In a first step the Dirac delta function δ(x − x ) is substituted by a kernel W (x − x , h). W has compact support, its value is depending on the particle distance and it is evaluated within the smoothing length h. The continuous kernel based function is subsequently transformed into a discrete weighted sum (Eq. (1) 2 ) of the field functions of neighbouring particles j, to determine the field function of a particle i, Here ρ j , m j and N represent density, mass and number of neighbouring particles, respectively. Divergence and gradient of vector fields or scalar fields can be derived and do not require the derivation of the field but only the derivation of the kernel function, which is computationally beneficial. Details about the discretisation of the Navier-Stokes equations with SPH can be found, e.g. in [3,5].
Computational Aspects: The solver is implemented on top of the software framework HOOMD-blue [1,2] since this allows for massively parallel CPU and GPU computations. The implementation includes the calculation of kernel, computation of field functions (ρ f , p,v f , v f ) and time integration whereas neighbour search and load balancing is handled by the HOOMDblue backend. A sufficient scaling behaviour as well as the numerical accuracy is demonstrated [4].
Accordingly, the total force comprises a volumetric force on a particle F B i = m i g as well as the pairwise dissipative viscous interaction forces F V ij and the pairwise conservative pressure interaction forces F P ij , whereby we distinguish between fluid-fluid and solid-fluid interactions [5].
Computation of permeability and Reynolds number: In classical continuum mechanics, coarse-grained values are used for the permeability. Here we use volume-averaged values to compute the Reynolds number and the intrinsic permeability, wherein || • || is the Euclidean Norm, κ s represents the intrinsic permeability and q is the filter velocity.

Results
For the following investigations regular three-dimensional porous structures (40 3 voxels) are used, which are similar to sandstones in their internal structure. A comparable study was carried out for two-dimensional artificial domains [6]. Due to the periodicity of pressure (boundary conditions), a body force g is set as the driving force and we vary the dynamic viscosity µ f to adjust Re.

Conclusions and future work
We have shown that SPH can be used to determine κ s (Re) relationships for regular porous structures and to generate velocity fields for non-artificial complex natural structures. To investigate the above relationships for actual reservoir rocks, larger and thus representative structures are necessary. In order to be able to represent the inherent structure of the porous material well as to map the flow fields with sufficient accuracy, the simulation for materials with low porosity and heterogeneous pore structures requires a high resolution. Future work will focus on massively parallel simulations with 1000 3 voxels and >10 9 degrees of freedom. Of course, besides the physical relationships, computational aspects will also be of interest.