Applied Mathematical Modelling 35 (2011) 1666–1675

Contents lists available at ScienceDirect

Applied Mathematical Modelling

journal homepage: www.elsevier.com/locate/apm

SPH modelling of fluid at the grain level in a porous medium

G.G. Pereira ⇑, M. Prakash, P.W. Cleary
CSIRO Mathematics, Informatics and Statistics, Private Bag 33, Clayton South, Vic 3169, Australia

a r t i c l e i n f o a b s t r a c t

Article history: Computational modelling of the flow of fluids in porous media has traditionally been at a
Received 30 August 2010 macroscopic level where the medium’s permeability and porosity are an input (from
Accepted 17 September 2010 experiments for example). In many cases this is difficult, especially if the porous medium
Available online 26 September 2010
changes its solid structure as a function of time. This situation occurs in reactive systems
such as ‘‘heap-leaching”, where biological and/or chemical solutions are introduced into
Keywords: the heap to dissolve or react with valuable materials. In this case, modelling fluid flow at
Smoothed particle hydrodynamics
the grain level is paramount and we show how this can be done with the SPH technique.
Porous media
Permeability
We present three-dimensional SPH simulations of fluid flow in an idealised porous med-
ium and show that the technique yields flows which are physically realistic. The perme-
ability of the medium is then predicted.
Crown Copyright Ó 2010 Published by Elsevier Inc. All rights reserved.

1. Introduction

Porous media flow is important in a number of areas such as oil recovery, groundwater flow, carbon sequestration, flow in
biological materials and recovery of valuable minerals (or metals) from ores via the process of heap-leaching. In each case,
one is interested in how the structure of the porous medium affects the fluid flow. Generally these flows have low Reynolds
numbers [1], where Re  qLv =l, where q is the fluid density, L is a typical length scale for the flow, v is the typical fluid speed
and l is the fluid viscosity. For these low Re, the flow is in the laminar regime and its macroscopic behaviour can be de-
scribed by Darcy’s law:
k
v ¼ rP; ð1Þ
l
where P is the pressure and k is the permeability of the medium.
In this study we model fluid flow in a representative porous medium at the grain-scale. By doing this we would like to be
able to calculate the permeability without making any a priori assumptions. To date there are only two other published
works which have used SPH methods to determine the permeability of a porous medium from first principles. They are
the studies by Zhu et al. [2] and by Jiang et al. [3].
The study by Zhu et al. [2] considered a two-dimensional (2D) porous medium made up of circular, solid particles either
arranged in a square or hexagonal array. This corresponds to fluid flow around long fibres. They were able to compute the
discharge velocity as a function of the applied pressure gradient and found that it followed a linear relationship. From the
gradient of this linear relationship they evaluated the medium’s permeability. They compared their SPH simulations with
Finite Element calculations and found good agreement.
The second work by Jiang et al. [3], considered flow through a slit, with the slit placed at a sequence of angles with respect
to the applied pressure gradient. This angle was varied from zero (parallel to pressure gradient) to 90° (orthogonal to the
pressure gradient). Once again this was a two-dimensional study. They were able to compute the average fluid velocity

⇑ Corresponding author.
E-mail addresses: Gerald.Pereira@csiro.au (G.G. Pereira), Mahesh.Prakash@csiro.au (M. Prakash), Paul.Cleary@csiro.au (P.W. Cleary).

0307-904X/$ - see front matter Crown Copyright Ó 2010 Published by Elsevier Inc. All rights reserved.
doi:10.1016/j.apm.2010.09.043
 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675 1667

as a function of the applied pressure gradient and found that this also followed a linear relationship. From this they deter-
mined the medium’s permeability. These results were then compared to the Blake–Kozeny–Carman [1,4] theoretical predic-
tions and good agreement was found.
Our study continues the development of the SPH method by considering flow in a three-dimensional (3D) porous struc-
ture. It is well known that permeability of 3D media is quite different to 2D because the extra dimension results in a vast
increase in the number of possible flow paths. We note here that a current technique popular for obtaining permeabilities
for porous media flow is the Lattice–Boltzmann (LB) method. See for example Chen and Doolen [5] or Heijs and Lowe [6].

2. SPH method

