7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Direct Numerical Simulation of Saturated Deformable Fibrous Porous Media
Irfan Khan and Cyrus K. Aidun
G.W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA
irfan.khan@gatech.edu and cyrus.aidun@me.gatech.edu
Keywords: Saturated deformable porous media, Direct numerical simulation, Fluidstructure interaction,
LatticeBoltzmann method, FiniteElement method
Abstract
Numerical techniques for modeling saturated deformable porous media have mainly been based on mixture theory
or homogenization techniques. However these techniques rely on phenomenological relationships for the constitutive
equations along with assumptions of homogeneous and isotopic material properties to obtain closure. Direct numerical
simulations of the multiphasic problem for flow in deformable porous media avoid such assumptions and thus can
provide significantly accurate understanding of the physics involved.
In this work, a hybrid method using Lattice Boltzmann Method (LBM) for fluid phase and Finite Element Method
(FEM) for solid phase is used for direct numerical simulation of saturated deformable porous media. The method
provides a number of unique features including scalability on distributed computing necessary for such a problem.
Additionally, the method has been used to understand the behaviour of model saturated porous media under
compressive loading conditions. In this paper, as a first approximation to fibrous porous media, regular arrangements
of cylinders are used as the model porous media. The compressional behaviour of the model fibrous porous media
match well with existing literature on saturated porous media. Thus it is found that macroscopic behaviour of a
generic porous media can be recovered with the use of model porous media constructed with simplified geometries.
As a result it is also seen that the LBMFEM method is an accurate and suitable method for further investigations on
deformable porous media.
Introduction
In the field of saturated deformable porous media the
poroelasticity theory developed by Biot (1941, 1956)
has been successfully used to predict the behaviour of
homogeneous isotropic porous media, such as soil, un
dergoing small deformation. A mathematically rigor
ous theory of porous media (TPM) developed by Bowen
(1980, 1982), de Boer R. (2005), Ehlers (1989); Ehlers
and Eipper (1999); Ehlers and Markert (2003), Bluhm
and de Boer (1997) has been used to model both small
and large deformation of homogeneous porous media.
Inspite of the rigour and success, these methods make
assumptions about the constitutive relationships, par
ticularly those that govern the behaviour of permeabil
ity and bulk modulus of the media. The assumption
of a smeared media used in these techniques requires
material properties of the bulk media, which are de
termined through either additional experimental work
or using phenomenological relations. Further, as men
tioned above these techniques are still limited to homo
generous media. However, in textile industry, there is
a need to understand the behaviour of nonhomogeneous
and anisotropic saturated deformable fibrous porous me
dia under various loading conditions.
Using direct numerical simulations (DNS), one can
capture all the pertinent physics without involving open
ended parameters. Due to the comprehensive nature of
DNS, there is no need for assuming the constitutive re
lationships. Thus providing significantly accurate infor
mation of any complex phenomena. DNS serve as a tool
to investigate the constitutive relationships in complex
geometries. DNS can be used to validate of existing
mixture theory models and determine their limitations.
Direct numerical simulations have been used to to in
vestigate fluid conductance in rigid fibrous porous media
(Sangani and Acrivos (1982a,b), Koch and Ladd (1997),
Sobera and Kleijn (2006), Zobel et al. (2007), Tahir and
Tafreshi (2009)). Boutt et al. (2007) used a discrete
element method coupled with LBM to simulate satu
rated porous media comprising of rigid movable parti
cles. However, the authors could not find any litera
ture or data on the use of DNS for understanding the
behaviour of saturated deformable porous media
A fully coupled and parallelized numerical method
using lattice Boltzmann method for fluid phase and finite
element method for solid phase is used for analysis of
saturated deformable porous media through DNS (Khan
and Aidun (2010b),Khan and Aidun (2010a). The lattice
Boltzmann method with its simplicity and local nature of
calculations is suitable to model complex flows (Aidun
and Clausen (2010)) including porous media. The finite
element method has been proven as a robust technique
for modeling elastic deformations.
In this study, the feasibility of the LBMFEM method
for direct numerical simulation of saturated deformable
fibrous porous media is evaluated. As an approximation,
regular arrangements of cylinders are used to develop
model fibrous porous media. In the next section, a brief
description of the methodology is provided including the
validation of the implemented contact model with Hertz
analytical solution for deformation of spherical surfaces
during contact. In the subsequent section the description
of the model porous media geometry, calculation of its
bulk properties and comparison with analytical models
are provided.
Method for modeling fluidstructure
interactions in complex geometries
In this work a parallel hybrid Lattice Boltzmann and fi
nite element method is used to model the fluidstructure
interaction in saturated porous media during deforma
tion. Details of the method and its validation can be
found in Khan and Aidun (2010b). Here only a brief
description of the method is provided.
The single relaxation D3Q19 implementation (Aidun
et al. (1998); Aidun and Lu (1995)) of lattice Boltzmann
method is used to model the fluid phase. The accuracy of
LBM has been compared to traditional numerical meth
ods in fluid dynamics such as finite difference and fi
nite volume methods (Bernsdorf et al. (1999); Sankara
narayanan et al. (2003)). Moreover, the local nature of
the calculations and the simple bounceback scheme for
noslip condition make LBM a favorable tool for model
ing singlephase flow through porous media (Noble et al.
(1996)), where the inherent parallel nature of LBM can
be leveraged through distributed computing.
In order to model the solid deformations, a linear elas
tic finite element method is used. The fournode tetrahe
dral element is used to discretize the geometry and solve
the weak form of the Cauchy's equation. The New
mark's scheme is used to perform time integration.
The transfer of momentum between the fluid and solid
phases is effected through the popular "linkbounce
back" scheme (Ladd (1994), Aidun and Lu (1995);
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Aidun et al. (1998)).
A near contact model developed by Ding and Aidun
(2003) and later modified for elastic finite element par
ticles by MacMeccan et al. (2009) is used in this work.
A linkwise repulsive force, based on the lubrication and
contact forces, is applied to the interacting surfaces. The
magnitude and direction of the linkwise force is given
as
6Fk =
q (1
+Acexp
 ^)Ua ek
