Improved PISO algorithms for modeling density varying flow in conjugate fluid–porous domains

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Improved

PISO

algorithms

for

modeling

density

varying

flow

in

conjugate

fluid–porous

domains

M. Nordlund

a,

∗

,

M. Stanic

b

,

A.K. Kuczaj

a,b

,

E.M.A. Frederix

b

,

B.J. Geurts

b,c a_{PhilipMorrisInternationalR&D,PhilipMorrisProductsS.A.,QuaiJeanrenaud5,2000Neuchâtel,Switzerland}

b_{MultiscaleModelingandSimulation,FacultyEEMCS,J.M.BurgersCenter,UniversityofTwente,P.O.Box217,7500AEEnschede,}

TheNetherlands

c_{AnisotropicTurbulence,FluidDynamicsLaboratory,FacultyofAppliedPhysics,EindhovenUniversityofTechnology,P.O.Box513,5600MB}

Eindhoven,TheNetherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory: Received11June2015

Receivedinrevisedform6November2015 Accepted12November2015

Availableonline21November2015

Keywords: Porous PISO Rhie–Chow Spuriousoscillations Collocated Segregated

Finitevolumemethod

Two modifiedsegregatedPISOalgorithms are proposed,whichare constructedto avoid the development of spurious oscillations in the computed flow near sharp interfaces ofconjugate fluid–porousdomains. Thenew collocatedfinitevolumealgorithms modify theRhie–Chowinterpolationtomaintainacorrectpressure–velocitycouplingwhenlarge discontinuous momentumsources associated with jumps inthe local permeability and porosityarepresent.TheRe-DistributedResistivity(RDR)algorithmisbasedonspreading flowresistivityoverthegridcellsneighboringadiscontinuityinmaterialpropertiesofthe porousmedium.TheFaceConsistentPressure(FCP)approachderivesanauxiliarypressure value atthefluid–porous interfaceusingmomentumbalancearound theinterface. Such derivedpressurecorrectionisdesignedtoavoidspuriousoscillationsaswouldotherwise arise with a strictly central discretization. The proposed algorithms are successfully compared against published data for the velocity and pressure for two reference cases ofviscous flow.Therobustnessoftheproposed algorithms isadditionallydemonstrated forstronglyreducedviscosity,i.e.,higherReynoldsnumberflowsandlowDarcynumbers, i.e.,low permeabilityofthe porousregions inthe domain, forwhichsolutions without unphysicaloscillationsarecomputed.BothRDRandFCParefoundtoaccuratelyrepresent porousmediaflowneardiscontinuitiesinmaterialpropertiesonstructuredgrids.

©2015ElsevierInc.All rights reserved.

1. Introduction

Fluidflowthroughporous mediaisfundamental tomanynaturalandindustrialprocesses,such asgroundwaterflows, filtration, and chemical andbiomass processing[1–5]. In order to efficiently simulate theseprocesses and predict their characteristicproperties,robustandaccuratenumericalmodelsareofhighimportance.Whiletheequationsgoverningthe flowandheattransferinporousmediaarereadilyspecifiedusingthemethodofvolumeaveraging

[6]

inconjunctionwith suitable closure models, it often remains challenging to obtain physically acceptable numerical solutions in the vicinity of fluid–porousinterfaces. Without specialcare in thealgorithm development, discontinuities inthe material properties, particularlythe resistivity,at theinterface mayyield spurious oscillationsinthe solution.Thisis especiallytrue forhigh

*

Correspondingauthor.

E-mailaddress:markus.nordlund@pmi.com(M. Nordlund). http://dx.doi.org/10.1016/j.jcp.2015.11.035

Reynolds(Re)numberflowsatlowDarcy(Da)number,forwhichthejumpinflowresistivityishigh,whiledampingeffects ofviscosityare ratherlow.The occurrenceofspurious oscillationsisparticularlypronouncedwhen segregatedalgorithms areapplied,inwhichthevelocityandpressureequationsaresolvedseparatelyandaniterativesolutionprocessisrequired. Numerical schemes that avoidspurious oscillations at sharpinterfaceswere proposed by [7]and[8]for bothstructured andunstructuredgridsusingacollocatedvariable,finitevolumeblock-coupledsolver.Inthiscasepressureandvelocityare solvedsimultaneously.Theserathercomplexschemes requirespecial,localtreatmentofthefluid–porousinterfaces.Inthis paper we turntosegregatedsolvers, which aremore commonlyadoptedinComputational FluidDynamics(CFD) studies andpresentalternativetreatmentsthatareequallywellsuitedtoalleviatetheproblemofspuriousoscillations.Wecompare theRDR

[9]

andFCPmethodsinthecontextofthePressureImplicitwithSplittingofOperators(PISO)algorithmandassess theirrelativeperformanceinsomedetailforseveraltestcasesofincreasingcomplexity.ItwillbeshownthatbothRDRand FCPapproachesareeffectiveineliminatingspuriousoscillations.

AcorrectiontotheRhie–Chowinterpolationincasethefluidissubjecttolargebodyforceswassuggestedby

[10]

anda force fielddiscretizationruleforthevolume-of-fluidmethodwassuggestedby

[11]

asremediesfortheoccurringspurious velocity oscillations.Thesecorrection schemesmaybe appliedforexplicittreatment oftheflow resistancetermspresent forporousmediaflow.However,duetothestiffnessoftheequationssuchexplicittreatmentsuffersfromverysmalltime stepswhichrenderthemethodlesseffective.

Inthispapertwomodified Rhie–Chow/PISOalgorithms areproposed, whichareconstructedtoavoidthedevelopment of spurious oscillations inthe solution when large discontinuous momentum sources are present. The RDR algorithm is based onaredistributionofthediscontinuous flowresistivityovergrid cellsadjacenttothediscontinuityinthematerial properties,ina similarwayaswasproposed by

[12]

fordiscontinuous bodyforce fields.Inordernottocreatetime step restrictionsforhighRe,lowDa numbers,forwhichtheresistanceishigh,thistermistreatedimplicitly.TheFCPalgorithm usesanapproachsimilartothatpresentedbyOxtobyetal.

[11]

,whereinterfacepressurevaluesarederivedbyensuringa force balancebasedonthemomentumequationsonbothsidesofthefluid–porousinterface. Thisisaimed ateliminating spurious oscillations [12]. Differently fromOxtobyetal.

[11]

the FCPmethod treatsall terms,except pressure,implicitly, thusrelaxingthetimestepconstraintimposedbyhighRe/lowDa numbers.

Throughout this paper we adopt a numerical formulation, which solves the equations directly forthe pressure field, rather than the densityfield asis commonfor compressible flows. The selected formulationis computationally cheaper than thedensity-based formulationfordensityvarying flows,anddoesnot requireresolutionof highfrequencypressure fluctuations(i.e.,soundpropagation)

[13,14]

.ThemethodsareimplementedintheopensourceOpenFOAM



environment, whichisbasedoncollocatedgriddiscretization,storingallvariablefieldsatcellcenters

[15]

.

The proposed algorithms are successfully compared with published data for the velocity and pressure fields for in-compressible, isothermal flow through a porous plug [7] and for flow parallel to a porous region, i.e., the so-called Beavers–Josephproblem

[16]

.Theseproblemsareusedasbenchmarks, sincetheir particulargeometryisknowntocause oscillatorysolutionsincasestandardinterpolationmethodsareusedfortheresistivityterm.Theso-calledporousresistance termin themomentum equation,which isthemajor sourceofthe discontinuityinthe domain,isa functionofvelocity andpermeability.ThetworelevantflowparametersareRe andDa numbers.Atfixedlength-scaleandfluidproperties,high velocitycorrespondstohighRe numberandlowpermeabilitycorrespondstolowDa number.Anincrease/decreaseofthese parameters isdirectlyinfluencingthemagnitudeoftheporousresistanceterm. Therobustness oftheproposedalgorithms isdemonstratedforbothhighRe numberflowsandlowDa numbers,whichconstitutesthemostchallengingsituation.RDR isfoundtosuccessfullyremovespuriousoscillationscomparedtothereferencemodelwhichadoptscentraldiscretization, forall testedvaluesofRe and Da numbers.FCPisfoundtobeadequate formostcasesinvolvinglower Re andhigherDa

