3B2v8: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