Ua ek
Ck
cTc
if g > g ,
c < g < Ck ,
0 < g < gc,
(1)
where Ua is the approach velocity of the interacting sur
faces, ek is the direction vector of the lattice link and
V', pF are the kinematic viscosity and density of the
fluid. ce is the length of a link, g is the linkwise gap be
tween surfaces, gc is the contact cutoff distance, A is the
local surface curvature and q is chosen to be 0.6 (Ding
and Aidun (2003)). The parameters gc and ac are depen
dent on the surface roughness of the geometry and are
determined a priori. A, is chosen such that the repul
sive contact force scales appropriately with the applied
stresses (a) i.e. A, z where ao is the contact area
a,
that varies with the applied stress. For nonconformal,
convex bodies, ao can be determined based on a and
surface curvature.
The contact model has been tested by analyzing the
deformation when two frictionless bodies are brought
into contact under an external force. Hertz solved the
problem within elastic limit (Timoshenko and Goodier
(1970)). The analytical expression for the relative dis
placement (a) of the centers of two elastic bodies with
convex surfaces under a compressive force is given as
a 9p2 8
[8R*E*2]
where P is the applied load, E* = ( where E is
the Young's modulus and v is the Poisson ratio. I is
equal to for spheres and R for cylinders; R being the
radius of curvature of the convex surface.
In the first case, two spheres of equal radii R are
aligned along the xaxis and a compressive load P is
applied. The diameter of the spheres was chosen to be
40 lattice units, which allows for resolving the transfer
of fluidstructure interaction forces and contact forces
accurately. Two different mesh resolutions were used
to discretize the geometry of spheres. The finite element
edge length of the coarse mesh is approximately 4 lattice
units which leads to approximately 10 elements across
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
a = 2R1
4p P
Load (P/HR2)
Figure 1: Comparison of the Hertz contact model
with the near contact model implemented in LBMFEM
method, a represents the deformation due to the com
pressive load P and R is the radius of the spheres. The
solid line represents the analytical model by Hertz and
the points represent the LBMFEM simulation results
using coarse mesh (x) and fine mesh (E).
the diameter of the sphere. The fine mesh has twice the
resolution of the coarse mesh, thus it has 20 elements
across the diameter of the sphere. In this study, since
the contact parameter Ac scales with the applied stress
load a = F the value of Ac increases as the load P
is increased.
Figure shows the comparison of nondimensional de
formation a/R between the implemented contact model
in LBMFEM with Hertz's contact model. As can be
seen, the predicted deformation matches well with the
analytical model, particularly near small deformation
(A < 0.04). The mesh resolution does not have a ma
jor effect on deformation when  < 0.04. However
for a > 0.04 the fine mesh gives better results for
the predicted deformation. Further, the simulation re
sults deviate from the analytical model as the deforma
tion increases. This deviation, also seen in finely meshed
spheres, is mainly due to lack of enough mesh resolu
tion. However it is also well known that for finite de
formations, Hertz's model overpredicts the value of a
Dintwa et al. (2008).
The second geometry used are cylinders of radii
R with o(hIh:in,.ill\ aligned axis. The cylinders are
brought into contact under a compressive force P. The
diameter of the cylinders is chosen to be 40 lattice units
to match those of the spheres. As in the case of spheres,
two mesh resolutions were tested. The coarse mesh has
10 elements across the diameter, while the fine mesh has
20 elements across the diameter of the cylinder. Fig
ure shows the comparison with Hertz's contact model.
A similar behaviour is seen, as in the case of spheres.
Figure 2: Comparison of the Hertz contact model
with the near contact model implemented in LBMFEM
method, a represents the deformation due to the com
pressive load P and R is the radius of the cylinders. The
solid line represents the analytical model by Hertz and
the points represent the LBMFEM simulation results
using coarse mesh (x) and fine mesh (E).
There is good match with the Hertz model for a < 0.04
and here the coarse and fine mesh results do not dif
fer a lot. The results deviate from the Hertz model as
S>= 0.04, again the fine mesh gives better results for
predicted deformation.
Ordered Orthogonally aligned Fibrous Porous
Media
de Boer et al. (1993) used the TPM to study wave prop
agation in an incompressible fluidsaturated porous me
dia. They solved the problems of onedimensional com
pression of a semiinfinite saturated porous media with
no body force terms. They assumed the strains to be in
finitisimally small, thus neglecting any variation in vol
ume fractions.
Subsequently, Diabels and Ehlers (1996) have used
the analytical solution to validate their nonlinear fi
nite element model based on the TPM to simulate bi
nary saturated porous media. More recently Hu et al.
(2009) used the analytical solution to validate their bi
nary porous media model based on Differential Quadra
ture Method. In the following the analytical solution
proposed by de Boer et al. (1993) is used to compared
with the results of LBMFEM direct numerical simula
tions using model porous media.
The problem of instantaneous compression of a semi
infinite medium as discussed above is duplicated using a
model porous media. The porous media is approximated
using uniform cylinders in an ordered layered arrange
ment where the axis of cylinder in one layer is orthog
a=2R1
Load (P/IR2)
i i h
(0,0,0)
Figure 3: 8 cylinders stacked vertically in a column of
fluid, representing an ordered orthogonal arrangement.
The dotted line represents the simulation domain.
onal to the axis of cylinders in the next layer. Figure
6 shows the ordered orthogonal arrangement. In such
an arrangement the geometric symmetry can be used to
simplify the simulation geometry. Thus, a column of
cylinders arranged vertically as shown in figure 3 is used
to represent the ordered orthogonal arrangement. The
column of cylinders are submerged in fluid, thus repre
senting a saturated porous media. This geometry rep
resents a first approximation of a more complex fibrous
porous media which form nonwovens.
The boundary condition at the bottom (z 0 ) is a
noslip, nodisplacement condition given as
v 0; u 0; at z =0; (3)
where v is the fluid velocity vector, u is solid displace
ment vector. At the top, an outflow condition is applied
on the fluid phase, and a load f(t) is applied on the solid
phase; given as,
9Vy 0v9
0; F(t) f(t)k; at z =h,
(4)
where F(t) is the load vector and h is the height of the
column. In this paper we only discuss the response of
using a constant load with time. On the four vertical
wall, symmetry conditions are applied, given as
v, 0; u, = 0; at x 0, 4R; (5)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
= 0; ,. 0; at y = 0, 4R. (6)
The initial conditions are given as
v(x,0) 0;
u(x, 0) 0; (8)
where x is the position vector.
The important nondimensional parameters that can
be used to define the problem of instantaneous compres
sion of a deformable porous media are p* the ratio of
solid to fluid densities ( ), n* the ratio of solid to fluid
volume fraction (), K* the fluid conductance or non
KFyLF /FP
dimensional permeability given as g where K
is the permeability defined in the next section, Dp is the
average diameter of particles that make up the porous
media and g is the acceleration due to gravity. Included
also are the Reynolds number (Re) and capillary number
(Ca) given as
UjF D
Re UavgDp Ca
1" .,p ,Ca
EaVF
Eav9
where UF av is the average flow velocity in the porous
media and vF is the fluid kinematic viscosity, 9 is the
strain rate and Eavg is the average bulk modulus of the
media.
In the following, the authors compare the deforma
tional behaviour of the model porous geometry through
DNS to the analytical model proposed by de Boer et al.
Since the analytical model is based on mixture theory,
which assumes the porous media to be a binary mixture
of superimposed but immiscible constituents, the model
makes use of bulk properties of the media. Two impor
tant bulk properties that appear in the analytical model
are the permeability (KF) and average bulk elastic mod
ulus (Eav,). In the following we describe the calculation
of these bulk properties for the chosen model porous me
dia. It should be noted that these parameters are evalu
ated only to compare the LBMFEM direct numerical
simulation results to the mixture theory based analytical
model. These parameters are not required and have not
been used in the direct numerical simulations.
Permeability of ordered orthogonally crossing
fibrous arrangement
As discussed in the section 1, there have been a lot of
studies, numerical, experimental and analytical in under
standing and developing models to predict permeability
of a porous media made up of cylindrical arrangements.
Most of these studies have been conducted to investi
gate flows in very low to low solid volume fractions
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1
(0,0,0) z
Figure 4: Setup and dimensions of the geometry used
to determine permeability of the ordered orthogonally
placed cylindrical arrangements. The color contours
represents pressure varation at steady state.
9e05
tice Boltzmann method has been used to simulate fluid
flow. The method has already been validated for accu
rate prediction of pressure drop in flow through porous
media Khan and Aidun (2010b).
The 8 cylinder arrangement as shown in figure 4 is
used to determine the permeability. Each cylinder has a
diameter of 40 lattice units and length of 80 lattice units.
A fluid domain of 720x80x80 lattice units is used in
these calculations. This length of the domain provides
enough space before and after the geometry to eliminate
entrance and outlet to effect the pressure drop calcula
tions. The fluid flows in the xdirection, the boundary
conditions are;
v = Vi; at x = 0;
9vy; 9vy O vz
9Ox Ox 9xv
0; at x = 1, (12)
0; at y = 0, 4R;
le05
00005 0001 00015 0002 00025 0003 00035
0004 0004E
Figure 5: Variation of pressure drop ( 1j) with aver
age flow velocity Vag. The slope of the line gives the
Darcy permeability (kF).
(ns < 0.3). However, the geometrical arrangement that
is intended to be used in these simulations has a solid
volume fraction of ns 0.3927. The only investiga
tions using orthogonally crossing cylinderical arrange
ments in high solid volume fraction geometries has been
performed by Sobera and Kleijn (2006). However, the
geometry used by them is a layered cylindrical arrange
ment where each layer consists of overlapping cylinders.
Further these layers are separated from each other by
cylinder half distance. Thus due to the lack of informa
tion available for permeability in orthogonally crossing
cylindrical arrangements where the flow is normal to the
cylinder axis, a permeability study has been carried out.
The Darcy's law relates the average fluid velocity to
the pressure gradient in the fluid and is given as
K
(v) VP, (10)
ft
where (v) is the average flow velocity, p is the fluid vis
cosity, VP is the resulting pressure gradient and K is
the hydraulic permeability tensor. For a homogeneous
isotropic system, K = kFI where I is the identity ten
sor.
In order to obtain the value of kF for the cylindri
cal arrangement, numerical simulations of fluid flow
through the arrangement have been performed. The lat
v, = 0; at z = 0, 4R; (14)
The inlet velocity is varied and the resulting pressure
drop is measured. In each case, the simulations are run
till a steady state distribution of pressure and velocity
values are attained throughout the domain.
Figure 5 shows the variation of pressure (1 dP) with
P1 dx
average flow velocity Vavg. The inverse of the average
slope of the line gives the Darcy permeability (kF). Thus
the average Darcy permeability from these simulations
is given as 4.94 x 10 7m2. The corresponding perme
ability KF is 1.2124 x 10 3m/s (table 1). This is the
value of permeability that is used in the analytical model
for comparison with LBMFEM simulations.
Bulk elastic properties of ordered orthogonally
crossing fibrous arrangement
Detailed derivation of the average bulk elastic modulus
of the cylindrical arrangement is provided in Khan and
Aidun (2010a). Here, only a brief description is pro
vided.
The simple ordered e(IIIh:i'n.Iill\ crossing homoge
neous and isotropic cylindrical arrangement is as shown
in figure 6. The two layers of cylinders are placed such
that the axis of one layer is orthogonal to the axis of the
other layer. Due to symmetry, the dotted region can be
taken as a unit cell which has the same elastic character
istics as the entire arrangement.
Following the work of Deresiewicz (1958), the instan
taneous bulk modulus of cylinders in an ordered orthog
onal arrangement is as
E*a
Eb 8
8R
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
I0
y
I0
1;
z
1;
Figure 6: Ordered orthogonally aligned axis, arrange
ment of two layers of cylinders. Due to symmetry, the
dotted region can be taken as a unit cell which has the
same elastic characteristics as the entire arrangement.
Table 1: Values of paramters used for comparing defor
mation of cylindrical arrangements with de Boer et al.
model
E
Eav9
noF
ns
pFR
PSR
As
Ps
KF
lFR
f(t)
46.08MPa
0.8973MPa
0.6164
0.3836
1000 kg/m3
2000 kg/m3
0.4803 MN/m2
11 :;I7IMN/m2
0.000738 m/s
10.0 KN/m2
0.01677 KN/m2
An average bulk modulus can be obtained by integrat
ing the above equation over the applied load. Thus the
average bulk modulus of the system is
3E*a, (16)
Ea g 32R (16)
Ear9 32R
where a, is the contact radius due to the applied load Po.
Table 1 lists the values of E, Ea, and Lame's constants
As, s.
Deformation of ordered orthogonally aligned
cylindrical arrangement
The analytical solution provide by de Boer et al. (1993)
gives the dynamic response of a semiinfinite saturated
porous medium under time dependent external load.
Since real geometries have the limitation of finite length,
geometries of three different lengths are chosen for com
0 20 40 60 80 100 120
time (tv/R2)
Figure 7: The average displacement of the top surface
nodes of the cylindrical arrangement with 8, 16 and
32 cylinders using LBMFEM direct numerical simu
lations. The displacement response is compared to the
analytical model proposed by de Boer et al.
prison with de Boer et al. model. Onedimensional or
dered orthogonal cylindrical arrangement of 8,16 and
32 cylinders are used in this study. All the arrange
ments contain cylinders of equal diameter and same ma
terial properties to model an isotropic and homogeneous
porous media.
The values of contact parameters correspond to those
chosen for a similar load to validate with Hertz model.
Thus no openended parameters are involved in this
study.
Further the effect of finite element mesh resolution on
the deformation response has been studied (Khan and
Aidun (2010a)). A coarse and a fine mesh have been
used to discretise the geometries. The coarse mesh has
10 elements across the diameter of a single geometry
(cylinder). While the fine mesh has 20 elements, twice
as many, across the diameter. It is found that the differ
ence in the deformational response of the coarse and fine
finite element mesh geometries is less than 3%. Thus in
the results presented here, all the simulations have been
run with the geometry where 10 elements have been
used across the cylinder diameter.
The values of nondimensional parameters are chosen
to be Re, Ca < 1, p* = 2.0, n* = 0.95 and K* =
0.003. The detailed list of all the parameters and their
values in SI units are listed in table 1.
Figure 7 shows the results of using a coarse mesh. The
average deformation responses of top surface nodes for
8, 16 and 32 cylinder geometries under impulse loading
using LBMFEM method are compared to de Boer et
al. model and found to match closely. An asymptotic
behaviour is seen in the deformation that varies with the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
LBMFEM(8cylinders) ..
LBMFEM(1cylnders)  
LBMFEM(32cylinders) **...
de Boer model 
de Boer et al 
LBMFEM(32 cyl) 
LBMFEM(16 cyl) 
LBMFEM(8 cyl) *1,
02
02 I
20 40 60 80 100 120
Figure 8: Comparison of the velocity response of the
model porous geometries made up of cylinder with the
analytical model.
domain size. The asymptotic behaviour marks the defor
mation reaching a steady state, when the internal stress
in the geometry balances the applied load. Concordant
to intuition, the onset of steady state is delayed with in
creasing domain size. This is because the deformation
corresponding to a given bulk strain is higher with in
creasing domain size. The 8 cylinder domain reaches
steady state quickly followed by the 16 cylinder geom
etry, while the 32 cylinder geometry does not reach a
steady state within the simulated time.
Figure 8 shows the average velocity response of the
top surface of nodes for 8, 16 and 32 cylinder geometries
and its comparison with de Boer et al. analytical model.
The velocity is large at the beginning of the compres
sion followed by a rapid decrease and then an asymp
totic behaviour to zero. This suggests rapid deformation
during the early part of compression followed by a slow
deformation phase. Again it is seen that the 8 cylinder
geometry seizes to deform first followed by 16 cylinder
geometry, whereas the 32 cylinder geometry follows de
Boer et al. model within the simulated time. The small
fluctuations in the velocity response are caused due to
minor instabilities arising from the contact model.
Figure 9 shows the normalized displacements along
the zaxis from the LBMFEM simulations of 8, 16 and
32 cylinder geometries, along with the de Boer et al.
model data. The zaxis is positive in the downward di
rection, such that z 0 indicates the top surface of
the geometries. The displacement values have been col
lected at time 20. At this time, the 8 cylinder
geometry has reached the final deformation state. Thus,
the 8 cylinder curve shows a linear behaviour, indict
ing the constant strain throughout the geometry. The 16
cylinder geometry also shows a similar behviour. How
04 035 03 025 02 015 01 005 0
Displacement (5z/R)
Figure 9: Comparison of displacement along the zaxis
from LBMFEM simulations of SC arrangement and the
de Boer et al. model at = 20. z 0 represents the
top surface and z is positive downwards.
ever, the 32 cylinder geometries has not reached the fi
nal deformation, hence its plot does not show a linear
behaviour. Further, the 32 cylinder geometry matches
with the analytical solution closely. The match is better
at the top and deviates towards the bottom of the geome
try. This shows that as the size of the geometry increases
there is a better match with the analytical solution.
Conclusions
A hybrid LBMFEM scheme is used to carry out direct
numerical simulations for the problem of instantaneous
compressive loading of a model saturated deformable
porous media. The contact model in the LBMFEM im
plemention has been tested with the Hertz model using
two different geometries i.e. spheres and cylinders. The
contact parameters are varied to obtain a good match
with the Hertz model.
Model porous media, approximated using regular ar
rangements of cylinders, are used to study the deforma
tion behaviour. Ordered o(i ilhi .'i n.ll. crossing cylindri
cal arrangement has been used with appropriate symme
try boundary conditions. Three different sizes of geome
tries involving 8, 16 and 32 cylinders were used. The
displacement and velocity responses are compared with
the analytical model based on the Theory of Porous Me
dia presented by de Boer et al. (1993).
To perform a comparison between the analytical
model and direct numerical simulations of model geom
etry, a correspondence was established between the bulk
elastic properties of the macroscopic porous media and
the single phase elastic properties of the solid. The per
meability of the media is also calculated for the compar
002
ison, through numerical simulations.
The LBMFEM simulations of compression of model
porous media compares well with the analytical model
proposed by de Boer et al. Based on the length of the ge
ometry the time for deformation to reach constant value
varies; being shorter for small length and longer for
larger length geometries. The velocity response is also
seen to match with the analytical model closely. Thus
as the length of the geometry increases, the match of the
response with the analytical model for semiinfinite ge
ometry, improves.
This study shows that the macroscopic behaviour of
a generic porous media can be recovered using regular
arrangements of simple geometries as long as the poros
ity, permeability and the bulk modulus is matched be
tween the porous media and the simplied geometry. The
study also shows that the average bulk elastic modulus
calculated based on small strain theory and Hertz contact
model provides a good estimate of the elastic property of
the model porous media in small strain limit.
The study also validates the accuracy and suitabil
ity of the hybrid LBMFEM method in modeling fluid
structure interaction in complex geometries like porous
media. It highlights the advantages of using this method
to understand the behaviour of real porous geometries,
as there is no need to calculate the bulk modulus, the per
meability or the porosity of the geometry. Further, using
the method, the study can be extended to the less in
vestigate nonhomogeneous and anisotropic deformable
porous media.
Acknowledgements
Simulation results presented in this paper make use of
computational resources provided trough TERAGRID.
Specifically the Steele cluster at purdue and Abe clus
ter at NCSA. Funding for I.K is provided through the
Institute of Paper Science and Engineering scholarship
program.
References
C.K. Aidun and J.R. Clausen. The latticeboltzmann
method for complex flows. Annual Review of Fluid Me
chanics, 42:439472, 2010.
C.K Aidun and Y. Lu. Lattice boltzmann simulation of
solid particles suspended in fluid. Journal of Statistical
Physics, 1:4961, 1995.
C.K Aidun, Y. Lu, and E.J Ding. Direct analysis of par
ticulate suspensions with inertia using the discrete boltz
mann equation. Journal of Fluid Mechanics, 373:287
311, 1998.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
J. Bernsdorf, F. Durst, and M. Schafer. Comparison of
celullar automata and finite volume techniques for sim
ulation of incompressible flows in complex geometries.
International Journal of Numerical Methods and Fluids,
29:251264, 1999.
M.A. Biot. General theory of three dimensional consol
idation. Journal of Applied Physics, 12:155164, 1941.
M.A. Biot. Theory of deformation of a porous viscoelas
tic anisotropic solid. Journal of Applied Physics, 27:
459467, 1956.
J. Bluhm and R. de Boer. The volume fraction concept
in the porous media theory. ZeitschriftfA r angewandte
mathematik und mechanik, 77:563, 1997.
D.F. Boutt, B.K. Cook, B.J.O.L. McPherson, and J.R.
Williams. Direct simulation of fluidsolid mechanics
of porous media using the discrete element and lattice
boltzmann methods. Journal of Geophysical Research,
112:B10209, 2007.
R.M. Bowen. Incompressible porous media models by
use of the theory of mixture. International Journal of
Engineering Science, 18:11291148, 1980.
R.M. Bowen. Compressible porous media models by
use of the theory of mixture. International Journal of
Engineering Science, 20(6):697735, 1982.
R. de Boer, W. Ehlers, and Z. Liu. Onedimensional
transient wave propagation in fluidsaturated incom
pressible porous media. Archive of Applied Mechanics,
63:5972, 1993.
de Boer R. Theoretical poroelasticity a new approach.
Chaos Soliton Fract, 25:861878, 2005.
H. Deresiewicz. Stressstrain relations for a simple
model of granular medium. Journal of Applied Mechan
ics, 25:402406, 1958.
S. Diabels and W. Ehlers. Dynamic analysis of fully sat
urated porous medium accounting for material and geo
metrical nonlinearities. International Journal ofNumer
ical Methods in Engineering, 39:8197, 1996.
E.J. Ding and C.K. Aidun. Extension of the lattice
boltzmann method for direct simulation of suspended
particles near contact. Journal of Statistical Physics, 112
(3):685708, 2003.
E. Dintwa, E. Tijskens, and H. Ramon. On the accuracy
of the hertz model to describe the normal contact of soft
elastic spheres. Granular Matter, 10:209221, 2008.
W. Ehlers. On thermodynamics of elastoplastic media.
Archives of Mechanics, 41(1):7393, 1989.
W. Ehlers and G. Eipper. Finite elastic deformations in
liquidsaturated and empty porous solids. Transport in
Porous Media, 34:179191, 1999.
W. Ehlers and B. Markert. On the viscoelastic behavior
of fluid saturated porous material. Granular Matter, 2:
153161, 2003.
Y. Hu, Y. Zhu, and C. Cheng. Dqm for dynamic response
of fluidsaturated viscoelastic porous media. Interna
tion Journal of Solids and Structures, 46:16671675,
2009.
I. Khan and C.K. Aidun. Modeling the macroscopic
behaviour of saturated deformable porous media using
direct numerical simulations. International journal of
multiphaseflow, page (In review), 2010a.
I. Khan and C.K. Aidun. Direct numerical simulation
of saturated deformable porous media using a parallel
hybrid latticeboltzmann finite element method. Inter
national journal of numerical methods in Engineering,
page (In review), 2010b.
D.L. Koch and A.J.C Ladd. Moderate reynolds num
ber flows through periodic and random arrays of aligned
cylinders. Journal offluid mechanics, 349:3166, 1997.
A.J.C Ladd. Numerical simulations of particulate sus
pensions via a discretized boltzmann equation. part 1.
theoritical foundation. Journal of Fluid Mechanics, 271:
285309, 1994.
R.M MacMeccan, J.R. Clausen, G.P. Neitzel, and C.K.
Aidun. Simulating deformable particle suspensions
using a coupled latticeboltzmann and finiteelement
method. Journal of Fluid Mechanics, 618:1339, 2009.
D.R. Noble, G.J. Georgiadis, and R.0. Buckius. Compar
ison of accuracy and performance for lattice boltzmann
and finite difference simulations of steady viscous flow.
Internation Journal ofNumerical Methods in Fluids, 23:
118, 1996.
A.S. Sangani and A. Acrivos. Slow flow through a pe
riodic array of spheres. International journal of multi
phase flow, 8(4):343360, 1982a.
A.S. Sangani and A. Acrivos. Slow flow past periodic
array of cylinders with application to heattransfer. In
ternational journal of multiphase flow, 8(3):193206,
1982b.
K. Sankaranarayanan, I.G. Kevrekidis, S. Sundaresan,
J. Lu, and G. Tryggvason. A comparative study of lattice
boltzmann and fronttracking finitedifference methods
for bubble simulations. International Journal of Multi
phase flow, 29:109, 2003.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
M. P. Sobera and C. R. Kleijn. Hydraulic permeability of
ordered and disordered singlelayer arrays of cylinders.
Physical Review E, 74(036301), 2006.
M. A. Tahir and H. Vahedi Tafreshi. Influence of fiber
orientation on the transverse permeability of fibrous me
dia. Physics offluids, 21(083604), 2009.
S.P. Timoshenko and J.N. Goodier. Theory ofElasticty.
McGrawHill Book Co., 1970. 3rd Edition.
S. Zobel, B. Maze, H. Vahedi Tafreshi, Q. Wang, and
B. Pourdeyhimi. Simulating permeability of 3d calen
dered fibrous structures. Chemical Engineering Science,
62:62856296. 2007.
