Multiscale extended finite element method for deformable fractured porous media

DOI

10.1016/j.jcp.2021.110287

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Xu, F., Hajibeygi, H., & Sluys, L. J. (2021). Multiscale extended finite element method for deformable

fractured porous media. Journal of Computational Physics, 436, 1-19. [110287].

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

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.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Multiscale

extended

finite

element

method

for

deformable

fractured

porous

media

Fanxiang Xu

a

,∗

,

Hadi Hajibeygi

b

,

Lambertus

J. Sluys

a

a_{Materials, Mechanics, Management and Design, Faculty of Civil Engineering and Geosciences, Delft University of Technology, P.O. Box 5048,} 2600 GA Delft, the Netherlands

b_{Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences, Delft University of Technology, P.O. Box 5048, 2600} GA Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Article history:

Availableonline19March2021 Keywords:

Fracturedporousmedia Extendedfiniteelement Multiscale

Geomechanics Scalableiterativesolver

Deformable fractured porous media appearin many geoscienceapplications. Whilethe extended finite element method (XFEM) has been successfully developed within the computationalmechanicscommunityforaccuratemodelingofdeformation,itsapplication in geoscientific applications is not straightforward. This is mainly due to the fact that subsurfaceformationsareheterogeneousandspanlargelengthscaleswithmanyfractures atdifferentscales.Toresolvethislimitation,inthiswork,weproposeanovelmultiscale formulation for XFEM, based on locally computed enriched basis functions. The local multiscale basis functions capture heterogeneity of th e porous rock properties, and discontinuities introduced bythe fractures. Inorder topreserveaccuracy ofthesebasis functions,reduced-dimensionalboundaryconditionsaresetaslocalizationcondition.Using these multiscale bases, a multiscale coarse-scale system is then governed algebraically and solved. The coarse scale system entails no enrichment due to the fractures. Such formulation allows for significant computational cost reduction, at the same time, it preservestheaccuracyofthediscretedisplacementvectorspace.Thecoarse-scalesolution is finally interpolated back to the fine scale system, using the same multiscale basis functions. The proposed multiscale XFEM (MS-XFEM) is also integrated within a two-stagealgebraiciterativesolver,throughwhicherrorreductiontoanydesiredlevelcanbe achieved.Severalproof-of-conceptnumericaltestsarepresentedtoassesstheperformance ofthedevelopedmethod.ItisshownthattheMS-XFEMisaccurate,whencomparedwith thefine-scalereferenceXFEMsolutions.Atthesametime,itissignificantlymoreefficient thanthe XFEMonfine-scaleresolution,as itsignificantlyreduces thesize ofthelinear systems. Assuch, it developsapromising scalable XFEMmethodforlarge-scale heavily fracturedporousmedia.

©2021TheAuthor(s).PublishedbyElsevierInc.Thisisanopenaccessarticleunderthe CCBYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Subsurfacegeologicalformationsareoftenhighlyheterogeneousandheavily fracturedatmultiplescales.Heterogeneity ofthedeformationproperties(e.g.elasticitycoefficients)canbeofseveralorders ofmagnitudewhichoccursatfinescale (cm)resolution.The reservoirs alsospanlarge scales,in theorderofkilometers. Numerical simulationofmechanical

de-*

Correspondingauthor.

E-mail addresses:f.xu-4@tudelft.nl(F. Xu),H.Hajibeygi@tudelft.nl(H. Hajibeygi),L.J.Sluijs@tudelft.nl(L.J. Sluys).

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

0021-9991/©2021TheAuthor(s).PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

hand, the immersedor embedded approach allows for independent grids formatrix block and fractures, by introducing enrichmentofthediscreteconnectivity(forflow)andshapefunctions(formechanics) [8–14].Theseenrichedformulations areaimedatrepresentingdiscontinuitieswithintheoverlappingmatrixcell,withoutanyadjustmentnorrefinementofthe grid[15].Theenrichmentstrategiesformodelingdeformationoffracturedmedia,usingfinite-elementschemes,arereferred toas‘extendedfiniteelement(XFEM)’methods.

XFEM enriches the partition of unity (PoU)[16] byintroducing additional degreesof freedom (DOFs) atthe existing element nodes.Thereexists setsofenriched functionstocapturethejump discontinuityinthedisplacement field,when thefracture elementcutsthroughtheentirecell,andthetipwhena fractureendswithin thedomainofanelement (i.e. itstipisinsideanelement) [17–22].Forthesejumpandtipscenarios,additionalshapefunctionsareintroducedwhichare multipliedbytheoriginalshapefunctionsandsupplementthediscretedisplacementapproximationspace.

Presenceofhighly heterogeneous coefficientswithhighresolutionwithin large-scale domainshas beensystematically

addressed in the computational geoscience community through the development ofmultiscale finite element andfinite

volume methods [23–26]. Recent developments also include mechanical deformation coupled with fluid pore pressure dynamics [27–31], in which localmultiscale basis functions were computed using either linear or reduced dimensional boundaryconditions.Inpresenceoffractures,speciallymanyfractures,however,complexityofthecomputationalmodel in-creasessignificantly.Assuch,developmentofarobustmultiscalestrategyfordeformationofheavilyfracturedporousmedia, whichalsoallowsforconvergentsystematicerrorreduction[32–34],isofhighinterestinthegeosciencecommunity.

Whenitcomestogeoscienceapplications,theclassicalXFEMisnotanattractivemethod,duetoitsexcessiveadditional degreesoffreedomtocapturethemanyfractures.

Whenseparationofscalescanbeappliedtotherelevantproperties,upscalingorhomogenizationapproachesarefound effective[35–39].Thisalsoincludestherich studiesrelevanttodeformation ofnonfractured[40–42] andfracturedmedia [43].Notethatgeologicalheterogeneitiesandtheirmultiscalefracturesdonotentailseparationofscales[44].

This paper develops a multiscale XFEM (referred to asMS-XFEM) which offers a scalableefficient strategy to model large-scalefractured systems.MS-XFEMimposesacoarse meshonthegivenfine-scalemesh.Themainnovelideabehind MS-XFEM istouseXFEMtocomputationallysolveforlocalcoarse-scale(multiscale)basisfunctions. Theseenrichedbasis functionscapturethefracturesandcoefficientheterogeneitywithineachcoarseelement.Thesolvingstrategyoftheselocal coarse-scale basisfunctionscan beeithergeometric oralgebraic [33,45,46]. Weprefer algebraic construction,asitallows forblack-boxintegrationofthemethodwithinanyexistingXFEMsimulator.Oncethebasisfunctionsaresolved,theywill beclusteredintheprolongationmatrix(P),whichmapsthecoarse-scalesolutiontothefine-scaleone.Notethattherewill benoadditionalmultiscalebasisfunctionsduetojumportips,andthatonly4multiscalebasisfunctionsperelementexist for2Dstructuredgrids(8in3D)ineachdirection(x,y,andz).

ThefinescaleXFEMsystemisthenmappedtothecoarsegridbyusingtherestriction(R)matrix,whichisdefinedbased on theFEM,asthe transpose ofP. The approximatefine-scalesolution isfinally obtainedaftermapping thecoarse-scale solutiontothefinescalegrid,byusingtheprolongationoperator.

TheapproximatesolutionofMS-XFEMisfoundacceptableformanyapplications,howevererrorcontrolandreductionto anydesiredlevelisnecessarytopreserveitsapplicabilityforchallenging cases.Assuch,theMS-XFEMisintegratedwithin thetwo-stageiterativesolverinwhichtheMS-XFEMispairedwithan efficientiterativesmoother(hereILU(0))toreduce the error[47,48]. One can alsouse theKrylov subspacemethods (e.g. GMRES) to enhance theconvergence, whichstays outsidethescopeofthispaper.