numbers,whileforhigherRe andlowerDa numberssomeoscillationscanstillbeseeninthesolution.

Spurious numerical oscillations ofvelocity oranyother variable can havea bigimpact on delicatephysical processes such asaerosolformationandevolution.Evenslightvariations inthermodynamicvariablescanhavealarge influenceon, e.g.,thenucleationrate,whichisatthecoreoftheformationofaerosoldroplets.Suchstrongdependenciesonrathersubtle changesinthelocalconditionsemphasizestheneedforanaccuratetreatmentoffluid–porousinterfaceswhereoscillations mayoccur.Thispaperprovidessimpleandeffectivenumericalmethodstosuppresstheaforementionedoscillationsalsoin caseonereliesonsegregatedsolverswithacollocatedarrangementofvariables.

The organizationofthispaperisasfollows.InSection2,thegoverningequationsforweaklydensityvaryingflowand heattransferinconjugatefluid–porousdomainsareintroduced.InSection3thefinitevolumediscretizationisdescribedand an analysisofthediscretizedequationsforlarge,discontinuoussourcetermsisgiven.Thereafter,Re-DistributedResistivity (RDR) andFace ConsistentPressure(FCP) algorithms are putforwardin Sections4and 5followed bytheir validationin Section6andconclusionsinSection7.

2. Governingequations

Fluid flow andconjugate heat transfer influid–porous domains can be described by the volume-averaged mass, mo-mentumandenergyconservationequations

[6]

.Imagineavolume V ,occupiedbytwoconstituents

α

(fluid)and

β

(solid), such that V

=

V α

+

Vβ. Then the extrinsic, orsuperficial, average andthe intrinsic volume average of a fluid property

ϕ

aredefinedas:



ϕ

α



=

1V



Vα

ϕ

αdV and



ϕ

α



α

₌

1 Vα



Vα

ϕ

αdV , respectively

[6]

.The



bracketsimplyvolume averaging overvolume V andforintrinsicpropertiestheexponentabovethebrackets



α or



β _{impliesvolumeaveragingonlyover}

volume V α or Vβ,respectively.The relationshipbetweenintrinsicandextrinsicvolumeaveraged propertiesisdefinedas



ϕ

α



α

= φ

ϕ

α



,where

φ

=

V α

/

V istheporosityoftheporousmedium.

Followingthederivationsin

[6]

andtheclosuremodelingin

[17]

and

[18]

,theequationsgoverningtheflowoffluidand heatinisotropic,heterogeneousporousmedia,consistingofafluidphase

α

anda solidphase

β

arepresented.Themass conservationforfluidflowis

∂

t

(φ



ρ

α



α

)

+ ∂

i

(



ρ

α



α



ui

) =

0 (1)

wheret is thetime,

ρ

the densityand ui the velocityin thexi direction. Because onlyweakdensityvariations are con-sidered, the mass dispersion term originatingfrom the volume averaging procedure is neglected. The operator

∂

t isthe temporalpartialderivativeand

∂

i isthepartialderivative withrespecttothespatialcoordinatexi.Weadoptamixed for-mulationinwhichtheextrinsicallyaveragedvelocityfield



ui



transportstheintrinsicallyaveragedfluidmassdensity



ρ

α



α

inthevolumeV .Themomentumequationforfluidflowis

∂

t

(



ρ

α



α



ui

) + ∂

j

(φ

−1



ρ

α



α



uj



ui

) = −φ∂

i



pα



α

+ ∂

j



τ

i j

 + φζ

i



α

− φ

D



ui



(2) withthevolume-averagedrateofstraintensor



τ

i j

 = 

μ

α



α

(∂

j



ui

 + ∂

i



uj

) − (

2 3



μ

α



α

_{− }

_κ

α



α

)δ

i j

∂

k



uk



(3)

InEqs.

(2)

and

(3)

p isthepressure,

ζ

iamomentumbodysourceincludingalsohigher-orderperturbationtermsoriginating from thevolume averaging procedure. The dynamic viscosity is denoted

μ

,

δ

i j is the Kronecker delta tensorand

κ

the dilatationalviscosity.The key termforthispaperisthe isotropicporousresistanceterm

φ

D



ui



,thelast terminEq.(2), whichtakesafinitevalueintheporousmediaandzerootherwise.Thisporousresistancetermisthesourceofdiscontinuity inthedomain,andthustheprimarycauseofspuriousoscillations,aslaterwillbedescribedindetail.Itisdefinedasascalar

D

= 

μ

α



α

(

K−1

+

K−1F

)

,whereK istheisotropicpermeabilityscalarandF

=

ρα

α

μααK

1/2

_|

_u

j

|

cE thenon-Darcyresistivity, inwhichcE istheformdragcoefficient,that togetherwiththepermeabilityaccountsforthemicroscopicstructureofthe porousmedia.

Thetemperatureequationsforthefluidandsolidconstituentsarerespectively:



cp,α



α



∂

t

(φ



ρ

α



α



Tα



α

)

+ ∂

i

(



ρ

α



α



ui



Tα



α

)



= ∂

i

(λ

eα

∂

i



Tα



α

)

+

hαβsαβ

(



Tβ



β

− 

Tα



α

)

+ ∂

t

(φ



pα



α

)

+ 

ui

∂

i



pα



α

+ ˆ

α (4)



cp,β



β



ρ

β



β

∂

t

((

1

− φ)

Tβ



β

)

= ∂

i

(λ

eβ

∂

i



Tβ



β

)

−

hαβsαβ

(



Tβ



β

− 

Tα



α

)

(5)

wherecp,misthespecificheatatconstantpressure,Tmthetemperature,

λ

emtheeffectivethermalconductivityforthephase m

∈ {

α

, β

}

,includingboththestagnanteffectiveconductivityandtortuosityeffects.Theinterfacialheattransfercoefficientis denotedbyhαβandsαβ isthespecificsurfacearea(surfaceavailableperunitofporousmediavolume).Thermaldispersion

isdenoted

ˆ

α .In Eqs.(2)to

(5)

the closuretermssuggestedin

[17,18]

wereadoptedandlocalthermal non-equilibrium

betweenthetwophaseswasassumed.

Thisis toillustrate that theso-calledone-domain approach [19] istaken, solving the described equationsforthe en-tire fluid–porous domain, rather than splitting the domain into parts which are governed by multiple sets ofequations and imposing interfacial boundary conditions in places where they merge. In the next section the discretization of the volume-averagedgoverningequationsforthefluidflowinconjugatefluid–porousdomainsispresented.

3. FinitevolumediscretizationandmodifiedRhie–Chowinterpolation

Application ofthe FiniteVolume (FV)method requiresflux evaluationat cellfaces, andinterpolation ofcell centered variables onto the cell faces. Segregated algorithms for FV methods, such asPISO [20], are based on either iterativeor stepwiseproceduresforobtainingthesolution.PISOconsistsoftwomainsteps:thepredictorandthecorrectorstep.Both ofthese steps requireevaluationof pressuregradients at cellfaces for calculatingpressure gradient values.An incorrect evaluationofthesepressuregradientsandpressurefacevaluesleadstopressure–velocitydecouplingorthe‘checker-board’ pressurepattern,whichisawell knownproblemoncollocatedgrids

[21,15]

.In1983, RhieandChowproposeda method whichmaintainsthepressure–velocitycouplingoncollocatedgrids evenwhenusingsegregatedsolvers

[21]

.Thismethod istodayknownastheRhie–Chowinterpolationandwill bedescribed insomedetailinthissectionto provideapointof referenceforunderstandingthemodificationsinRDRandFCPaspresentedmomentarily.

The Rhie–Chow interpolation was developed assuming governing equations without discontinuous body forces acting on the flow. As shownin Section 2,this is not the case when porous media are present in part ofthe domain. When body force discontinuities are present, in particular in the porousresistance term

φ

D



ui



,Rhie–Chow interpolation fails tomaintainthepressure–velocity couplingresultingin‘velocityripples’,aswasdiscussedearlierbyMencingeretal.

[12]

. Theexact mechanismbehindits failurewill alsobepresentedhere. Inthissectionwe willpresentamodified version of the Rhie–Chowinterpolation which doespreserve pressure–velocity couplingeven when large body force discontinuities are present inthe flow. In orderto discretize the governing equations using a collocated FV method,the equations are

