Application of control volume based finite element method for solving the black-oil fluid equations

Pet.Sci.(2013)10：361-372 361DOI lO.1007/s12182-013.0284.3Application of control volume based finiteelement method for solving the black-oil fluidequationsGHoREISHIAN AM IRI S A件，SADRNEJAD S A ，GHASEM ZADEH Hand M 0NTAZERI G HFaculty of Civil Engineering，K.N.Toosi University ofTechnology，Tehran，IranResearch and Development Dept.，Iranian Central Oil Field Co.，Tehran，Iran◎China University ofPetroleum(Beijing)and Springer-Verlag Berlin Heidelberg 2013Abstract：This is the second PaDer of a series where we introduce a control volume based finite elementmethod(CVFEM1 to simulate multiphase flow in porous media.This is a fully conservative methodable to deal with unstructured grids which can be used for representing any complexity of reservoirgeometry and its geologica1 obiects in an accurate and efficient manner.In order to deal with the inherentheterogeneity of the reservoirs.all operations related to discretization are performed at the elementlevel in a manner similar to classical finite element method(FEM).Moreover，the proposed method caneffectively reduce the so-called grid orientation effects.In the first Paper of this series，we presented thismethod and its application for incompressible and immiscible two.phase flow simulation in homogeneousand heterogeneous porous media.In this paper，we evaluate me capability of the method in the solutionof highly nonlinear and coupled partial diferential equations by simulating hydrocarbon reservoirs usingthe black-oi1 mode1.Furthermore.the efrect of grid orientation is investigated by simulating a benchmarkwate ooding problem.The numerical results show that the formulation presented here is efficient andaccurate for solving the bubble point and three.phase coning problems。

Key words：Control volume based finite element，black-oil model，grid orientation，porous medial IntrOducti0nPrediction of performance for primary and secondaryoil recovery processes has been one of the main concernsof reservoir engineers through the history of the petroleumindustry.With the advent of high speed computing，reservoirsimulators have proven to be invaluable tools to this end。

In this respect，various flow models are employed by thereservoir simulators.These models range from simplesingle.phase flow models(Aronofsky and Jenkins.1 954)tosophisticated multiphase.multicomponent compositionalflow models(Chen et a1.2006).Among these.the black-oilmodel is a standard three-phase flow model which is mostoften used by petroleum reservoir simulators.This is mainlybecause the black.oi1 model not only provides a reasonablygeneral representation of the multicomponent.multiphaseflow，but also avoids the necessity of using complicated phaseequilibrium models。

The black-oil models consist of a set of partial diferentialequations describing the conservation of mass for thewater，oil and gas components that generally coexist in aCorresponding author.email：sa ghoreishian###dena.kntu.ac.irReceived May 21，2012hydrocarbon reservoir.Darcys law and equations of state。

Since the equations are strongly nonlinear and coupled，theirnumerical solution is still a challenging task for reservoirengineers(Bergamaschi et al，1 998：Li et al，2003；2005：Naderan et al，2007；Lee et al，2008)even with the continualprogress made in both computational algorithms andcomputer hardware。

A reliable numerical solution should be able to take intoaccount the complexity of a real reservoir.The irregulargeological and geometrical morphology of hydrocarbonreservoirs affect the computational domain.The reservoirpermeability and porosity fields may experience very largelpeal variation up to 8 or 1 0 orders of magnitude(Durlofskyet aL 1 992)which results in highly discontinuous terms in thediscretized forln of the equations.This may lead the mode1to solve the equations inaccurately if the solution methodis not appropriate.M oreover，the strongly nonlinear natureof the equations can produce highly diffusive non.physicaloscillations at saturation fronts.The key solution for theseissues is to develop a conservative numerical scheme whichable to employ unstructured grids for spatial discretization。

Various numerical methods have been developed to modelmultiphase fluid flow through homogenous and heterogeneousporous media.The finite diference method(FDM1 is thePet.Sci.(2013)10：361-372 363the FVM to take the advantage of the fact that fluid velocitieswhich are discontinuous between tw o elements of the mesh。

are continuous through the control volume faces of the FVgrid(Durlofsky，1 993)。

The numerical method used by Forsyth f l 990).Fung et al(1992)，Gottardi and Dall Olip(1992)and Verlna f1996)inpetroleum reservoir simulation literature is called CVFEM .Inmost CVFEM employed in the field of reservoir simulation，the discrete equations of the multiphase flow system areobtained by integrating the equations of single-phase flowmodels.and then extending them by inoducing the mobilityterms.By considering this strategVthe method encounterssome grid orientation problems which restricts the meshangles to be equal or less than a right angle(Cordazzo et a1。

2004a；2004b).This restriction on mesh generation of thephysical domain can be di伍 cult to follow for most of thereservoirs due to their complex geometries.In addition.incommon control volume finite element approaches the porousmedium properties such as the absolute permeability and theporosity are stored in the center ofthe contro1 volumes(Verrna。

1 996).If this strategy is considered for heterogeneous porousmedia.inter-nodal permeability evaluation will be requiredsince the integration point lies on the interface of dif-ferentmaterials.The problems related to the inter-nodal permeabilityevaluation are investigated by Romeu and Noetinger r1 995)and Cordazzo et al(2003).E伍cient treatment of a numericalscheme when it faces with discontinues material properties isa very important issue in modeling heterogeneous reservoirssuch as fractured formations re.g.Nick and MatthN，20 ll：Eikemo et al，2009；Reichenberger et al，2006)。