Severalproof-of-concept numericaltestsare presentedto assesstheaccuracy ofthepresented MS-XFEMwithout and withiterativeimprovements.Thetestcasesincludelargedeformationswhichmaynotberealisticingeoscienceapplications, butimportanttobe studiedinordertoquantifythe errorsinlargedeformation scenarios.From theseresultsitbecomes clearthattheMS-XFEM,despiteusingnoenrichedbasisfunctionsatcoarsescale,presentsanefficientandaccurate formu-lationtostudydeformationoffracturedgeologicalmedia.

The structureofthispaperissetasthefollowing.Next,thegoverningequationsandthefinescale XFEMmethodare introduced.Then,theMS-XFEMmethodispresentedindetail,withemphasisontheconstructionoflocalmultiscalebasis function andthe approximate fine scale solution.Then, different numericaltest cases are presented. Finally, concluding remarksarediscussed.

Fig. 1. An illustration of fractured domain setup.

2. Governingequationsandfine-scaleXFEMsystem

Considerthedomain



boundedby



asshowninFig.1.PrescribeddisplacementsorDirichletboundaryconditionare imposedon



u,whiletractionsareimposedon



t.Thecracksurface



c (linesin2-Dandsurfacesin3-D)isassumedtobe traction-free.

Themomentumbalanceequationsandboundaryconditionsread

∇ ·

σ

+

f

=

0 in



(1)

σ

· −

→

n

= ¯

t on



t (2)

σ

· −

→

n

=

0 on



c (3)

u

= ¯

u on



u (4)

where

σ

is the stress tensor andu is the displacement field over the whole domain.

−

→

n is the normal vector pointing outsidethedomain[49,50].

Theconstitutivelawwithlinearelasticityassumptionreads

σ

=

C

:

ε

=

C

: ∇

su

,

(5)

where,

∇

s_{denotesthesymmetricaloperatorandC isthepropertytensordefinedas}

C

=

⎡

⎣

λ

+

λ

2

μ

λ

+

λ

2

μ

00 0 0

μ

⎤

⎦ ,

with

λ

and

μ

denotingtheLame’sparameters[5,51]. Thestraintensor

ε

isexpressedas

ε

= ∇

s_u

₌

1

2

(

∇

u

+ ∇

T_u

_),

₍₆₎

where,

∇

denotesthegradientoperator.

SubstitutingEqs. (5) and(6) inthegoverningEq.(1) resultsina2nd orderPartialDifferentialEquation(PDE) for dis-placementfieldu

∇ · (

C

: ∇

su

)

+

f

=

0

.

(7)

Eq. (7) is then solved for computational domains withcracks (representing faults andfractures). Thisis done by the extendedfiniteelement(XFEM)method,whichisbrieflyrevisitedinthenextsection.

2.1. ExtendedFiniteElementMethod(XFEM)

FEMemployssmoothshapefunctionsNiforeachnodei

∈

I,whereI is thesetofallnodes,toprovideanapproximate numericalsolutiontoEq.(7) fordisplacementunknowns,i.e.,

u

=



i∈I

Fig. 2. Illustration of polar coordinates(r, θ )towards the fracture, for tip enrichments.

Fig. 3. Enrichment mechanism: node I and J will be enriched using tip and jump functions.

This smooth FEM approximation is insufficient to capture discontinuities imposed by the existence of the fractures and faults. As such, the XFEM method introduces two sets of enrichment to the original FEM to capture the discontinuities without adapting the grid. Theseenrichment sets are associated withthe body andtip ofthe fracturesand faults. The body is enriched by jumpfunctions, andthe tip by tipenrichment functions[17]. Below briefdescriptions ofthesetwo enrichmentfunctionsareprovided.

2.1.1. Jumpenrichment

Thejumpenrichmentrepresentsthediscontinuityinvolvedinthedisplacementfield acrossthefractureandfaultmain body.ThejumpenrichmentisoftenchosenasthesteporHeavisidefunction,whichcanbeexpressedas

H

(

x

)

=



1 on



+

−

1 on



−

.

Notethat



+and



−zonesaredeterminedbasedonthenormalvectorpointingoutofthefracturecurve.Forlinefractures, thedirectioncanbeanyside,aslongasalldiscreteelementsusethesame

+

and

−

sidesforafracture.

2.1.2. Tipenrichment

Thetipenrichmentrepresentsthediscontinuityofthedisplacementfieldnearthefracturetip.Thistypeofenrichment function,denotedby Fl,isbasedontheauxiliarydisplacementfieldnearthefracturetipandcontainsfourfunctions,i.e.,

Fl

(

r

, θ )

= {

√

rsin

(

θ

2

),

√

rcos

(

θ

2

),

√

rsin

(

θ

2

)

sin

(θ ),

√

rcos

(

θ

2

)

sin

(

θ

2

)

}.

Here,

θ

rangesfrom

[−

π

,

π

]

,whichisdefinedbasedonthelocalpolarcoordinatesofapointwithrespecttothefracture tip.AnillustrationisgiveninFig.2.

2.1.3. Enrichmentmechanism

Todecidewhetherthenodeisenrichedornot,thenodelocationrelatedtothefractureisthekeyfactor.Thesketchof theenrichmentmechanismisshowninFig.3.Moreprecisely,inthisfigure,thetipandjumpenrichednodesarehighlighted incircularpointandsquarepoint,respectively.

2.2. XFEMlinearsystem

TheXFEMapproximatesthedisplacementfieldu atfine-scaleresolutionh byuh whichisdefinedas

u

≈

uh

=



i∈h uiNi

+



j∈J ajNjH

(

x

)

+



k∈K Nk

,



_

4 l=1 Fl

(

x

)

blk

,

(9)

whereN, H andFl represent,respectively,theclassicalFEMshapefunctions,theHeavisidefunctionandthetipenrichment functions. The fine-scale mesh has



h nodes. Moreover, u denotes the standard DOFs associated to the classical finite

element method.Furthermore, a denotes the extra DOFsassociated to the jump enriched node. For the jump enriched

nodes,inthe2Ddomains,eachnodewouldcontain2extraDOFs.Similarly,b indicatestheextraDOFsassociatedtothetip enrichment,whichaddsfourextraDOFsperdirection(8intotalina2Ddomain)foreachtipinsideanelement.

Thefirsttermintheright-hand-side(RHS)ofEq.(9) isthecontributionoftheclassicalfiniteelementmethod.Thisterm captures thesmooth deformation, usingclassical shape functions. Thesecond term, however, representsthe contribution of the jumpenrichment. Note that the jumpenrichment is modeled by the weighted Heaviside functions, with weights beingtheclassicalshapefunctions.Therewillbeasmanyjumpenrichmentfunctionsasthenumberoffracturesinsidean element.Finally,thethirdtermintheRHSisthecontributionofthefracturetips.Notethatifseveralfracturetipsendup inanelement,therewillbe4additionalDOFspertipperdirectioninthatelement.

Theresultinglinearsystementailsthenodaldisplacementunknownsu,aswellasthejumplevela andtipweightb per

fracture(andfault).TheaugmentedXFEMlinearsystemKh_dh

₌

_fh_{,therefore,reads}

⎡

⎣

KKuuau KKuaaa KKubab Kbu Kba Kbb

⎤

⎦



Kh

⎡

⎣

ua b

⎤

⎦



dh

=

⎡

⎣

ffu_a f_b

⎤

⎦



fh

.

(10)