The SPH method is reviewed in detail by Monaghan [7]. The method starts with the interpolation of any function A at any
position r using SPH smoothing which is given by:
X Ab
AðrÞ ¼ mb Wðr  rb ; hÞ; ð2Þ
b
qb
where mb and rb are the mass and density of particle b and the sum is over all particles b within a radius 2h of r. Here W(r, h)
is a C2 spline based interpolation or smoothing kernel with radius 2h that approximates the shape of a Gaussian function. The
gradient of the function A is given by differentiating the interpolation Eq. (2) to give:
X Ab
rAðrÞ ¼ mb rWðr  rb ; hÞ: ð3Þ
b
qb
Using these interpolation formulae and suitable Taylor series expansions for second order derivatives, any parabolic partial
differential equations can be converted into ordinary differential equations for the motion of the particles and the rates of
change of their properties.
The SPH continuity equation, taken from Monaghan [7,8] is:
dqa X
¼ mb ðv a  v b Þ  rW ab ; ð4Þ
dt b

where a is the density of particle a with velocity va and mb is the mass of particle b. We denote the position vector from par-
ticle b to particle a by rab = ra  rb and let Wab = W(rab, h) be the interpolation kernel with smoothing length h evaluated at
distance |rab|. This form of the continuity equation is Galilean invariant (since the positions and velocities appear only as dif-
ferences), has good numerical conservation properties and is not affected by free surfaces or density discontinuities.
The SPH method used here is quasi-compressible and an equation of state is used to specify the relationship between par-
ticle density and fluid pressure. A form suitable for incompressible fluids is:
 c 
q
P ¼ P0 1 ; ð5Þ
q0
where P0 is the magnitude of the pressure and q0 is the reference density [7,8]. For water and similar fluids c = 7.
The momentum equation, when converted to SPH form, becomes an acceleration for each particle a, i.e.,
X   
dv a Pb Pa n 4la lb v ab rab
¼g mb þ  rW ; ð6Þ
dt b
qb
2 q2a qa qb ðla þ lb Þ r2ab þ g2 a ab
where Pa and la are pressure and viscosity of particle a and vab = va  vb. Here n is a factor associated with the viscous term
[9], g is a small parameter used to smooth out the singularity at rab = 0 and g is gravity.
The pressure scale factor P0 is given by:
gP0
¼ 100V 2 ¼ c2 ; ð7Þ
q0
where V is the characteristic or maximum fluid velocity and c is the speed of sound. This means that the sound speed is ten
times the characteristic speed. This ensures that the density variation is less than 1% and that the flow is close to incompress-
ible. The computational speed of sound is much lower than the actual (physical) speed of sound so as to increase the numer-
ical time-step.
As two particles approach each other, their relative velocity is negative (as is the gradient of the kernel) so that there is a
positive contribution to dqa/dt. If this rate of change is positive then the density of particle a rises. This leads to a positive pres-
sure that pushes the particles apart again. If two particles move apart, then their densities decrease creating a negative pres-
sure that pulls the particles back towards each other. This interplay of velocity and density/pressure ensures that the particle
remains ‘‘on average” equally spaced and that the density is close to uniform, so that the fluid is close to incompressible.
To simulate confined or partially confined fluid flow, such as is typically found in industrial and geophysical fluid flows,
the modelling of physical boundaries is important. Boundaries in SPH can be modelled in a range of ways, but here we use
1668 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675

approximately evenly spaced SPH particles at which a Lennard-Jones type force (a very steep polynomial potential relating
the force with the displacement into the boundary) is applied in the normal direction. In the tangential direction, the bound-
ary particles are included in the summation for the shear stress using Eq. (6) (with the pressure gradient term set to zero) to
give non-slip boundary conditions on the walls.

3. SPH porous media configuration

In this work we seek to model flow around individual (pore-sized) grains. To do this we model the solid grain as a sphere
made up of stationary SPH particles. The sphere has diameter D and is placed in a cube of fluid with side-length L. Periodic

Fig. 1. Fluid and solid particle distribution at three different times: (a) t = 0 s, (b) at t = 0.2 s, and (c) t = 2.0 s for 2D simulation with L = 14 mm. Particles are
coloured by velocity from blue (stationary) to red (0.11 m/s).
 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675 1669