Fig. 1. ExamplecellofacollocatedgridforFVanalysis.TheownercellisdenotedbyP ,theneighborcellbynb.Thevectordi,Pnbisthedistancevector

betweenthenodeP andthenodenb.ThecommonfaceofcellsP andnb isdenotedby f .Thesurfacenormalvectorofthefacef isdenotedSi,f;this

vectorhasamagnitudeequaltothesurfaceareaofthefacef .ThevectorbetweenthecellcenterP andfacecenterisdenotedwithdi,Pf.Theindexi

standsfortheithcoordinatedirection.

integratedoveracontrolvolumeVP centeredaroundanode P andboundedby N faceswithsurface areavector Si,f and facecenters f ,asshownin

Fig. 1

.

Applying Gauss’ divergence theorem and using Einstein’s tensor notation, the semi-discretized form of the volume-averagedmass,momentumandenergyequationsforacollocatedvariablediscretizationyields:

VP

∂

t

(φ

P

ρ

P

)

+



f

(

ρ

ui

)

fSi,f

=

0 (6) VP

∂

t

(

ρ

Pui,P

)

+



f

(φ

−1

ρ

uj

)

fui,fSj,f

= −φ

P



f pfSi,f

+



f

τ

i jSj,f

+ φ

P

ζ

i,PVP

− φ

PDPui,PVP (7) VPcp,P

[∂

t

(φ

P

ρ

PTP

)

+



f

(

ρ

ui

)

fTfSi,f

]

=



f

λ

e_f

(∂

iT

)

fSi,f

+

VPhαβ,Psαβ,P

(

Tβ,P

−

TP

)

+

VP

∂

t

(φ

PpP

)

+

ui,P



f pfSi,f

+

VP

ˆ

P (8) VP

(

cp,β

ρ

β

)

P

∂

t

((

1

− φ

P

)

Tβ,P

)

=



f

λ

e_β,f

(∂

iTβ

)

fSi,f

−

VPhαβ,Psαβ,P

(

Tβ,P

−

TP

)

(9) Here andinthesequelthe spatialaveraging operators



α and



andthesubscript

α

forthefluid phasearedropped forbrevity, suchthat, i.e.,



Tα



αP

=

TP and



Tβ



βP

=

Tβ,P.Notethat ui refers tothesuperficial velocity



ui



andtheother variablesareintrinsic.

Inthediscretizedequationsabove,someofthepropertiesandvariablesarerequiredatthecellcentersandsomeatthe facecenters.Inacollocatedvariablerepresentation,thepropertiesandvariablesareallstoredinthecell-centersandneed to beinterpolatedby interpolationschemestothe facecenters.Applyingdiscretization andinterpolationschemes forthe termsandvariablesin

(7)

,dividingby VP andreplacingthediscretizedpressuregradientwithitsnon-discretizedformfor easeofnotation,thesemi-discretizedmomentumequationforthenode P canbewrittenas:

APui,P

=

Hi,P

− φ

P

(∂

ip

)

P (10)

AP isa scalarcontaining allthe centralcoefficients, includingimplicittime derivative contributions,and Hi,P isa vector containing contributionsfromtheneighboringcells,aswell astimederivative contributions. AP andHi,P areleft unspec-ifiedastheirexact discretizedformdependson thechoiceofdiscretization.Moredetails onthediscretizationof AP and Hi,P can befound in

[20,22]

.Thechoicesofthediscretization schemesofthisparticularwork arespecifiedinSection 6. Dividingby AP resultsinthefollowingequation:

ui,P

=

A−_P1Hi,P

−

A−_P1

φ

P

(∂

ip

)

P (11)

giving usthe expressionfor thevelocity at cell P .Eq. (11)isused in boththe predictor and corrector stepsof the PISO algorithm

[20]

.Inthepredictorstepitisusedtofindtheinitialsolutionofthevelocitypriortoenteringthecorrectorloop, wherethesameequationisusedtoupdatethevelocityfieldwiththeimprovedpressurefield.

3.1. DiscretizedpressureequationandRhie–Chowinterpolation

Inthechosenformulationthereisnoexplicitequationforthedensityintheconservationequations.Theprimaryreason forchoosingthisformulationiscomputationalefficiency,asdensity-basedsolverstendtorequiresignificantlysmallertime stepsandresolvepressurefielddetails,suchassoundwavepropagation

[14]

,thatarenotofimportancetoourapplication. Inthe selectedformulation,thepressureequation isobtainedby combiningmomentum andmassconservation.Oncethe pressure field issolved for, weobtain the densityfield throughthe useofthe equation-of-state forwhichwe adoptthe idealgaslaw.

Toderive theequation forthepressureEq.

(11)

ismultipliedby

ρ

andthereafterinsertedinto themassconservation equation

(6)

.Furthermore,replacing

ρ

inthetimederivativewiththeequationofstate:

where

ψ

isthecompressibilityfactor,resultsinthefollowingpressureequation: VP

∂

t

(φ

P

ψ

PpP

)

+



f

(

ρ

A−1Hi

)

fSi,f

=



f

(

ρ

A−1

φ∂

ip

)

fSi,f (13)

InEq.

(13)

,

(

A−1

φ∂

ip

)

f needs tobe specifiedusingan appropriate interpolationrulebased oncell-centeredvaluesof thevariables.Onewaytodetermine

(

ρ

A−1

_φ∂

ip

)

f istointerpolateA−P1

φ

P

(∂

ip

)

P anditsneighborvalue A_nb−1

φ

nb

(∂

ip

)

nbonto thecellface f usinglinearinterpolationdefinedas:



f

= (

1

−

r

)

P

+

r



nb (14)

wherer

= |

di,Pf

|/|

di,Pnb

|

,and



istheproperty tobeinterpolated.However,ifthiswouldbeadoptedonauniform, collo-catedgridthepressure–velocitydecouplingwouldbeinevitableasthecontributionofthecentralpressurevalue(node P )

wouldcancelout inthepressuregradientcalculationforcell P [15].It isatthispointthatRhieandChow

[21]

proposed that insteadofdirectlyinterpolatethistermonto thefacesit isapproximatedas

(

ρ

A−1

_φ∂

ip

)

f

= (

ρ

A−1

φ)

f

(∂

ip

)

f,where

(∂

ip

)

f iscalculateddirectlyas:

(∂

ip

)

f

= ˆ

ni,f

pnb

−

pP

|

di,Pnb

|

(15)

withn

ˆ

i,f beingtheunitnormalvector. AdoptingthemethodofRhieandChow,thepressureequation takesthefollowing form: VP

∂

t

(φ

P

ψ

PpP

)

+



f

(

ρ

A−1Hi

)

fSi,f

=



f

(

ρ

A−1

φ)

f

(∂

ip

)

fSi,f (16)

As it is well known,this approach ensures pressure–velocity coupling andavoids the unphysical‘checker-board’ flow pattern.AnticipatingfurthermodificationsinthenexttwosectionsweremarkthatthisRhie–Chowapproximationisonly practicalincasethegoverningequationsdonotcontain discontinuousbodyforces.Modificationsthat copealsowith dis-continuousbodyforceswillbeformulatedintheRDRandFCPmethods.

In thediscretization of thegoverning equations,the mass flux over thecell faces plays an importantrole andneeds toberegularly updatedwithin thealgorithm.The expressionfortheupdate ofmassflux followsfromthepressure equa-tion(16):

(

ρ

ui

)

fSi,f

= (

ρ

A−1Hi

)

fSi,f

− (

ρ

A−1

φ)

f

(∂

ip

)

fSi,f (17) where, again, the pressure gradient at the face is directly computed using the Rhie–Chow assumption from Eq.

(15).

3.2. ModifiedRhie–Chowinterpolationfordiscontinuousbodyforces

Thissection presentsamodifiedRhie–Chowinterpolationmethodmakingitsuitablefordomainsthatincorporatebody forcediscontinuities.Weshowtheproposedmodifications,startingfromthemomentumequationandderivenovelpressure andcorrectorstepequations.Foreasiernotation,we revertthediscretizedpressure term



_fpfSi,f fromEq.

(7)

backto itsnon-discretizedform

(∂

ip

)

P.