The formulation presented here is derived directly fromthe multiphase flow equations in order to reduce the so-calledgrid orientation effects.Unlike the usual formulation ontriangular and tetrahedral elements.any attempt of adaptingthe discretized equations to the conventional form s of FVMis discarded.Consequently，the concept of transmissibility iscompletely abandoned in the form ulation。

In the first Paper of this series(Sadrneiad et al，2012)，wedeveloped a control volume based finite element method tosolve the governing equations of incompressible two.phasefluid flow in heterogeneous porous media.The capability ofthe method to handle discontinuous material properties andits efficiency for capturing saturation fronts with minimumnumerical dispersion and diffusion errors are evaluated byseveral numerical examples.In this paper，the proposedmethod is adopted for numerical solution of the black.oil fluidequations.This model iS able to consider the compressibilityand the mass transfer effects between the phases.Thenumerical results for the benchmark problems of the firstand second SPE comparative solution projects fOdeh，1981：Weinstein et a1.1 986)are presented and compared with thereported solutions to evaluate the stability and convergence ofthe form ulation.Moreovethe effects Of grids orientation areinvestigated by a benchmark waterflooding problem。

2 Governing equationsIn this section，the basic equations describing the black-oil model for reservoir simulation are derived based onthe classical continuum theory of mixtures fGoodman andGowin，1 972).The reservoir fluid is considered as a mixtureof water，oil and gas phases.It is assumed that the only massexchange occurs between oil and gas phases.and no masstransfer occurs between water and the other two phases。

Furthermore，tw o distinct zones are considered in the porousmedium.which are a dominant water-oil zone and a dominantoil.gas zone.The system in the water-oil zone is consideredto be water-wet.while in the oil.gas zone is oil-wet。

The mass conservation equations are presented asv.( )：Mw Otfor the water phase，掣v( )-m，ogtJtforthe oilphase，andV.(ng，OgWg)：-砌。g(1)(2)(3)for the gas phase，where月 ，P ，w ，and represent thevolume fraction，the mass density，the relative velocity andthe source/sink terms of phase c[c：water(w)and oil(o)and gas(g)]，respectively；and m。g stands for the exchangeof mass between the oil and gas phases.It is worth notingthe source/sink term s( )are calculated based on the welmodel presented by Wan(2002)。