ComparedtotheclassicalFEM,thereexistseveraladditionalblocksinvolvedinthestiffnessmatrix,duetotheexistenceof thediscontinuities.TheadvantageofXFEMisthatitdoesnot relyoncomplexmeshgeometry,instead,itallowsfractures tooverlapwiththematrixelements.Ontheotherhand,forgeoscientificfracturedsystems,theadditionalDOFsduetothe enrichmentprocedureresultsinexcessivecomputationalcosts.ThisimposesasignificantchallengefortheXFEMapplication ingeoscienceapplications.Inthispaper,wedevelopascalablemultiscaleprocedurewhichconstructsacoarse-scalesystem basedonlocallysupportedbasisfunctions.Themethodisdescribedinthenextsection.

3. MultiscaleExtendedFiniteElementMethod(MS-XFEM)

Amultiscaleformulationprovidesanapproximatesolutionuh tothefine-scaleXFEMdeformationuh through uh

≈

uh

=



i∈H

NH_i u_iH

,

(11)

whereN_iH arethecoarse-scale(multiscale)basisfunctionsanduH_i arethecoarse-scalenodaldisplacementsatcoarsemesh



H. Note that thismultiscale formulationdoes not include any enrichment functions. Instead, all enrichment functions are incorporatedintheconstructionofaccuratecoarse-scale basisfunctions NH.Thisallowsforsignificantcomputational complexityreduction,andmakestheentireformulationattractiveforfield-scalegeoscientificapplications.

Next,constructionofthecoarse-scalesystemandthebasisfunctionswillbepresented.

3.1. Coarsescalelinearsystem

MS-XFEM solves the lineardeformation systemon a coarse mesh,imposed on a givenfine-scale mesh, asshown in

Fig.4.Thecoarseningratioisdefinedastheratiobetweenthecoarsemeshsizeandfine-scalemeshsize. Themultiscaleformula(11) canbealgebraicallyexpressedas

uh

≈

uh

=

P dH

,

(12)

whereP isthematrixofbasisfunctions(i.e.,prolongationoperator)anddH isthecoarse-scalealgebraicdeformationvector, correspondingtothephysicalnode-baseduH _{unknowns.Algebraicformulationallowsforconvenientimplementationofthe} proposedMS-XFEM,anditsintegrationasablack-boxtoolforanygivenclassicalXFEMsolver.Therefore,theremainderof thearticlewillbedevotedtotheformulation.

The coarse-scalesolutiondH _{needs tobefound bysolving acoarse-scale system.Toconstructthecoarse-scale system} andsolvefordH_{,onehastorestrict(map)thefine-scalelinearsystem(K}h_dh

₌

_fh_{)tothecoarse-scale,i.e.,}

Fig. 4. Illustration of the multiscale mesh imposed on the given fine-scale mesh, with a coarsening ratio of 3×3.

(

R KhP

)



KH

dH

=

R fh

.

(13)

Here,R istherestrictionoperatorwiththesizeof



H

× 

h+j+t,where



h+j+t isthesizeofthefine-scaleenrichedXFEM systemincludingjumpandtipenrichment.ProlongationoperatorP hasthedimensionof



h+j+t

× 

H.Thisresultsinthe coarse-scalesystemmatrixKH sizeof



H

× 

H.

Thefinite-element-basedrestrictionfunctionisintroducedasthetransposeoftheprolongationmatrix,i.e.,

R

=

PT

.

(14)

Therefore, thecoarse-scale matrix KH _{issymmetric-positive-definite (SPD), if K}h _{is SPD. Oncethe coarse-scalesystem is} solvedon



H spacefordH,onecanfindtheapproximatefine-scalesolutionusingEq.(12).Overall,themultiscaleprocedure canbesummarized asfindinganapproximatesolutiondh accordingto

dh

≈

dh

=

PdH

=

P

(

R KhP

)

−1R fh

.

(15)

Next,theprolongationoperatorP,i.e.,thebasisfunctionsareexplainedindetail.OnceP isknown,alltermsinEq.(15) aredefined.

3.2. Constructionofmultiscalebasisfunctions

Toobtainthebasisfunctions,thegoverningEq.(7) withoutanysourcetermusingXFEM,i.e.,Eq.(10) needstobesolved ineachcoarseelement

ω

H_{.Thiscanbeexpressedassolving}

∇ · (

C

: (∇

SN_iH

))

=

0 in



H

,

(16)

subjecttolocalboundaryconditions.Here,wedevelopareduced-dimensionalequilibriumequationtosolveforthe bound-arycells[25,52],i.e.,

∇

· (Cr

: (∇

SNiH

))

=

0 on



H

.

(17)

Here,



H denotestheboundarycellsofthecoarseelement



H.Inaddition,

∇

S denotesthereduceddimensionaldivergence andsymmetricalgradientoperators,whichact paralleltothedirectionofthelocaldomainboundary.Moreover, Cr isthe reduced-dimensional(here,1D)averageelasticitytensoralongtheboundaryofthecoarsecell.Moreprecisely, for bound-aries parallelto x-direction,theonlynonzero entriesare Cx_r

(

1

,

1

)

= ¯λ

andC_rx

(

3

,

3

)

= ¯

μ

;whileits nonzeroentries forthe boundariesparalleltoy-directionareCry

(

2

,

2

)

= ¯λ

andCry

(

3

,

3

)

= ¯

μ

.Here,the

¯λ

and

μ

¯

areaveragedelasticitypropertiesof theadjacent2Dcellsbelongingtotheboundary1Dedges.For2Dgeometries,thereduced-dimensionalboundarycondition representsthe1D(rod)deformationmodelalongthecoarseelementedges.Thelocaldeformationinresponsetohorizontal andverticalunitynodaldisplacementswillbeingeneraltwo-dimensional.Therefore,theprolongationmatrixP willhave off-block-diagonalentries(i.e., Pxy

=

0 and Pyx

=

0)andreads

Fig. 5. Illustration of the multiscale local basis functions, constructed using XFEM for the node H in x direction (a) and y direction (b).

Fig. 6. Illustrationofabasisfunctionthatcaptures thediscontinuityofafracture.Modelsizeis10×10 m2_{.Yellowsegmentrepresentsthefracture} extendingfrom(3,3)to(7,3)coordinates.Theshownbasisfunctionbelongstothecoarsenodelocatedatthecoordinate(4,4).(Forinterpretationofthe colorsinthefigure(s),thereaderisreferredtothewebversionofthisarticle.)

P

=



Pxx Pxy Pyx Py y



.

(18)

Fig.5showsanexampleofalocalsystemtobesolvedforbasisfunctionsbelongingtothehighlightednode H inx and y directions.NotethattheDirichletvalueof1issetatH foreachdirectionalbasisfunctions,whileallother3coarsemesh nodesare setto0.Furthermore,asshowninFig.5,theboundaryproblemissolvedforbothedges whichhavethenode

H atoneoftheirendvalues.Moreprecisely, e.g.,tofindthebasisfunctioninx-directionfornode H ,we setthevalue of

ux

(

H

)

=

1 atthelocationH .Thiscausesextensionofthehorizontalboundarycellsandbendingoftheverticalboundary. Thesub-blockmatricesPxxandPxyrefertothedeformationalongx-axisandy-axis,respectively,inresponsetotheunit horizontaldisplacementatthecoarsenode H ,asshowninFig.5(a).Similarly, Py y andPyxrefertothedeformationalong y-axisandx-axis,respectively,inresponsetotheunitverticaldisplacementatthecoarsenode H ,asshowninFig.5(a).

