Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 12.6.2 - Direct Numerical Simulation of Saturated Deformable Fibrous Porous Media
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00311
 Material Information
Title: 12.6.2 - Direct Numerical Simulation of Saturated Deformable Fibrous Porous Media Fluid Structure Interactions
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Khan, I.
Aidun, C.K.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: saturated deformable porous media
DNS
fluid-structure interaction
Lattice-Boltzmann method
finite-element method
 Notes
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 LBM-FEM method is an accurate and suitable method for further investigations on deformable porous media.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00311
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1262-Khan-ICMF2010.pdf

Full Text



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, Fluid-structure interaction,
Lattice-Boltzmann method, Finite-Element 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 LBM-FEM 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 LBM-FEM 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 fluid-structure
interactions in complex geometries

In this work a parallel hybrid Lattice Boltzmann and fi-
nite element method is used to model the fluid-structure
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 bounce-back scheme for
no-slip condition make LBM a favorable tool for model-
ing single-phase 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 four-node 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 "link-bounce-
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 link-wise repulsive force, based on the lubrication and
contact forces, is applied to the interacting surfaces. The
magnitude and direction of the link-wise 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 link-wise 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 non-conformal,
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 x-axis 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 fluid-structure 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 = 2R-1

4-p P


Load (P/HR2)


Figure 1: Comparison of the Hertz contact model
with the near contact model implemented in LBM-FEM
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 LBM-FEM 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 non-dimensional de-
formation a/R between the implemented contact model
in LBM-FEM 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 over-predicts 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 LBM-FEM
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 LBM-FEM 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 fluid-saturated porous me-
dia. They solved the problems of one-dimensional com-
pression of a semi-infinite 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 non-linear 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 LBM-FEM 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=2R-1


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 non-wovens.
The boundary condition at the bottom (z 0 ) is a
no-slip, no-displacement 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 non-dimensional 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 LBM-FEM 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.

9e-05


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 x-direction, 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;


le-05
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 LBM-FEM 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 :;I7-IMN/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 semi-infinite 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 LBM-FEM 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. One-dimensional 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 open-ended 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 non-dimensional 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 LBM-FEM 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


LBM-FEM(8cylinders) -----..
LBM-FEM(1cylnders) - -
LBM-FEM(32cylinders) **...
de Boer model -


de Boer et al -
LBM-FEM(32 cyl) -----
LBM-FEM(16 cyl) --
LBM-FEM(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 on-set 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 z-axis from the LBM-FEM simulations of 8, 16 and
32 cylinder geometries, along with the de Boer et al.
model data. The z-axis 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 z-axis
from LBM-FEM 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 LBM-FEM 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 LBM-FEM 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 LBM-FEM 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 semi-infinite 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 LBM-FEM 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 non-homogeneous 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 lattice-boltzmann
method for complex flows. Annual Review of Fluid Me-
chanics, 42:439-472, 2010.

C.K Aidun and Y. Lu. Lattice boltzmann simulation of
solid particles suspended in fluid. Journal of Statistical
Physics, 1:49-61, 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:251-264, 1999.

M.A. Biot. General theory of three dimensional consol-
idation. Journal of Applied Physics, 12:155-164, 1941.

M.A. Biot. Theory of deformation of a porous viscoelas-
tic anisotropic solid. Journal of Applied Physics, 27:
459-467, 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 fluid-solid 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:1129-1148, 1980.

R.M. Bowen. Compressible porous media models by
use of the theory of mixture. International Journal of
Engineering Science, 20(6):697-735, 1982.

R. de Boer, W. Ehlers, and Z. Liu. One-dimensional
transient wave propagation in fluid-saturated incom-
pressible porous media. Archive of Applied Mechanics,
63:59-72, 1993.

de Boer R. Theoretical poroelasticity a new approach.
Chaos Soliton Fract, 25:861-878, 2005.

H. Deresiewicz. Stress-strain relations for a simple
model of granular medium. Journal of Applied Mechan-
ics, 25:402-406, 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:81-97, 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):685-708, 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:209-221, 2008.

W. Ehlers. On thermodynamics of elasto-plastic media.
Archives of Mechanics, 41(1):73-93, 1989.











W. Ehlers and G. Eipper. Finite elastic deformations in
liquid-saturated and empty porous solids. Transport in
Porous Media, 34:179-191, 1999.

W. Ehlers and B. Markert. On the viscoelastic behavior
of fluid saturated porous material. Granular Matter, 2:
153-161, 2003.

Y. Hu, Y. Zhu, and C. Cheng. Dqm for dynamic response
of fluid-saturated visco-elastic porous media. Interna-
tion Journal of Solids and Structures, 46:1667-1675,
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 lattice-boltzmann 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:31-66, 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:
285-309, 1994.

R.M MacMeccan, J.R. Clausen, G.P. Neitzel, and C.K.
Aidun. Simulating deformable particle suspensions
using a coupled lattice-boltzmann and finite-element
method. Journal of Fluid Mechanics, 618:13-39, 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:
1-18, 1996.

A.S. Sangani and A. Acrivos. Slow flow through a pe-
riodic array of spheres. International journal of multi-
phase flow, 8(4):343-360, 1982a.

A.S. Sangani and A. Acrivos. Slow flow past periodic
array of cylinders with application to heat-transfer. In-
ternational journal of multiphase flow, 8(3):193-206,
1982b.

K. Sankaranarayanan, I.G. Kevrekidis, S. Sundaresan,
J. Lu, and G. Tryggvason. A comparative study of lattice
boltzmann and front-tracking finite-difference 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 single-layer 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.
McGraw-Hill Book Co., 1970. 3rd Edition.

S. Zobel, B. Maze, H. Vahedi Tafreshi, Q. Wang, and
B. Pourdeyhimi. Simulating permeability of 3-d calen-
dered fibrous structures. Chemical Engineering Science,
62:6285-6296. 2007.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs