A computational continuum model of poroelastic beds

2.5 MB·

by  Uǧis Lācis

A computational continuum
model of poroelastic beds
rspa.royalsocietypublishing.org
1 2 1
U. La¯cis , G. A. Zampogna , S. Bagheri
1
LinnéFlowCentre,DepartmentofMechanicsKTH,
Research
SE-10044Stockholm,Sweden
2
DICCA,ScuolaPolitecnica,UniversitàdiGenova,via
Articlesubmittedtojournal Montallegro1,16145Genova,Italy
7 Despite the ubiquity of fluid flows interacting with
1 SubjectAreas:
porous and elastic materials, we lack a validated
0 Biomechanics,computational
non-empirical macroscale method for characterizing
2
mechanics,fluidmechanics the flow over and through a poroelastic medium.
n
We propose a computational tool to describe such
a
J Keywords: configurationsbyderivingandvalidatingacontinuum
3 Poroelasticity,connected-structures, modelfortheporoelasticbedanditsinterfacewiththe
1 numericalsimulation,anisotropy abovefreefluid.Weshowthat,usingstresscontinuity
conditionandslipvelocityconditionattheinterface,

] the effective model captures the effects of small
n
y Authorforcorrespondence: changes in the microstructure anisotropy correctly
d ShervinBagheri and predicts the overall behaviour in a physically
u- e-mail:[email protected] consistent and controllable manner. Moreover, we
show that the performance of the effective model is
l
f
accuratebyvalidatingwithfullymicroscopicresolved
.
s
simulations. The proposed computational tool can
c
i be used in investigations in a wide range of fields,
s
y including mechanical engineering, bio-engineering
h andgeophysics.
p
[

1
1. Introduction
v
6
Recent advances in surface micro- and nano-fabrication
9
5 techniques[1–3]areprovidingnewtechnologicaloppor-
3 tunities of enormous potential. However, we lack high-
0 fidelitymodelscapturinghowtheunderlyingsmall-scale
.
1 physicochemical processes interact with the large-scale
0 flow and heat- and mass-transport phenomena. The
7
reasonisthevastrangeofscalesinbothtimeandspace
1
thatneedberesolvedinordertocapturethefullphysical
:
v picture,whichrendersfull-scalenumericalinvestigations
i
X extremely costly[4]. The development of multi-scale
models with reduced complexity is a necessary enabler
r
a in an increasing number of applications where the
control of events at the small length scales determines
the properties of the flow over much larger space and
time scales. Examples include, but are not limited to,
(cid:13)c TheAuthors. PublishedbytheRoyalSocietyunderthetermsofthe
CreativeCommonsAttributionLicensehttp://creativecommons.org/licenses/
by/4.0/, whichpermitsunrestricteduse, providedtheoriginalauthorand
sourcearecredited.
thedesignofnovelsurfacestocontrolheattransfer,skinfrictionorpressuredrag,acousticnoise
2
andfluidmixing.
Thereisthusaclearneedtodeveloptechniquesthatcombinecomputationalfluidmechanics
with mesoscopic or microscopic models of materials. In this work, we take a step in this ....rspa
dBiyre“cmtioesnosbcyopmico”dewlleinmgetahnesitnrtuecrtaucrteiosnthoaftmcaensobsecodpesiccrsiuberfdacbeytceoxntutirneusuwmithmemthaocrdoss,cboupticarfleoswtisll. ......roya
significantlysmallerthanthelarge-scaleflowphenomena.Whentheseconstructions,whichcan ...lso
.c
taketheformofhigh-aspectratiostructures(suchasfibres)orgranular-typestructures(suchas ..ie
particles),aresaturatedwithafluidandcanabsorbthefluidstressandpossiblydeform,wehave ...typ
aporoelasticmaterial. ...ub
The description of poroelastic medium as a continuum has an extensive and broad history, ...lish
bfoergminunliantigonwoitfhtotthalestermespsirtiecnaslormoofddelesfoormf apbolreoupsoromuesdmiaedbiya bDyarBciyo[t5[6],a7n].dTofodlaloy,wtehderebyexitshtes ........ing.org
asignificantamountofwork[8–29]onporoelasticmediaadoptingtop-downapproachesbased .
.P
empiricalmacroscalestress/deformationtests,bottom-upapproachesusinghomogenizationand ....roc
volumeaveragingtechniques,physicalapproachesusinganalyticalandmechanicalmodels,or ..R
.
othermethods. .S
Inthiswork,weemploythemethodofhomogenizationviamulti-scaleexpansiontomodel ...oc
..A
flows through and over poroelastic surfaces. Within the general multi-scale analysis platform ..0
.0
(MAP) classification proposed by Scheibe et al.[30], the current method falls into the category ...00
of methods for which the different scales can be completely decoupled from each other (e.g. .0
“motif B” in [30]). In other words, our work is “formal upscaling”, which means that after the ...00
.
macroscale governing equationshave been derived, the parametersof themodel aregoverned
by fully decoupled microscale problems (from now on, we use the term “microscale” instead
of “mesoscale” to conform with the terminology within two-scale expansion approach). This
methodlimitsourinvestigationstoset-upsforwhichtheunderlyingmicroscaleclosureproblems
are linear. For example, the fluid phase in the poroelastic material must be slow enough such
thatinertialeffectscanbeneglectedormodelledthroughsomekindoflinearisation.Essentially
all previous upscaling works on deformable porous media has been in the same motif B
category. Examples include the method of volume averaging by Whitaker[31] and the method
ofhomogenization,asemployed,forexample,byMei&Vernescu[12].
There exists now a number of studies [32, and references within] on development and
applicationsofthemotifBmulti-scalemethodsforporoelasticmedia.Itis,however,theauthors’
opinionthatwestilllackacomputationalframeworkthatlaysoutthesequentialstepsneeded
to be taken in order to obtain the physical parameters describing the poroelastic medium not
only in the interior of the material, but also at the interface with freely moving fluid. Ideally,
such a framework, when provided a particular microscopic structure in terms of its geometry
and properties (skeleton elasticity, connectivity, etc.), will provide the anisotropic macroscopic
material properties (permeability, elasticity, etc.) of the effective continuum fields without any
fittingparametersobtainedfromexperiments.
To the best of the authors’ knowledge, the theoretical frameworks of Whitaker[16] and Mei
& Vernescu[12] have not been validated from microscopic point of view for the flow over a
surface, which is coated with porous and elastic material. Validations of these methods most
oftenconsideronlymacroscopicmeasures.Morespecifically,whilemicroscaleproblemsinunit
cellsneededintheupscalingprocedureshavebeenpresentedandsolvedpreviously[12,16,33],
thereisnocomparisonbetweenglobalmacroscalesimulations–usingpropertiesobtainedfrom
thosemicroscalesolutions–andcorrespondingfullyresolvedsimulations.Inthiswork,wedo
not only compute effective material tensors by numerically solving microscale problems, but
alsocomparetheobtainedeffectivecontinuumdescriptionwithfullyresolveddirectnumerical
simulation(DNS)ofthefluidflowinsideandabovetheporoelasticsurface.Inthisway,wecan
assesstheaccuracyoftheeffectivemodelquantitativelybothintheinteriorandneartheinterface
ofthemediumwithafree-flowingfluid.
A particularly important contribution is the insight it provides on how the velocity and the
3
stress interface conditions between the poroelastic medium and the free fluid region model
the microscale effects in averaged manner. Gopinath & Mahadevan[29] as well as Minale[34]
point out that the effective interface condition for a poroelastic region may in fact be more ....rspa
stotrabiaglhatnfocerwtahredfltuoidmasttcrhes,scofrmopmartehdetofreteheflruigididwpiothrouthsemsaotleidriaslt,resisnsceofththereesiusrafanceatumraaltewriaayl. ......roya
However,thereexistsnovalidation–wheretheinterfaceboundaryconditionsfordisplacement, ...lso
.c
fluid velocity and pressure are coupled to the Stokes equations above – for a non-trivial flow, ..ie
where there is transport across the interface. Most works treating the boundary conditions are ...typ
empirical[35–38]. Those contributions which have recently treated the interface problems from ...ub
first-principleshavefocusedonrigidporousmediaandone-dimensionalproblems,suchasthe ...lish
ltahmeminaatrecrhiaalnannedlflthoewfr[e3e9,fl4u0i]do),nolyri(nfofirltwrahtiiochntflhoewreoisnlnyo[4t1r]a.nsferofmassormomentumbetween ........ing.org
Insummary,theobjectivesofthecurrentworkareto(i)presentaframework,derivedusing .
.P
multi-scaleexpansion,suitabletomodelflowthroughandoverporoelasticmaterials,(ii)validate ....roc
the framework with respect to the fully resolved direct numerical simulations, and (iii) gain ..R
.
insight on stress transfer near the interface and evaluate the accuracy of interface boundary .S
conditions.Ourworkispartiallythenumericalcounterpartoftheanalyticalandasymptoticstudy ...oc
..A
byGopinath&Mahadevan[29];weaimatcomputing(insteadofphysicallymodelling)physical ..0
.0
materialpropertiesofanisotropicporoelasticmaterialsaswellastheinterfacewiththefreeflow. ...00
WhileGopinath&Mahadevan[29]considerparticularlybiologicallyrelevantmicrostructuresof .0
ordered or disordered filaments, we consider connected materials consisting of linked spheres ...00
.
andellipsoids.
This paper is organized as follows. In section 2, we introduce both the microscale and the
macroscale/effectiveequationsgoverningourproblem,whichconsistsofaporoelasticmaterial
attheinterfacewithamovingfluid.Amethodtocomputetheeffectivepropertiesofporoelastic
materials is presented in section 3. In the same section, we provide effective tensor results for
cubic-symmetricandmonoclinic-symmetricporoelasticmaterialsandanalysethem.Alid-driven
cavity problem to investigate poroelastic material response to a steady two-dimensional flow
vortexisproposedandsolvedusingthehomogenizedequationsinsection4.Inthesamesection,
wereportresultsobtainedfromresolveddirectnumericalsimulationsandexplaintheshearstress
transferbetweenfreefluidandporoelasticmedium.Insection5,wediscussthelimitationsofthe
presentedtheory.Finally,insection6weconcludethisworkandoutlinefuturedirections.
2. Micro-andmacroscaleequationsdescribingaporoelasticbed
Foradenseporoelasticmediumexposedtoafreeflow(Fig.1),onemaydefineatleasttwolength
scales; a microscopic (pore) scale l, which characterizes the size of voids in the material and
a macroscopic (global) scale H, which characterizes the size of the large-scale processes in the
mediumornearby.Thelargescaleflowischaracterizedalsobythecreatedpressuredifference
∆P. In this section, we present the microscale equations that resolve every scale in the full
physical domain, and a set of macroscale equations that model the effective average behavior
ofthe poroelasticbedunder agivenset ofassumptions. Althoughthe effective fieldequations
fordescribingafluid-saturatedporoelasticmaterialareknown[6,29],were-derivethemfrom
first principles in the supplementary appendices. The main reason is to uncover the detailed
microscale problems in unit cells that are required for determining the physical coefficients
appearinginthemacroscalesystems.
(a) Microscale governing equations
Themicroscaleequationsresolvingthefluid-structurephysicsatscale∼laretheconservationof
massandmomentum.ForaNewtonianfluidwithdensityρf andviscosityµ,themomentumin
Stillfluid
Flow 4
µ,ρ vortex
f
∼H
∆P ....rspa
.....roy
.a
..ls
.o
.c
..ie
..ty
.p
...ub
l ...lish
..in
Lsoilniedasrkeellaesttoicn,Csk,ρs ......g.org
.
Rigidwall .P
....roc
..R
Figure1. Illustrationofafreefluidvortexinteractionwithaporoelasticmaterial.Thematerialiscomposedofmultiple .
.S
instancesofunit-cellskeletongeometry–spherewithcircularcylinderconnectionsinalldirections.Werepresentaslice ...oc
onepore-structurethick.Theunit-cellisacubewithsidelengthl.Thesolidskeletonischaracterizedbylinearelasticity ..A
TtehnesosrolCidsskkealnedtodneinssditeyfoρrsm.eTdheunfldoewrtvhoertienxfluheanscaeloefntghtehfsrecaelefluHidvaonrdtecxa.usesacharacteristicpressuredifference∆P. ......0000
.0
...00
.
thefreefluidregionandwithintheporousskeletonisgovernedbytheincompressibleNavier-
Stokesequations
ρf[∂tu+(u·∇)u]=−∇p+µ∆u, (2.1)
∇·u=0, (2.2)
whereuandparevelocityandpressurefields,respectively.
The material is defined by the solid skeleton density ρs and linear elasticity tensor Csk.
Assuming an isotropic material the elasticity tensor is defined by Young’s modulus E and
Poisson’sratioν.Thesolidskeletonmomentumisgovernedbyabalancebetweensolidinertia
andstress,obtainedusingthelinearstress-strainrelationship
(cid:26) 1(cid:104) (cid:105)(cid:27)
ρs∂t2v=∇· Csk:2 ∇v+(∇v)T , (2.3)
wherevisthedisplacementfieldofthesolidskeleton.Inordertocouplethefluidandstructure
problems, the no-slip condition and continuity of stress is prescribed at the boundary between
solidskeletonandsurroundingfluid,i.e.,
(cid:110) (cid:104) (cid:105)(cid:111) (cid:26) 1(cid:104) (cid:105)(cid:27)
u=∂tv and −pδ+µ ∇u+(∇u)T ·nˆ= Csk: ∇v+(∇v)T ·nˆ. (2.4)
2
Here,δ isthesecond-rankidentitytensorandnˆ isunit-normalvectorattheboundary.Solving
thegoverningequationseverywhereattheporescaleiscomputationallyveryexpensivedueto
agloballylargedomainandduetorequirementoffineresolutionnearthepores.Thismotivates
the development of an alternative continuum description, where the pore fluid and solid are
considered as one composite, which given appropriate equations and constitutive relations
presentstheaveragebehavioroftheactualporoelasticbed.
Artificialboundarybetween
poroelasticmaterial Flow 5
vortex
andfreefluid
∼H
∆P ....rspa
.....roy
.a
..ls
.o
Homogenizedporoelasticmaterialwith .c
secffaelcetsiveepaerlaatsitoicni,tyε,≪C,1s,opliodrodseintysi,tθy,,pρes,rmeability,K, .....ietyp
porepressurecontribution,α ...ub
..lis
.h
..in
......g.org
.
Rigidwall .....Proc
..R
Figure 2. Illustration of a free fluid vortex interaction with a poroelastic material. The material is described in a .
.S
homogenizedtwo-domainsettingandischaracterisedbythescaleseparationparameter(cid:15)andporosityθ,aswellas ...oc
effectivepermeabilityK,effectiveelasticityC,densityρsandporepressurecontributionα.Thehomogenizeddomain ....A0
isnotvisiblydeformedundertheinfluenceoffreefluidvortex.Thereisanartificialboundaryintroducedbetweenthe .0
homogenizedporoelasticmaterialandfreefluid. ...00
.0
...00
.
(b) Effective field equations
Wedividethephysicaldomainintotwoparts;onecontainingonlythefreefluid,andtheother
containingthefluidandsolidskeleton(Fig.2).Thefree-fluidregionisgovernedbytheNavier-
Stokesequations(2.1–2.2).
Thecontinuumdescriptionofthecomposite(solidandfluid)isbasedonaseparationbetween
thepore-scaleandthesystemscale.Mathematicallythiscanbeformulatedbyascaleseparation
parameter (cid:15)=l/H(cid:28)1. The effective displacement field v of the homogenized poroelastic
material(seeFig.2)isgovernedby,
(cid:20) 1(cid:16) (cid:17) (cid:21)
(1−θ)ρs∂t2v=∇· C:2 ∇v+(∇v)T −αp . (2.5)
Here, θ=Vf/V is the porosity (which in general is a function of space, but in this work is a
constant),whereVf isthefluidphasevolumeandV isthetotalvolumeofthecompositeinone
unitcell(definedlater).Moreover,Cisthefourth-rankeffectiveelasticitytensorofthematerial
andαisthecoefficientforthecontributionoftheporepressure(p)inthetotalstress.Whereas
themicroscaleelasticityCskonlydependsonEandν,theeffectivetensorCingeneraldepends
ontheporosity,poregeometry(asmodelledinthecurrentwork)andontheparticulartypeof
boundaryconditionimposedneartheinterfacewiththefreefluid,ortheimpermeablewalls[29].
Thetensorαdoesnothaveacorrespondenceinthemicroscaleandisaneffectofsolidskeleton
deformationduetoseepageflowthroughthepores.WecharacterizeindetailbothCandαfor
twodifferentporoelasticmediainsection3.
Moreover,forsufficientlydenseporoelasticmaterialthefluidflowbetweentheporesisslow,
suchthatinertialeffectsarenegligible.Thereforetheporepressureisthedominantcontribution
inthefluidflowandtheleading-orderequationisthe(relative)Darcy’slaw,
K
u−θ∂tv=− ·∇p. (2.6)
µ
Thisexpressionrelatesthegradientoftheporepressure,thesolidvelocity∂tvandtheflowfieldu
intheporoelasticmediumtoeachother.HereKistheinteriorpermeabilitytensor.Thetermθ∂tv
arisesduetoviscousfrictionbetweenthesolidstructureandtheporefluid;ifthereisamotionof
6
thesolidskeleton,thesurroundingporefluidisdraggedalongthroughboundarycondition(2.4).
Finally,theconservationofmassrequiresthat
∇·u=D:21∂t(cid:16)∇v+(∇v)T(cid:17)−E∂tp, (2.7) ........rspa.roy
where D=θδ−α is a dimensionless second-rank tensor. It determines how the strain of the ...als
.o
displacement(causedeitherbytheflowthroughtheporesorbyaboundarycondition)modifies .c
..ie
thesolidstructurevolumewithinonepore,thussqueezingtheporefluidinoroutofthepore. ..ty
ThescalarE characterizesthechangeofsolidstructurevolumewithinoneporewithrespectto ....pub
timevaryingpressure,which,similarlytosolidstrain,cancausechangeinporefluidcontentand ..lis
.h
consequentlyintroduceapparentcompressibilityoftheflowfield. ..in
comTbhienisnygst(e2m.6)owfitehqu(2a.t7io),ntshe(2s.5y-s2t.e7m) dceatnerbmeinreedsutcheedsteovfeonurunuknnkonwownsns(u(v,vanadndp).p)T,hbeufltubidy ......g.org
.
velocityucanthenbecomputedasapostprocessingsteponcepressureanddisplacementfields .P
areknown.Followingthisapproach,theequationfortheporepressureis ....roc
E∂tp−∇·(cid:18)Kµ ·∇p(cid:19)=−α:12∂t(cid:16)∇v+(∇v)T(cid:17). (2.8) .......RSoc
..A
Notethattheeffectivesystemabovefordescribingaporoelasticsystemcanbederivedusing ..0
.0
multi-scale expansion (see section 3 and supplementary appendices), motivated using mixture ...00
tahpeporroyp,roiraptehybsoiucanldlyarmyocdoenldleidti.oTnhseamt ainjoterrcfhacaellsenwgiethisstooluidsewthalilssd,efrveeeloflpueiddesffoerctoivtheesryssttreumctwuritehs .....000
to describe problems arising from various applications. The aim of this work is to provide a
frameworkfromwhichonecanformafullyclosedeffectivesystemtodescribetheresponseand
interactionofaporoelasticmaterialwithasurroundingfreefluid.Inthenextsection,wetherefore
providetheneededboundaryconditionsbetweenaporoelasticmaterialandarigidwall,andalso
afreefluid.
(c) Effective interface conditions
(i)Conditionsfortheporoelasticbed
To solve the governing equations for the poroelastic material, one needs to impose boundary
conditionsforboththepore-pressure(equation2.8)andthedisplacementfield(equation2.5).On
rigidwalls,weimpose(similarto[29])zerodisplacementandzerotranspiration(normalfluid
velocity),whichthroughtherelativeDarcy’sequation(2.6)leadsto
(cid:12)
∂p(cid:12)
v|rw=0 and ∂nˆ(cid:12)(cid:12) =0. (2.9)
rw
Here, “rw” means “rigid wall” and nˆ is the unit-normal vector at the wall. Physically, the no-
slip condition should be satisfied at the wall, but this is not compatible with the leading-order
presentation of a poroelastic media based on Darcy’s law. The Darcy’s law only describes the
directproportionalitybetweenthepore-pressuregradientandthevelocity,anddoesnotinclude
anymacroscopicdiffusioneffects.
Attheartificialinterfacewiththefreefluid,showninFig.2,weimposeapressurecontinuity
condition
p−=p. (2.10)
Here,theporepressurebelowtheinterfaceisdenotedbyp−andpressureofthefreefluidabove
theinterfacewithp.Notethatthechoiceofconditionforpressuredependsontheassumptions
made about the flow as well as the material geometry. For example, La¯cis & Bagheri[42] have
shownthat,iftheinterfacevelocitycausedbytheshearstressisofthesameorderasthevelocity
induced by the pore pressure gradient, the pressure continuity is the leading order boundary
condition for any pore geometry. On the other hand, Mikelic´ et al.[43, 44] have shown that, if
velocitycontributionfromtheshearstressattheinterfaceisoneorderhigherthancontribution
7
fromthepressuregradient,thereisapressurejumpforanisotropicporegeometry.However,if
theporegeometryisisotropic/cubic-symmetric,thepressurecontinuityisstillapplicable.Finally,
forthematerialdisplacement,weimposecontinuityofstressesattheinterface ....rspa
(cid:20)C:21(cid:16)∇v+(∇v)T(cid:17)−αp(cid:21)·nˆ=(cid:104)−pδ+µ(cid:16)∇u+(∇u)T(cid:17)(cid:105)·nˆ. (2.11) ........royals
.o
.c
We thus assume, similarly to [29], that the total stress of the free fluid is transferred to the ..ie
totaleffectivestressoftheinteriorporoelasticmedium.Ingeneral,theeffectiveelasticityofthe ..ty
.p
compositeneartheinterfacecouldhoweverbedifferentfromitsvalueintheinterior(itisargued ...ub
in general to depend on boudary conditions[29]) and, to arrive to a more accurate boundary ..lis
.h
condition,onecouldconstructaninterface-cellwithanelasticityproblem,similarlyasdonefor ..in
tshheeavreslotrceistysibsoturanndsafreyrrceodntdoittihoen[s4o2l]id.Ostnreesosb(jfiecrtsitvaecorofstshitshepainpteerrfiascteotuhnendeinrsstiadnedthheowbedth)eaflnduiidf ......g.org
.
.P
tehffeeicntitveerfsatcreescsowrrietchtitohneiisnnteercioesrspaaryra;minesteecrtsioisna4b,lewteowcaiplltsuhroewthtehtartaanpsfperrooxfimstraetsiosnreoafstohneaibnltyerwfaeclle. ......rocR
.
.S
(ii)Conditionsforthefreefluid ...oc
..A
To solve for the Navier-Stokes equations in the free fluid domain, we need to impose velocity ..0
.0
–inprinciple,onecouldalsohavestresscondition,butitisalreadyusedforthedisplacements ...00
.0
of the poroelastic material – boundary conditions at the interface. In this work, we extend the ...00
velocity boundary condition derived by La¯cis & Bagheri[42] for a rigid porous bed to include .
poroelasticity1.Theconditionforthetangentialinterfacevelocityis
(cid:32) (cid:33)
Kif (cid:104) (cid:105)
u·τˆ=∂tv·τˆ+ − ·∇p−+L: ∇u+(∇u)T ·τˆ, (2.12)
µ
wheretheunitvectorτˆdenotesbothtangentialdirectionsofthesurface.Notethatthepressure
gradientisthepore-pressuregradientfromporoelasticmaterialsideoftheinterface(hencethe
minus superscript), whereas the flow velocity field u is on the free-fluid side. The interface
velocity has two distinct contributions; i) the no-slip contribution governed by the movement
of the solid structure; ii) the slip contribution, which is caused by the porosity of the solid
structureanddependsbothonporepressuregradientandfree-fluidshear.Theslipcontributionis
characterizedbythesecond-rankinterfacepermeabilitytensorKif andthethird-ranksliplength
tensorL.Foradensematerial,thefirstsliptermscalesasl2andissignificantlysmallerthanthe
secondterm,whichscalesasl.Themotionoftheporoelasticmaterialortheno-slipcontribution
depends not only on the pore-length l, but also on the elasticity of the material and the flow
regime.Theconditionforthenormal“penetration”interfacevelocitycomponentissetbymass
conservation,i.e.
(cid:18) (cid:19)
K
u·nˆ=∂tv·nˆ− ·∇p− ·nˆ, (2.13)
µ
whereKistheinteriorpermeabilitytensordefinedinequation(2.6)oftheporousmediumand
∇p−denotesthepressuregradientwhenapproachingtheinterfacefromthebed.Similarlyasfor
thetangentialcomponent,thevelocityhastwoparts–theno-slippartandthe“slip”part.The
slip inthe normaldirection isessentially thefluid masstransport inand outof theporoelastic
material,whichhastobeequaltotherelativevelocityfromtheinterior.Thissameconditionalso
arisesfromtheinterfacecell[42],becausethetensorLcomponentscorrespondingtopenetration
interfacevelocityarezero.
Inthenextsection,wewillcomputethephysicalparameterswhichcharacterizetheinterior
poroelasticmedium(C,α,KandE)aswellastheparameterscharacterizingtheinterfacewith
the free fluid (Kif,L) for two specific pore-scale geometries. In section 4, we will investigate
inmoredetailsinterfaceconditionsatboundarybetweenporoelasticmaterialandfreefluidby
1ThediscussioncanbefoundinsupplementaryAppendixB,afterequation(B19).
comparingtheeffectivemodelagainstthemicroscopically-resolvedsystemusingthelid-driven
8
cavityproblem.
....rspa
3. Properties of poroelastic material exposed to free fluid .....roy
.a
Wenowturntodeterminingtheeffectivepropertiesofaporoelasticmaterialexposedtothefree ..ls
.o
fluid,byintroducingandsolvingasetofparticularmicroscaleproblemsinaninteriorunit-cell ...cie
andinaninterface-cell.Wefocusonaperiodicsolidstructurethatcanbeobtainedbyduplicating ..ty
.p
a single pore structure in all directions. The sketch in Fig. 3a shows how the material and its ...ub
interfacewithfree-fluidisdividedintocubicinteriorcellsandelongatedrectangularinterface- ..lis
.h
cells.Theeffectiveparametersarecomputedbysolvingtwoelasticityproblems(toobtainC,α ..in
ianntderEfa)cea-ncdellon(teofloubitdainprKobilfe,mL).(tTohoebsetaimnicKro)sicnalteheprionbtelermiorscaerlelfianrsdt tilwluostflrautieddpfroorbtlehmesinitnertihoer ......g.org
.
.P
dmoimcroaisncaulesipnrgobalewmesaaktlythaenidntearfsatrcoenwgiltyhathneisofrtereopfliucidm.icrostructure, before we move on to the ....roc
..R
.
.S
...oc
(a) Interior of the poroelastic material ..A
..0
.0
(i)Microscaleproblemsforelasticityparameters ...00
.0
The effective poroelastic bed obtained from the unit-cell approach can not have fully isotropic ...00
.
elasticpropertiesduetotheboundariesbetweencubicunitcells,asshowninFig.3a.Theresulting
effective elasticity C will at most exhibit a cubic symmetry; a symmetry, which is defined by
planes–paralleltotheunit-cellsidesanddiagonalacrosstheunit-cellinalldirections–going
throughthecenterofthecube.
Inordertopreservecubicsymmetry,theunit-cellcanbefilledwithanystructurethatitself
iscubicsymmetric,suchasWigner-Seitzgrains[12,33]orcubewithsphericalholes[45].Inthis
work,weuseasphereatthecenterofunit-cellwithradiusR,whichisconnectedtoneighbouring
cellsviacircularcylindersofradiusr,asshowninFig.3b.Thesolidskeletonstructureisbuilt
fromisotropicelasticmaterialwithelasticitytensor
(cid:26) (cid:27)
Csk =E ν δ δ + 1 (cid:0)δ δ +δ δ (cid:1) .
ijkl (1+ν)(1−2ν) ij kl 2(1+ν) ik jl il jk
Theeffectiveelasticitytensorcanbedetermined(seeAppendixBforderivation)usingathird-
rankdisplacementtesttensorχas
(cid:28)1(cid:104) (cid:105)(cid:29)
C=(1−θ)Csk+Csk: ∇χ+(∇χ)T , (3.1)
2
where we define transpose of the fourth-rank tensor ∇χ acting on the first two indices
(∇χ)T =(∇χ) .Theexpressionaboveprovidesadirectlinearlinkbetweenthemicro-and
ijkl jikl
macroscaleelasticitytensorsthroughtheknownporosityθandthepore-scalegeometrycaptured
byχthatwillbecomputedbysolvingequations(3.3–3.4).Here,itisassumedthattheskeleton
elasticitytensorCskisconstantinspace.Thebracketsdenotethevolumeaverageovertheinterior
unit-cellvolume,whichisfilledeitherbysolidorbyfluid
(cid:90)
1
(cid:104)f(cid:105)= fdV, (3.2)
l3
Vσ
where Vσ is the volume of either the solid or fluid phase. In general, the effective properties
shouldbere-evaluatedastheporositychanges.However,weassumesmalldisplacementsofthe
effectivesystemsothattheporosityisroughlyconstantandthereforeitisenoughtocomputethe
averagesonlyonce.
(a) (b)
9
Freefluid
Icnetlelrface l ........rspa.roy
.a
..ls
.o
.c
..ie
l ..ty
.p
2r l ...ub
z y z y .....lishin
x Interiorcell x ......g.org
.
.P
Figure3. Inframe(a)weshowaconstructedporoelasticmaterialconsistingof84unit-cellswithvolumel3,whichcan ....roc
exhibitonlycubic-symmetryduetoboundariesbetweentheunit-cells.Abovethematerial,thereisafreefluid.Interior ...R
.S
unit-cellisdepictedwithgreen,interface-cellisdepictedwithred.Inframe(b)weshowacubic-symmetricmicro-structure ...oc
placed in a single interior unit-cell. The structure is built using a sphere with radius R=0.3l and cross-connected ..A
cylinderswithradiusr=0.1l. ...00
...00
.0
...00
.
Themicroscaledisplacementtensorχisdefinedinthesolidphase(havingvolumeVs)ofthe
unit-cellandisthesolutionofthefollowingproblem,
(cid:20) 1(cid:110) (cid:111)(cid:21)
∇· Csk: ∇χ+(∇χ)T =0, (3.3)
2
(cid:20) 1(cid:110) (cid:111)(cid:21) (cid:104) (cid:105)
Csk: ∇χ+(∇χ)T ·nˆ= Csk:δ(4) ·nˆ. (3.4)
2
Here, the fourth-rank identity tensor appearing in the boundary condition is defined as
(cid:16) (cid:17)
δ(4) =δikδjl. This equation corresponds to a standard steady linear-elasticity problem
ijkl
generalisedtothird-ranktensorsubjectedtodifferentcombinationsofsurfaceloading(appliedon
theinterfacebetweenthesolidskeletonandtheporefluid)inordertocharacterisetheresponse
ofthestructuretoallpossiblesurfaceloadingscenarios.Tocompletetheformulation,periodic
boundary conditions are applied to solid surfaces, which are in contact with the boundary of
theunit-cell.Inordertorenderthesolutionunique,weimposeconstraintsonaveragevaluesof
displacement[12]usingpenaltytermsinweakformulation[46].
The unit-cell domain has been meshed2 using GMSH software[47] and equations (3.3 – 3.4)
hasbeensolvedusingFreeFEM++[48].ThesphereradiusisR=0.3landthecylinderradiusis
r=0.1l, which results in porosity θ=0.85. Assuming the Poisson’s ratio to be ν=0.3 (due to
thelinearityoftheproblem,thesolutionisvalidforanyYoung’smodulusvalueE),theeffective
elasticitytensorforthecubic-symmetricporoelasticmediuminVoigtnotation[49,p.136]is
 