boundary conditions are imposed in all three directions at the edges of the fluid cube. The remainder of the region is filled
with fluid particles. The density of the fluid used is 1000 kg/m3 and the viscosity is 1 Pa s. In all the cases reported the solid
particle has a diameter of 4 mm. We estimate the Reynolds number for the flow to be of the order of 0.001–0.1, given the
average flow velocities are in the range of 0.001–0.1 m/s and a characteristic length scale of 1–10 mm. The important
length-scale here is the distance between particles, rather than particle size, since it is this distance which has the dominant
effect on the flow. No-slip boundary conditions are implemented at the solid sphere boundary. Gravity (or a constant body
force) is imposed downward along the z direction and this drives the fluid flow. This set-up was also implemented for 2D
flow, by replacing a solid circle (of diameter 4 mm) for a sphere and a square of fluid for the cube.
The two-dimensional simulations represent flow around a pack of cylinders, for example flow around long fibres while
the three-dimensional simulations represent flow around a periodic cubic array of spherical grains. Initially, simulations
were run in 2D with a number of different fluid particle sizes to test resolution dependence. We then monitored quantities
such as kinetic energy of particles, fluid density and particle velocities to determine the largest fluid particle size at which
these quantities did not show significant changes. By doing this we are able to achieve accurate results in a reasonable time.

4. Fluid flow results

We first consider 2D simulations primarily to decide on the resolution for further simulations. In all simulations reported
in this paper the speed of sound is set at 10 m/s. We use a computational speed of sound of 10 m/s, which is much lower than
the actual speed of sound for this medium. Because characteristic fluid velocities are very small (0.5 m/s) using such a low
speed of sound does not affect the accuracy, but greatly increases the computational efficiency. Fig. 1 shows the fluid velocity
at different times for a driving acceleration of 9.8 m/s2 downwards. The square has side length of L = 14 mm giving a solid
volume fraction of 0.064 and the fluid particle size is 0.25 mm. Initially all fluid particles are at rest and over time maximum
flow velocities appear in the gap (red coloured particles) between the solid particle and its periodic image while above the
solid particle the fluid velocities are small. The time taken to reach a stable velocity profile (which is determined by the com-
petition between the gravitational force and viscous drag) is 0.3 s.

Fig. 2. (a) Fluid and solid particle velocities for a 2D simulation (L = 14 mm) with fluid particle size of 0.1 mm. The maximum velocity (red) is 0.15 m/s. (b)
Average steady-state kinetic energy as a function of fluid particle size in 2D simulations.
1670 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675

Fig. 3. (a) Steady-state fluid velocity distribution in 3D for a solid fraction of 0.012, and (b) average flow velocity through two planes which are orthogonal
to the direction of gravity. The first plane goes through the centre of the solid grain while the second plane is just below the solid grain.

To make a suitable resolution study we ran a number of simulations for different fluid particle sizes from 0.28 mm down
to 0.1 mm. Fig. 2a shows the fluid velocity for the case of a resolution of 0.1 mm. The steady-state limit is around 0.13 m/s,
which is slightly larger than for the much coarser 0.25 mm resolution. To make a quantitative comparison of the effect of
fluid resolution we monitored the average kinetic energy over all the particles as a function of time. Here we observed that
by about 0.2 s this reaches a plateau value (steady-state conditions) which was approximately 0.55 mJ.
Fig. 2b shows the steady state value of the average kinetic energy for different resolutions. It is almost constant over the
entire range, so a resolution of 0.25 mm is suitable for giving sufficiently accurate results. However, it is important to note
that as the separation between solid grains decreases the fluid particle size also needs to decrease so that fluid particles can
access the decreasingly small crevices between solid grains. In all subsequent simulations we monitored relevant quantities
(such as average particle velocities, kinetic energies, etc.) to ensure resolution was adequate.

4.1. Three dimensional simulations

In 3D, we initially consider the case where the solid volume fraction is 0.012 (which is a small solid fraction but results
here are typical of all simulations). The fluid particle size for all simulations in this section is 0.25 mm. We impose a body
force (per unit mass) of 9.8 m/s2 downwards. The density variation throughout this simulation was less than 0.1% indicating
the fluid is virtually incompressible. Fig. 3a shows the particle and velocity distribution after the system has reached steady-
state conditions. One slice is displayed which is a fluid particle spacing thick in the x–z plane and located at the centre of the
periodic cell. The colour of particles represents their velocity from blue (stationary) to red (0.4 m/s).1 The flow is similar to
that described for the 2D case, with high velocities in the gap between solid grains and lower velocities just above and below
the solid grain.
For purposes of determining a medium’s permeability we are interested in the average steady-state fluid velocity. To cal-
culate this we consider a (virtual) plane which is perpendicular to the gravity vector and calculate the velocity, averaged over

