Multiscale reconstruction in physics for compositional simulation

Multiscale reconstruction in physics for compositional simulation

Ganapathy, Chandrashekar; Voskov, Denis

DOI

10.1016/j.jcp.2018.08.046

Publication date

2018

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Ganapathy, C., & Voskov, D. (2018). Multiscale reconstruction in physics for compositional simulation.

Journal of Computational Physics, 375, 747-762. https://doi.org/10.1016/j.jcp.2018.08.046

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

'You share, we take care!' - Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Multiscale

reconstruction

in

physics

for

compositional

simulation

Chandrashekar Ganapathy

a

,

Denis Voskov

a

,

b

,

∗

a_{Department of Geoscience and Engineering, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, the Netherlands} b_{Department of Energy Resources Engineering, Stanford University, 367 Panama Street, 065 Stanford, CA 94305, United States}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Article history:

Received 16 March 2018

Received in revised form 10 July 2018 Accepted 24 August 2018

Available online 29 August 2018 Keywords:

Multiscale transport Upscaling Fractional flow Gas injection

Acompositionalformulationisareliableoptionforunderstandingthecomplexsubsurface processes and the associated physical changes. However, this type of model has a great computational cost, since the number of equations that needs to be solved in each grid block increases proportionally with the number of components employed. To address this issue, we herewith propose a multiscale reconstruction in physics for compositional simulation. The ideology consists of two stages, wherein two different sets of restriction and prolongation operators are defined based on the dynamics of compositional transport. Inthe first stage, an operator restricting the arbitrarynumber of components to onlytwo equations for flowand transport is implementedwith the objectiveofaccuratelyreconstructingthemultiphaseboundariesinspace.Theprediction of multiphase front propagation is the most critical aspect of the approach, as they involve alot ofuncertainties. Once the position oftwo-phase boundaries is identified, the fullconservative solutioninthe single-phaseregioncan beaccuratelyreconstructed based on the prolongation interpolation operator. Subsequently, in the second stage, the solution for the multicomponent problem (full system)in the two-phase regionis reconstructedbysolvingjusttwotransportequations withtheaidofrestrictionoperator definedbasedonaninvariantthermodynamicpath.Theproposedreconstructionstrategy resultsincoarseningofthecompositionalproblemintermsofthephysicalrepresentation (numberofequations),therebyappreciablyreducingthesimulationtimebyseveralfolds withoutsignificantlossintheaccuracy.Wedemonstratetheapplicabilityoftheproposed multiscalestrategyforseveralchallenginggasinjectionproblems.

©2018ElsevierInc.Allrightsreserved.

1. Introduction

Acompositionalmodelisan inevitableandimportanttoolinunderstandingthesubsurfaceprocesses.The solutionfor thismulti-component multiphaseflowproblemisdeterminedbysolvingtheassociatednonlineargoverningequations[1]. Anequationsofstate(EoS)modelisgenerallyemployedtodescribethephasebehaviorofthesystem,andtheyareinturn resolvedintwo stages:phasestabilitytest [2] andflashcalculation[3]. Anonlinearsystemofthermodynamicconstraints is used to represent the instantaneous thermodynamic equilibriumof the mixture. The aforesaid procedure is generally referredtoasthestandardEoSbasedapproachforsolvingcompositionalproblem[4].Ontheotherhand,thesesimulations

*

Corresponding author at: Department of Geoscience and Engineering, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, the Netherlands. E-mail address:[email protected](D. Voskov).

https://doi.org/10.1016/j.jcp.2018.08.046

arecomputationallyintensiveanddemandingsincethenumberofequationsthatneedstobesolvedineachgridblockand ineachiterationisremarkablyhigh.

Onepossibilitytoimprovetheperformanceofreservoirsimulationistoapplyamulti-scaleapproach.Ithelpsto recon-structafine-scalesolutionclosetothecostofacoarse-scalesolutionwiththecontrollederrorinaccuracy,therebymaking it anattractiveproxymodel.Variousmulti-scalemethods havebeendevelopedfortheapplicationofreservoirsimulation [5–7].Amongst them,Multi-ScaleFinite Volumemethods havebecome astate ofthe artinmodern reservoirsimulation [8–10],astheyimprovedtheperformancebyseveralfolds.Butmostofthemulti-scalemethodsarelimitedtotheelliptic solution(globalpressureorconductivetemperature),withoutimprovingthehyperbolicpart(saturationorcomposition).The onlyexceptionistheworkby[11],wheretheyuseddifferentoperatorsforsaturationreconstructionintheheuristic-driven fronttrackingand,[12],wheretheyusedanadaptivemeshrefinementbasedonanalgebraicmulti-scale.

Significant progresshas been madein the past to enhance the performance of compositional simulators.Specifically, severaleffortswere madetodescribetheentiresystemonthebasisofpseudo-ternaryrepresentation[13,14].Such repre-sentationswerefoundtobevalidforgeneralcompositionalproblem[15].Itwasalsoshownthatthesolutionrouteenters andexits thetwo-phase region by ashockandthat theseshockscoincide withthetie-line anditsextension [15]. Later, an insightfuldemonstrationonthebasisofaunifiedMethodofCharacteristics(MOC)frameworkwascarriedoutonagas injectionprocesses[16,17].Thetheorywaslimitedtoadispersion-free1Ddisplacementproblemanditwasshownthatthe entiredisplacementpathcanbeconstructedbasedonthekeytie-linesofthesystem.Itwasalsoproventhatthestructure andpropertiesofthetie-linedictatethesolutionrouteinthecompositionalspace[18].

Basedontheseinsights,ageneralapproachforsolvingacompositionalproblemwiththeassumptionofconstant parti-tioning coefficientswasproposedin[19].Thisapproachusedthefactthatthestructureofthesolutionintie-linespaceis independent ofchanges inhydrodynamicpropertiesofthesystem.However, thisfact wasbasedona generalkey tie-line theory andwas not appropriately proven. Anumericalframework describing similarideology (based onthe tie-line con-cept) was later implementedon moregeneralcase[20], wherethe original compositionalproblemwas projected onto a tie-linespace,andlater,theparameterspertainingtoatwo-phaseregionofEoSbasedsystemaredeterminedbypolynomial approximations.