4.792 0.239 0.239 0 0 0
 0.239 4.794 0.239 0 0 0 
 
 0.239 0.239 4.792 0 0 0 
C= ·10−2E. (3.5)
 0 0 0 0.246 0 0 
 
 
 0 0 0 0 0.246 0 
0 0 0 0 0 0.246
2Weusemeshspacing∆s1=0.025latthesolidskeletonand∆s2=0.10lattheboundariesoftheunit-cell.Wehave
alsocarriedoutthesimulationoncoarsermeshes(∆s1=0.05land∆s2=0.20l;∆s1=0.10land∆s2=0.40l)and,by
observingtheconvergenceofresults,wesetentriescontainingonlynumericalnoisetozero.
Note that the entries of the effective elasticity tensor are significantly smaller (up to 300 times)
10
compared to these of the skeleton elasticity. This is because most of the solid mass has been
removed,leavingonly15%solidvolumefractionifcomparedtoacompletelyfilledcase,hence
theresultingmaterialismuchsofter.Anotherobservationisthattheshearcoefficients(diagonal ....rspa
ewleamveencotsefifinci3en×ts3[50b,optt.o2m2]r(idgihatgomnaatlreixlembleonckts)ianre3×no3wtompulecfhtmsmatarlilxerblcoocmk)p.aTrheidsetoffepctrecsasnurbee- ......roya
attributedtoconnectingnetofcylindersbetweenspheres;thecylindersaremucheasiertobend, ...lso
.c
comparedtobulkmaterial,i.e.thereisnocontinuoussupportfromthesides,astherewouldbe ..ie
inthecontinuousmaterialcase.Theeffectivetensorcanbecharacterizedbythreeindependent ...typ
parameters3andhastheformofamaterialwithcubicsymmetry[49,p.99]. ...ub
The other two elasticity-related effective parameters α and E can be computed through ...lish
eexqps.re2s.4s,i2o.n7s,24.8in].vAolltveirnngattihveeliyn,voenresecaonfsthoelvfeofuorrtha-nraandkdietliaosntaiclittyesttednissoprlaecfefemcetinvtefieelladstηicigtyovCern[2e9d, ........ing.org
bythelinearelasticityproblem .
.P
∇·(cid:34)CEsk :12(cid:110)∇η+(∇η)T(cid:111)(cid:35)=0, (3.6) .......rocR
.S
(cid:34)CEsk :12(cid:110)∇η+(∇η)T(cid:111)(cid:35)·nˆ=δ·nˆ. (3.7) .......ocA0
.0
Notethatthisproblemisdrivenbydifferentsurfaceforcingcomparedtotheχproblem,butis ....000
otherwise solved under the same conditions. Given η, the pore-pressure contribution tensor is ...00
.
computedfrom,
(cid:42) (cid:43)
Csk 1(cid:104) (cid:105)
α=θδ+ : ∇η+(∇η)T , (3.8)
E 2
andthecoefficientdescribingtheporoelasticmaterialresponsetothetime-variationofthepore
pressureiscomputedfrom
(cid:28) (cid:29)
∇·η
E= .
E
Bysolving(3.6–3.7)fortheparticularmicrostructurewithcubicsymmetryshowninFig.3b,we
obtain
α=0.9789δ and E=0.1571E−1.
The tensor α is not unity, therefore one can conclude that the solid phase of the poroelastic
materialisnotincompressible.Thisconclusionagreeswiththeusedelasticityparameter(weuse
ν=0.3,whileincompressiblesolidshaveν=0.5).However,itisnon-trivialtousetensorαasa
measureofcompressibility,becauseitisalsoafunctionofporosityθ.
(ii)Microscaleproblemforpermeabilitytensor
Turning attention to the fluid flow, the poroelastic material is characterized by a permeability
[5],whichinourframeworkisobtainedbysolvingthefollowingStokesproblemintheinterior
unit-cell(Fig.3b),
−∇A+∇2K=−δ, (3.9)
∇·K=0. (3.10)
Here,Kisasecond-ranktensorfield;Kij istheithvelocitycomponentofjthvectorialStokes
problemassociatedwiththepressurefieldAj.ThethreevectorialStokesproblemscharacterize
howtheporeflowrespondstovolumeforcinginonespatialdirectionatatime.Thisequation
systemiscomplementedwithno-slipboundaryconditionattheinterfacewiththesolidskeleton,
3Thedifferenceof2·10−5betweenthefirstandsecondseconddiagonalelementisattributedtonumericalerrorthatcould
becausedby,forexample,imperfectionsinthegeneratedmesh.
4Forexplicitrelationshipsobtainedusingmulti-scaleexpansion,seesupplementaryappendixC.