1
For interpretation of colour in Figs. 1–3, 6 and 7, the reader is referred to the web version of this article.
 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675 1671

all particles passing through this plane at any particular instant, as a function of time. We chose two such planes, the first one
placed through the centre of the solid grain and the second one just below the solid grain. The velocity through these planes
is shown in Fig. 3b and tends to equilibrate after about 0.2 s. The light blue curve is for the plane that goes through the centre
of the solid grain while the pink curve corresponds to the plane just below the solid grain. Since fluid particles have a smaller
cross-sectional area to pass through when flowing around the solid grain, to conserve fluid incompressibility they must
move correspondingly faster. Hence the average velocity of particles passing through this plane is slightly larger than for
the plane placed just below the solid grain. The difference in average fluid velocity between these two planes is, however,
quite small and we can take an average of these quantities for the overall (average) fluid velocity. In Fig. 4 we show the
iso-surfaces corresponding to velocities in the z-direction for this case. The lowest velocity iso-surface, shown in Fig. 4a, clo-
sely follows the outer surface of the solid particle, since we have no-slip boundary conditions on this surface. The highest
velocity iso-surfaces are located in the corners of the periodic cube since these are furthest away from the particle (and hence
furthest away from any obstructions). The velocity iso-surface (in the mid velocity range) shows four-fold symmetry due to
the use of a cubic cell.

Fig. 4. Iso-surfaces (at steady state conditions) of the z-velocity for the 3D simulation at a solid fraction of 0.012 (same case as Fig. 3). (a) Iso-surface at
vz = 0.13 m/s. (b) Iso-surface at vz = 0.26 m/s. (c) Iso-surface at vz = 0.29 m/s. Note, the solid particle in the centre is not shown.
1672 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675

Simulations were also run for body forces per unit mass of 4.9, 2.45, 1.23 and 0.61 N/kg giving flow velocities in the same
way as was shown in Fig. 3a. Fig. 5a shows these steady-state flow velocities as a function of imposed body force. In addition
to the points from simulations we also know that with a zero body force there is no flow. Since the flow is laminar we expect
to find a linear relationship. The best fit line that goes through these points has gradient 0.034 s. Note that the permeability is
related to the imposed body force per unit mass, Fb, and average velocity, hui, via

l
k¼ hui: ð8Þ
qF b
Since our fluid viscosity is 1 Pa s and density is 1000 kg/m3 this yields a permeability of k = 3.4  105 m2.
Simulations were performed in 3D for a number of different fluid volumes, i.e. decrease the edge length of the periodic
cube in regular steps (of 2 mm) from L = 14 mm to 8 mm. In each case, a series of different body forces were used and the
average steady-state flow velocity was determined. These are plotted in Fig. 5a along with a line of best fit. The average stea-
dy-state velocity is clearly a linear function of the imposed body force, which is in agreement with Darcy’s Law (Eq. (1)).
Moreover, the gradient of the graph is related to the permeability of the medium and hence we have demonstrated that
we can determine, from first principles, the permeability of a 3D medium using SPH simulations. The numerical values of
the permeabilities derived from these simulations are plotted in Fig. 6.

4.2. Three dimensional simulations – large solid fraction

The solid volume fractions used so far have been quite small (less than 0.1) while in most porous media one would expect
this to be much larger. In the present set-up the maximum solid fraction we can reach is 0.5236 for the case where the
spheres touch. However, simulating at this level requires us to have a much finer particle resolution otherwise fluid particles
cannot access regions where solid grains touch. In each of the following simulations we decrease the fluid particle size until
we are sure we are obtaining accurate results. We do this by monitoring the average fluid velocities and kinetic energies and

Fig. 5. Average flow velocity after equilibrium steady-state conditions have been reached versus imposed body force per unit mass. (a) For 3D cases at
different solid volume fractions of 0.012 (circles), 0.019 (squares), 0.034 (diamonds), and 0.065 (triangles). (b) Average steady state flow velocity for 3D
simulation with solid volume fraction of 0.155 (circles), 0.393 (squares) and 0.5236 (triangles).
 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675 1673