The water phase densitv is determ ined by等(1Cw(Pw-Pi) (4)where Pwis the density of water at standard conditions；Pis the water phase pressure；Bwl is the water form ation volumefactor at the initial formation pressure ofpi；and Cw representsthe water compressibility factor.The gas phase density iscalculated by；where Pg is the density of gas at standard conditions withpressure p。and temperature ；B is the gas formation volumefactor；Z is the deviation factor；T is the reservoir temperature；and Pg is the gas phase pressure.The oil phase density shouldbe determined with respect to the fact that it consists of oiland gas components-Rsop PoPo - - (6)where p0s represents the density of the oil phase at the standardconditions；R 。is the gas solubility in the oil phase；and B。isthe oi1 formation volume factorThe mass exchange term(唬g)in the mass balanceequations，could be calculated as鲁 no (7)The relative velocity of each phase(regarding the solidphase)could be calculated by Darcys law-饥-364 Pet.Sci.(2013)10：361-372-kr，o，x[p，g- n。 (8) 丢(dpg-dpo)where K is the absolute permeability tensor； and are and from Eqs.(1 5)and(9)we havethe relative perm eability and the dynamic viscosity ofphase 0c，respectively. dn。 dn1-dIn addition，there is a constraint for the fluid volumefractions dngdn-dnlw g (91 substituting Eq. (1 6)into Eq.(1 9)，one obtainswhere F/represents the rock porosity，and it is assumed tohave the following form for the slightly compressible systemsn"i(1CR(P-Pi)) (10)where Hi state the porosity at the reference pressure i)；CRrepresents the rock compressibility；and P is the volumeaveraged pore pressure which has been defined by Pap et al(2001)PSwPwSop0Sgpgwhere S represents the saturation of phase .Replacing thephase saturations by the phase volume fractions，Eq.(1 1)canbewriten aspl (nwpw ) (12)In order to specify the interacting motion of each phaseon the other phases，constitutive equations are required whichcan link the fluid phase pressures to their volume fractionsAccording to Hassanizadeh and Gray(1 993)，the mostpractical method for considering this interacting motion iSto use empirical correlations relating the capillary pressure( )to the phase volume fractions.The capillary pressure iSdefined as the pressure difference of two immiscible fluidsacross their interface.For the water.oi1 zone.capillary(18)(19)(20)d，z0. (呶 -dpo)- (dp。-dpw) (21) 印 。

for the evolution of oil volume fraction.Substituting Eqs. (1 7)and(2 1)into Eq.(1 6)，and then substituting the result intoEq.(20)，one obtains the folowing relation for the evolutionof gas volume fractiond ：- ×n-niCR(1P )- c n'wPow) (州c -(22)The final form of the flow equations can be simplydeveloped by substituting Eqs.(4)·(8)，(17)，(21)and(22)into Eqs.(1)-(3)as below("pw)O p o PwKw(Pwg-Vp )-Mw0pressure is indicated by p～(P。-P )，and for the oil. for the water phase，in whichgas zone，capillary pressureFollowing this，we can writen F(p 。 ) F( )where r/1 representsVOlume fraction：nw nois indicated by Pcg。(Pg-P。)： 0w (13)) the total liquid phase(water and oil) similarly，for the oil phase(15)Partially differentiating Eqs.(10)，(13)and(14)，oneobtainsd 丽nag ( w 。

(16)F/gdpg dF/wPodnoPg dng)dn - wopc0(dp。-dP ) (17))划 ( 。 - ，-c 警 ( )VT Epo，o(；og- )]- 0inwhich(23)(24)(25)上 (26)(27)塾堕Pet.Sci.(2013)10：361-372 365，- n S0 - Cp。

勰 cCpoK。：-kro-K/toand finaly，for the gas phase"iQ- CR(1Pg) Hwn警 ×。- 警V IpgKg(&g- )] 0inwhich(28)(29)(30)- Po)l/(31)(32)(33)Eqs.(23)，(26)and(3 1)represent a system of highlynonlinear and coupled equations describing the black。

oi1 now in a hydrocarbon reservoir.The major sources ofnOn1inearities in these equations.i.e.the phase volumefraction(，2 )，and relative permeability( )are stronglydependent on the primary unknown variables and thereforeshould be continuously updated during the solution procedure。

However，effects of the weak nonlinearities.i.e.the formationvolume factor( )，viscosity( )，gas solubility( 。)，andporosity 1 could not be neglected.Treatment of these twotypes of nonlinearities is carefully described by Setari andAziz f1 9751.In order to complete the description of thegoverning equations，it is necessary to denne appropriateinitial and boundary conditions.The initial conditions specifythe fu1 field of phase pressures at time tOpO in Q and on F (34)where Q is the domain of interest and F is its boundary.Theboundary condition can be of two types or a combinationof these：the Dirichlet boundary condition，in which thephase pressures on the boundaries(rp)are known，and theNeumann boundary condition，in which the values of phasefluxes at the boundaries(I1q)are imposed(where FFp U rq)p on Fpl P ( g- 。)I.，l 广 1 J on Fq(35)(36)where n denotes the outward unit normal vector on theboundary and孽 is the imposed mass flux which is normal tothe boundary。

3 Numerical solutionDiscretization of the governing equations can be nowexpressed by the use of the CVFEM developed by Sadrneiadet al(20 1 2)in terms of the nodal phase pressures(i.e. )which are selected as the primary variables.In this method，the physical domain is discretized using hexahedral elements(Fig.1(a))，and further subdivision of the elements intocontrol volumes is performed in the transform ed space(Fig。

1(b11.Furthermore，in order to represent the discretizedequations at the element level(instead of the control volumeleve1).the control volumes are also divided into sub.controlvolumes.Each of these sub.control volumes belongs to aspecific element which is in association with the given controlvolume(Fig.1(c))。

Fig.1 (a)System of finite element mesh in the physical space；(b)representation of control volume around a node in the transformed space；(e)ilustrating a sub-control volume belongs to the eliminated element(afterSadmad et al，2012)The values of phase pressures(p )at any point within anelement are approximated by the folowing expressionP ( ，7， )N(4，7，Ob (37)where N(4，7， )is the vector of standard finite elementshape functions。

CVFE discretization procedure，as presented by Sadmejadet al(20 1 2)，when applied to Eqs.(23)，(26)and(3 1)alongwith the boundary condition(36)，they yield366 Pet.Sci.r2013)10：361-372只 0 l0 。

l 0 0 LWOjd FW1。 l J l厂g J -IFS.C.V-Fq ( N)T.n) dI1。 -IFS，C.V，-Fq ( N)T./， dF -IFs.c.v.-Fq ( ) ·即)dI1 ， 、l 。 l 。 jⅣdQfn。 - - e" pgs]ⅣdQ。

。 - 等]ⅣdQCw。 ( 'w，Ow)NdfCo ( "Po)NdnCo (，z[Po)Nd I ( ，z p)jⅣdQ。 等 丽niCg 叫c 刊]ⅣdQ 。MwW dQ-fr dr-fr。 [ ( g) . ] dr dQ-j.r dr-j.r [ ( 。g) . ] dr dQ-r dr-r 。[ ( g T.hiW dr(38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48)(49)(50)(51)(52)where Qs c v and Fs c v denote the domain of a sub.control volume and its faces，respectively；and W is the vector of thewei tini nctions.The weighting functions are chosen such that the ith weighting function of an element takes a constantva1ue of unity over the sub.control volume belonging to node i，and zero elsewhere in the element，1·e·- l inthe sub。controlW Volumebe1ongs o node ⅢkJ.-，) ， l 0 otherwiseThe tempora1 discretization of Eq.(38)is performed by the fuly implicit first order accurate finite diference che-..................L 1， J 儿 0 0 -Pet.Sci.(2013)10：361-372 367l △ Cw。

乳0 Jt l Iwhere At is the time step length and 1p 1-p 。

Eq.(54)represents a very large system of coupled andhighly nonlinear equations onAp )f1 which should besolved by an appropriate iterative method.In the presentwork.the globa1 inexact afine invariant Newton techniquerGIANT1 developed by Nowak and Weimann f 1 990)isemployed to solve the system of Eq.f54).In this method，the global inexact Newton technique(Deuflhard.1 990)iscombined with the fast secant method rDeuflhard et al，1 990)，as an iterative linear solver，to obtain an efficient and robustnumerical solution of a very large scale highly nonlinearsystem of equations.An algorithmic overview of the varipusparts ofthe calculations is given as a pseudo code in 1'ab1e 2。