Subsequently,furtherdevelopmentsweremadeindeterminingthesolutionofthecompositionalprobleminthetie-line space. Oneapproachwas theuseofmulti-linearinterpolationtechniquesto representthethermodynamicphasebehavior ofthesystemwhichwas referredtoastheCompositionalSpaceParameterization(CSP)[21]. Also,itwas proventhat the projection ofthesolutionto thetie-linespaceisinvariant tothehydrodynamic parameters oftheproblem[22].The CSP ideology eventually led to the development of CompositionalSpace Adaptive Tabulation (CSAT) method[21], where the tie-linesarestoredinatableandlaterlookeduptoobtainthephasestateofthemixture.ThekeytoCSATimplementation is based on the fact that the phase state identification in a gas injection processes involves only a limited number of tie-lines. The proposed techniquesnot onlyacceleratedthephase behavior computations,butalsoprovided thebasis for thedevelopmentoffullyEOS-freecompositionalsimulationapproach[23].

Another way to reduce the computational burden of compositional simulation is to reduce the number of equations solved percontrolvolume.Themoststraightforwardwayisbyadoptingan ImplicitPressureExplicitComposition(IMPEC) reduction[24,25].Inthisapproach,theonlyglobalequationthatneedstobesolvedineachgrid-blockisthepressure equa-tion. Theframework isconstructed usingan IMPECreductionof fullconservationequations withexplicitly approximated composition. Whilethisapproachhelpstoreducethenumberofglobalequations,itintroduces asevererestrictiontothe computationaltimestepduetotheexplicittreatmentofhyperbolicvariables[26].

ThislimitationcanbeovercomebyapplyingAdaptiveImplicitMethod(AIM)[27,28].Inthismethod,theIMPESreduction isappliedtoonlythosegridcellswhereCourant–Friedrich–Lewy(CFL)condition[29] isbelowone,andtherebyremoving the time step restriction of the IMPES method.The AIM approach can significantly reduce the number ofequations, in the casewhena significantnumberof control-volumesstayscloseto theoriginal state [30].However, theAIMapproach introduces an error atthe interface between explicitand implicitblock which can affectnonlinear convergence [31]. In addition,theAIMstategystronglyreliesonnaturalvariablesformulation.

Inthiswork,weproposeamultiscalestrategyforthesolutionofacompositionalproblemusingmolarformulationwith significant changesincomposition. Weintroduce twodifferentsetsofmultiscaleoperators basedonthe dynamicsofthe hyperbolic problem. Initially, afine scalerestrictionoperator basedon MOCtheory [18] isappliedto reduce thenumber of conservationequationsto asingle transport equation, whichhelpsto reconstructboundariesofthe two-phase region. Oncetheseboundariesareidentified,thefullsolutioninsingle-phaseregionsisreconstructedusinganinterpolation-based prolongation operator. Next,the solutionfor thefull compositional problemin thetwo-phase region is reconstructedby solvingjusttwotransportequationswithrestrictionandprolongationoperatorsbasedonCompositionalSpace Parameteri-zationtechnique.

The paperis organized asfollows: abrief overviewofthestandard EOSbased simulationapproach isdescribed first. Then the methodology and framework pertaining to the multiscale reconstruction strategy of a compositional problem are explained.Later,their applicationtomulticomponent gasinjectionproblems(includingvaporizingandcondensinggas drives)alongwiththesupportingresultsarediscussed.Theconclusionandscopeforfurtherdevelopmentsaregiveninthe lastsection.

2. Conventionalcompositionalapproach

Thesolutionframeworkforthecompositionalproblemcanbebroadlyseparatedintotwosections,namelythe thermo-dynamicpartandthehydrodynamicpart.

2.1. Hydrodynamicmodel

The hydrodynamic framework attributes to the solution of flow and transport equations, thereby determining the pressure (p) and the compositional changes (z) of the system. For an isothermal compositional system consistingof nc

componentsandtwohydrocarbonphases(vaporandliquid)[23],theconservationequationsofmolarmasscanbewritten as

∂

∂

t

φ



(

xi

ρ

oSo

)

+ (

yi

ρ

gSg

)



+ ∇ ·



(

xiuo

ρ

o

)

+ (

yiug

ρ

g

)



+



(

xi

ρ

oqo

)

+ (

yi

ρ

gqg

)



=

0

,

i

=

1

, . . . ,

nc

.

(1)

Thevelocityofthefluidphase,neglectinggravityandcapillarity,isdescribedbasedonDarcy’slawas uj

= −

K

kr j

μ

j

∇

p

,

j

=

1

,

2

.

(2)

Fora1Dincompressibleisothermalproblem,thegoverningmassconservationequationcanbesimplified[21] asshown in

∂

zi

∂

t

+

ut

φ

∂

Fi

∂

x

=

0

,

i

=

1

, . . . ,

nc

.

(3)

Wewillbeusingthisformforsimplicityofdescriptionwithoutlimitingageneralapplicabilityoftheproposed technique. Hereziisdefinedas

zi

=

xi

ρ

oSo

+

yi

ρ

gSg

ρ

oSo

+

ρ

gSg

,

(4)

andcanbefoundfromthevaporfraction

v

=

ρ

gSg

ρ

oSo

+

ρ

gSg

,

(5)

zi

=

xi

(

1

−

v

)

+

yiv

.

(6)

The‘Fi’isfurtherexpandedbasedonthelinearrelationshipasshown,

Fi

=

xi

(

1

−

fg

)

+

yifg

,

i

=

1

, . . . , (

nc

−

1

).

(7)

Finally,theauxiliaryrelationsusedinclosingthesystemofequationsaregivenby

nc



i=1 xi

=

1 & nc



i=1 yi

=

1

,

(8) So

+

Sg

=

1

.

(9) 2.2. Thermodynamicmodel

Thethermodynamicframeworkenablesustodeterminethephasebehaviorofthesystemwiththeaidofphasestability test[2] –toidentifythenumberofphaseinaparticulargridcellandflashcalculation[3] –todeterminethesplitfraction of components amongst the phases present. The flash calculation can be carried out either based on equilibrium ratios (alsoknownasKvalues)orbasedonequationsofstate (EoS).Inthisresearchwork, wehavefollowedtheK-valuebased methodologytoperformflashcalculations.Thisassumptiondoesnotlimittheapplicabilityoftheproposedtechniquetoa generalEoS-basedsystem.

InaK-valuebasedsystem,thegeneralphasesplitprocedureisreducedtothesolutionofthefollowingequation: h

(

v

)

=

nc



i=1



zi