Fig. 6. Comparison of permeability from the two methods described in this study as a function of solid fraction. Black circles correspond to SPH simulations.
The red curve corresponds to the theoretical predictions of Zick and Homsy [11]. In (a) we show the entire range of solid fractions simulated and in (b) we
zoom into the higher solid fraction region.

make sure they do not fluctuate appreciably (after reaching steady state). Thus for the simulations reported in this section
we use a fluid particle size of 0.1 mm.
Fig. 7 shows vertical slices of the flow field after steady state has been reached for the case with a driving force (per unit
mass) of 9.8 N/kg. We verified steady-state conditions had been achieved by monitoring the average fluid velocity and
observing that it had reached a plateau. Note, even though there are some red coloured particles near the solid boundary
in Fig. 7, the magnitude of velocities is very small (0.00114 m/s) and the no-slip condition is still satisfied. In Fig. 7a, which
corresponds to a vertical plane through the middle of the solid grain, fluid is mostly stationary (blue), since in this plane the
solid grains touch. Here fluid cannot move at any appreciable speed in the direction of the gravity vector. The next slice
(Fig. 7b) cuts through part of the solid grain and hence fluid particles in this plane still have a small velocity. Finally in
Fig. 7c we show a slice at the extremity of the solid grain. In this plane fluid has a comparatively unrestricted flow path
and so the fluid velocity is maximal, although still quite small (of the order of 0.001 m/s) compared to the previous cases
reported.
Once again we predicted the permeability for this solid fraction (by running simulations at a variety of body forces) and
found the permeability to be 1.3  107 m2 (see Fig. 5b). Hence, we have shown the SPH method is applicable to real porous
media with realistic solid fractions.

5. Theoretical predictions

At low solid fractions, the fluid flow becomes more like an infinite fluid flowing around an isolated solid grain. This is the
Stokes limit where the fluid drag is 3plDv. We are not quite at this infinite limit in our simulations even at low solid fraction.
For us the fluid drag on a solid grain is modified by the flow and pressure effects of the surrounding particles in the lattice.
The geometry of our simulations is that of a periodic array of spheres in a simple-cubic packing. Theoretical predictions from
the permeability for this geometrical construction have been worked out numerically by Hasimoto [10] at low solid fractions
and by Zick and Homsy [11] for high volume fractions. Note, that in this set-up the highest obtainable solid fraction is 0.5236.
The Blake–Kozeny–Carman relationships, which assume a porous medium is composed of a bunch of capillary tubes of equal
length [1], is appropriate for a random array of spheres rather than a regular/periodic array of spheres.
1674 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675

Fig. 7. Snap-shots (in the x–z plane) of the fluid and solid particle distribution at steady-state conditions at three different slices for maximum solid fraction
of 0.5236 and a body force of 9.8 m/s2: At (a) y = 0, (b) y = 1 mm and (c) y = 2 mm. The colour of particles represents their velocity from blue (stationary) to
red (large velocity). Note, even though there are some red-coloured particles near the solid boundary in (c) these velocities are very small and the no-slip
condition is still valid.

Hasimoto [10] obtained a perturbation solution (in terms of Fourier series) for Stokes flow of an incompressible viscous
fluid past an array of regularly ordered spheres. The expansion was made in terms of fractional powers of the solid fraction
and is thus limited to relatively small solid fractions. In fact Hasimoto [10] gave his results in terms of the fluid drag on each
sphere and wrote his results as

F drag ðfs Þ ¼ 3plDv C D ðfs Þ; ð9Þ

 G.G. Pereira et al. / Applied Mathematical Modelling 35 (2011) 1666–1675 1675

and CD is called the drag coefficient and is a function of the solid fraction fs. Thus CD is 1 in the limit of an infinitely large fluid
volume and increases as fs increases. To obtain a permeability from this we need to obtain the force per unit mass of fluid.
The mass of fluid in our simulations is pqD3(1  fs)/(6fs). Thus the fluid drag (per unit mass of fluid) is
18C D fs lv
F drag ¼ : ð10Þ
D2 ð1  fs Þ q
Comparing this equation with Eq. (8) gives us the permeability as
ð1  fs Þ
ktheory ¼ D2 : ð11Þ
18C D fs
Zick and Homsy [11] used a numerical method, based on reducing the Stokes partial differential equations to a set of integral
equations and then using a Galerkin method to solve these equations. Their results are also expressed in terms of a fluid drag
coefficient (Eq. (9)) and hence the permeability can also be expressed as Eq. (11), but with a different value of the drag coef-
ficient (compared to Hasimoto [10]). The numerical result of Zick and Homsy [11] is essentially the same as Hasimto’s [10]
result for solid fractions less than 0.065, but shows increasing deviations at higher solid fractions.
This permeability is a monotonic decreasing function of solid fraction and for small solid fractions is inversely propor-
tional to fs. Permeabilities calculated by this method are compared with SPH in Fig. 6. We see good overall agreement
(Fig. 6a). Zooming into the region of higher solid fraction we see there is a small deviation around solid fractions of 0.1.

