Table Of Content3B2v8:06a=wðDec52003Þ:51c PHYSA : 9959 Prod:Type:FTP ED:JyothiG
þmodel pp:1215ðcol:fig::NILÞ PAGN:Uday SCAN:
ARTICLE IN PRESS
1
3
PhysicaA](]]]])]]]–]]]
www.elsevier.com/locate/physa
5
7
Effective medium equations for fractional Fick’s law in porous
9 media
11
(cid:1)
Francisco J. Valdes-Parada , J. Alberto Ochoa-Tapia, Jose Alvarez-Ramirez
13
DepartamentodeIngenierı´adeProcesoseHidrau´lica,UniversidadAuto´nomaMetropolitana-Iztapalapa,AparatadoPostal55-534,Mexico
D.F.09340,Mexico
15
Received8February2006;receivedinrevisedform26May2006
17
F
19 O
Abstract
21 O
This paper studies reaction–diffusion phenomena in disordered porous media with non-Fickian diffusion effects. The
aimistoobtainaneffectivemediumdescriptionoftheconcentrationdynamicshavingafractionalFick’slawdescription
23 for the particles flux. Since the methodology is based on a volume averaRging approach, a fractional spatial averaging
theoremisdevelopedto interchangeaveragingintegration andfractionaldifferentiation. Modelstructure simplifications
P
25 aremadeonthebasisofanorderofmagnitudeanalysisfromphysicalinsights.Theclosureproblemassociatedwiththe
effective diffusivity definition is also developed, showing that th e macroscale diffusion parameter is affected by (i) the
27 scaling from mesoscalesto macroscales, and(ii) bythe disordeDred structure of the porous medium.
r2006PublishedbyElsevier B.V.
E
29
Keywords: Non-Fickian diffusion; Fractional calculus; Reaction–diffusion; Porous media; Effective medium equations; Pseudo-
T
homogeneousequations
31
C
33 E
1. Introduction
35 R
Chemicalreactioniscommonlycoupledtotransportphenomenainmanynaturalandindustrialsystems.At
R
37 certain length scales, an equation widely used to describe the process is the traditional reaction–diffusion
equation section
O
39
qc
¼D=2cþCRðcÞ,
qt
41
where c is a gNiven particle concentration, D is the diffusion coefficient, and RðcÞ is a reaction rate. The
underlying hypothesis behind this reaction–diffusion model isthat the transport mechanisms are constitutive
43
propertiesUequivalent to the average of unmeasurable transport properties at arbitrarily small scales, and so
the parameters of the model are scale independent. That is, the transport mechanisms are invariant under
45
spatialscaling,implyingthatthemodelisabletodescribethetransportprocessatanyspatialandtimescales.
This corresponds to Fickian behavior, which is valid when the particle jump size (i) is uncorrelated in time,
47
(cid:1)
Correspondingauthor.Tel.:+525858044648;fax:+525558044900.
49
E-mailaddress:[email protected](F.J.Valdes-Parada).
51 0378-4371/$-seefrontmatterr2006PublishedbyElsevierB.V.
doi:10.1016/j.physa.2006.06.007
PHYSA : 9959
ARTICLE IN PRESS
2 F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]]
1 and(ii)hasfinitemeanandvariance[1–3].Fornon-reactivesolutes,(i.e.,RðcÞ¼0),thefundamentalsolutions
to the Fickian diffusion over time will be Gaussian densities with (finite) means and variances based on the
3 transportcoefficientD.However,fieldandexperimentalmeasurementshaveshownthatthestandardFickian
diffusion equation typically underestimates concentrations in the leading and/or trailing edges tracer plumes
5 [5–7]. Application of the standard framework to describe field data reveals an apparent scale-dependence of
dispersity, not consistent with Fickian transport, complicating the prediction of plume evolution in time and
7 space[3].Sincenon-Gaussianbreakthroughcurvesareoftenobserved,atleastoneoftheconditionsimplying
Fickian transportisfailing [3,4].Violation of theuncorrelated jump assumption leads toenhanced diffusion,
9 which is faster than Gaussian analytical solution predictions [8]. Most non-Fickian transport theories are
based on the effects of long-range temporal correlations due, for instance, to solute sorption or preferential
11 pathways. Recently, it has been suggested that non-Gaussian plumes can be explained by a violation of the
finite-varianceassumption[4,9,10].Schumeretal.[3]haveshownthatnon-Gaussiandistributionswithheavy
13 leading edges can be the result of the infinite-variance particle jump distributions that arise during the
transport in disordered porous media. They demonstrated that a fractional Fick’s law, as suggested by
15 Schumer et al. [3,4], is a governing equation for solute transport in porous media in cases where temporally
correlatedvelocityfieldsdonotdominatethetransportprocess.Inthisway,Fickiandispersioncanonlyoccur
17 in homogeneous media.
F
Scaling problems are encountered when Fickian models are applied to non-Fickian processes. In fact, the
19 infinite-varianceassumptionleadingtonon-FickiantransportimpliesthattheOFickiantransportpropertiesat
a given spatial scale are not necessarily the same at macroscales, leading to scale-dependent transport
21 parameters. For instance, from a Fickian perspective, a large particleOthat jumps at mesoscales could be
interpretedasnormal(Fickian)particlejumpsatmacroscalesbyvirtueofanaverageeffect.Hence,thereisa
23 scaling problem for this kind of particle transport. This is particRularly important for interpreting lab-scale
experimentalmeasurements.Infact,experimentsarecarriedoutatrelativelysmallspatialandtemporalscales,
P
25 and measurements arising from such should recover the transport behavior at real scales. Wrong
interpretation of transport coefficients can lead to, e.g. m alfunctioning of industrial equipments or failures
27 inthepreservationofaquifersperturbedbycontaminanDts.Toaddressthisscalingproblem,consideraprocess
involving diffusion and surface reaction in a heterogeneous (porous) medium at characteristic length scales
E
29 described in Fig. 1. Suppose that there are only two phases present in the particle; namely, the solid (k) and
fluid phase (g). The solid phase is assumed toTbe impermeable and bounded by a surface where a chemical
31 reaction takes place. Therefore, the chemical reaction occurs only at the surface of the solid phase. The
C
conservation equation that governs the transport process in the g-phase is given by
33 qc E
gþ=(cid:1)N ¼0 in the g(cid:2)phase, (1)
qt g
35 R
wherec andN aretheconcentrationandthemolarfluxofthereactantintheg-phase,respectively.Assume
g g
37 heterogeneous,firstorder,irRreversiblechemicalreactiononthecatalystsurface.Thecorrespondinginterfacial
O
39
C
41
N
43
U
45
47
49
51
Fig.1. Porousmediumandaveragingvolumewithcharacteristiclengths.
PHYSA : 9959
ARTICLE IN PRESS
F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]] 3
1 reaction–diffusion condition is given by
n (cid:1)N ¼kc at the gk(cid:2)interface, (2)
gk g g
3
where n is the normal surface vector directed from the g-phase toward the k-phase. For dilute solution
gk
5 conditions, the Fick’s constitutive equation
N ¼(cid:2)D =c (3)
g g g
7
has been used to describe dynamical molar flux at all spatial and time scales. This is an acceptable approach
for,e.g.,highporosityporousmediawherethetransportdynamicsaredominatedbythediffusionprocessin
9
the g-phase. Referred to Fig. 1, the Fick’s constitutive equation (3) describes the diffusion process at local
(pore)scalel .However,thissituationmaynolongerholdwhenoneconsidersscaleslargerthanl .Infact,as
11 g g
discussed above, the Fick’s law is implied from the assumption that the motion of diffusing particles in the
porousmediumisgeneratedbyasequenceofrandomjumpswithfinitevariance.However,theFick’slawcan
13
be unsuitable for describing particle paths within non-local scales at which experimental measurements are
made. In the case of the system described in Fig. 1, some non-Fickian diffusion phenomena due to particle
15
transport obstruction can arise when the disordered porous structure is accounted at spatial scales r bl .
0 g
Extendedirreversiblethermodynamicswasemployedtomodelthetransportinporousmedia[12].Thebasic
17
ideawastoraisethediffusionfluxtothestatusofindependentvariablesandtoviFewtheporousmediumasa
binary mixture formed by a perfectly rigid solid and a fluid. Following similar ideas, del Rio and Lopez de
19 O
Haro [13] obtained approximate time evolution equations for fluxes predicting a finite propagation velocity
within the porous medium. A tool, that is candidate to describe anomalous diffusion is based on fractional
21 O
calculus in that it provides a methodology to take into account these spatial correlations. Some fractional
23 extensionstotheFick’slawshavebeenproposedrecently.GorenfloRandMainardi[14]obtainedageneralized
diffusion equation in which the second-order spatial derivative is replaced by a pseudo-differential operator.
Chaves[15]showedthatfractionalFick’slawcanrecoverLe´vyPstatistics.Schumeretal.[3,4]derivedasimple
25
fractionalFick’slawinthecontextofhydrologicalmodels.Paradisietal.[16]obtainedafractionalFick’slaw
27 generating the Le´vy–Feller statistics. Solutions to advecDtion–dispersion equations with fractional Fick’s law
were discussed by Schumer et al. [3,4], showing that the model solution is able to resemble solute plumes in
granular aquifers. E
29
T
1.1. Problem statement
31
C
A vectorial version of the fractional Fick’s equation is given by Refs. [3,4,16]:
33 E
N ¼(cid:2)D =ac ; a2ð0;1(cid:3), (4)
g a;g g
35 R
where D is the fractional diffusion coefficient, and the fractional derivative operator is defined in the
a;g
Riemann–Liouville’ssense[17R].Inthelimita!1,theclassicalFick’sequationN ¼(cid:2)D =c isrecovered.A
37 g g g
combination of Eqs. (1) and (4) yields
O
39 qcg ¼=(cid:1)ðD =ac Þ in the g(cid:2)phase. (5)
qt a;Cg g
41 To obtain the full description of the reaction–diffusion problem at the spatial scale r bl , the following
0 g
N
interfacialconditionresultingfromthecombinationofEq.(2)andthefractionalFick’sconstitutivelaw(4)is
43 considered:
U
(cid:2)n (cid:1)ðD =ac Þ¼kc at the gk(cid:2)interface. (6)
45 gk a;g g g
It should be stressed that the fractional derivative in the diffusion equation leads to a super-diffusion
47 phenomena. It is also observed that the kinetics constant for the surface reaction leads to exponential
retardation(i.e.,notfractalorpower-lawscaling)sinceitisnotaffectedbythefractionalderivatives.Thatis,
49 surface chemical reaction is independent of fractal particle transport.
As it stands, the reaction–diffusion model given by Eqs. (5) and (6) is not useful for most practical
51 applications.Infact,thenature of catalyst reactivesystems consisting ofdisorderedporousstructures makes
almostimpossiblethederivationofasolutionfor,e.g.,analyzingordesigningcontaminantprocesses.Froma
PHYSA : 9959
ARTICLE IN PRESS
4 F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]]
1 practical standpoint, macroscopic models intended to describe changes at length scales Lbr , where L is a
0
given reacting system characteristic length at which typical field and industrial measurements are commonly
3 made,andr istheradiusoftheaveragingvolume.Inthisway,acorrectinterpretationofexperimentaland/or
0
field measurements require the model structure to be known for macroscopic scales. In principle, this should
5 also lead to estimations of transport parameters that are consistent with the scale considered.
Commonly, macroscopic models for reactive systems are derived from a mass-balance given by
7
qhci
þ=(cid:1)hNi¼k hci, (7)
qt eff
9
wherek andhciare,respectively,aneffectivekineticsconstantandacertainaverageconcentrationwithina
eff
11 spatial macroscale L. Notice that, while in the local model described above the chemical reaction is
heterogeneous, here the chemical reaction considered as pseudo-homogeneous. In a similar way, consider a
13 fractional Fick’s law equation for macroscale conditions given by
hNi¼(cid:2)D =ahci; a2ð0;1(cid:3), (8)
15 a;eff
whereD isaneffectivefractionaldiffusioncoefficientforthereactiveparticle.ThecombinationofEqs.(7)
a;eff
17 and (8) yields
F
qhci
19 qt ¼=(cid:1)ðDa;eff=ahciÞþkeffhci, O (9)
21 which, together with some suitable macroscopic boundary and initial coOnditions, govern the dynamics of the
averageconcentrationhci.Importantquestionsariseregardingthemacroscopicmodel(9):(i)Isthestructure
23 of the model (9) consistent with the structure of the detailed heteRrogeneous model (6), (7)? Notice that the
intuitiveidea behind thepseudo-homogeneousmodel (9) isthat it shouldretain, upto some spatialand time
P
25 scales,themaindynamicalcharacteristicsoftheheterogeneousmodel.(ii)Howthemacroscopictransportand
reactionparameters(Da;eff andkeff)arerelatedtotheloca lones(Da;g andk)?Theimportanceofthisquestion
27 relies on the fact that experimental and field data reveaDl scale-dependence of transport parameters. Focusing
on these questions, the aim of this paper is threefold:
E
29
(cid:1) ToderiveafractionalspatialaveragingtheoTrem(fSAT)asatooltoobtaintheeffectivemediumdynamics,
31 where the spatial scale r will be used as the characteristic length for volume averaging. Eventually, the
0 C
spatial scale r will become the scale at which averaging should be made. For convenience and clarity in
0
33 presentation, this issue is addressedEin the appendix.
(cid:1) To pose closure problems for the computation of the effective transport parameter D .
a;eff
35 (cid:1) To use spatial averaging methRods to demonstrate the consistency of the pseudo-homogeneous model (9)
relativetotheheterogeneousmodel(6),(7),wherethefractionalFick’sconstitutiveequationisassumedto
R
37 hold at spatial scales r bl , where r is significantly smaller than the macroscopic characteristic length L.
0 g 0
SufficientconditionsOonthespatialandtimescalesarederivedinorderthattheeffectivemediumdynamics
39 be reduced to the pseudo-homogeneous model (9).
C
41
Since transport models with fractional derivatives are being widely used in many areas of physics and
N
engineering, it is expected that our results will find applications in many different natural and application
43
areas. U
45
2. Effective medium dynamics
47
2.1. Preliminaries
49
Thelocalconcentrationisafunctionoftimetandpositionx,i.e.,c ðx;tÞ.Inthesequel,foreaseinnotation,
g
51 positionandtimefunctionalitieswillbeassumedimplicitly,sothatthesimplesymbolc willbeusedtodenote
g
the actual functionality c ðx;tÞ.
g
PHYSA : 9959
ARTICLE IN PRESS
F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]] 5
1 Consider an averaging volume V for the gk-system (see Fig. 1). Let V be the volume of the g-phase
g
contained within V.In this form,the volume fraction e ¼V =V isidentified asthe porosity.For any given
g g
3 spatial concentration distribution c , the application of the operator
g
1 Z
5 hcgi¼ cgdV
V
Vg
7 leads to the so-called superficial average concentration. Analogously, the intrinsic averaging operator is
defined by
9 1 Z
hc ig ¼ c dV.
g V g
11 g Vg
Inthisway,hc iandhc ig becometheaverageconcentrationassociatedtotheaveragingvolumeVandtheg-
g g
13 phase volume, respectively. Hence, the above operators are related in the following form:
hc i¼e hc ig.
15 g g g
When the averaging operator hc i is applied to the local conservation equation (5), it will be required to
g
17 interchange integration and differentiation in order to express the diffusive flux in terms of hc ig. For integer
g
derivatives, this is done by means of the so-called spatial averaging theorem (SATF) [18]:
19 1 Z O
h=c i¼=hc iþ c n dS,
g g g gk
V
21 Agk O
where A is the area of the gk-interface contained within V. However, in this paper we will need to
gk
23 interchange integration and fractional differentiation. In the AppeRndix, for a2ð0;1(cid:3), we have developed a
fSAT given by
P
25 1 Z
h=acgi¼a(cid:2)1=ahcgiþ Jð1(cid:2)aÞðcgÞngkdS,
V
27 Agk D
where J FðwÞ is the fractional integration operator [19]. Notice that the standard SAT is recovered in the
ð1(cid:2)aÞ E
29 limit as a!1.
T
31 2.2. Volume averaging of g-phase dynamics
C
33 A spatial average over a control volEume element will be carried out to obtain a smooth description of the
local non-uniformities. To this end, the following assumptions are made:
35 R
(cid:1) The catalyst particle has a characteristic length L (see Fig. 1).
R
37 (cid:1) The volume averaging procedure is applied over a volume of radius r 5L.
0
(cid:1) The dominant porousOcharacteristic length is equal to l 5r .
g 0
39 (cid:1) Onthebasisofthisdisparityoflengthscales,spatialvariationsofthediffusion parameterD ,relaxation
a;g
time-constant taCand kinetics parameter k, are neglected in the averaging volume.
r
41
N
The application of the superficial averaging operator to Eq. (5) gives
43
(cid:1)qcU(cid:2)
g ¼h=(cid:1)ðD =ac Þi.
45 qt a;g g
SinceV isnotafunctionoftime,theaveragingandtimederivativeoperatorscanbeinterchangedtoreduce
g
47
the above equation into
qhc i
49 g ¼h=(cid:1)ðD =ac Þi. (10)
qt a;g g
51 Inordertointerchangeintegrationandintegerdifferentiationintheright-handsideofEq.(10),thestandard
SAT is applied to obtain
PHYSA : 9959
ARTICLE IN PRESS
6 F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]]
1 qhcgi¼=(cid:1)ðD h=ac iÞþ 1 Z n (cid:1)ðD =ac ÞdS. (11)
qt a;g g V gk a;g g
Agk
3
The use of the gk-interfacial condition (see Eq. (6)) in the integral term of the above equation yields
5 qhc i 1 Z
g ¼=(cid:1)ðD h=ac iÞ(cid:2) kc dS
qt a;g g V g
Agk
7
¼=(cid:1)ðD h=ac iÞ(cid:2)ka hc i ,
a;g g v g gk
9 where hc i is the area averaged concentration defined as
g gk
1 Z
11 hcgigk ¼ cgdS,
V
Agk
13 and a is the area per unit volume given by
v
A
15 av ¼ gk.
V
At this point, we use the fSAT to interchange integration and fractional differentiation to obtain
17
" !# F
qhc ig 1 1 Z
19 eg qtg ¼=(cid:1) Da;g aðeg=ahcgigþhcgig=aegÞþV ngkJð1(cid:2)aÞðcgÞdS O (cid:2)kavhcgigk, (12)
Agk
21 where we have used the fact that hc i¼e hc ig. In this description, thOere are two different concentrations,
g g g
namelytheintrinsicaveragehc ig andthepointc concentration.Thepresenceofthepointconcentrationc in
g g g
23 theintegraltermR n J ðc ÞdAdoespresent aproblemsinceRit isonlyavailable intermsofthesolution
ofthepointboundAagrkygvkalðu1(cid:2)eapÞrogblemgivenbyEqs.(5)and(6).Sincetheintrinsicaverageconcentrationhc ig
P g
25 isthevariableofpracticalinterest,inasimilarwaytotemporaldecompositionusedinthestudyofturbulence,
the point concentration is decomposed as the sum of th e average concentration and the spatial deviations.
27 That is, D
c ¼hc igþc , E (13)
29 g g eg
where ecg denotes the spatial concentration dTeviation. Eq. (13) represents a decomposition of length scale
31 concentrationvariations.Whenthelengthconstraintl 5r 5L issatisfied,theaverageconcentrationcanbe
C g 0 c
handled constant into the averaging volume scale r and undergoes significant changes only over a certain
0
33 catalyst particle lengthscale Lc.OnthEeotherhand,thedeviationsecg aredominatedbythesmall lengthscale
R
l . When Eq. (13) is used in the integral term n J ðc ÞdS, the following expression is obtained:
35 g R Agk gk!ð1(cid:2)aÞ g
1 Z 1 Z 1 Z
n J ðc ÞdS ¼ n J ð1ÞdS hc igþ n J ðc ÞdS,
gk ð1(cid:2)aÞ g R gk ð1(cid:2)aÞ g gk ð1(cid:2)aÞ eg
37 V Agk V Agk V Agk
where we have used theOproperty R n J ðhc igÞdS ¼ðR n J ð1ÞdSÞhc ig. This approximation can
39 be taken on the basis that l 5r anAdgkr2g5k Lð1L(cid:2)a.Þ Ong the otherAhgkangdk, tðh1(cid:2)eaÞfSAT impglies that
C g 0 0 e
41 1 Z n J ð1ÞdS ¼h=a1i(cid:2)1=ah1i.
V Ngk ð1(cid:2)aÞ a
Agk
43 It is notUcomplicated to show that =a1¼0 and h1i¼ð1=VÞR 1dV ¼e . In this way, one can write
Vg g
45 1 Z 1
n J ð1ÞdS ¼(cid:2) =ae .
V gk ð1(cid:2)aÞ a g
Agk
47
Therefore, Eq. (12) is reduced into the following one:
" !#
49 qhc ig e 1 Z
eg qtg ¼=(cid:1) Da;g ag=ahcgigþV ngkJð1(cid:2)aÞðecgÞdS (cid:2)kavhcgigk. (14)
51 Agk
For the moment, let us take the approximation
PHYSA : 9959
ARTICLE IN PRESS
F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]] 7
1 hc i (cid:4)hc ig (15)
g gk g
inaneighborhoodofthegk-interface.Ajustificationofthisapproximationwillbemadelater.Eq.(15)allows
3
to write Eq. (14) as follows:
5 qhc ig (cid:3)D (cid:4) (cid:3)D (cid:4)
g ¼=(cid:1) a;g=ahc ig þe(cid:2)1=e (cid:1) a;g=ahc ig
qt a g g g a g
7 !
Z
D
þe(cid:2)1=(cid:1) a;g n J ðc ÞdS (cid:2)e(cid:2)1ka hc ig. ð16Þ
g V gk ð1(cid:2)aÞ eg g v g
9 Agk
Summingup,theoriginalboundary-valueproblem(5)–(6)hasbeenreducedtotheaveragedonerepresented
11
by Eq. (14). However, two different concentration representations are present; namely, the intrinsic volume
averaged concentration hc ig, and the spatial deviation concentration c . As hc ig is the concentration of
g eg g
13 practical interest, one should provide a procedure for representing the averaged dynamics in terms of hc ig
g
solely.
15
2.3. Closure problem
17
F
Inprinciple,thec -fieldcanbeobtainedbysolvingsimultaneouslytheboundaryvalueproblems(5),(6)and
19 eg O
(16). Since catalyst particles are disordered media, the exact determination of the deviation c is not an easy
eg
task. In order to develop a procedure to estimate the effects of the deviation dynamics on the average
21 O
concentration, it is convenient to substract Eq. (16) from the heterogeneous model Eq. (5) to obtain an
23 equation governing the dynamics of the deviations ecg in the g-phasRe; namely,
(cid:3) (cid:3) (cid:4) (cid:4) (cid:3) (cid:4)
qc 1 D
25 qetg ¼=(cid:1)ðDa;g=aecgÞþ=(cid:1) Da;g 1(cid:2)a =ahcgig (cid:2)e(cid:2)g1=Peg(cid:1) aa;g=ahcgig
" #
Z
27 (cid:2)e(cid:2)1=(cid:1) Dg n J ðc ÞdA þe(cid:2)1ka hDc ig; in the g(cid:2)phase. ð17Þ
g V gk ð1(cid:2)aÞ eg g v g
Agk E
29
Fromanorderofmagnitudeanalysis[11],somesimplificationsofEq.(17)arepossibleoncelengthconstraints
T
l 5r 5L are accepted. In fact, the order of magnitude estimates of the first and fourth terms in the RHS of
31 g 0
the above equation are C
33 D c ! E
=(cid:1)ðD =ac Þ¼O a;geg (18)
a;g eg laþ1
35 g R
and
R
37
" # !
D Z D c
e(cid:2)1=(cid:1) g n JO ðc ÞdA ¼O a;geg , (19)
39 g V gk ð1(cid:2)aÞ eg Lla
Agk g
C
respectively. Therefore, if the following length-scale constraint is satisfied
41
N
l
g51, (20)
43
L U
45 tdhifefunsoivne-ltoecraml d=if(cid:1)fuðsDive=tearcmÞ,ea(cid:2)gn1d=E(cid:1)q½ð.D(1g=7V) cÞaRnAgbkengskiJmð1p(cid:2)liafiÞðeecdgÞidnAto(cid:3) tchaenfboellonweignlegctoende:with respect to the local
a;g eg
(cid:3) (cid:3) (cid:4) (cid:4) (cid:3) (cid:4)
47 qc 1 D
eg ¼=(cid:1)ðD =ac Þþ=(cid:1) D 1(cid:2) =ahc ig (cid:2)e(cid:2)1=e (cid:1) a;g=ahc ig
qt a;g eg a;g a g g g a g
49 þe(cid:2)1ka hc ig; in the g(cid:2)phase. ð21Þ
g v g
51 The interfacial boundary condition for this equation is obtained by using the decomposition c ¼hc igþc
g g eg
into Eq. (6):
PHYSA : 9959
ARTICLE IN PRESS
8 F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]]
1 (cid:2)n (cid:1)ðD =ahc igÞ(cid:2)n (cid:1)ðD =ac Þ¼khc igþkc ; at the gk(cid:2)interface. (22)
gk a;g g gk a;g eg g eg
After performing an order of magnitude estimate on Eq. (22), one obtains
3
"la=Laþkla=D #
c ¼O g g a;g hc ig, (23)
5 eg 1þkla=D g
g a;g
7 in a neighborhood of the gk-interface.
Thetermkla=D canbeexpressedintermsoftheso-calledThielemodulustoshowthatinmostpractical
g a;g
9 cases klag=Da;g51 [11]. If in addition ðlg=LÞa51, which is implied by the inequality (20), one obtains that
c 5hc ig; in a neighborhood of the gk(cid:2)interface.
eg g
11
This simplifies the boundary condition (22) into the following one:
13 (cid:2)n (cid:1)ðD =ahc igÞ(cid:2)n (cid:1)ðD =ac Þ¼khc ig; at the gk(cid:2)interface. (24)
gk a;g g gk a;g eg g
Notice that the fact that c 5hc ig allows Eq. (15) to be valid. The order of magnitude of the first term in the
15 eg g
LHS of the above equation gives
17 1 Z n (cid:1)ðD =ahc igÞdS ¼O(cid:3)Da;ghcgig(cid:4).
V gk a;g g l La F
Agk g
19 O
On the other hand, the order of magnitude analysis of the diffusive terms in Eq. (21) gives the following
relationships:
21 O
1 Z (cid:3) (cid:3) 1(cid:4) (cid:4) (cid:3)D hc ig(cid:4)
=(cid:1) D 1(cid:2) =ahc ig dV ¼O a;g g (25)
23 V a;g a g Laþ1 R
Vg
and P
25
1 Z e(cid:2)1=e (cid:1)(cid:3)Da;g=ahc ig(cid:4)¼O(cid:3)Da;ghcgig(cid:4), (26)
27 V g g a g L La D
Vg e
where L is the characteristic length scale of the gE-phase. In addition, one can treat the closure problem as
29 e
quasi-steady state [11] when
T
31 qc
eg5=(cid:1)ðD =ac Þ. C
qt a;g eg
33 This condition will be satisfied wheneEver the characteristic time, say t , is large enough so that the following
c
constraint is valid:
35 R
D t
a;g cb1.
37 l1þa R
g
When the macroscopicOprocess is transient, the characteristic time t should be constrained by
c
39
D t
a;g c ¼Oð1ÞC.
41 L1þa
However, byNhypothesis Lbl , which implies D t =l1þa51 can be assumed. Therefore, if
g a;g c g
43
l U
g 51, (27)
L
45 e
andtheconstraintgivenbyEq.(20)ismet,thenEq.(21)thatgovernstheconcentrationspatialdeviationcan
47 be approximated as follows:
ka
49 =aþ1ecg ¼ (cid:2)e Dv hcgig (28)
g a;g
|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}
51 Volume reactive source
with boundary condition
PHYSA : 9959
ARTICLE IN PRESS
F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]] 9
1 (cid:2)n (cid:1)ðD =ac Þ¼n (cid:1)ðD =ahc igÞþ khc ig ; at the gk(cid:2)interface. (29)
gk a;g eg gk a;g g g
|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflffl{zfflffl}
Diffusive surface source Reactive surface source
3
In this form, the boundary value problem for the deviation c is given by Eqs. (28) and (29). The following
eg
5 comments are in order: (i) The idea of constructing the above model is not to solve it over the entire
macroscopicregion(seeFig.1).Instead,onewantstosolveforc insomerepresentativeregioninthevolume
eg
7 averaginglengthscaler (seeFig.2).Thisisreasonablesincethescaleofthec -sources(hc ig and=ahc ig)are
0 eg g g
oftheorderofthesmallcharacteristiclengthl .Letusassumethattheporousmediumisdisorderedandany
g
9 sample at the length scale r is representative of the porous medium structure. This allows us to consider the
0
porousmediumasspatiallyperiodic(seeFig.2).(ii)Forhc ig ¼0,theuniquesolutionoftheboundaryvalue
g
11 problem (28), (29) is c ¼0. As the average concentration departs from zero, a non-trivial concentration
eg
deviation is generated. In this form, the ‘‘sources’’ in the right-hand side of Eqs. (28) and (29) (i.e., hc ig and
g
13 =ahc ig) determine the behavior of c . These observations, together with the linear nature of the boundary
g eg
value problem, suggest the following superposition solution:
15
c ¼b (cid:1)=ahc igþs hc ig, (30)
eg g g g g
17
wherethevectorb andthescalars arereferredasclosurevariables.SubstitutionFoftheproposedsolutionin
g g
19 (28)and(29)alongwiththeperiodicityconditiongivethefollowingboundaryOvalueproblemsfortheclosure
variables:
21 O
23 (cid:1) Boundary value problem for bg: R
=(cid:1)ðD =ab Þ¼0; in the g(cid:2)phase
a;g g P
25
(cid:2)n (cid:1)=ab ¼n ; at the gk(cid:2)interface
gk g gk
27 bgðrþliÞ¼bgðrÞ;i¼1;2;3ðperiodicityÞ. D
E
29 (cid:1) Boundary value problem for s :
g
T
31 =(cid:1)ðDa;g=asgÞ¼(cid:2)e(cid:2)g1kav; in the g(cid:2)phase
C
(cid:2)n (cid:1)ðD =as Þ¼k; at the gk(cid:2)interface
gk a;g g
33 s ðrþlÞ¼s ðrÞ;i¼1;2;3ðperEiodicityÞ.
g i g
35 R
R
37
O
39
C
41
N
43
U
45
47
49
51
Fig.2. Periodicunitcellforarepresentativeregionoftheporousmedium.
PHYSA : 9959
ARTICLE IN PRESS
10 F.J.Valdes-Paradaetal./PhysicaA](]]]])]]]–]]]
1 In the above descriptions, r represents the center of the representative porous medium cell, and the l’s
i
represent vector that are required to represent a spatially periodic porous medium. This means that the
3 geometry of the representative cell is invariant to a transformation of the type r!rþl.
i
Theclosureproblemissolvedafterobtainingthesolutionoftheabovedescribedboundaryvalueproblems
5 for the closure variables. Substitution of Eq. (30) into Eq. (16) leads to a closed form for the governing
equation for the averaged concentration hc ig:
g
7 " # !
qhc ig 1 1 Z
9 eg qtg ¼=(cid:1) egDa;g aIþVg AgkngkJð1(cid:2)aÞðbgÞdS (cid:1)=ahcgig
( ! )
1 Z
11 þ=(cid:1) Da;g ngkJð1(cid:2)aÞðsgÞdS hcgig (cid:2)kavhcgig. ð31Þ
V
Agk
13
Notice that in the above dynamics the only variable is the averaged concentration hc ig.
g
15
3. Discussion
17
Thebehavioroftheaverageconcentrationhc ig isaffectedbythesolutionofthFeclosurevariables.Inturn,
g
one obtains the pseudo-homogeneous model (30), which is more complex than the expected pseudo-
19 O
homogeneous model given by Eq. (9). To obtain a model with a structure similar to that in Eq. (9), the
following simplifications should be made:
21 O
23 (cid:1) Theeffectivekineticsparameterisgivenbykeff ¼avk.NoticethaRtthekineticscoefficientisaffectedonlyby
the averaging process, and not by the fractional flux constitutive law.
25 (cid:1) The term =(cid:1)fDa;gðð1=VÞRAgkngkJð1(cid:2)aÞðsgÞdSÞhcgigg shouPld be neglected. In fact, from an order of
magnitude analysis, it is obtained that when the length scale constraint
27 l D
g51
L E
29
is met, the following inequality is satisfied:
T
" ! #
31 1 Z (cid:6)k (cid:7)
=(cid:1) Da;g V ngkJð1(cid:2)aÞðsgÞdS Chcgig 5O Lhcgig .
33 Agk E
In this case, the boundary value problem for the closure variable s needs not to be solved.
g
35 (cid:1) An effective diffusion parametRer tensor can be defined as
" !#
1 R1 Z
37 Da;eff ¼Da;g aIþ V ngkJð1(cid:2)aÞðbgÞdS . (32)
O g Agk
39
The above considerations lead to the following pseudo-homogeneous model for the reaction–diffusion
C
dynamics:
41
e qhcNgig ¼=(cid:1)ðe D (cid:1)=ahc igÞ(cid:2)k hc ig. (33)
43 g qt g a;eff g eff g
U
Thestandardvolumeaveragedreaction–diffusionmodel[11]isrecoveredinthelimitasa!1.Noticethat
45 D is composed by two parts: (i) a(cid:2)1D I, which is the contribution of the local transport mechanisms.
a;eff a;g
The factor a(cid:2)1 arises as a consequence of the infinite-variance of the velocity field (represented in the
47 fractional Fick’s equation) which, when averaged for macroscale conditions, produces the effects of an
enhanced fractional diffusion coefficient. Notice that such effect is not present in Fickian diffusion where
R
49 a¼1,suchthatlocaltransportpropertiesareinvariantunderscaling.(ii)Da;gðð1=VgÞ AgkngkJð1(cid:2)aÞðbgÞdSÞ,
which reflects the effects of the disordered porous media. Notice that, contrary to Fickian diffusion, the
51 fractalityoftheporousmediahasanimportanteffectontheeffectivefractionaltransporttensorD .In
a;eff
the case of an isotropic porous medium (i.e., D ¼D I) with =e (cid:4)0 in the volume averaging scale,
a;eff a;eff g