(

Ki

−

1

)

v

(

Ki

−

1

)

+

1



=

0

.

(10)

The nonlinear equ. (10) is also referred as the Rachford–Rice equation [32], which can be solved by a combination of bisectionandNewton’smethodtodeterminethevaporphasefraction(v).Itisalsoworthnoticingthattheconceptofflash calculationisnotlimitedtotwo-phasemixturesaloneastheycanalsobeperformedonsingle-phasecompositions[32,33].

Such a phasecomputation is generallyreferredto asthenegative-flashcalculation,where theequ. (10) issolved forthe vaporphasefraction(v)ofsingle-phasecompositions.Theideabehindthenegativeflashcalculationisthatasingle-phase mixtureispresentasalinearcombinationofphasecompositionanditalwayslies ontheextensionofthetie-lines.Inthis study,wehaveadoptedthenegativeflashbasedframeworktodeterminethephasebehaviorofthesystem.

Inordertoclosethesystemofequations,additionalconstraintsbetweenthecomponentsandthephasesexistunderthe statement of instantaneousthermodynamic equilibrium–which isinturn expressedasequalityinfugacity ofthevapor andliquidphases[34],referequ. (11) below

fi,g

(

p

,

T

,

yi

)

=

fi,o

(

p

,

T

,

xi

).

(11)

Thefugacityisafunctionofpressure(p),temperature(T )andthephasecompositions(xi,j)anditisgenerallydetermined

basedontheequationsofstate(EoS).

A sequentialframework, explicitin time, isused insolving thesystem ofequationsequ. (3), (9) withvarious closing relations,therebydeterminingthesolutionfortheunknownvariables(1pressure,two-phasesaturationsandnc

×

np phase

compositions)ineachgridblock.

3. Multiscalereconstructioninphysics

In ordertogain betterinsight intotheproposed multiscaletechnique, thebasicspertaining totheconstruction of so-lution is initially discussed. Ina compositional formulation, thephysical systemis represented asa mixture of different componentswithdifferentvolumefractions.Forinstance,ina3-component system,thesemixturesare referredtoasthe light,intermediateandheavyhydrocarboncomponent.Atgiventhermodynamicconditions(pressureandtemperature),the light componentispresentina gasphaseandtheheavy componentispresentina liquidphase.Theintermediate com-ponent isgenerallypresentin eitherofthesetwo phases(liquidorgas). Thephase behavior ofsuch systems(at afixed pressure andtemperature)is directlydrivenby the compositionalchanges anditisgenerallyrepresented ona phase di-agram. Fig. 1showsthe phasediagram, the fractional flowcurve andthe compositionalprofile oflight componentfora ternarygasinjectionprocess.

Thedewpointandbubblepointcurves(blueandgreenlineinFig.1(a)),whicharecollectivelyreferredasthebinodal curve,distinguishthesingle-phasecompositionswiththoseexistingintwophases.Thetwo-phasecompositionslieinsideof thebinodalcurveandtheyarefurthersplitintoliquidandvaporphasefraction.Thesephasefractionsinsidethetwo-phase region are eventually connected by a tie-line (blackdotted linein Fig.1(a)). Both the binodal curve andthe equilibrium tie-linesarestrongfunctionsofpressureandtemperaturewithweakdependencyoncomposition.

Further, thetrajectory (or path) that describesthecompositional changes betweentheinitial oil compositionand the injected gascomposition isreferredtoasthecompositionalpath.Thesalientprinciplesthat govern theconstructionofthe solution (compositional path) are the conservation andentropy conditions, resulting in a series ofself-sharpening waves (shocks) andspreading waves(rarefaction)asshowninFig. 1.Thesechangesare generallyconsistent withtheabove de-scribedphysicalprinciples.

The fractionalflowtheory givesinsightonthefluidmobilityinthepresence ofother fluidphases. Fig.1(b)showsthe plot of fractional flow as a function of composition. The two curves, shownin the plot, refer to fractional flow on the injection(redcurve)andinitialtie-line (bluecurve)fora three-componentsystem, andtheconstruction ofsolutionpath issubsequentlyexplained.Thevelocity(a.k.a. slopeofthefractionalflowcurve)ofeachcompositionfollowsasimplerule where thecharacteristic velocityshould monotonicallyincrease frominjectionto initial composition.This correspondsto thefactthatinthestablesolution,theslowerwavescannotpropogateaheadofthefasterwaves.Iftheslowerwavesfrom compositionsclosetotheinitialconditionsoriginateaheadoffasterwaves,ashockwillformasthefasterwavesovertake the slower waves.Hence, by satisfyingthis condition, theconsistent solution contains two shocksconnecting initial and injectioncompositionsanditisgenerallyindependentofthenumberofcomponentspresentinthesystem[16,17]. These twoshocksarereferredtoastheleading andtrailingshocks.

Fromaphysicalperspective,theleadingshockisformedbecauseoftherapidmovementoftheinjectedfluidwhen com-paredwiththeinitialoilcomposition,andthetrailingshockisformedduetoevaporationoftheintermediatecomponents intotheunsaturatedinjectedfluid.Thevelocityofthesetwoshocksdependsontheinjectionandtheinitialcompositions. Betweenthesetwoshocks,thefanofcharacteristicsassociatedwiththecontinuousvariationofcompositionforms,which isalsoknownasthespreading orrarefaction waves[15].Thiscontinuousvariationinthecompositionscanoccuralongthe tie-line(known astie-linerarefactionpath),andaswellasbetweenthetie-lines(knownasnon-tie-linerarefactionpath), referFig.1.

FurtherinsightintotheshocksolutionisdescribedbasedontheRankine–Hugoniotcondition[15],whichisgivenbythe equ. (12):



=

dFi dzi

=

Fi,2

−

Fi,1 zi,2

−

zi,1 i

=

1

, . . . ,

nc

,

(12)

where, Fi,1 & zi,1 representsthe fractional flow andcomposition ofcomponent‘i’ onthe downstreamside ofthe shock

Fig. 1. General