Table 2 Algorithmic outline of the various parts of calculationsin the present modelPre-processing and initializingD0 F0R EACH time stepUPDATE variables and parametersDo FoR EACH Newton iterationD0 FoR EACH elementUPDATE the Jacobian matrix ofEq.(54)ASSEM BLE to the global Jacobian matrixEND D0D0FoR EACH elementUPDATE the residual of Eq.(54)ASSEM BLE ofthe global residua1 vectorEND D0DoLinear solver with fast secant methodIF(Linear solver Conv.TRuE.、ExITEND D0IF(Newton Conv.TRUE.)THENEXITELSEUPDATE variablesEND IFEND D0Post-processingEND D0END4 Numerical experimentsThe benchmark problems of the first and secondcomparative solution projects(CSP)of the SPE are usedto evaluate the validity of the presented formulation andstability and convergence of the method to deal with abubble-point and a three-phase coning problem.Furthermore，grid orientation effects on the results obtained by the modelare investigated by a benchmark waterflooding example。

4.1 Gas displacementThis simulation problem is adopted from the first caseof the benchmark problem of the first CSP f0deh，1 98 1 1。

This benchmark problem is a challenging case.and it wasdesigned to evaluate the stability of black.oil reservoirsimulators to deal with strong nonlinearity of the govemingequations and transition of the reservoir condition from anundersaturated to a saturated state.The state of a reservoiris caled undersaturated when it initially exists at a pressurehigher than its bubble.point pressure.At the undersaturatedcondition，oil and water are the only fluid phases present inthe reservoir.W hereas，at the saturated state of the reservoir，the free gas exists and its relatively high compressibilityproduces a strong source of nonlinearity in the governingequations.It is worth noting that the state of a location in areservoir can change from saturated to undersaturated state，or vice versa，during the solution.The complete details of theproblem can be found in Odeh f19811。

The results obtained by 7 organizations which participatedin the solution project have been reported by Odeh r1 98 l 1.All description of the simulators used by these participantscan be found in the reference.Generally，all the models weredeveloped based on the traditional finite diference or finitevolume methods.In this paper，the results obtained by thepresent model are compared with those obtained from twocompanies，namely，Shell Development Co.and IntercompResource Development and Engineering Inc.。

Figs.2 and 3 show the oil production rate and the gas-oil ratio obtained by the present model and those reportedby Odeh(1981).As seen in the figure，the models showsimilar trends in the results.however a slight discrepancy isobserved during 2.5 to 6.0 years.The maximum ratio ofthesedi日erences is about l 0% which is observed in the calculatedoil production rates at 2.95 years。

The averaged pressure values in the blocks containingproduction and injection wells are compared with the resultsobtained by Shell Development Co.and Intercomp Resource2 3 4 5 6 7 8 9Time，YearFig.2 Oil production rate versus time1O加 倡 侣 他 们 8 6 4 2 0凸、∞J∞乏 000 1..m .1 co-l。了 0J 01 O variation ofgas.。il rati。versuS嘛 C ：- J -370 Pet.Sci.(2013)10：361-372elements is in good agreement with the former.Regardingto these results，one could conclude that sensitivity of themethod to the grid orientation is very low.These figures alsodemonstrate the convergence of numerical solution scheme。

1 00 80 60 40 200 0 5 1 0 1 5 2 0Injection volume，PVFig-12 Maximum and minimum ofwater cut at the production welsSolid Iine：307 elementsDash Iine：1 96 elementsFig-1 3 Water saturation contour after 2 PV ofwater injection250000 0 2 0 4 0.6 0 8 1 0 1 2 1 4 1 6 1 8 2 0Injection volume，PVFig.14 Variation ofpressure at the injection blockFig.1 5 shows the distribution of local mass balance errorfor the present test case.As it iS seen in the figure. for theinjection rate of 0.8 m /day，the local mass balance eror isin the range of 1 0 m /day.Fig.1 6 describes the cumulativeglobal mass balance error over injection time.Dividing theerror by the totaflow rate Of 0.8 m /day，we get less than1 0 % erors.These figures show that the present model iS afully conservative scheme。

-1 E-07 5E-08 2E.07Fig.1 5 Distribution of local mass balance eror2 3 4 5 6 7 8 9Time.YearFig.1 6 Cumulative global mass balance ell'or over injection time5 Conclusions10In this paper.a fully coupled control volume finiteelement model iS developed to simulate the black.oi1 flOWin hydrocarbon reservoirs.The method iS fully conservativeand able to deal with unstructured grid systems.It combinesthe mesh flexibility of FEM with the loca1 conservativecharacteristic Of FDM and consequently overcomes thedeficiencies Of FDM in dealing with geometrically complexreservoirs，or inability of FEM to conserve mass locally。

Loca1 and global conservative characteristic of the methodlo Ja ∞主-o/0.co BJru ∞∞∞O∞ ∞ ∞ ∞加 侣 ∞ 5.aJ了∞∞ .1372 Pet.Sci.(2013)10：361-372Lu Q.Malgorzata P and Gai X.Implicit black-oil mode1 in IPARSframework IPARSv2，ModelI，keyword BLAC .Version 1.1.TexasInstitute for Computational and Applied Mathematics，The UniversityofTexas atAustin.2001Mathai S K.Geiger S，Roberts S G，et a1.Numerical simulation ofmulti-phase fluid flow in structurally complex reservoirs.GeologicalSociety，London.Specia1 Publications.2007.292：405.429Monteagudo J and Firoozabadi A.Control-volume model for simulationof water injection in fractured media：incorporating matrixheterogeneity and reservoir wetability efects.SPE J.2007.12：355-366Naderan H.Manzari M T and Hannani S K.Application and performancecomparison of high-resolution central schemes for black oil mode1。

Int.J.Num.Meth.for Heat&Fluid Flow.2007.1 7f7)：736-753Nick H M and Matthi S K.Comparison of three FE-FV numericalschemes for single--and two.phase flow simulation of fracturedporous media.Transp.Porous Med.20l1.90：421-444Nowak U and Weimann L.GIANT：A software package for thenumerical solution of very large systems of highly nonlinearequations.Konrra.Zuse-Zentrum fuer Informationstechnik Berlin。