Oncetheboundaryvaluesarefound,theinternalcellsaresolvedsubjectedtoDirichletvaluesfortheboundarycells.An illustrationofa basisfunction obtainedusingthisalgorithm ispresentedinFig.6.It isworth tobe emphasized thatthe localboundary,whencrossedbyfractures,willbeenriched byjumpenrichmentfunctions.Specially,thereexistsachance that fracture tipsend exactlyatlocal boundaries. Inthat case, since the reduced-dimensionalboundary condition solves 1D problems alongthe boundary, thetip enrichmentis replaced bythe jumpenrichment, andconsequentlythe original fine-scaleXFEMstiffnessmatrixisalsomodifiedinthesamemanner.

Note that the illustrated basis function captures the fractures, because of theXFEM enrichment procedure. The basis function NH

i willbestoredinthecolumni oftheprolongationoperatorP.Onceallbasisfunctionsarefound,theoperator P isalsoknownandonecanproceedwiththemultiscaleprocedureasexplainedbefore.

Fig. 7. Illustration of the 3 categories of Internal, Edge, and Vertex cells, corresponding to the position of each fine cell within a coarse element.

Next,weexplainhowthebasisfunctionscanbealgebraicallycomputedbasedonthegivenXFEMfine-scalesystem.This crucialstepallowsforconvenientintegrationofourmultiscalemethodintoagivenXFEMsimulator.

3.3. Algebraicconstructionofmultiscalebasisfunctions

Thebasisfunctionformulation(16) subjectedtothelocalboundarycondition(17) canbeconstructedandsolvedpurely algebraically.Thisisimportant,sinceitallowsforconvenientintegrationofthedevisedmultiscalemethodintoanyexisting XFEMsimulator.

Consider the coarse cell (local domain) asshownin Fig. 7.The cells are split into3 categories ofinternal, edge and vertex(node),dependingontheirlocations[33].

Notethat the vertexnodesarein factthe coarsemesh nodes,wherethecoarse-scale solutionwill becomputed. The basisfunctionsareneededtointerpolatethesolutionbetweenthevertexcellsthroughtheedgeandinternalcells.

Todevelopthebasisfunctions,firstthefine-scalestiffnessmatrixKh ispermuted,suchthatthetermsforinternal,then edge,andfinallythevertexcellsappear.ThepermutationoperatorT assuchreorders Kh_{into K}v _suchthat

Kv

=

TKhTT

=

T

⎡

⎣

KKauuu KKaaua KKubab Kbu Kba Kbb

⎤

⎦

TT

=

⎡

⎣

KKE II I KE EKI E KE VKI V KV I KV E KV V

⎤

⎦ .

(19)

Here, I representstheinternalnodes,E representstheedgenodesandV representsthevertexnodes.Thepermutedlinear system,therefore,reads

⎡

⎣

KKI IE I KKE EI E KKE VI V KV I KV E KV V

⎤

⎦

⎡

⎣

ddIE dV

⎤

⎦ =

⎡

⎣

ffIE fV

⎤

⎦ .

(20)

NotethatthepermutedsystemcollectsallentriesoftheXFEMdiscretesystembelongingto I,E,andV cells.Therefore, theXFEMenrichmententriesduetotipsandjumpsarewithintheircorresponding I,E,andV entries.

Thereduced-dimensionalboundaryconditionisnowbeingimposedby replacingthe2D equationfor E bya1D XFEM discretesystem.ThisleadstheentryK

¯

E I tovanish,astherewillbenoconnectivitybetweentheedgeandinternalcellsfor theedgecells.This1Dedgeequationscanthenbeexpressedas

K_{E E}R dE

+

KE VR dV

=

0

.

(21)

Knowing thatthe solutionatvertexcells willbe obtainedfromthe coarse-scalesystem, thereorder fine-scalesystem matrixcannowbereducedto

⎡

⎣

KI I0 KKI E_{E E}R KKI V_{E V}R 0 0 KH

⎤

⎦

⎡

⎣

d  I dE dV

⎤

⎦ =

⎡

⎣

00 fH

⎤

⎦ .

(22)

Note that theequations forbasis functionsdonot consider anysource termsin their right-hand-side, andthat we have replacedtheequationsforthevertexcellsbythecoarse-scalesystem(13).

Oncethe coarse-scale solutionsd_V are found,the upper-triangular matrixofEq.(22) can be convenientlyinvertedto obtain theprolongationoperator. Moreprecisely, giventhecoarsenodes solutionsd_V,onecan obtainthesolution atthe edgevia

dE

= −(

KE ER

)

−1KRE VdV

=

PEdV

,

(23) similarly,thesolutionattheinternalcellsreads

dI

= −

K−I I1

(

KI EdE

+

KI VdV

)

= −

K−_{I I}1

(

−

KI E

(

K_{E E}R

)

−1KRE V

+

KI V

)

dV

=

PIdV

.

(24) Notethat PE andPI arethesub-matricesoftheprolongationoperator,i.e.,

d

=

⎡

⎣

d  I d_E d_V

⎤

⎦ =

⎡

⎣

−K

−1 I I

(

−KI E

(

KE ER

)

−1KRE V

+

KI V

)

−(

K_{E E}R

)

−1KRE V IV V

⎤

⎦



P d_V

.

(25)

Here, IV V isthediagonalidentitymatrixequaltothesizeofthevertexnodes.

After defining theprolongationoperatoralgebraically, basedonthe entriesofthe2D XFEM(forinternal cells)and1D XFEM(foredgecells),onecanfindthemultiscalesolution.

3.4. Iterativemultiscaleprocedure(iMS-XFEM)

ThemultiscalesolutionswiththeaccurateXFEMbasisfunctionscanbeusedtoprovideanapproximateefficientsolution formanypracticalapplications.However,itisimportanttocontroltheerrorandreduceittoanydesiredtolerance[32] if needed. Assuch, theMS-XFEM ispaired withafine-scale smoother(here,ILU(0)[47])to allowforerrorreduction. Note thatthisiterativeprocedurecanalsobeusedwithina GMRESiterativeloop[53] to enhanceconvergencerates.Thestudy ofthemostefficientiterativestrategytoreducetheerrorisoutsidethescopeofthecurrentpaper.Theiterativeprocedure reads

•

ConstructtheP andR operators

•

Iterateuntil

rν+1

₌

_fh

₋

_Kh_dν+1

_<

_e r

MS-XFEMstage:

[1a]

δ

dν+1/2

=

P

(

RKhP

)

−1Rrν

[1b]updateresidualrν+1/2

Smoothingstage:iteratefrom j

=

1 till j

=

nstimes [2a]ApplyILU(0):

δ

dν+1/2+1/2(j/ns)

₌

_M−1

ILU(0)rν+1/2+1/2(j−1)/ns [2b]Updateresidualrν+1/2+1/2(j/ns)

Here,theapproximateILU(0)smootheroperatorMILU(0)isappliedfornstimesontheupdatedresidualvector. 4. Numericaltestcases

Inthissectionseveraltestcasesare consideredtoinvestigatetheperformanceofMS-XFEM bothasapproximatesolver andintegratedwithintheiterativeerrorreductionprocedure.

4.1. Testcase1:singlefractureinaheterogeneousdomain

Inthefirsttest,asquare2DdomainofL

×

L withL

=

10 [m]isconsidered,whichcontainsasinglehorizontalfracture inits center,asshowninFig.8(a).The fine-scalemeshconsistsof40

×

40 cells,whiletheMS-XFEM containsonly 5

×

5 coarsegrids.Thisresultsinacoarseningratioof8

×

8.TheheterogeneousYoung’smodulusdistributionisshowninFig.8(b), whilethePossion’sratioisassumedtobeconstant0.2everywhere.ThefracturetipcoordinatesareshowninFig.8(a).The Dirichlet boundary condition is set at the south face, while the north boundary is under distributed upward load with

q

=

5

×

105 [N/m]magnitude.Theeastandwestboundariesaresetasstress-free.