structure of compositional solution: (a) ternary phase diagram, (b) fractional flow curve and (c) solution profile (light component). (For

interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Hugoniotconditionisanintegralversionoftheconservationequationshowninequ. (3) anditillustratesthefactthatthe volumeisconservedacrosstheshock.

Now, the outset ofmultiscale based reconstruction technique is attributedto the saturation history(Fig. 1(c)) ofthe compositional problem, since the profilecan be classified intothree distinct regions based on thecharacteristic changes [11]. The region1 corresponds toa placid zone, wherethe gassaturation value remains (almost)constant. The region2 correspondstoshocksolution,wherethereisapronouncedjumpinthesaturationvalues.Thecharacteristicsofthisregion are defined based on the Rankine–Hugoniotcondition, see equ. (12). Finally, the region 3 refers to a fan characteristics (rarefactions),wherethesaturationchangesaremorecontinuous.Theregion1andregion3arereferredtoas“single-phase region”and“two-phase region”respectively,whileregion2isa“transitionzone”.Onthebasisofthisclassification,itcan be observed that the shockseparates two contrasting zonesofsaturation, withthese zoneslocated on eitherside ofit. Thisapprehending factforms thebasis ofthe proposed multiscaletechnique, wherespecially constructedrestriction and prolongationoperatorsareimplementedonthedifferentregionsofsaturationtoreconstructthesolutionofthehyperbolic problem(fullsystem).

Initially, we reconstruct theposition of shock inspace as they hold the key, and it is carriedout using a fine scale restrictionoperator. Oncethepositionofshocks isidentified,theinterpolation-basedprolongation operatorisapplied for the reconstruction of the solution in single-phaseregions similar to the approach suggested in [11]. Later, the solution for nc components lying in between the leading andtrailing shock is reconstructed with the aid ofthe restriction and

prolongationoperatorsdefinedbasedonCSP technique[21].Both vaporizingandcondensinggasdrives(referTable1& 2 forcompositions)havebeenimplementedbasedontheproposed approachandthecorrespondingresultsarediscussedin thesubsequentsections.

Table 1

Vaporizing gas drive.

Comp. C1 CO2 C4 C10

Initial 0.10 0.18 0.37 0.35 Injection 0.97 0.01 0.01 0.01 K-value 2.5 1.5 0.5 0.05

Table 2

Condensing gas drive.

Comp. C1 CO2 C4 C10

Initial 0.01 0.01 0.39 0.59 Injection 0.59 0.39 0.01 0.01 K-value 2.5 1.5 0.5 0.05

Fig. 2. Analytically constructed fractional flow curve for light component (CH4) in a vaporizing gas drive.

4. Multiscalereconstructionofsingle-phaseregion

Thetrackingoffrontpropagationforacompositionaltransportproblemiscarriedoutbasedonanalyticallyconstructed fractionalflowcurves,whichinturnplaytheroleofrestrictionoperator.Initially,theliquid(xi)andthevapor( yi)fraction

values ofthecompositionslyingon theextension ofinjectionandinitial tie-lines(single-phaseregion)are selectedfrom the negative flash calculation.Further, by keepingthisxi and yi valuesconstant, thefractional flow curve is analytically

constructedinthetwo-phaseregionbasedonequ. (7) [20].

The fractional flow plot is an S-shaped curve within the two-phase region. In the single-phaseregion, the fractional flow value ofa particularcomponentvarieslinearly withthemolar fractionofthe corresponding component. Hence the fractionalflowcurveisjustastraightlineinthesingle-phaseregion(referFig.1(a)).Twofractionalflowcurvesforachosen component I correspondingtoinitialandinjectiontie-linecanbefoundfrom:

Fini_I

=

xini_I

(

1

−

fg

)

+

yiniI fg

,

FinjI

=

x inj

I

(

1

−

fg

)

+

yinjI fg

.

(13)

Once thesecurves areanalyticallyconstructed basedonthe keytie-lines ofthesystem, they areeventually combined together,andanequivalentsolutioncurveisconstructedbytakingaconvexhullontheunionofbothcurves:

FR

=

conv

(

FinjI

∪

FiniI

).

(14)

Thisprocedureresultsinasingleanalyticcurve,thatswitchesbetweentheinitialandinjectionfractionalflowcurveatthe pointoftheirintersection.

The analyticallyconstructed curve shownin Fig. 2 represents thepseudo-fractional-flow curve for thereconstruction ofleading andtrailingshocks.Later,thesevaluesaretabulated anddefinedasrestrictionoperatorduring thecourse ofa simulation,therebyenablingustoaccuratelylocatetheshockpositionsinspace.Thesimulationframeworkisanalogousto afullyimplicitsolverexceptforthefactthatthedefinedrestrictionoperatorisappliedfortherestrictedtransportequation:

∂

zR

∂

t

+

ut

φ

∂

FR

∂

x

=

0

.

(15)

4.1. Applicationfor1Dtransportproblem

Fig. 3showsthe reconstructedshocksolutionfora 4-component system(vaporizing gas drive).The resultsofthe re-constructedsolutionareevaluatedbycomparingitwithareferencefull-physicssolution.Itcanbeobservedthattheshock

Fig. 3. Analytically

constructed fractional flow curves (on left), and the corresponding reconstructed solution of transport problem (on right) for (a), (b)

light component C1, (c), (d) intermediate component C4and (e), (f) heavy component C10.

positionsareaccuratelyreconstructedtofine-scaleaccuracywiththeuseoftherestrictionoperator.Ontheotherhand,the solutionbetweenthe shocks(in thetwo-phase region)differs.Thoughthere isa validreasonbehind thismismatch(see discussioninAppendixA.2),thisdifferencecanbeneglectedsinceourobjectiveislimitedtotheaccuracyofshock

recon-struction. Inthisstage, ourmultiscaleapproachreducesthecomplexity (numberofequationssolved per controlvolume) from

(

nc

)

conservationequationstotwo:

[

nc

×

nb

⇒

2

×

nb

]

,wherenb isthenumberofgridblocksandnc isthenumber

ofcomponentsinthesystem.Further,theproposedmethodologyremainsvalidforanychangesintheinjectionandinitial streams. Undersuch circumstances, the negative flash calculation (equ. (10)) is performedon the new (injection/initial) composition,andequ. (13) –equ. (15) iseventuallysolvedtodeterminethecorrespondingphaseboundariesofthesystem. Subsequently, thesolution inthesingle-phaseregion isreconstructed usingtheprolongation operatorbased on inter-polationofentiresolutionbetweeninjectionandinitialcomposition,andbyusingtheconcentrationofpseudo-component fromthesolutionoftherestrictedproblem:

κ

(

zR

)



R

1

_{=⇒ R}

nc−1



_:

_z

₌

_I

{zini_,zinj_}

(

zR

).

(16)

Here, ‘

κ

’ istheinterpolation-prolongationoperator, ‘zR’isthe reconstructedshocksolutionand‘I’isthepiecewise linear

interpolation function. Thisprolongation will matcha conservativesolution insingle-phase regions dueto theweak de-pendency betweeninitial/injection compositions andcompositions atshocks, controlled by Rankine–Hugoniot conditions showninequ. (12).Noticethat,conceptuallythisapproachissimilar totheprolongationoperationforreconstructionofa fine-scalespatialsolutionproposedin[11].

4.2. Applicationfor2Dtransportproblem

The proposed firststage reconstruction was alsotestedforflowandtransport compositionalproblemina highly het-erogeneousreservoir.Forthistest,thetoplayerofSPE 10modelwithheterogeneousporosityandpermeabilitymapswere used [35]. Themodelis builtona regular Cartesiangrid (60

×

220 cells),withonegas injectionwellin thecenterand four productionwellsatcorners. The implementationofthe firststage reconstruction was carriedout based onthe OBL technique described in[36,37].Here, we translate aconstructed pseudo-fractional flowcurve to convectionoperator and tabulateitasafunctionofnonlinearparameters(p andzR).A pseudo-binarysystemissubsequentlysolved,thereby

yield-ingthesolutionforpressureandindicatorcomposition.Later,thesolutionforthefullsystem(nc components)iseventually

reconstructed by linearlyinterpolating the full setof compositions betweenthe initial andinjection mixtures using the solutionofindicatorcompositionzR.Oncethefullcompositionisreconstructedusing(16),thenegativeflashprocedureis

appliedtopredictphasebehaviorandlocatetheboundariesofthemultiphaseregion.Insummary,irrespectiveofthesize of thesystem, a pseudo binary modelis sufficientenough to determinethe phase boundariesforthecorresponding full system.

Fig. 4showsthereconstructed shocksolutionofSPE10model attwodifferentsteps, forvaporizinggas drivesystem (mainparametersaregiveninTable5andFig.10).Thesubplotinthefigureshowsthe(1)compositionalchanges,(2)the phasesstateineachgridblockand(3)thedifferenceinphasestateidentification.Thelimitedmismatchinphaseboundaries betweenthefullconventionalapproachandtheproposedmulti-scalereconstructiondemonstratestheprofoundnessofthe proposedfirststagereconstruction.

5. Multiscalereconstructionoftwo-phaseregions

Thesolutioninthetwo-phaseregionisdeterminedbasedonCompositionalSpaceParameterizationtechnique [21].The CSP approachwasproposed tospeedupthephasebehaviorcalculationsbyreplacingtheflashcalculationwith interpola-tions intheparameterspaceoftheproblem.Thephasebehaviorofgas-injectionprocessesispredominantlycontrolledby the propertiesoftwo keytie-lines that extendthroughthe initialandinjectioncompositions, andhenceitis convenient to parameterize the problem basedon these two tie-lines. It has alsobeen proven that the projection of compositional solutions withdifferenthydrodynamicpropertiesontothetie-linespaceisinvariantandonlydependsonthermodynamic propertiesofsolution[22].

On the basis of parameterization technique, an advanced simulation framework has been proposed for the solution of a hydrodynamic problem, thereby reconstructing the solution in two-phase region [23]. This reconstruction is based on the invarianceof a compositional solutionin tie-linespace. Notice that this technique requireseither thesolution of auxiliarythermodynamicproblem[20],orthesolutionofafullproblemanditsprojectionontothetie-linespace[21].Once constructed,thissolutionandits projectionprovidethebasisformultiscalereconstructionindependentofanychanges in thehydrodynamicproperties.

Initially, the existing tie-lines of the solution are parameterized based on ‘

γ

’ variables. By this, each tie-line in the system isuniquely identified. It isworth noticing that theparameterizing variables (

γ

) can be definedarbitrarily in the compositionalspace,exceptforthefactthatthisdefinitionshouldbeunique.Forinstance,theycanbedefinedasatie-line traceonthefaceofcompositionaldiagram[20],orataplanerepresentingthetie-linecenters[21].Also,the

γ

isconstant along atie-line,andchanges onlywhenthesolutionis switchedbetweentie-lines.Once,thetie-linesare parameterized, the

γ

variable(vectorquantity)isplottedontoatie-linespace(

γ

-space),andtheresulting

γ

-pathpresentdiscontinuities (shocksandrarefactions)associatedwiththekeytie-linesofthesolution.InFig.5,the

γ1

and

γ2

refertothephasefractions ofintermediatecomponents(CO2& C4forthiscase)inafour-componentsystemintroducedabove.

Fig. 4. First

stage reconstruction in SPE-10 vaporizing gas drive for two different timesteps:

(upper) composition of light component, (middle) phase

identification and (lower) blocks with misidentified state.

Fig. 5. Parameterized space (γ-path): (a) vaporizing gas drive and (b) condensing gas drive.

Oncethe



-pathisconstructed,anewvariable‘

σ

’isintroduced[38].Theobjectiveofthis‘

σ

’variableistoparameterize the

γ

-pathand the sameis accomplishedby orthogonally projecting the ‘

γ

’ valuesonto ‘

σ

’such that itranges from 0 and 1(Fig.5).Inthisrepresentation,allphasefractionsbecomenonlinearwith yi

(

σ

)

andxi

(

σ

)

.

Followingthisprocedure,thecompositionalvariationoftheentireproblemisparameterizedintie-linespace.Noticethat thissolutionis invariant ofchangesin hydrodynamicproperties oftheproblemandonly dependenton thermodynamics and injection/initialconditions [20,22]. Now thesolution for the two-phase region (area within the leading andtrailing shocks)isdeterminedbysolvingthetransportproblemintheparameterizedspace.First,werepresentoverallcomposition atnewtimestepthroughthenonlinearrelationinvolvingphasefractions:

zn_i+1

=

yi

(

σ

n+1

)

·

vn+1

+

xi

(

σ

n+1

)

· (

1

−

vn+1

).

(17)

Now,we cansubstitutethisrelationintotwoconservationequationsfrom(3) indiscretizedform. Theresultingrestricted systemtakesthefollowingresidualform:

Ri

=

yi

(

σ

n+1

)

·

vn+1

+

xi

(

σ

n+1

)

· (

1

−

vn+1

)

−

zn

−

dt

φ

·

dx

{(

Fi

·

Ui+1

)

− (

Fi−1

·

Ui

)

}

n _i

₌

₁

_,

₂

_.

₍₁₈₎

Any two componentscanbe arbitrarilychosen tosolve equ. (18),herewe selectlight andtheheavy component. Fur-ther, anexplicit fluxapproximationtechnique hasbeenused tosolvethe non-linearsystem.Notice, thatthe typeofflux approximation doesnot limittheapplicationofthe proposed approach.Theequ. (18) isnonlinear innature,evenforan explicitfluxapproximation,anditisinducedbythetie-lineparameters(xiand yi).Asaresult,Newton’smethodhasbeen

appliedtoresolvethesystemandobtainthesolutioniteratively.Throughthistechnique,thenc

−

1 transportequationsof

the originalproblemare replacedby justtwo equationsdependingon twovariables only(

σ

and v). Inother words,the equ. (18) issolvedwithparameterizedvariable

σ

andvaporfractionv asitsprimary unknownsinsteadofdirectlysolving for thecompositional changes(zi), thereby resultinginthe restrictionof afull system. The final systemofequations(in

algebraicform)isgivenbyequ. (19).



_∂R 1 ∂σ ∂ R1 ∂v ∂R2 ∂σ ∂ R2 ∂v

 

δ

σ

δ

v



=



−

R1

−

R2



.

(19)

In the previous stage, the shock positions are located in space based on the solution with restriction operator and compositions insingle-phaseregionare reconstructedusingprolongationoperator basedoninterpolation.Now, the solu-tion lies in thetwo-phase regionis restored by solving two restrictedconservation equations(18) in the parameterized space.

Oncethesolutionisobtained,theprolongationoperatorisappliedagainbasedontheprojectionof

σ

solution(points ontheredcurve)backto

γ

-path(pointsonthebluecurve),seeFig.5fordetails.Thisoperatorcanbeexpressedas

ψ(

v

,

σ

)

R

2

_{=⇒ R}

nc−1

_:

_z

₌

_y

₍

_σ

₎

_·

_v

₊

_x

₍

_σ

₎

_{· (}

₁

₋

_v

_).

₍₂₀₎

Notice that

γ

-parametersuniquelyidentifytie-linesinthecompositionalspace,whichdefinesthecorresponding valuesof

x and y.Oncetie-lineisreconstructedbasedon

σ

,thecorrespondingoverall compositionsz arerecoveredusingsolution in‘v’and‘equ. (20)’.

Fig. 6andFig.7showsthesolutionofa vaporizingandcondensinggasdrive problem.It canbe seenthat the recon-structedtwo-phaseregion(’Z-CSP’curve)accuratelymatcheswiththereferencefinescalesolution.Itisalsoworthnoticing thatthereconstructedsolutionoftheentirecompositionalsystemisobtainedbysolvingtherestrictedsystem (19) instead ofthefullsetofconservationequations.

This makes thetotalcost offull multi-scalereconstruction

[

nb

×

2

+

nb,2

×

2

]

insteadof

[

nb

×

nc

]

,wherenb andnb,2

correspondtothetotalnumberofblocksandthenumberofblocksintwo-phaseregionrespectively.Onthewhole,thecost offullreconstructionfornc-componentcompositionalsolutionbecomesonlyslightlyhigherthanthesolutionofthe

black-oilproblem.Also,thedecouplingbetweenthethermodynamicandhydrodynamicproblemsprovidesasimplerplatformto carry outamulti-dimensionalanalysisandsensitivitystudies,sincethe



-pathisindependentofflowpropertiesandonly dependsonthermodynamicparametersofthesystem.

Further,basedonthedescriptionoftheparameterizationtechnique,itisalsoevident thattheconstructionof

γ

-space is of paramount importance – which requires unique identification and parameterizing of the tie-lines involved in the conservativesystem.This,inturn,impliesthefactthatwithanychangesinthethermodynamicsofthesystem(initialand injectioncompositions),thetielinesinvolvedintheconstructionofconservativesolutionalsochanges,therebyresultingina different

γ

-path.Fig.8showsthechangesinthe

γ

-pathforasystemwithvaryinginitialoilandinjectiongascompositions showninTables3,4.Thesecompositionalpaths canbeparameterizedindependentlyandusedintheOBLframework[36, 37],therebydirectlyincorporatingthecompositionalchangesandthusmakingtheproposedtechniquerobust.

Also, the parameterization technique remains valid irrespective of the size of system, and it is enough to solve the same restrictedsystemofequations withv and

σ

asprimary unknowns. Theonly difference isinthe dimensionof the gamma-space,whereinvariantpathisconstructed.Inotherwords,forasystemofnc components,thecompositionalspace

parameterization technique results ina tie-linespace (

γ

-path) having (nc

−

2) spatial dimensions[20]. Forexample, the

Fig.9belowshowsthegammapathfora5-component systemandthecorrespondingresultsofreconstructionareshown inAppendixA.3.

6. Conclusion

Through thisresearch work, a multiscale technique has been proposed to enhance the performance ofcompositional simulation.The methodologyisbasedonthelocalcompositionalprofileofthesystemandeventually,atwo-stage recon-structionofthecompositional problemiscarriedout.Inthisstudy,we presenttheapplicationoftheproposed approach foranidealizedcompositionalproblemwithconstantpartitioncoefficients.

Fig. 6. Full

reconstruction of four-component problem for vaporizing gas drive (a) light component CH

4(b) first intermediate component CO2(c) second

intermediate component C4(d) heavy component C10.

First,we presentthetechnique forpredicting thefront (shock)positions inspaceusingtherestriction operatorbased onfractionalflowtheory.Theapproachwastestedforseveral1Dgasinjectionproblemsofpracticalinterest,includingthe vaporizingandcondensinggasdrives,anditaccuratelypredictsboundariesofthetwo-phaseregion.Aninterpolation-based prolongationoperator was applied forthereconstruction ofthe single-phaseregion,thereby precisely matchingthe con-servative solution of thefull problem inthis region. Subsequently, theapproach was alsoextended to multidimensional problemswithhighlyheterogeneouspropertiesusingOperator-BasedLinearization(OBL)formalism.

The second stageof theproposed multiscalestrategy isapplied toreconstruct thecompositional solutioninthe two-phaseregion.ThisreconstructionisperformedusingtherestrictionoperatorbasedonCompositionalSpaceParametrization (CSP)technique.Therestrictionoperatorreducesthefullsystemofconservationequationstotwoequationspereach con-trolvolumeinthetwo-phasestate.Based onthesolutionofthisreducedsystem,aprolongationinterpolationoperatoris appliedinthetie-linespace,whichiscapableofreconstructingaconservativefinescalesolutionofthefullsystem.

Whiletheproposedmultiscaleapproachisfoundtoberewarding,thereareseveralpossibilitiesforitsfurther improve-ments. The first and theforemost advancement wouldbe to implement thistechnique in a fullyimplicit compositional frameworkinvolvingcomplexnonlinearphysics.Theequationsofstatecanalsobedirectlyincorporatedforthephase be-havior of the system. In addition, the approach could also be extended to problems with a strong initial compositional gradientorchangesininjection/initialstreamandthesamecanbearchivedbythedirectapplicationofOBLmethodology. Whilethecurrentimplementationislimitedtoimmisciblegasinjectionprocess,theextensionofthemethodologyto near-miscibleregimesisstraightforwardandisthefocusofourcurrentwork.Finally,asimilarapproachcanalsobedeveloped forthermalproblemswithanarbitrarynumberofphases.

Fig. 7. Full

reconstruction of four-component problem for condensing gas drive: (a) light component CH

4, (b) first intermediate component CO2, (c) second

intermediate component C4and (d) heavy component C10.

Table 3

Variations in initial oil composition.

Ini. comp. C1 CO2 C4 C10 Case 1 0.10 0.18 0.37 0.35 Case 2 0.10 0.15 0.50 0.25 Case 3 0.30 0.30 0.25 0.15

Table 4

Variations in injection gas composition.

Inj. comp. C1 CO2 C4 C10 Case 1 0.97 0.01 0.01 0.01 Case 2 0.85 0.13 0.01 0.01 Case 3 0.80 0.18 0.01 0.01

Fig. 9. Parameterized space (γ-path) for five-component system.

Nomenclature

Symbol Description Units

k Absolute permeability Darcy

μ

g Viscosity of gas Pa S

μ

o Viscosity of oil Pa S

ρ

g Density of gas Kg/m3

ρ

o Density of oil Kg/m3

φ

Porosity –

Ki Equilibrium ratio of component i –

v Mole fraction of vapor phase Mol Fr

L Mole fraction of liquid phase Mol Fr

γ

Tie-line parameter –

fg Fractional flow of gas phase –

Fi Fractional flow of component i –

Sg Saturation of gas –

So Saturation of oil –

λ

j Mobility of phase j –

uj Darcy velocity of phase j m/s



Shock velocity m/s

xi Composition of component i in liquid phase Molar Fraction

yi Composition of component i in vapor phase Molar Fraction

zi Overall composition of component i Moles

p Pressure N/m2

σ

Parameterized variable –

nc Number of components –

np Number of phases –

Acknowledgement

Table 5

Reservoir fluid properties.

Properties Gas Oil

Residual saturation (Sjr) 0.0 0.0

End point relative perm. (Krej) 1.0 1.0

Saturation exponent (nj) 2.0 2.0

Phase viscosity in cP (μj) 1.0 5.0

Fig. 10. Upper layer of SPE10 test case: (a) porosity and (b) log of permeability.

Fig. 11. Analytical fractional flow curve based on non-tie-line rarefaction.

Appendix A

A.1. Simulationparameters

InthisAppendix,wepresentsimulationparametersforallperformedsimulationsinadditiontotheparametersalready describedintheresultsection.

A.2. Nontie-linerarefaction

Theconstructionofafractionalflowsolutioncurveissimilartotheconstructionofacompositionalpath,whichinturn isgenerallyrepresentedbya seriesofshocksandrarefaction(bothtie-lineandnon-tie-line).Basedonthis, thefractional flow solutioncurve (Frac.SolcurveinFig.11)isconstructedincompliancewiththeconservationlawsasshowninbelow. Thecardinaldifferencebetweenthetracedsolutioncurve,implementedintheconventionalsimulation,andtheoneshown hereisrelatedtothepathofnon-tie-linerarefaction(indicatedbythestraight lineinbetweenthetangentdots).Though thetracedsolutioncurveindicatedaboveensuresconservation,ontheotherhand,locatingthesepointsonthecurvemakes theapproachcumbersome.Hencethesolutioncurveisconstructedbytracingtheinjectionandtheinitialflowcurveswith a switch at their point of intersection, without tracing the non-tie-line path (sinceit is sufficient enough to accurately

Fig. 12. Reconstructed

two-phase region of 5C system, vaporizing gas drive (a) light component CH

4, (b) first intermediate component CO2, (c) second

construct or trace the fractional flow solution curve with respect to the shockregions). As a result, the solution profile withinthetwo-phaseregionissporadic.

A.3. Reconstructionof5Csystem

The modified transport equationbased ontheCSP technique issolved over the3D parameterizedspace(as shownin Fig.9) andthecorrespondingresultsofreconstructionareshowninFig.12.Basedontheplots,wecanconcludethatthe proposed technique yields exactresultsasthat ofaconventional fine-scalesimulation.Also, theresultsreinstate thefact thatthedevelopedmultiscaleapproachresultsincoarseningofsystemrepresentation,therebyenhancingtheperformance.

References

[1]K.Aziz,T.Settari,PetroleumReservoirSimulation,AppliedSciencePublishers,1979.

[2]M.Michelsen,Theisothermalflashproblem.PartI.Stability,FluidPhaseEquilib.9 (1)(1982)1–19.

[3]M.Michelsen,Theisothermalflashproblem.PartII.Phase-splitcalculation,FluidPhaseEquilib.9 (1)(1982)21–40.

[4]K.Coats,Anequationofstatecompositionalmodel,SPEJ.20 (05)(1980)363–376.

[5] T. Hou, X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1) (1997) 169–189, https://doi.org/10.1006/jcph.1997.5682.