5.1. Comparisons between theoretical predictions and SPH

We have been able to estimate the permeability of idealized three-dimensional porous media from first principles using
the SPH method. In principle this method can be extended to a medium with many more grains although, clearly, the method
is limited by computational considerations. Looking at Fig. 4 it is clear the average fluid velocity is a linear function of the
imposed body force. This is as one would expect in the laminar case, i.e. it follows Darcy’s Law. This is reassuring and tends to
lend support to the validity of the results of the SPH method. All predictions agree qualitatively in that the permeability de-
creases with increasing solid fractions. The overall agreement between SPH permeability predictions and the theoretical pre-
dictions is relatively good and validates the implementation of SPH for obtaining permeabilities in simple porous media
geometries.

6. Conclusions

We have modelled flow through an idealized three-dimensional porous medium using the SPH technique. As this is a first
principle type study we have used a very simple model of the porous medium – flow around a single solid particle in a cubic
cell with periodic boundary conditions at the edges. We have found the flow rate is a linear function of the imposed body
force, which obeys Darcy’s law for laminar flow, and from this we have been able to estimate the medium’s permeability.
The SPH estimates of the permeability, are in qualitative and relatively good quantitative agreement with theoretical calcu-
lations based on solving Stokes equation in the presence of a periodic array of solid spheres.
We have made no direct comparison insofar as computational efficiency of our SPH method versus LB methods for eval-
uating the permeability of this porous medium. This is because we have not optimized our method and we envisage the SPH
method as an alternative to LB methods, especially in cases where one has irregularly shaped particles or dissolving particles.
The present SPH method can easily be used for solid particles of arbitrary shape and size. All that requires to be done is to
construct digital images of the solid particle surface which would then be imported into our SPH code. As such, we can obtain
the local permeability of an arbitrarily complex porous medium. Another important development of this procedure would be
to use in cases where the solid particles change in shape and size. This is particularly relevant to applications such as heap-
leaching. In principle, the modelling of a liquid which flows through a porous medium which dynamically changes due to the
flow of the liquid can be implemented with an SPH technique. In fact SPH would be highly suitable for this process because
each fluid and/or solid particle (individually) can carry additional information such as its acid concentration or whether it is
moving.

References

[1] J. Bear, Dynamics of Fluids in Porous Media, first ed., American Elsevier Publishing, New York, 1972.
[2] Y. Zhu, P.J. Fox, J.P. Morris, A pore-scale numerical model for flow through porous media, Int. J. Numer. Anal. Methods Geomech. 23 (1999) 881–904.
[3] F.M. Jiang, M.S.A. Oliveiria, A.C.M. Sousa, Mesoscale SPH modelling of fluid flow in isotropic media, Comput. Phys. Commun. 176 (2007) 471–480.
[4] T.E. Faber, Fluid Dynamics for Physicists, first ed., Cambridge University Press, Cambridge, 1995.
[5] S. Chen, G.D. Doolen, Lattice Boltzmann methods for fluid flow, Ann. Rev. Fluid Mech. 30 (1998) 329–364.
[6] A.W.J. Heijs, C.P. Lowe, Numerical evaluation of permeability and Kozeny constant for two types of porous media, Phys. Rev. E 51 (1995) 4346–4352.
[7] J.J. Monaghan, Smoothed particle hydrodynamics, Rep. Prog. Phys. 68 (2005) 1703–1759.
[8] J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994) 399–406.
[9] P.W. Cleary, Modelling confined multi-material heat and mass flows using SPH, Appl. Math. Model. 22 (1998).
[10] H. Hasimoto, On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres, J. Fluid
Mech. 5 (1959) 317–328.
[11] A.A. Zick, G.M. Homsy, Stokes flow through periodic arrays of spheres, J. Fluid Mech. 115 (1982) 13–26.