ResultsareshowninFig.9.ItisclearthattheresultsofMS-XFEMononly5 gridcellsisinreasonableagreementwith thatofthefine-scaleXFEMsolverusinga40

×

40 mesh.NotethatnoenrichmentfortheMS-XFEMisused,andthebasis functionsarecomputedusingtheXFEMmethodonlocaldomains.

The errorat each node i is defined asthe difference between fine-scale XFEM (i.e., ui,f) and MS-XFEM (i.e., ui,M S) solutions,i.e.,

Ei

=

ui,M S

−

ui,f

.

Fig.10illustratestheerrordistributionforTestCase1.Notethatthezonesofhigherrorsaremainlynearthefracture. Thestressfield

σy y

,obtainedfromthedisplacementfield,isalsoplottedinFig.11(a)and(b),respectively,forfine-scale

Fig. 8. Testcase1:(a)illustrationofthemodelsetup,(b)heterogeneousYoung’smodulusdistribution.NotetheunitsareSI.Theblacklinesrepresentthe coarsemesh.

Fig. 9. TestCase1:displacementfieldforaheterogeneousfracturedreservoirusing(a)finescaleXFEMand(b)MS-XFEM.Blacklinesrepresentthecoarse scalemesh.

Fig. 10. TestCase1:displacementdifferencebetweenfinescaleXFEMandMS-XFEMin(a)

x direction (b) y direction.

Blacklinesrepresentcoarsescale mesh.

Fig. 11. TestCase1:stressfieldσy ycomparisonbetween(a)fine-scaleXFEMand(b)MS-XFEM(c)iMS-XFEMafter5iterationswith5smoothingloopsper

iteration.

Fig. 12. Testcase:basisfunctionsofsinglefracturetestcase.Singlediscontinuityiscapturedbyaxialequilibriumandtransverseequilibriumsolutions.This basisfunctioniscenteredatcoordinate(4,6).

AbasisfunctionforafracturedlocaldomainisillustratedinFig.12.Notethatthediscontinuityiscapturedbythebasis functions,sinceXFEMisusedtosolveforit.

Coarsening ratioindicates the numberof fine-cells ineach coarse cell. The effectofthe coarseningratio isshownin Fig.13.Todothisstudy,thefine-scalemeshiskeptconstant,whilethecoarseningratioischanged.Thisresultsindifferent coarse-scalesystemsizestosolveforthesamefine-scalesystem.Thenormalizederrore inthisfigureiscomputedusing

ei

=

ui,M S

−

ui,f

2

N

,

∀

i

∈

x

,

y

,

(26)

where N is the number of fine-scale mesh nodes.ui,M S and ui,f denote the MS-XFEM solution displacement field and fine-scalesolutiondisplacementfield,respectively.

The MS-XFEM errorsare dueto thelocalboundary conditionsusedto calculatethe basis functions,andalso because noadditionalenrichmentfunctionsareimposedatcoarsescale.Notethismeansthatforheterogeneousdomainstherecan be a finerresolution forcoarse cells atwhich thelocal boundaryconditions imposemore errorscompared withcoarser resolutions.Inspiteofthis,Fig.13shows,forthisexample,adecayingtrendoftheerrorwithrespecttothefinercoarse meshisobserved.

As discussed in section 3.4, one can pairthe MS-XFEM in an iterative strategy in whichthe error is reducedto any desiredlevel[33].FromFig.14thatwith5timesoffinescalesmoothersappliedinthesecondstageafter5iterationsthe multiscalesolutionerrorisdecreasedsignificantly.AlsoasshowninFig.11(a)and(c)after5iterationsthestressfieldsof iterativeMS-XFEM(iMS-XFEM)andfine-scaleXFEMmatchwell.ConvergencehistoryoftheiMS-XFEMisshowninFig.15. Thebluecirclerepresentsthenormalizederroroftheresultwithouttheiterativestrategy, whiletheredcirclerepresents thenormalizederrorafter5 iterationswithin each5smoothingoperationsare applied.These2circlescorrespondto the resultsshowninFig.14andFig.10.Differentsmoothingstepsperiterationvaluesnscanbeused.NotethatneitherGMRES [53] noranyotheriterativeconvergenceenhancingprocedureisusedhere.Clearly,onecanreducethemultiscaleerrorsto themachineaccuracybyapplyingiMS-XFEMiterations.Inparticular,forpracticalapplications,onecanstopiterationsafter

Fig. 13. Test case 1: change of normalized errors with different coarsening ratios.

Fig. 14. TestCase1:displacementdifferencebetweenfinescaleXFEMandMS-XFEMin(a)

x direction (b)

y direction after 5iterations,withinwhich5 smoothingstepswereemployed.Solidblacklinesrepresentthecoarsemesh.

afewiterationcounts,oncetheerrornormfallsbelowacertaintolerance

τ

,consideringthelevelofuncertaintywithinthe parametersoftheproblem,i.e.,

ei



τ

,

i

=

x

,

y

.

(27)

Here,

τ

ischosenas10−10toconfirmconvergencecanbeachievedtomachineaccuracy. Fig.15showstheconvergencehistoryfordifferentsmoothingstepcountsnsperiteration. 4.2. Testcase2:heterogeneousreservoirwithmultiplefractures

The second test caseis setto modeldeformation in aheterogeneous reservoir withmorefractures. The size andthe heterogeneous properties of thistest caseare the same as those intest case1. Here, more fracturesare considered. In addition,comparedtotestcase1,theeastandwestboundariesarealsoobservingdistributedloads,asshowninFig.16.

Simulation resultsfor both fine-scaleXFEM andMS-XFEM are shown inFig. 17.Multiscale solutions are obtainedby 8

×

8 coarseningratio,appliedtothefine-scalemeshwith40

×

40 finecells.Thisleadstosignificantcostreduction,asthe multiscale coarsesystemhasonly 5

×

5 cells withnoadditionaldegreesoffreedom duetoenrichments. Theblacklines on Fig.17(b)representthe coarsescalemesh.It isclearthat theMS-XFEM(without anyiterations)resultissignificantly differentthanthefine-scaleoneinthezonesoflargedeformations.

Stressfieldsforallthreecasesoffine-scaleXFEM,MS-XFEMandiMS-XFEMareshowninFig.18.TheiMS-XFEMresults areshownafter10iterationswithns

=

5.Notethatafteriterations,multiscalestressdistributionsagreewiththefine-scale one.Additionally,onecanemployadaptivecoarseningratios,sotominimizetheneedtouseaniterativestrategy[54].

Fig. 15. Testcase1:iterationhistoryofiMS-XFEMprocedurewithdifferentnumberofsmoothingperstep.Errorsfordisplacementin

x (a)

and y (b) directionsareshown.TheblueandtheredcirclesrepresentthenormalizederrorscorrespondingtotheresultsshowninFig.10andFig.14,respectively.

Fig. 16. Test case 2: multiple fractures within a heterogeneous reservoir under tension stress across three boundaries.

Fig. 17. TestCase2:displacementfieldforaheterogeneousmediawithmultiple fracturesfor(a)finescaleXFEMand(b)MS-XFEMwithoutiterative improvements.Notethe40×40 fine-scaleXFEMsystemissolvedbyapplyingonly5×5 coarsecells,usingMS-XFEM.Blacklinesrepresentthecoarse scalemeshlines.

An example oftwo basis functionsfor thistest caseis shownin Fig. 19. It illustrates, in2 differentviews, how two fracturesarecapturedbythelocalbasisfunctions.

The iMS-XFEM procedure, as explained before, is now applied to reduce the multiscale errors to machine precision. Fig.20(c)and(d)presentimprovedresultscomparedto20(a)and(b),after10iterationswith5timesfine-scalesmoother appliedineach iteration.Thelargedeformation ofthemostleft fracturestill requiresmoreiterationshoweverthe