1990Odeh A.Comparison of solutions to a three.dimensional black.oilreservoir problem.J.Pet.11ech.1 98 1.1 3.25Pao W K.Lewis R W and Masters I.A fully coupled hydro-thermo。

mechanical model for black oil reservoir simulation.Int.J.Numer。

Ana1.Meth.Geomech.2001.25：l229.1256Pruess K，Oldenburg C and Moridis G.TOUGH2 users guide，version2.Earth Science Division.Lawrence Berkeley National Laboratory。

2012Reichenberger、Jakobs H，Bastian P.et a1.A mixed.dimensional finitevolume method for two.phase flow in fractured porous media.Adv。

Water Resour.2006.29：1020.1036Romeu R and Noetinger B.Calculation of internodal transmissibilitiesin finite diference models of flow in heterogeneous porous media。

water Resour.Res.1 995.3 l f4)：943-959Sadmejad S A.Ghasemzadeh H Ghoreishian Amiri S A.et a1.A coiltrolvolume based finite element method for simulating incompressibletwo-phase flow in heterogeneous porous media and its application toreservoir engineering.Petro1.Sci.2012.9f4 ：485-495Setari A and Aziz K.Treatment of nonlinear terms in the numericalsolution ofpartial differential equations for multiphase flow in porousmedia.Int.J.Multiphase Flow.1975.1：817.844Simbnek J.Huang K.and Van Genuchten M Th.The SWMS.3D codefor simulating water flow and solute transport in three-dimensionalvariably.saturated media.U.S.Salinity Laboratory，Agricultura1Research Service，U.S.Department of Agriculture，Riverside，Califomia.1 995Verma S.Flexible grids for reservoir simulation.Ph.D.Thesis。

Department ofPetroleum Engineering.University ofStanford.1996肋 n J Stabilized finite element methods for coupled geomechanics andmultiphase flow.Ph.D.Thesis.Stanford University.2002W lang W ，Rutqvist J，G6rke U J，et a1.Non-isotherm a1 flow in lowperm eable porous media：a comparison Of Richardsan d two-phaseflow approaches.Environ.Earth.Sci.20 l1.f621：l1 97.1 207W instein H.Chappelear J and Nolen J.Second comparative solutionproject：A付1ree-phase coning study.J.Pet.Tech.1986.345.353White M D and Oostrom M.STOMP：Subsurface Transport OverMultiple Phases，version 2，theory guide.U.S.Department ofEnergy，2000Young L.An efficient finite element method for reservoir simulation。

Proceedings of the 53rd SPE Annual Technica1 Conference andExhibition，Houston。1978(paper SPE 74131(Edited by Sun Yanhua)