[6] P. Jenny, S. Lee, H. Tchelepi, Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, J. Comput. Phys. 187 (1) (2003) 47–67, https://doi.org/10.1016/S0021-9991(03)00075-5.

[7] Y. Efendiev, V. Ginting, T. Hou, R. Ewing, Accurate multiscale finite element methods for two-phase flow simulations, J. Comput. Phys. 220 (1) (2006) 155–174, https://doi.org/10.1016/j.jcp.2006.05.015.

[8] H. Hajibeygi, G. Bonfigli, M. Hesse, P. Jenny, Iterative multiscale finite-volume method, J. Comput. Phys. 227 (19) (2008) 8604–8621, https://doi.org/10. 1016/j.jcp.2008.06.013.

[9]O.Møyner,K.-A.Lie,Amultiscalerestriction-smoothedbasismethodforhighcontrastporousmediarepresentedonunstructuredgrids,J.Comput. Phys.304(2016)46–71.

[10]S.Shah,O.Møyner,M.Tene,K.-A.Lie,H.Hajibeygi,Themultiscalerestrictionsmoothedbasismethodforfracturedporousmedia(f-msrsb),J.Comput. Phys.318(2016)36–57.

[11]H.Zhou,S.Lee,H.Tchelepi,Multiscalefinite-volumeformulationforsaturationequations,SPEJ.17 (01)(2012)198–211.

[12] M. Cusini, B. Fryer, C. van Kruijsdijk, H. Hajibeygi, Algebraic dynamic multilevel method for compositional flow in heterogeneous porous media, J. Com-put. Phys. 354 (2018) 593–612, https://doi.org/10.1016/j.jcp.2017.10.052.

[13]G.I.Barenblatt,V.M.Entov,V.M.Ryzhik,TheoryofFluidsFlowThroughNaturalRocks,SpringerScienceandBusinessMediaDordrecht,1990,Series ISSN0924-6118,HardcoverISBN978-0-7923-0167-7,SoftcoverISBN978-90-481-4042-8.

[14]D.Tang,A.Zick,Anewlimitedcompositionalreservoirsimulator,in:SPESymposiumonReservoirSimulation,SocietyofPetroleumEngineers,1993.

[15]F.Helfferich,Theoryofmulticomponent,multiphasedisplacementinporousmedia,SPEJ.21 (01)(1981)51–62.

[16]W.Monroe,M.Silva,L.Larson,F.J.Orr,Compositionpathsinfour-componentsystems:effectofdissolvedmethaneon1DCO2floodperformance,SPE Reserv.Eng.5 (03)(1990)423–432.

[17]R.Johns,B.Dindoruk,F.J.Orr,Analyticaltheoryofcombinedcondensing/vaporizinggasdrives,SPEAdv.Technol.Ser.1 (02)(1993)7–16.