differ-Fig. 18. TestCase2:stressfieldσy ycomparisonbetween(a)finescaleXFEMand(b)MS-XFEM(c)MS-XFEMafter10iterationswith5roundsoffinescale

smootherappliedinthesecondstage.

Fig. 19. Illustration of the basis functions Py ywhere two discontinuities are captured in a local coarse cell.

encearoundotherthreefractureshavebeenimprovedsignificantly.NotethatnoGMRESnoranunconditionally-convergent smootheris used,butthe ILU(0)isused asthesecond stage smoother.ILU(0) provides anapproximate (incomplete) de-compositionofthefinescalesystemformaintainingtheefficiency[53].AsshowninFig.21,ahighernumberofsmoothing steps wouldresultinfasterconvergence.However,it comeswithadditionalcomputationalcosts. Similarastheextensive studiesformultiscalesimulationofflowinporousmedia(seee.g.[55–57]),fordeformationsimulation,aCPU-basedstudy shouldbeperformedtoobtaintheoptimumcombinationofcoarse-gridresolution,typeandcountofsmoothingsteps.

4.3. Computationalefficiencyanalysis

AppropriateCPU-basedspeedupanalysesisoutofthescopeofthispaper.However,toprovideanestimateofthe com-putationalefficiencyoftheproposedMS-XFEM strategy,comparedwithafine-scaleXFEM,anoperator-cost-basedanalysis is presented below. Then, a scalability test with respect to the performance of error reduction strategy with increasing densityofthefracturesisalsostudiedandpresented.

4.3.1. Approximatespeedupbasedonanoperation-basedanalysis

Theapproximateoperation-basedapproachbuildsontheassumptionthatsolvingalinearsystemofsize

ζ

costs O

(ζ

m

)

operations, i.e.,C

=

α

ζ

m,whereC isthecost function.As such,theapproximatecostfunction forfine-scaleXFEMsystem (i.e.,Cf)withNf nodesis

Cf

=

α

(

2Nf

+ ψ)

m

.

(28)

Here,

ψ

representsthenumber ofextradegreesoffreedom duetoXFEM enrichments.also sincethereare x and y

dis-placementdirectionsforeachnode,afactor2isused.

MS-XFEM (withandwithoutiterations)provides anapproximatesolution tothefine-scalesystemina procedurethat includesthreemainsteps:(1)calculatebasisfunction,(2)solvecoarsescalesystem,and(3)iterativelyimprovetheresults. Considerthecoarseningratioof



,thus,thecoarsegirdssizeNc is

Fig. 20. TestCase2:displacementdifferencebetweenfinescaleXFEMandMS-XFEM(a)

x direction (b) y direction (c) x direction after

10iterations(d)

y

direction after10iterations.5roundsoffinescalesmootherareusedinthesecondstage.Blacklinesrepresentthecoarsemesh.

Fig. 21. IterationhistoryofiMS-XFEMprocedurewithdifferentnumberofsmoothingperstep.Errorsfordisplacementin

x (a)

and y (b) directionare shown.

Nc

=

Nf



.

Assuming that

γ

enrichmentsper eachbasis functionis considered,basis functionsforeach coarsenode ofa coarsecell costs

Fig. 22. Speedupfactorfordifferentcoarseningratiosanddifferentiterationnumbersforfine-scalegirdof

N

f=108(a)and

N

f=104(b).TheextraDOFs

duetoXFEMenrichmentisψ=0.1Nf is107.Forbothcasesthenumberofsmoothingstepsperiterationissetto

n

s=8.

Cbasis

=

2

×

α

(

+

γ

)

m

,

(29)

wherethefactor2representsdisplacementsinx and y directions.Noteall Nc coarsenodeswillhavebasisfunctions. Thecoarsescalesystementailsnoextradegreeoffreedomduetoenrichments.Assuch,itscostCc canbeapproximated as

Cc

=

α

(

2Nc

)

m

=

α

(

2

×

Nf



)

m

_.

₍₃₀₎

Againthefactor2isforthe2D (x

,

y) unknowns.ThecostofthelinearsmootherI LU

(

0

)

isapproximatedbysettingm

=

1 andaconstantis

β

,i.e.,

Csmooth

=

ns

× (β × (

2Nf

))

1

,

(31)

wherethensisthenumberofsmoothingstepsperiMS-XFEMiterations.

Ifnoiterativestrategyisapplied,then,theMS-XFEMcomputationalcostCM S canbeapproximatedas CM S

=

Nc



2

α

(

+

γ

)

m



+

α

(

2Nc

)

m

,

(32) whereNc

=

Nf  .

IftheiterativestrategyisappliedforNit times,thentheiMS-XFEMcomputationalcostCiM S canbeapproximatedas CiM S

=

CM S

+

Nit



α

(

2Nc

)

m



Cc

+

ns

(β

× (

2Nf

))



Csmooth



.

(33)

NotethatifNit

=

0,thecostfunctionCiM S becomesidenticalwithCM S.

Thespeedupfactor

φ

isnowintroducedasthecostratiosbetweeniMS-XFEMandfine-scaleXFEM,i.e.,

φ

=

CiM S

Cf

.

(34)

ThisfunctionresultsinthespeedupforMS-XFEM(noiterations)whenNit

=

0.

Toprovideabettersenseaboutthespeedup,twotestswithdifferentfinescalegridnumbersofNf

=

104andNf

=

108 areconsidered.ForthecasewithNf

=

104,coarseningratiosof



∈ {

2

×

2

,

5

×

5

,

10

×

10

,

20

×

20

,

50

×

50

}

areconsidered, while for the casewith Nf

=

108,



∈ {

10

×

10

,

20

×

20

,

50

×

50

,

100

×

100

,

1000

×

1000

]

are employed. It is assumed that

α

/β

=

1,meaningthescalabilityfactorofthesmootherisclosetoalinearsolver.Inaddition,wesetm

=

1

.

3,asthe exponent oflinear solver computational complexity. The samefactor was used in theliterature [58]. Furthermore, extra unknownsduetoXFEMenrichmentsareconsideredto be10% offine-scalegridcells, i.e.,

ψ

=

0

.

1Nf.The basisfunctions arethenconsideredtohavetheseextraDOFsequallydistributed,whichleadsto

γ

≈

ψ

Nc.Withthesesettings,onecanstudy

theoverallapproximatespeedupoftheproposedMS-XFEM,basedonanoperational-basedanalyses.ForthatMS-XFEM(no iterations),iMS-XFEMwithNit

=

2 and Nit

=

5 areconsidered.ForbothiMS-XFEMcases,ns

=

8.Thespeedupfactor

φ

for alltestsisplottedinFig.22.

FromFig.22,itisclearthattheMS-XFEM providesasignificantspeedupcomparedwithfine-scaleXFEM.Thespeedup becomesbigger whentheproblemsizeisbigger.Also,fromthepresentedoperational-basedanalyses,onecanconsideran

Fig. 23. IllustrationoftheincreasingfracturedensityfromphaseF1toF7,withadditionaldegreesoffreedomof94,132,192,262,300,366,404, respec-tively,forF1toF7.

Fig. 24. IterationhistoryofiMS-XFEMprocedurewithdifferentnumbersofextraDOFsperstepfordisplacementin

x (a)

and y (b) directions.F1toF7 phasesrepresenttheincreasingfracturedensityasshowinFig.23.

optimumcoarseningratio,inwhichthemultiscaleprocedurewouldperformoptimal. Notethat,similarasformultiphase flow simulations[57,59],aCPU-based analysesisneededtodrawmoreconclusiveremarks abouttherealspeedupofthe MS-XFEM.