To avoid interpolation of the discontinuous

φ

from the cell-centers to the face centers in the Laplacian term (



_f

(

ρ

A−1

φ)

f

(∂

ip

)

fSi,f) of thepressure equation

(13)

, the discretized momentumequation (16)is divided by

φ

P. Fur-therdivisionby VP yields:

1

φ

P

∂

t

(

ρ

Pui,P

)

+

1

φ

PVP



f

(φ

−1

ρ

uj

)

fui,fSj,f

= −(∂

ip

)

P

+

1

φ

PVP



f

τ

i j,fSj,f

+ ζ

i,P

−

DPui,P (18) Insteadof discretizingandinserting DPui,P asacoefficient into AP and Hi,P,aswas done in(10), DPui,P iskept as a separate term, independent of AP and Hi,P. The rest ofthe terms, except the pressure gradient term, are discretized accordingto: APui,P

−

Hi,P

=

1

φ

P

∂

t

(

ρ

Pui,P

)

+

1

φ

2_PVP



f

(

ρ

uj

)

fui,fSj,f

−

1

φ

PVP



f

τ

i j,fSj,f

− ζ

i,P (19) Notethatthereciprocalporosityintheconvectivetermof

(19)

hasbeenmovedoutsideofthesumtoavoidinterpolation ofthediscontinuousporositytothecell-faces.Thisapproximationisrelevantonlyattheinterface,andisgenerallysmallin charactersinceporosityvaluesareboundedbetween0and1.Theresultingsemi-discretizedmomentumequationtakesthe form:

Fig. 2. Examplecellinthevicinityofthefluid–porousinterface.TheownercellisdenotedbyP ,theeastneighborcellbyE andthewestneighborbyW . Thelocalspatialstepx isassumeduniformforeaseofanalysis.Cellfacesfromwesttoeastaredenotedbyw,e andee,respectively.Thehatchedareas symbolizetheporousregion.

Dividingitfirstby AP andthencollectingthevelocityterms,thecell-centeredvelocitybecomes:

ui,P

=

BPA−_P1

(

Hi,P

− (∂

ip

)

P

)

(21)

where BP

= (

1

+

DPA−P1

)

−1.WeemphasizethatonlythematrixBP containstheporousresistivityDP.Thisisimportantto maintainastrongpressure–velocitycouplingthroughoutthealgorithm.Inordertoavoidoscillationsnearthediscontinuity andtointerpolatethediscontinuousvariablesconsistentlytotheinterfaceinthefacemassfluxexpression,B isinterpolated separatelytothefacesresultinginthefollowingmodifiedRhie–Chowinterpolation:

(

ρ

ui

)

fSi,f

=

Bf

[(

ρ

A−1Hi

)

fSi,f

− (

ρ

A−1

)

f

(∂

ip

)

fSi,f

]

(22) where Bf

= (

1

+

Df

(

A−1

)

f

)

−1.Theterms

(

ρ

A−1Hi

)

f,

(

ρ

A−1

)

f,Df and

(

A−1

)

f arefoundusingthelinearinterpolationin

(14),and

(∂

ip

)

f isdiscretizedusing

(15)

.NotethattheexpressioninthesquarebracketsinEq.

(22)

isinfacttheoriginal Rhie–Chow

[21]

expressionformassflux,Eq.

(17)

.Wekeeptheoriginal Rhie–Chowmethodforthecontinuousorslowly varyingvariables(densityandvelocity),andisolatethecontributionfromthelargediscontinuityinthematrixBf.

Inordertoobtainthemodifiedpressureequation,weinsertEq.

(22)

intothediscretizedmassconservationequation

(6)

andreplace

ρ

P inthetimederivativewith

(12)

tofind:

VP

∂

t

(φ

P

ψ

PpP

)

+



f Bf

(

ρ

A−1Hi

)

fSi,f

=



f Bf

(

ρ

A−1

)

f

(∂

ip

)

fSi,f (23)

The modified velocitydefinitionEq.

(21)

andthemodified Rhie–Chowinterpolation forthemassfluxare keysteps in accountingforthediscontinuousresistivitycontribution.Thismodified Rhie–Chowinterpolationconstitutesanew simula-tion methodthat in principle incorporatesdiscontinuities in D. As willbe shown innumericalillustrations inSection 6, this particularmodification by itself doesnot eliminate velocity oscillations to a sufficientdegree. This drives ustoward the constructionoffurtheradaptationsthat moreeffectivelyhandlesuchspurious oscillations,which,nevertheless, follow thepatternaslaiddowninEq.

(21)

.Thefirstoptionforfurthermodificationsaddressesthedefinitionofthematrix BP – thisoptionwillbereferred toastheRDRapproach.Alternatively,asecondoptionpresentsitselfthroughanadaptationof thepressuregradientinEq.

(21)

–thisoptionwillbereferredtoastheFCPapproach,closelyfollowingOxtobyetal.[11]. In thenext twosections wewill elaboratetheRDR andFCPoptions respectivelybeforecoming toa detailedassessment of the original Rhie–Chow approach Eq. (11), the modified Rhie–Chow approach Eq. (21) andthe further RDR and FCP modificationsinSection6.

3.3. Pressure–velocitydecouplingfordiscontinuousbodyforces

TheRhie–ChowinterpolationandthePISOalgorithmwereproposedforsingle-phasefluidflows,forwhichthematerial properties are smooth andcontinuously varying in response to temperatureand densitychanges. The Rhie–Chow inter-polationasdescribed previously hasbeenshownto resultinspuriouspressure andvelocity oscillationsinthevicinity of discontinuities inporouspropertiesoringenerallargesourcetermsinthegoverningequations[7,8,10,12].Therearetwo reasonswhyRhie–Chowinterpolationfailsandwewilldedicatesomespacetoclarifythese.Forthepurposeofclarification, weconcentrateon1Dexamplegridinthevicinityofafluid–porousinterface,shownin

Fig. 2

.

The first reason for the failure of Rhie–Chow interpolation in the presence of discontinuities is that the Rhie–Chow approximation

(

ρ

A−1

_φ∂

ip

)

f

= (

ρ

A−1

φ)

f

(∂

ip

)

f,isinaccuratenearfaceswherelarge bodyforce discontinuities arise.This canbeeasilyverifiedbywritingouttheaboveapproximationinitsdiscreteformforthegridshownin

Fig. 2

:

1 2



ρ

PA−_P1

φ

P pe

−

pw



x

+

ρ

EA −1 E

φ

E pee

−

pe



x



=

1 2



ρ

PA−_P1

φ

P

+

ρ

EA−_E1

φ

E



_pE

₋

_pP



x

(24)

withthe localspatial step



x

=

const, andfacevaluesofpressure pw, pe and pee beinginterpolatedusing Eq.(14). Ex-pandingthelefthandsideofEq.

(24)

usingtheinterpolationEq.

(14)

andsomealgebrawegetthat

1 2



ρ

PA−P1

φ

P pe

−

pw



x

+

ρ

EA −1 E

φ

E pee

−

pe



x



=

1 2



ρ

PA−_P1

φ

P

+

ρ

EA−_E1

φ

E



_pE

₋

_pP



x

+

1 2



ρ

PA−_P1

φ

P 1 2

_pP

₋

_pW



x

−

pE

−

pP



x

−

ρ

EA−E1

φ

E 1 2

_pE

₋

_pP



x

−

pEE

−

pE



x



(25)

Fig. 3. Examplegridshowingvisuallythevalueoftheporousresistanceacrosstheinterface.TheownercellisdenotedbyP andtheneighborcellbyE. Cellfacesaredenotedw,e andee.x isthelocalspatialstepofthegrid,assumedconstantforsimplicity.Thetexturizedgreyareaisrepresentingthe valueoftheporousresistancetermφDuifromEq.(2).Thesolidlinesmarkthestateofporousresistancevaluepriortousinganyadditionalinterface

treatmentmethods,whilethedashedlinescombinedwithright-slantedhatchingpatternareporousresistancevaluesaftertheapplicationoftheRDR method.

wherethelast twotermsontherighthandsideofEq.

(25)

canbe thoughtofas‘correction’termsrelativetotheoriginal approximation Eq.(24).In fact,the firstterm onthe righthandside corresponds directly toEq.

(24)