[18]F.M.OrrJr.,TheoryofGasInjectionProcesses,Tie-LinePublications,2007,ISBN87-989961-2-5.

[19]P.Bedrikovetsky,M.Chumak,Riemannproblemfortwo-phasefour-andmorecomponentdisplacement(idealmixtures),in:ECMORIII-3rdEuropean ConferenceontheMathematicsofOilRecovery,1992.

[20]D.Voskov,V.Entov,Problemofoildisplacementbygasmixtures,FluidDyn.36 (2)(2001)269–278.

[21]D.Voskov,H.Tchelepi,Compositionalspaceparameterization:theoryandapplicationforimmiscibledisplacements,SPEJ.14 (03)(2009)431–440.

[22]A.Pires,P.Bedrikovetsky,A.Shapiro,Asplittingtechniqueforanalyticalmodellingoftwo-phasemulticomponentflowinporousmedia,J.Pet.Sci.Eng. 51 (1)(2006)54–67.

[23]R.Zaydullin,D.Voskov,H.Tchelepi,Nonlinearformulationbasedonanequation-of-statefreemethodforcompositionalflowsimulation,SPEJ.18 (02) (2012)264–273.

[24]H.Kazemi,C.R.Vestal,G.Shank,Efficientmulticomponentnumericalsimulator,Soc.Pet.Eng.AIMEJ.18 (5)(1978)355–368.