4.3.2. Scalabilitytest

Inthistestthe effectofthe densityoffracturesonconvergencerateoftheiMS-XFEMisinvestigated.Thisisspecially relevantto geoscientificapplications.Afixeddomain sizeof10

×

10 m2 _{isconsidered. Fineandcoarsegrids are 40}

_×

₄₀ and 5, respectively. As shownin Fig. 23, starting from phase F 1 tophase F 7 one fracture is addedin the domain. For convenience, homogeneous Young’s modulus of E

=

107 Pa andthe Possion’sratio of 0.2 are considered. The boundary conditionsarethesameasintestcase2.

As itcan be seeninFig. 24,theconvergence rateofthedevelopediMS-XFEM isnot muchaffected by theincreasing fracture density.Thisproof-of-theconcept analysisindicates theapplicabilityofthedevisedmethodforreal-field applica-tions.

5. Conclusion

A multiscale procedure forXFEM is proposed to model deformation ofgeological heterogeneousfractured fields.The method resolves the discontinuities through localmultiscale basis functions, which are computed using XFEMsubjected to local boundaryconditions. The coarse-scalesystem isobtained by usingthe basis functions, algebraically, which does not have anyadditional enrichmentfunctions, incontrast to thelocal basis function systems.This procedure makes the

whenlocalpropertieschange.Nevertheless,therealspeedupwillbefoundwhenthemethodisimplementedincompilable language.

CRediTauthorshipcontributionstatement

Fanxiang Xu: Conceptualization, Methodology, Software, Writing – original draft. HadiHajibeygi: Conceptualization, Methodology, Supervision, Writing – review & editing. LambertusJ.Sluys: Conceptualization, Methodology, Supervision, Writing–review&editing.

Declarationofcompetinginterest

The authors declare that they haveno known competingfinancial interests or personal relationships that could have appearedtoinfluencetheworkreportedinthispaper.

Acknowledgements

FanxiangXuissponsoredbytheChineseScholarshipCouncil(CSC),GrantNumber201807720039(CSC).HadiHajibeygi wassponsoredbyDutchNationalScienceFoundation(NWO)underVidiProject“ADMIRE”,GrantNumber17509(NWO). Au-thorsacknowledgeallmembersoftheDelftAdvancedReservoirSimulation(DARSim)andADMIREresearchgroup,specially YaoluLiu,forfruitfuldiscussions.

References

[1]C.A.Barton,M.D.Zoback,D.Moos,Fluidflowalongpotentiallyactivefaultsincrystallinerock,Geology23(1995)683–686,08.

[2]E.L.Majer,R.Baria,M.Stark,S.Oates,J.Bommer,B.Smith,H.Asanuma,Inducedseismicityassociatedwithenhancedgeothermalsystems,Geothermics 36 (3)(2007)185–222.

[3]M.Rashid,Thearbitrarylocalmeshreplacementmethod:analternativetoremeshingforcrackpropagationanalysis,Comput.MethodsAppl.Mech. Eng.154(Feb.1998)133–150.

[4]T.Bittencourt,P.Wawrzynek,A.Ingraffea,J.Sousa,Quasi-automaticsimulationofcrackpropagation for2dlefmproblems,Eng.Fract.Mech.55(1996) 321–334,09.

[5]R.Cook,FiniteElementModelingforStressAnalysis,Wiley,1995.

[6]J.Jiang,R.M.Younis,Hybridcoupleddiscrete-fracture/matrixandmulticontinuummodelsforunconventional-reservoirsimulation,SPEJ.21(2016) 1009–1027.

[7]T.T.Garipov,P.Tomin,R.Rin,D.V.Voskov,H.A.Tchelepi,Unifiedthermo-compositional-mechanicalframeworkforreservoirsimulation,Comput.Geosci. 22 (4)(2018)1039–1057.

[8]G.Wells,L.Sluys,Three-dimensionalembeddeddiscontinuitymodelforbrittlefracture,Int.J.SolidsStruct.38(Feb.2001)897–913.

[9]H.Hajibeygi,D.Karvounis, P.Jenny, Ahierarchicalfracturemodelforthe iterativemultiscalefinitevolumemethod, J.Comput.Phys.230(2011) 8729–8743.

[10]M.Tene,S.B.Bosma,M.S.A.Kobaisi,H.Hajibeygi,Projection-basedembeddeddiscretefracturemodel(pEDFM),Adv.WaterResour.105(Jul.2017) 205–216.

[11]A.R.Khoei,M.Vahab,E.Haghighat,S.Moallemi,Amesh-independentfiniteelementformulationformodelingcrackgrowthinsaturatedporousmedia basedonanenriched-FEMtechnique,Int.J.Fract.188(Jun.2014)79–108.

[12]Y.Efendiev,J.Galvis,G.Li,M.Presho,Generalizedmultiscalefiniteelementmethods:oversamplingstrategies,Int.J.MultiscaleComput.Eng.12 (6) (2014)465–484.

[13]J.-Y.Wu,F.-B.Li,AnimprovedstableXFEM(is-XFEM)withanovelenrichmentfunctionforthecomputationalmodelingofcohesivecracks,Comput. MethodsAppl.Mech.Eng.295(Oct.2015)77–107.

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

[15]T.Belytschko,Y.Y.Lu,L.Gu,Element-freeGalerkinmethods,Int.J.Numer.MethodsEng.37(Jan.1994)229–256.

[16]J.Melenk,I.Babuška,Thepartitionofunityfiniteelementmethod:basictheoryandapplications,Comput.MethodsAppl.Mech.Eng.139 (1)(1996) 289–314.

[17]N.Moes,J.Dolbow,T.Belytschko,Afiniteelementmethodforcrackgrowthwithoutremeshing,Int.J.Numer.MethodsEng.46(Sep.1999)131–150.

[18]T.Belytschko,T.Black,Elasticcrackgrowthinfiniteelementswithminimalremeshing,Int.J.Numer.MethodsEng.45 (5)(1999)601–620.

[19]A.M.Aragón,A.Simone,Thediscontinuity-enrichedfiniteelementmethod,Int.J.Numer.MethodsEng.112(Sep.2017)1589–1613.

[20]F.P.Meer,L.J.Sluys,Aphantomnodeformulationwithmixedmodecohesivelawforsplittinginlaminates,Int.J.Fract.158(May2009)107–124.

[21]G.N.Wells,L.J.Sluys,Anewmethodformodellingcohesivecracksusingfiniteelements,Int.J.Numer.MethodsEng.50 (12)(2001)2667–2682.

[23]H.Hajibeygi,P.Jenny,Multiscalefinite-volumemethodforparabolicproblemsarisingfromcompressiblemultiphaseflowinporousmedia,J.Comput. Phys.228(Aug.2009)5129–5147.

[24]N.Castelletto,H.Hajibeygi,H.A.Tchelepi,Multiscalefinite-elementmethodforlinearelasticgeomechanics,J.Comput.Phys.331(Feb.2017)337–356.

[25]I.Sokolova,M.G.Bastisya,H.Hajibeygi,Multiscalefinitevolumemethodforfinite-volume-basedsimulationofporoelasticity,J.Comput.Phys.379(Feb. 2019)309–324.

[26]P.Jenny,S.Lee,H.Tchelepi,Multi-scalefinite-volumemethodforellipticproblemsinsubsurfaceflowsimulation,J.Comput.Phys.187(May2003) 47–67.

[27]R.Deb,P.Jenny,Finitevolume-basedmodelingofflow-inducedshearfailurealongfracturemanifolds,Int.J.Numer.Anal.MethodsGeomech.41(Jul. 2017)1922–1942.