whilethe last two termslistadditionalcontributions.Theseadditionalcontributionswerewrittenintermsofthedifferencebetweentheleft andrightsided derivative ofthepressure atthe cell faces,which can be showntoreduce to zeroin caseofa differen-tiable pressurefield. However, in thecaseof discontinuous body forces,i.e., a kinkin thepressure field indicative ofits non-differentiability atthe interface, thesecorrection termsare not necessarilysmalland alternativeapproaches suchas RDRandFCPareessential.

ThesecondreasonforthefailureoftheRhie–ChowversionofPISOiswhatMencingeretal.

[12]

refertoasimbalance of

thecellpressuregradientandcellbodyforceterm. Wewilldescribethisimbalance in somemoredetail,therebyderiving thebasisfortheRDRmethod.

Near an interface ofdiscontinuitywe equate thepressure gradientterm

φ∂

ip and the porousresistance term

φ

D



ui



fromthediscretizedmomentumequation,Eq.(7),forcell P in

Fig. 3

φ

P



f

pfSi,f

= −φ

PDPui,PVP (26)

Inthechosen1DsettingVP

= 

x.Dividingby

φ

P and



x,andcarryingoutthesumwearriveat

pe

−

pw



x

= −

DPui,P (27)

Weuselinearinterpolationtoobtainvaluesofthepressureatthefaces w ande.

pE

+

pP

−

pP

−

pW

2



x

=

pE

−

pW

2



x

= −

DPui,P (28)

ExpandingEq.

(28)

relativetothefacevalueatfacee wefind

(

pE

−

pe

)

+ (

pe

−

pW

)

2



x

= −

DPui,P (29)

whichcansuggestivelyberewrittenintermsoftheweightedsumoftwopressurederivatives

1 2



₁ 2 2

(

pE

−

pe

)



x

+

3 2 2

(

pe

−

pW

)

3



x



= −

DPui,P (30) orsimplified 1 4



₍

_pE

₋

_pe

₎

1 2



x



+

3 4



₍

_pe

₋

_pW

₎

3 2



x



= −

DPui,P (31)

Theformulationisintermsoftheintrinsicvelocityforwhichweapproximateui,P

=

ui,W

=

ui,E

=

ui.Foreach ofthe pressuregradientsin

(31)

wecanidentifyacorrespondingresistivitycontribution.Thiscanbeexpressedas

−

1 4DEeui

−

3

4DWeui

= −

DPui (32)

whereDEeandDWeareporousresistivityvaluesassociatedwiththeregionsbetweenpointE andfacee andbetweenpoint W andfacee, respectively.Inthe settingaspresented inFig. 3we observethat DWe

=

0, sincethisregionis purefluid, while DEe

=

DE sincethisregion islocated within theuniformporousmedium. Bythe samerationale DP

=

0.Inserting thesevaluesintoEq.

(32)

wenotice

−

1 4DEui

−

3

ThisiswhatMencingeretal.

[12]

refertoasimbalance betweenthepressuregradientandthebodyforce.Thereisaneasy remedytothisproblemwhichconsistsofre-defining DP forwhichweproceedasfollows.Althoughphysicallytheporous resistance incell P iszero,Eq.

(33)

showsthat aforce balance wouldbe maintainedifwe wouldassign DP a newvalue Drd_P

=

1₄DE.Withasimilarargumentappliedtothecellcenteredaround E itcanbeconcludedthatthevalueofDE should bemodifiedto Drd

E

=

3 4DE.

Thisdiscretizationresultsinaredistribution ofporousresistancebetweencellE andcell P ,whichiswhythemethodis calledRe-DistributedResistivity (RDR).Agraphicaldepictionofporousresistancevaluesassociatedwiththeutilizationof theRDRmethodisdenotedbydashedlinesandright-slantedhatchingin

Fig. 3

.

Apartfromtheporousresistanceredistribution,wherelinearinterpolationofthepressurefromcellcentersto facesis appliedanddiscontinuitiesarehandledbyredistributingtheresistivity,analternative‘accounting’ofdiscontinuitiesinthe bodyforcetermcanbeachievedbyassigningalternativevaluestothepressureattheinterface pe.Thishasdirectinfluence onthepressuregradientsincells P andE.Derivingsuchanauxiliaryvalueforpe istheessenceoftheFCPmethod,which isobtainedcloselyfollowingOxtobyetal.

[11]

.Technical detailsoftheimplementationofbothRDRandFCPmethodswill becoveredindetailinSections4and

5

.

4. Re-distributedresistivityalgorithm

Due tothe presenceofthe bodyforce discontinuity,represented bythe porousresistanceterm

φ

D



ui



inEq.

(2)

,the pressureatthefluid–porousinterfacebasedonthestandardcentraldiscretizationissuchthatthevelocitysolutionadjusts withvelocityoscillations.Theseoscillationshavebeenexploredbeforein

[12,10]

.Inordertorestorethebalancebetween thepressure gradientandthebody force,weproposeto adapttheresistivity inthecells adjacenttotheinterface.Thisis thebasicprinciplebehindthenovelRe-DistributedResistivity(RDR)algorithm

[9]

whichwillbepresentedbelow.

4.1. Flowresistivityredistribution

Toavoidthegenerationofspuriousvelocityoscillationsinthevicinityofthediscontinuity,whensolvingthediscretized momentum equation, Eq.(20) or correcting the velocity with Eq.(21), a balance betweenthe discontinuous porous re-sistance field and the cell-centered pressure gradient is required [12]. Because of the body force discontinuity at the fluid–porous interface,theinterface pressurevalue pf iserroneouslydeterminedbylinearinterpolation usingEq.

(14)

as was showninSection 3.3.Thesepressuregradientsinthetwocellsneighboringthediscontinuityarenolongerinbalance withthecell-centered flowresistanceinthesecells.Because ofthis, thecell-centered resistancesaremodified tobalance thecell-centeredpressuregradientscalculatedfromthelinearlyinterpolatedfacepressures,asshowninSection3.3.

Inordertodeterminethecellsrequiringredistributionofporousresistivity,acellindicatorfunctionisdefinedas:



P

=



f

|(φ)

f

|

max

(



_f

|(φ)

f

|, ξ)

(34)

where

(φ)

f

= φ

nb

− φ

P and

ξ

isasmallnumberontheorderof10−15 toavoiddivisionbyzeroif



f

|(φ)

f

|

=

0.



P hasthevalue1 forcellsrequiringredistributionofresistivityand0 otherwise.Anotherindicatorfunction

θ

f isdefinedas:

θ

f

=

|((φ))

f

|

max

(

|((φ))

f

|, ξ)

(35)

where

(φ)

f

= 

nb

φ

nb

− 

P

φ

P, in orderto determine the faces required forthe calculation of the redistributedflow resistivity inthe cells where



P

=

1.

θ

f hasthe value 1 for the fluid–porous interface faces and the faces of the cells adjacenttoadiscontinuityand0 otherwise.TheredistributedresistivityDrd_P iscomputedusinginversedistanceweighting fortherequiredfacesaccordingto:

Drd_P

= (

1

− 

P

)

DP

+



P



f

ω

f

θ

f



f

ω

f

θ

fDf (36)

where

ω

f

=

1

/

|

di,Pf

|

andDf isinterpolatedtothefacecentersbythelinearinterpolationin

(14)

.Thegraphicaldepiction ofthisporousresistance redistributionprocess isshowninFig. 4.Theredistributedresistivity

(36)

is basedona‘double’ interpolation incells withdiscontinuities, i.e.,firstlinearinterpolationisusedtoestimate thefacevaluesDf,afterwhich inverse distanceweightingisused tocompletethe re-distribution.Thisapproachissimilar tothat proposed byRaeini et al.(2012)

[23]

forsmoothingtheinterface representationinavolume-of-fluidmethod,byusing‘doubleinterpolation’and Gaussiansmoothingtodefinethenewinterfacenormals.

Themodifieddiscretizedmomentumequationwiththeredistributedresistivityyields:

1

φ

P

∂

t

(

ρ

Pui,P

)

+

1

φ

2_PVP



f

(

ρ

uj

)

fui,fSj,f

= −

1 VP



f pfSi,f

+

1

φ

PVP



f

τ

i j,fSj,f

+ ζ

i,P