[25]L.Nghiem,D.Fong,K.Aziz,Compositionalmodelingwithanequationofstate,SPEJ.21 (6)(1981)687–698.

[26] K. Coats, Impes stability: selection of stable timesteps, SPE J. 8 (2) (2003) 181–187, https://doi.org/10.2118/84924-PA.

[27] G. Thomas, D. Thurnau, Reservoir simulation using an adaptive implicit method, SPE J. 23 (5) (1983) 759–768, https://doi.org/10.2118/10120-PA. [28]T.F.Russell,StabilityAnalysisandSwitchingCriteriaforAdaptiveImplicitMethodsBasedontheCFLCondition,SocietyofPetroleumEngineers,

Jan-uary 1,1989.

[29] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen differenzengleichungen der mathematischen physik, Math. Ann. 100 (1) (1928) 32–74, https:// doi.org/10.1007/BF01448839.

[30]H.Cao,DevelopmentofTechniquesforGeneralPurposeSimulators,PhDThesis,StanfordUniversity,2002.

[31] R. de Loubens, A. Riaz, H.A. Tchelepi, Error analysis of an adaptive implicit scheme for hyperbolic conservation laws, SIAM J. Sci. Comput. 31 (4) (2009) 2890–2914, https://doi.org/10.1137/080724502.

[32]C.Whitson,M.Michelsen,Thenegativeflash,FluidPhaseEquilib.53(1989)51–71.

[33]A.Iranshahr,D.Voskov,H.Tchelepi,Generalizednegative-flashmethodformultiphasemulticomponentsystems,FluidPhaseEquilib.299 (2)(2010) 272–284.

[34]D.Voskov,Anextendednaturalvariableformulationforcompositionalsimulationbasedontie-lineparameterization,Transp.PorousMedia92 (3) (2012)541–557.

[35]M.Christie,M.Blunt,TenthSPEcomparativesolutionproject:acomparisonofupscalingtechniques,SPEReserv.Eval.Eng.4 (4)(2001)308–316.

[36] D. Voskov, Operator-based linearization approach for modeling of multiphase multi-component flow in porous media, J. Comput. Phys. 337 (2017) 275–288, https://doi.org/10.1016/j.jcp.2017.02.041.

[37]M.Khait,D.Voskov,Operator-basedlinearizationforgeneralpurposereservoirsimulation,J.Pet.Sci.Eng.157(2017)990–998.

[38]V.Entov,F.Turetskaya,D.Voskov,Onapproximationofphaseequilibriaofmulticomponenthydrocarbonmixturesandpredictionofoildisplacement bygasinjection,in:ECMORVIII-8thEuropeanConferenceontheMathematicsofOilRecovery,2002.