[28]G.Ren,J.Jiang,R.M.Younis,AfullycoupledXFEM-EDFMmodelformultiphaseflowandgeomechanicsinfracturedtightgasreservoirs,Proc.Comput. Sci.80(2016)1404–1415.

[29]N.Castelletto,H.Hajibeygi,H.Tchelepi,Hybridmultiscaleformulationforcoupledflowandgeomechanics,in:ECMORXV- 15thEuropeanConference ontheMathematicsofOilRecovery,EAGEPublicationsBV,Aug.2016.

[30]A.Fumagalli,A.Scotti,AnefficientXFEMapproximationofDarcyflowsinarbitrarilyfracturedporousmedia,OilGasSci.Technol.–Revued’IFPEnergies nouvelles69 (4)(2014)555–564.

[31]B.Giovanardi,L.Formaggia,A.Scotti,P.Zunino,UnfittedFEMformodellingtheinteractionofmultiplefracturesinaporoelasticmedium,in:Lecture NotesinComputationalScienceandEngineering,SpringerInternationalPublishing,2017,pp. 331–352.

[32]H.Hajibeygi,G.Bonfigli,M.A.Hesse,P.Jenny,Iterativemultiscalefinite-volumemethod,J.Comput.Phys.227(Oct.2008)8604–8621.

[33]Y.Wang,H.Hajibeygi,H.A.Tchelepi,Algebraicmultiscalesolverforflowinheterogeneousporousmedia,J.Comput.Phys.259(2014)284–303.

[34]E.T.Chung,Y.Efendiev,G.Li,Anadaptivegmsfemforhigh-contrastflowproblems,J.Comput.Phys.273(2014)54–76.

[35]G.Allaire,Homogenizationandtwo-scaleconvergence,SIAMJ.Math.Anal.23 (6)(1992)1482–1518.

[36]A.Abdulle,E.Weinan,Finitedifferenceheterogeneousmulti-scalemethodforhomogenizationproblems,J.Comput.Phys.191 (1)(2003)18–39.

[37]Y.Amanbek,G.Singh,M.F.Wheeler,H.vanDuijn,Adaptivenumericalhomogenizationforupscalingsinglephaseflowandtransport,J.Comput.Phys. 387(2019)117–133.

[38]U.Hornung,HomogenizationandPorousMedia,vol.6,SpringerScience&BusinessMedia,1997.

[39]K.Kumar,T.vanNoorden,I.S.Pop,Upscalingofreactiveflowsindomainswithmovingoscillatingboundaries,DiscreteContin.Dyn.Syst.7(2014) 623–644.

[40]H.-W.Zhang,J.-K.Wu,J.Lü,Z.-D.Fu,Extendedmultiscalefiniteelementmethodformechanicalanalysisofheterogeneousmaterials,ActaMech.Sin. 26(Dec.2010)899–920.

[41]R.Patil,B.Mishra,I.Singh,AnewmultiscaleXFEMfortheelasticpropertiesevaluationofheterogeneousmaterials,Int.J.Mech.Sci.122(Mar.2017) 277–287.

[42]A.Khoei, M.Hajiabadi,Fullycoupledhydromechanicalmultiscale modelwith microdynamiceffects,Int.J.Numer.MethodsEng.115(Apr.2018) 293–327.

[43]M.R.Hajiabadi,A.R.Khoei,Abridgebetweendualporosityandmultiscalemodelsofheterogeneousdeformableporousmedia,Int.J.Numer.Anal. MethodsGeomech.43(Oct.2018)212–238.

[44]K.Bisdom,B.Gauthier,G.Bertotti,N.Hardebol,Calibratingdiscretefracture-networkmodelswith acarbonatethree-dimensionaloutcropfracture network:implicationsfornaturallyfracturedreservoirmodeling,Am.Assoc.Pet.Geol.Bull.98(Jul.2014)1351–1376.

[45]M.HosseiniMehr,C.Vuik,H.Hajibeygi,Adaptivedynamicmultilevelsimulationoffracturedgeothermalreservoirs,J.Comput.Phys.X(May 2020) 100061.

[46]M.HosseiniMehr,M.Cusini,C.Vuik,H.Hajibeygi,Algebraicdynamicmultilevelmethodforembeddeddiscretefracturemodel(f-ADM),J.Comput. Phys.373(Nov.2018)324–345.

[47]E.Chow,Y.Saad,ExperimentalstudyofILUpreconditionersforindefinitematrices,J.Comput.Appl.Math.86(Dec.1997)387–414.

[48]H.Zhou,H.A.Tchelepi,Two-stagealgebraicmultiscalelinearsolverforhighlyheterogeneousreservoirmodels,SPEJ.17(Jun.2012)523–539.

[49]J.T.Camargo,J.A.White,R.I.Borja,Amacroelementstabilizationformixedfiniteelement/finitevolumediscretizationsofmultiphaseporomechanics, Comput.Geosci.(2020).

[50]K.M.Terekhov,Cell-centeredfinite-volumemethodforheterogeneousanisotropicporomechanicsproblem,J.Comput.Appl.Math.365(2020)112357.

[51]H.Wang,TheoryofLinearPoroelasticitywithApplicationstoGeomechanicsandHydrogeology,PrincetonSeriesinGeophysics,PrincetonUniversity Press,2017.

[52]N.Castelletto,S.Klevtsov,H.Hajibeygi,H.A.Tchelepi,Multiscaletwo-stagesolverforbiot’sporoelasticityequationsinsubsurfacemedia,Comput. Geosci.23 (2)(2019)207–224.

[53]Y.Saad,M.H.Schultz,GMRES:ageneralizedminimalresidualalgorithmforsolvingnonsymmetriclinear systems,SIAMJ.Sci.Stat.Comput.7 (3) (1986)856–869.

[54]J.E.Aarnes,S.Krogstad,K.-A.Lie,Ahierarchicalmultiscalemethodfortwo-phaseflowbaseduponmixedfiniteelementsandnonuniformcoarsegrids, MultiscaleModel.Simul.5 (2)(2006)337–363.

[55]M.Tene,Y.Wang,H.Hajibeygi,Adaptivealgebraicmultiscalesolverforcompressibleflowinheterogeneousporousmedia,J.Comput.Phys.300(2015) 679–694.

[56]M.Cusini,A.Lukyanov,J.Natvig,H.Hajibeygi,Constrainedpressureresidualmultiscale(cpr-ms)methodforfullyimplicitsimulationofmultiphase flowinporousmedia,J.Comput.Phys.299(2015)472–486.

[57]M.Tene,M.S.AlKobaisi,H.Hajibeygi,Algebraicmultiscalemethodforflowinheterogeneousporousmediawithembeddeddiscretefractures(f-ams), J.Comput.Phys.321(2016)819–845.

[58]R.J.deMoraes,J.R.Rodrigues,H.Hajibeygi,J.D.Jansen,Multiscalegradientcomputationforflowinheterogeneousporousmedia,J.Comput.Phys.336 (2017)644–663.

[59]A.Kozlova,Z.Li,J.R.Natvig,S.Watanabe,Y.Zhou,K.Bratvedt,S.H.Lee,Areal-fieldmultiscaleblack-oilreservoirsimulator,SPEJ.21(2016)2049–2061.

[60]K.-A.Lie,O.Moyner,J.R.Natvig,Useofmultiplemultiscaleoperatorstoacceleratesimulationofcomplexgeomodels,SPEJ.22(2017)1929–1945.

[61]O.S.Klemetsdal, K.-A.Lie,Dynamiccoarseningand localreorderednonlinear solversfor simulatingtransportinporous media,SPEJ. 25(2020) 2021–2040.