−

DrdPui,P (37)

Fig. 4. Graphicaldepictionshowingthestateofvariablesaftertheredistributionofthe porousresistancehasbeen carriedout.Thenewstateofthe porousresistancefieldcompensatesforoverorunderestimatesofthepressuregradientsinthecellsadjacenttotheinterfaceandforcebalance isrestored, retainingthecorrectvelocityfield.

andtheexpressionforthevelocitycorrectionreads:

ui,P

=

BrdPA−P1

(

Hi,P

− (∂

ip

)

P

)

(38)

where Brd_P

= (

1

+

Drd_PA−_P1

)

−1. The face velocity ui,f in the convective term can be discretized by any suitable scheme (e.g.linear) andthe massflux

(

ρ

ui

)

fSi,f iscalculated fromEq.(22).We emphasizethat relative tothe first Rhie–Chow modificationEq.

(21)

,theRDRapproachcanbeinterpretedbyachangeinBP onlythrougharedistributionoftheresistivity D inthedirectvicinityofthefluid–porousinterfaceorporousdiscontinuity.

4.2. Re-distributedresistivityalgorithm

TheoriginalPISOalgorithm

[24]

isastepwisetechniqueforthesolutionoftheimplicitlydiscretizedtime-dependentflow equationsinasegregatedmanner.Thesolutionisaccomplishedateachtimestepthroughasequentialpredictor–corrector processbywhichthedifferentdependentvariablesareupdated individually.Followingthenotation in

[24]

,variableswith asterisksareintermediatestepsolutions.TheRDRalgorithmschemeisshownin

Fig. 5

.

The RDR algorithm differs from the original PISO algorithm [24] in that (i) the equations are differently formulated,

(ii) the modified Rhie–Chowinterpolation isused and(iii) an extrastep prior tothe predictorstepis introduced, where weexecutetheporousresistanceredistributioninaccordancewithEq.

(36)

.Nosignificantcomputationalcostwasnoticed becauseofthisextrastep.Subsequently,we continue theuseofthe Brd_P matrixinthefirstandsecond correctorstepsin ordertokeepthemethodconsistent.TheoriginalPISOalgorithm

[24]

wasintendedtohaveonlyasinglepassthroughthe predictorand two correctorsteps. In orderto improvethe precision of thealgorithm andto control conservationerrors coming fromthesegregated approach,we have addeda convergencecondition atthe endofthe algorithm asproposed by[25].Thisconditionisbasedontheresidualvalueobtainedbyinsertingthelatestpressuresolutionp∗∗intothepressure equation ofthe first correctorstep. Oncethe residualvalue goes downbelow apre-set value, e.g. 10−12,the solution is considered asconverged.Theconservationerrors cantherefore becontrolled to an arbitrarilyprecision by specifyingthe pre-setvalue.

InthenextsectionwewillpresenttheFCPmethod,whichachievessimilarresultsasRDRbymodifyingthevalueofthe interfacepressureinsteadoftheresistivity.

5. FaceConsistentPressurePISOalgorithm

IntheprevioussectionwepresentedtheRDRmethodwhichaimstoestablishabalancebetweenthepressuregradient and the body force term by redistributing the body force termbetween the cells adjacent to the interface. Considering Eq.

(21)

an alternative approach presents itself: if thepressure at the interface would be adapted, the balance between thepressure gradientandthe bodyforce terms couldbe achievedaswell. Thisapproachwas previously explored by [7]

and

[11]

,wherespecialfacepressureinterpolationschemeswereintroduced toaccuratelyestimatetheinterfacepressure. Sincesuchinterfacepressurevaluesarederivedtoaccountfortheinfluenceofthelocaljumpinthebodyforceacrossthe interface,thepressureattheinterfaceisreferredtoasconsistent withtheresistivityjump.Schemesin

[7]

and

[11]

require eitherexplicit correctionsbased onthe mass flux orextrapolations ofthe cell-centered pressurefrom eitherside ofthe interface.

Inthissectionwe presenttheFaceConsistentPressure(FCP)algorithm,whichutilizesauxiliarypressure valuesatthe interface in order to prevent spurious velocity oscillations.The method builds on thework ofOxtoby etal. [11], butis differentinthreekeyways:FCPsolvesforextrinsicvelocityfield,FCPtreatsall terms,exceptpressuregradient,implicitly, andFCPincorporatesthemodifiedRhie–Chowinterpolation,Eq.

(21)

,presentedinSection 3.2.TheFCPmethodcalculates theinterface pressurevalue whichis consistentwiththeadjacentflow resistance.Thisresultsin apressurefield that by defaultsatisfiesthebalancebetweenthediscontinuousflowresistancefieldandthecell-centeredpressuregradient

[12]

in cellsadjacenttotheinterfaceatsteady-stateconditions.Thederivationoftheequationsandthefunctioningofthemodified algorithmarepresentedbelow.

Fig. 5. RDRalgorithm.ThemajordifferencewithrespecttotheoriginalPISO[24]isthatthetimestepbeginswithredistributingtheporousresistance values.Theconvergenceconditionisbasedontheresidualvalueobtainedbyinsertingthelatestpressuresolutionp∗∗intothepressureequationofthe firstcorrectorstep[25].

Fig. 6. Controlvolumearoundtheinterface(dashedline)whichisusedforderivationofthepressureinterpolationattheinterface.Themajordifference withrespecttothegenericcellrepresentationofFig. 1isthatwenowconcentrateontheimmediatevicinityofthefluid–porousinterface.Todifferentiate theinterfacefacesfromtheotherfacesdenotedwith f ,wedenotetheinterfacevalueswithindexε.Thecontrolvolumeusedforderivationofauxiliary interfacepressurevalueismarkedwithdashedlines.Thevectorconnectingthetwocellcentersisdenoteddi,PEandthedistancefractionsr and(1−r)

arealsomarked.

5.1. Derivationoftheconsistentinterfacepressure

Letusanalyzethemomentumbalanceofafictionalcontrolvolumearoundtheinterfaceasshownin

Fig. 6

.Cellcenter P isinthepurefluiddomainandcellcenter E isintheporousdomain.Inordertodifferentiatebetweenthespeciallyderived interface pressureandthe other face values(calculated by linearinterpolation), we denote theinterface values withthe index

ε

.Theinterface

ε

ispositionedbetweenpoints P and E atdistancefractionr awayfrompoint P and

(

1

−

r

)

away frompointE.Thedistancebetweenthetwopointsisdenotedby

|

di,PE

|

.

Followingthecombinednotation of

Fig. 6

andSection 3,we firstwrite outthemassconservationexpressions foreach halfofthecontrolvolume:

(∂

t

(φ

ρ

))

P

+

ˆ

ni

(

ρ

ui

)

ε

− ˆ

ni

(

ρ

ui

)

P r

|

di,PE

|

=

0 (39) and

(∂

t

(φ

ρ

))

E

+

ˆ

ni

(

ρ

ui

)

E

− ˆ

ni

(

ρ

ui

)

ε

(

1

−

r

)

|

di,PE

|

=

0 (40)

BycombiningEqs.

(39)

and

(40)

weget:

r

|

di,PE

|(∂

t

(φ

ρ

))

P

+ (

1

−

r

)

|

di,PE

|(∂

t

(φ

ρ

))

E

= ˆ

ni

(

ρ

ui

)

P

− ˆ

ni

(

ρ

ui

)

E (41)

ThemomentumfluxesfortheleftandrighthalfofthefictionalcontrolvolumecanbeexpressedonthebasisofEq.

(21)

(

ρ

ui

)

P

= (

B

ρ

A−1Hi

)

P

− (

B

ρ

A−1

)

P pε

−

pP r

|

di,PE

|

ˆ

ni (42) and

(

ρ

ui

)

E

= (

B

ρ

A−1Hi

)

E

− (

B

ρ

A−1

)

E pE

−

pε

(

1

−

r

)

|

di,PE

|

ˆ

ni (43)

wherethepressuregradientisapproximatedusingtheauxiliarypressurevalue pε.InsertingEqs.

(42)

and

(43)

intoEq.(41), wecansolvefortheconsistentpressureattheinterface:

pε

= (

B

ρ

A−1

)

−ε1

[(

B

ρ

A−1p

)

ε

−

r

(

1

−

r

)

|

di,PE

|

2

((∂

i

(

B

ρ

A−1Hi

))

ε

+ (∂

t

(φ

ρ

))

ε

)

]

(44)

withinterpolatedtermsbeing:

(

B

ρ

A−1

)

ε

= (

1

−

r

)(

B

ρ

A−1

)

P

+

r

(

B

ρ

A−1

)

E (45)

(∂

i

(

B

ρ

A−1Hi

))

ε

= |

di,PE

|

−1

(

n

ˆ

i

(

B

ρ

A−1Hi

)

P

− ˆ

ni

(

B

ρ

A−1Hi

)

E

)

(46)

(∂

t

(φ

ρ

))

ε

=

r

(∂

t

(φ

ρ

))

P

+ (

1

−

r

)(∂

t

(φ

ρ

))

E (47)

and

(

B

ρ

A−1p

)

ε

= (

1

−

r

)(

B

ρ

A−1p

)

P

+

r

(

B

ρ

A−1p

)

E (48) We refer to pε as consistent pressurebecause its value istaking intoconsideration the local jumpin porous resistance throughthecoefficientmatrixB.

In orderto obtain the expression for thefinal pressure field, we utilizethe same maskingfunction

θ

f asdefined in Section4.1toidentifytheporousfluidinterface.Theexpressionofthenewfacepressurefielddenoted



f thenbecomes:



f

= (

1

− θ

f

)

pf

+ θ

fpε (49)

Notethatbecauseoftheuseofthemaskingfunction,thenewly devisedfaceconsistentpressurecorrection(secondterm inEq.

(49)

)will beappliedonly attheinterface.The other facevaluesofpressure(pf) inthedomain willbecalculated accordingto the interpolationscheme ofchoice (e.g.linear, upwindetc.). Thiscan be alteredso that the faceconsistent pressureinterpolationwouldbeappliedthroughoutthedomain



f

=

pε,regardlessoftheinterfacepresence.Wereferto thisapproachastheGlobal-FCPmethod,sincetheinterpolationofEq.

(44)

isappliedtoallthefacesinthedomain.Wedo notpursuethisapproachinthispaperandleaveitsanalysisforfuturework.

Byusingthenewlydevisedpressure



f wecannowformapressuregradientterminthemomentumequation(Eq.(18)) toget: 1

φ

P

∂

t

(

ρ

Pui,P

)

+

1

φ

PVP



f

(φ

−1

ρ

uj

)

fui,fSj,f

= −



f



fSi,f

+

1

φ

PVP



f

τ

i j,fSj,f

+ ζ

i,P

−

DPui,P (50) Bycreatingthisconsistentvalue offacepressure



f attheinterface, weaimatavoidingtheforceimbalanceinthecells adjacenttotheinterface.

Acorrectionofthepressurevalueattheinterfacehasbeenderivedbyanalyzingthemomentumbalanceinthefictitious control volume around the interface. This corrected pressure incorporates the influence of the porous resistance values adjacenttothefluid–porousinterface.Therefore,thepressuregradientsatpoints P andE satisfytheforcebalancebetween thediscontinuousflowresistancefieldandthecell-centeredpressuregradient

[12]

.Thisforcebalanceaimstoremovethe spuriousoscillationsinthevelocityfield.

Fig. 7. FCPalgorithm.ThedifferencewithrespecttotheRDRalgorithmisthatinsteadofredistributingporousresistivity,thefacepressurefieldpf is

updatedtof usingEq.(49)throughoutthealgorithm.

5.2. FaceConsistentPressurealgorithm

The FCPalgorithm isverysimilar totheoriginal PISOalgorithm proposedby Issaetal.[24].The maindifferencesare

(i) the differentlyformulated equations,(ii) themodified Rhie–Chowinterpolationand(iii) theuseofamodified pressure (Eq.(49))inthemomentumpredictorandcorrectorstepsthroughoutthealgorithm.Following

[24]

,theFCPalgorithmtakes theformshownin

Fig. 7

.

TheFCPalgorithmusesthesamestepsastheRDRalgorithm,exceptthatinsteadofupdating BP to BrdP inthepredictor andcorrectorsteps,we updatethefacepressurefield pf to



f.The convergencecondition attheendisthesameasfor theRDRalgorithmtocontrolconservationerrorscomingfromthesegregatedapproach.

TheproposedRDR,FCPandtheoriginalPISOalgorithmsdescribedinSections4,

5

and

3

,respectively,havebeen imple-mentedusingtheOpenFOAM



opensourcecomputationalfluiddynamicsC++library(version2.2.0).

6. PorousplugandBeavers–Josephtestcases

Thecaseforflowperpendiculartoaporousregion

[26]

andtheBeavers–Josephproblem

[16]

forflowparalleltoaporous region are considered, in order to demonstrate the accuracy and robustness of the proposed algorithms. The respective geometriesofthesetestcasesareshownin

Fig. 8

.Intheseillustrationsweassumeincompressibleandisothermalflowand

Fig. 8. Geometries for the (a) porous plug and (b) Beavers–Joseph cases.

an isotropic porous medium even though the algorithms are implemented ina compressible formulation. Explorationof compressibilityandheattransferareconsideredinfuturework.

Thefirstcaseissetina2Dchannel,ofheighth andlength8h,inwhichanisotropicporousplughasbeeninstalled be-tween3h and5h[26].Theinletvelocityboundaryconditionisafullydevelopedparabolicprofile,withmeanvelocityU .The ReynoldsnumbertakesthevaluesRe

= {

1

,

100

,

1000

}

andisdefinedasRe

=

ρ_{μ .}U h TheDarcynumberoftheflowisdefined asDa

=

K

/

h2_{,where K isthe porousmedium permeability,andfor testingpurposeswe chose Da}

_{= {}

₁₀−2

_,

₁₀−3

_,

₁₀−7

_}

_. ThesevaluesofRe andDa numberscorrespondtolowandmoderateamplitudesoftheporousresistanceterms.The poros-ityoftheporousplugiskeptconstantat

φ

=

0

.

7.Theoutletboundaryconditionforvelocityiszerogradientorn

ˆ

j

∂

jui

=

0i wheren

ˆ

jisthenormalvectortotheoutletand0i isazerovector.Thepressurevalueisspecifiedasconstantattheoutlet, whileattheinletwerequirethezerogradientvalueofpressureorn

ˆ

i

∂

ip

=

0.Forchannelwalls,no-slipboundaryconditions (ui

=

0) wereapplied.Thegridisuniform, structuredanditsresolutionisequivalenttothegridresolutionusedin

[7]

for bettercomparison.Thegeometryispresentedin

Fig. 8

(a).

ThesecondcaseweconsiderwasoriginallyproposedbyBeaversandJoseph

[16]

,andiscommonlyusedintheliterature asavalidationcase.Weadoptaslightlymodified versionfollowing

[26]

and

[7]

.Thedomainisa2Dchannelofheight 2h andlength8h.Halfofthechannelheighth isfilledwithisotropicporousmaterial.Fortheinletvelocityboundarycondition auniformvelocityofmagnitudeU isadoptedatRe

=

1 sothattheflowhastimetofullydevelopbeforereachingtheend ofthe domain.The testedDarcynumbers are Da

= {

10−2

_,

₁₀−3

_}

_{, andtheporosity is keptconstant forall cases}

_φ

₌

₀

_.

_7. Theoutletboundarycondition forvelocityiszerogradientorn

ˆ

j

∂

jui

=

0i wheren

ˆ

j isthenormalvectortotheoutlet.The pressure valueis specifiedasconstant attheoutlet, while attheinlet we requirethezerogradient value ofpressureor

ˆ

ni

∂

ip

=

0.Thechoiceofthegridresolutionwasmotivatedbytheworkof

[7]

forbettercomparison.Wallsofthechannels haveanimposedno-slipboundaryconditionforthevelocityorui

=

0.Thegeometryisshownin

Fig. 8

(b).

Thetimederivative isdiscretizedbyasecond-orderimplicitbackwarddifferencingschemeandtheconvectivetermby the second-order linear upwind differencing (LUD) scheme [27].The discretized equations foreach coordinate direction are solved usingstandard solvers andutilitiesavailable in OpenFOAM



such as:thesmooth solver witha Gauss–Seidel smootherdowntoatolerance10−11_{forvelocitycomponents,andthepressureequationsaresolvedwithaPreconditioned} ConjugateGradient (PCG)solveradoptingthe FasterDiagonal Incomplete-Cholesky (FDIC)preconditionerdown toa toler-anceof10−12.

Theproposed RDR andFCPmethods differfromastandardcentral discretizationwhereverthereare variationsor dis-continuities in the porosity andconsequently inthe resistivity term. Central discretizationmethods have alreadyproven theirsuitabilitytosimulateflowinporousmediawithsmoothlyvaryingporosity

[28]

.Therefore,werestricthereto discon-tinuitiesintheporositynearfluid–porousinterfaces,wherespuriousoscillations mayarise.Thetwo testcasesconsidered aredesignedtoaddressthesechallengesclosetofluid–porous interfaces,wherediscontinuities aremostpronounced.The proposedalgorithms arehoweveralsosuitable forgeneralheterogeneousfluid–porousmedia,e.g.,wheretheporosityand resistivityaresmooth functionsofthepositioninthedomain.Wevalidatethenewmethodforsegregatedsolversagainst numericalresultsfromliteraturebasedonfullyblock-couplednumericalmethods

[8]

.Validationisdonefor2Dgeometries, on structured Cartesian grids, becauseno standard validationcase basedon 3D geometries are available. The robustness andapplicabilityofthealgorithmsforcomplex3Dgeometriesandgeneralunstructuredmesheswillbepursuedfurtherin futureresearch.

6.1. Porousplugcase

In

Fig. 9

,thevelocitymagnitudealongthecenterline oftheporousplug geometrylocatedatheight0

.

5h ispresented. Since thevelocity inlet is specifiedas fullydeveloped parabolicprofile withmeansuperficial velocity U ,the peak value at thecenterline is1

.

5U . Oncethe flow reaches the porous plug, we registera dipin the centerline velocity. This may appearcounter-intuitive,astheporousplugconstrictstheflowandaccelerationisexpected.However,whathappensisthat thevelocityprofileinthe y-directionchangesoncethefluidenterstheporousplug,goingfromparabolictoaflatprofile, thereby maintaining conservationofmass even though thecenterline velocity is much reduced. Afterthe fluidexits the porousplug,thevelocityprofilein y-directionisseentorapidlyrecoverbacktothePoiseuilleprofile.

Fig. 9. Centerline velocity for porous plug problem for Re=1 and Da=10−3_.

Fig. 10. Figure(a)showingcenterlinepressurefortheporousplugproblematRe=1 andDa=10−3_{.Figure(b)iszoomedinaroundtheinterface,where}

thepressuregradientisdiscontinuous.

It canbe seen inFig. 9(a)that boththe RDRandFCP algorithmsprovide anoscillation-free solutionforthereference case(Re

=

1,Da

=

10−3), whilethesolutionobtainedby theoriginalPISO algorithmexhibitssubstantial oscillations.This is emphasizedin

Fig. 9

(b),wherewesee azoomed-inarea aroundtheinterface. Wealsoseethat boththeRDR andFCP algorithmsareingoodagreementwiththereferencedata

[26]

.AthigherRe and/orlowerDa thespuriousoscillationswere foundtobecomeevenmorepronounced.

In

Fig. 10

(a)we presentthepressure fieldalongthe centerlineofthegeometry,correspondingto thevelocitysolution shownin Fig. 9.Thepressure gradientinthe fluidarea upstream oftheporousplug islow astheflow doesnotrequire a high pressuregradient to maintain the requiredmass flowrate. Once the fluid reachesthe porousplug at 3h,we see a discontinuityin the pressure gradient and a much stronger forcing ofthe flow, i.e.a much steeper negative pressure gradient, is neededin theporous region to keep the mass-flowrate constant.The pressure gradient rapidlyrecovers its valuethatwasalsoobservedaheadoftheplugafterexitingtheporousplugat5h.

Allsolutionsseemtomatchthereferencedataof

[26]

.Whenwezoomaroundtheinterface,

Fig. 10

(b),weseethatthe originalPISOalgorithmisproducingsmallpressurevariations,whicharelinkedtothevelocityoscillations.Wealsoseethat RDR andFCPslightlyunderestimatethevalue ofthepressurewithrespecttothereferencecase, butshow nounphysical oscillations. The slightunderestimation by about 0.5% is a measure for the remaining discretization errorat the spatial resolutionthatwasadopted.

Next,wecomparetheRDRandFCPalgorithmsfordifferentflowandporosityconditionsandcomparevelocitysolutions forfour differentcasesatRe

= {

100

,

1000

}

andDa

= {

10−3

_,

₁₀−7

_}

_.

_{Fig. 11}

_{showscenterline velocity solutionsoftheRDR} andFCP methodsforthesefour cases.Thesesolutionslook identical,withthe onlynoticeabledifference beinga slightly smoother velocity profileof the FCPsolution atthe porous plug inlet(h

=

3) at Da

=

10−7. The original PISO algorithm failedatproducingaccurateresultsforthesedemandingconditions.Thediscrepanciesin

Fig. 11

becomeobviousoncewe concentrateontheareaaroundtheinterface,whichareshownin

Figs. 11(c)

and11(d).

Fig. 11. CenterlinevelocitiesfortheporousplugproblematRe= {100,1000}andDa= {10−3_,10−7_{}.TheFigures(a)and(c)showtheRDRsolutions,while}

theFCPsolutionsarepresentedintheFigures(b)and(d).Noticethesmall-amplitudeoscillationsintheFCPsolutionsinFigure(d),particularlyforcases withDa=10−3_.

From

Figs. 11(c)

and11(d)itcanbeimmediatelyconcludedthatRDRproducessmoothersolutions.

Theconvergenceofbothalgorithmshasbeentestedtoconfirmcorrectimplementationofthemethodsandtoquantify thetruncationerrorcomingfromspatialdiscretization.

Fig. 12

(a)showsthel2 _{normplotofboththeRDRandFCPsolutions} fortheporous plugcase. The referenceresolution casehada fourtimeshigherresolution thanthehighestplottedpoint in

Fig. 12

(a).

Fig. 12

(b)qualitativelydemonstratestheconvergence,whilealsoshowinggoodagreementwiththereference dataofCostaetal.

[26]

.Theconvergencerateforbothalgorithmsisroughlyofthefirstorder.

Asa remark,wementionthat inadditiontotheoriginal PISO,RDR andFCP algorithms,we alsoinvestigatedthe pos-sibilityofusingonlya modifiedRhie–ChowPISO algorithm,withoutthe RDRorFCPextensions.Thisversion ofthecode producedsmoothvelocitysolutionfortheRe

=

1 andDa

=

10−3 case,butatmoreaggressiveparameters,e.g.Re

>

100,the spuriousoscillationsoccurred. Thisconfirms thatmakingRhie–Chowinterpolationconsistentisnotenoughtoremovethe spuriousoscillationsingeneral,assuggestedinSection3.2.

Weusedaporousplugcase

[26]

tocomparetheperformanceofthenewlydevelopedRDRandFCPalgorithms.Wefind thatforlowRe/highDa numbers,bothRDRandFCPalgorithmsgivesatisfactorysolutions.OncepushedtohigherRe/lower Da numbers,theRDRalgorithmmaintainsthesmoothnessofthevelocitysolution,whiletheFCPsolutionsexperiencesome oscillationsaroundtheinterface.

6.2. Beavers–Josephcase

FortheBeavers–Josephcase[26,7,16]wepresenttheresultsatRe

=

1 andDa

= {

10−2

_,

₁₀−3

_}

_{forwhicha comparison} withareferencesolutionfromliteratureisshownin

Fig. 13

.

All algorithms are seen toperform equally well (nearly overlapping),compared tothe referencedata [26]. As can be seen,theoriginalPISOalgorithmalsoprovidesoscillationfreeresults.Thisisbecausethepressuregradientovertheporous interfaceisclosetozero,leadingtoanegligibleflowinthe y-direction,i.e.,normaltotheinterface betweenthefluidand theporousmedium.Asadirectconsequence,thespuriousoscillationsarenottobeexpected.Thereforewelimitourselves topresentresultsonlyforafewselectedparameters.