Multiscale simulation of inelastic creep deformation for geological rocks

Multiscale simulation of inelastic creep deformation for geological rocks

Ramesh Kumar, Kishan; Hajibeygi, Hadi

DOI

10.1016/j.jcp.2021.110439

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Ramesh Kumar, K., & Hajibeygi, H. (2021). Multiscale simulation of inelastic creep deformation for

geological rocks. Journal of Computational Physics, 440, 1-19. [110439].

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

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.

www.elsevier.com/locate/jcp

Multiscale

simulation

of

inelastic

creep

deformation

for

geological

rocks

Kishan Ramesh

Kumar

∗

,

Hadi Hajibeygi

FacultyofCivilEngineeringandGeosciences,DepartmentofGeoscience andEngineering,DelftUniversityofTechnology,Stevinweg1,2628CV, Delft,theNetherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Availableonline18May2021 Keywords:

Subsurfaceenergystorage Geomechanics

Nonlinearmaterialdeformation Multiscalebasisfunctions Algebraicmultiscalemethod

Scalablephysics-basednonlinearsimulation

Subsurfacegeologicalformationsprovidegiantcapacitiesforlarge-scale(TWh)storageof renewableenergy,once thisenergy(e.g.fromsolarandwindpowerplants)isconverted to greengases, e.g.hydrogen. Thecritical aspects ofdevelopingthistechnologyto full-scale will involve estimation ofstorage capacity, safety, and efficiency of a subsurface formation. Geological formations are oftenhighly heterogeneous and,when utilized for cyclicenergystorage,entailcomplexnonlinearrockdeformationphysics.Inthiswork,we presentanovelcomputationalframeworktostudyrockdeformationundercyclicloading, inpresenceofnonlineartime-dependentcreepphysics.Bothclassicalandrelaxationcreep methodologiesareemployedtoanalyzethevariationofthetotalstraininthespecimen overtime. Implicittime-integration schemeisemployed topreservenumerical stability, duetothenonlinearprocess.Oncethe computationalframeworkis consistentlydefined using finite element method on the fine scale, a multiscale strategy is developed to represent the nonlinear deformation not only at fine but also coarser scales. This is achieved by developinglocally computed finiteelement basisfunctions atcoarse scale. Thedeveloped multiscalemethodalsoallowsforiterativeerrorreductiontoanydesired level, afterbeingpairedwithafine-scalesmoother. Numerical testcases arestudied to investigatevariousaspectsofthedevelopedcomputationalworkflow,frombenchmarking with experiments to analysing the impact of nonlinear deformation for a field-scale relevantenvironment.Resultsindicatetheapplicabilityofthedevelopedmultiscalemethod inordertoemploynonlinearphysicsintheirlaboratory-basedscaleofrelevance(i.e.,fine scale), yetperformfield-relevant simulations. Thedeveloped simulatoris madepublicly availableathttps://gitlab.tudelft.nl/ADMIRE_Public/mechanics.

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

1. Introduction

A successful transition towards cleaner energy relies on the availability of large scale renewable energy storage. To achieve this, renewable energy can be converted to green fuels such as hydrogen[1] or green methaneto be stored in large-scale quantities in geological subsurface [2]. Compressed air energystorage [3], energy storedin the form ofheat [4,5] are otheralternativestostoreenergyinthesubsurface,whicharecyclicinnaturetoo.Dependingonthesupplyand demand, theenergystoredcan varyin therangeofdifferent time scalesanddifferent capacitiesofstorage volumes[6].

*

Correspondingauthor.

E-mailaddresses:K.Rameshkumar-2@tudelft.nl(K. RameshKumar),H.Hajibeygi@tudelft.nl(H. Hajibeygi). https://doi.org/10.1016/j.jcp.2021.110439

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

Table 1

Constitutiverelationspresentedintheliteraturetoexpresscreepstrainrateandcreepstrainfordifferent mate-rials.

Reference Formulation Material

Konstantin Naumenko[21] ε˙cr=32aσ

n−1

v M s Rock salt

Bérest and Brouard[8] εcr= β0σntme−Q/R Tt+A0(δσ)ne(−Q/R T)t Rock salt Betten[22] ε˙cr= ˙ε0σ_σ0e f f n e( Q R T0− Q R T) _{Rock salt} Tian et al.[23] dεc i j= λc∂G cr ∂σi j Sandstone

Xu et al.[24] ε˙cr=32AnSi jσne−1tm−1exp

_−Q R T  Sandstone Jia et al.[25] εcr=C3C+11(σ¯) C2_tC3+1 _ε˙_cr=_C1(σ¯₎C2_tC3 _{Clayed rock} Yang et al.[26] εcr=s−Dbs cr 1−Db  1 2G1+ 1 2G2(1−e −G2t η1 )+ t 2η2  Coal

SubsurfacestoragetechnologyhasbeenimplementedworldwideinmanycountriesincludingtheNetherlands[7].Tosafely utilize the underground storage capacities (e.g. in porous rocks andsalt caverns[8–10]),it is critical to understand the underlyingphysicsinvolvedwhenacyclicloadisimposed.

Cyclicloadinghasbeenstudiedinminingsciencestostudyhydraulicfracturing,rockcutting,forecastingvolcanichazard, designingtunnelsagainstearthquakes(seee.g.[11–14]).Veryfewresearchershavestudiedtheapplicationsofcyclicloading onoilandgasstoragebothnumericallyandexperimentally[15–17].Relevanthystericrockandfluidconstitutivegoverning equationscanbeimplementedtocapturetheimportantphysicsinthecyclicsubsurfacestoragetechnology[18–20]. How-ever,there isa significantneed forabetter understandingofthe complexitiesinvolvedinthe subsurfacestoragephysics toimplementitsafely,efficiently, andreliably onalongtermbasis.Especiallythisbecomesacruciallyimportantfieldof studywhenstoragecitiesbecomeclosertotheurbanareas.

Ofparticularinterestistodevelopascalable(i.e.,multiscale)mechanicaldeformation modelling concepttoincludenot only theelastic butalso time-dependent nonlinearcreep physics. As geologicaldomains span large lengthscales, having inelastic(nonlinear)andheterogeneousproperties,developmentofamultiscalestrategywhichallowstocapturethese fine-scalephysicsandparametersinthecoarse-scalesystemsiscrucial.Suchamultiscalestrategywouldalsoallowforaccurate quantificationsoftheconsequencesofthepore-pressuredynamicsinthesubsurfacegeologicaldomains.Thismotivatesthe developmentofthecurrentstudy.

Complexthermal,geo-mechanical,andchemicaleffectsmightoccursimultaneouslymakingitdifficulttopredictthe fea-sibilityoftheenergystorage[27].Thisworkisfocused onstudyingtheinfluenceofin-elasticgeo-mechanicaleffects.These effectsincludecreep,plasticityandviscoplasticity.Fewresearchershavestudiedplasticityandvisco-plasticityextensivelyin porousmedia[28–30].Inthisstudywefocusonnonlineartime-dependentcreepdeformation.

Creepdeformation dependsonthe typeofmaterial andthetype ofloadinginvolved. Examplescanbe found indeep resource mining [31] andradioactive nuclear waste disposal [25], in geological reservoirs such as shalerocks [32,33] to name afew.Several constitutivelawshave beenintroduced inthe past tostudy creep andthisstill is an active field of research intherockphysicslabsandgeophysicalfieldstudies.Table1show differentcreepformulationsemployed inthe pastfordifferentmaterials.Inthiswork,theconstitutiveequationformodelling inelasticcreepdeformationischosenbased ontheexistingliteraturewhichiselaboratedinsection3.

The presentwork introduces a multiscale 3D non-linear finite element method (FEM) to study the creep mechanics

that isinvolvedinthecompressionofrocks undernormalandcyclicloading.Classicalcreepandrelaxationcreeparethe two methodologiesthatare usedtoanalyzeandquantifythecreep. Bothare consideredinthiswork. Thecomputational frameworkallowsusingbothexplicitandimplicittime integrationscheme.However,implicittimeintegrationisemployed to guarantee convergence irrespective of the size of the time steps. For field-scale applications, one has to employ the physicalconstitutivelawswhicharevalidatedinthelaboratorysettings(cm-scale)only.

Applicationsoftheseexperimentallawsforlarge(severaltens ofmeters) gridblocks,havingheterogeneouselasticand nonelastic properties,are questionable. To avoid excessive upscaling of these nonlinear terms and highly heterogeneous elasticitycoefficients, a multiscaleprocedure isalso introduced.The multiscale methodisbuilt onfiniteelement (hence, MSFE)toobtaincoarse-scalesystemsusinglocally-supportedfinite-elementcoarse-scalebasisfunctions.Themultiscale for-mulation ofthistype, i.e.local basis function-based,hasbeen applied anddevelopedby various researchersin thepast. ([34–43]).Researchershaveemployed multiscalestrategytostudyflowinporousmedia[44–46],geo-mechanicsinporous media[47–51] andothernon-linearproblems[52–54].Notethattherichdomainofmultiscalemodelling for inelastic de-formationsalsoincludesnon-basis-function-basedapproaches,suchasFE2 _{forcompositematerials[55],FEM-DEM(Discrete}

ElementMethod) forgranular materials [56] andhigh-porosityrocks[57]. Themultiscalestrategy developedinthiswork aimsto (1)develop localcoarse-scalebasis functions,(2)resolvefine-scaleheterogeneitiesincoarsesystem, (3)allow for iterativeerrorreductiontothemachineaccuracy,(4)allowforalgebraicformulationandimplementation.Intheseaspects, itstaysnovelcomparedtotheexistingliterature ofmultiscaleFEMforinelasticdeformationmodelling.Notethatinorder toguaranteeconvergenceoftheiterativemultiscaleerrorreductionstrategy, thedevelopedmultiscaleprocedureispaired withafine-scalesmoother[58,59].

Different creep mechanisms, asinelastic deformation, are then introduced in section 3. Constitutive equationsemployed to modelcreepare presented.The algorithmsemployed to modelcreepare described indetail.Next,the algebraicMSFE procedurefornonlineardeformationmechanismisdevelopedtoallowforfield-scaleanalyseswithoutrelyingonexcessive upscaling ofthe nonlinearcreep physicsandheterogeneous materialpropertiesin section4. MSFEformulationis further extendedbyemployingiterativeMSFE.Numericalresultsarethenpresentedinsection5startingwith(i)Uniaxial compres-sion testcaseforrocksaltforboth fine-scaleandalgebraic MSFE. Next,(ii)triaxial loadingofcoal isstudied usingboth fine-scale anditerativemultiscalestrategy. Next,(iii)tri-axialcyclicloading testcaseis studiedforsandstoneusing fine-scaleanditerativemultiscalestrategy.Fortheabovetestcases,thenumericalresultsarecomparedwithexperimentaldata byusingfittingparameterstounderstandthevariationofdeformationovertime.Next,theinfluenceofcreeponsubsidence isstudiedinthenorthernItalianfield[60] byassumingplainstrain.Finally,insection6,concludingremarksarepresented.

2. Linearelasticity

Themomentumbalancefordrainedsolidporousmediumisgivenby

∇ · (

σ

)

=

f

,

(1)

whereporepressureisconstant.Thestressisgivenby

σ

=

C

:

ε

,inwhichtheelasticitycoefficientmatrixandstraintensor are denotedby C and

ε

,respectively. Here f is thevolumetric forceexpressed as[N/m3_{]. Thegoverning equation when}

expressedintermsofdeformationreads

∇ · (

C

: ∇

s

·

u

)

=

f

.

(2) Here,

∇

s_u

₌

1 2



∇

u

+ ∇

T_u



isthegradientsymmetricofdisplacementvector.TheelasticitytensorisexpressedusingLamé parameters

λ

and

μ

,i.e.,

C

=

⎡

⎢

⎢

⎢

⎢

⎢

⎣

(λ

+

2

μ

)

λ

λ

0 0 0

λ

(λ

+

2

μ

)

λ

0 0 0

λ

λ

(λ

+

2

μ

)

0 0 0 0 0 0

μ

0 0 0 0 0 0

μ

0 0 0 0 0 0

μ

⎤

⎥

⎥

⎥

⎥

⎥

⎦

.

(3) Here

λ

=

_(1+νν₎₍E_{1−2ν and}

μ

=

E

2(1+ν), where E is Youngs modulus and

ν

is Poisson’s ratio. Using Voigt’s notation and transposeofdivergenceoperator,thefinalsetofequationsin3Dare

∂

∂

x



(λ

+

2

μ

)

∂

ux

∂

x

+ λ

∂

uy

∂

y

+ λ

∂

uz

∂

z



+

∂

∂

y



μ

∂

ux

∂

y

+

μ

∂

uy

∂

x



+

∂

∂

z



μ

∂

ux

∂

z

+

μ

∂

uz

∂

x



=

fx (4a)

∂

∂

x



μ

∂

ux

∂

y

+

μ

∂

uy

∂

x



+

∂

∂

y



λ

∂

ux

∂

x

+ (λ +

2

μ

)

∂

uy

∂

y

+ λ

∂

uz

∂

z



+

∂

∂

z



μ

∂

uy

∂

z

+

μ

∂

uz

∂

y



=

fy (4b)

∂

∂

x



μ

∂

ux

∂

z

+

μ

∂

uz

∂

x



+

∂

∂

y



μ

∂

uy

∂

z

+

μ

∂

uz

∂

y



+

∂

∂

z



λ

∂

ux

∂

x

+ λ

∂

uy

∂

y

+ (λ +

2

μ

)

∂

uz

∂

z



=

fz

.

(4c)

2.1. FEMformulationofthesystem

Tosolve equations(4a)-(4c), finite elementformulationis followed.Using minimization ofweighted residual Galerkin approach, thedisplacement ofeach element is interpolatedwithin an element usingfine-scale shape functions. The dis-placementinsideanelement

(

u

˜

e

)

canbeapproximatedfromnodaldisplacements

(

u

ˆ

e,i

)

as

˜

ue

≈

Nnodes



,e

i=1

Ni,eu

ˆ

e,i

.

(5)

Here, u

ˆ

e,i can be written asdisplacements

(

ux

,

uy

,

uz

)

foreach node i. Ni,e are theshape functionsused to interpolate

Fig. 1. Typical creep curve shows the variation of total strain with respect to time when the temperature is constant [61,22].

Ni,e

≈

1

8

(

1

± ξ

1

)(

1

± ξ

2

)(

1

± ξ

3

)

−

1

≤ (ξ

1

, ξ

2

, ξ

3

)

≤

1

.

(6)

UsingtheGaussdivergencetheorem,thefinaldiscreteequationforanelementreads



e DeTC Deu

˜

edV

=



e NeTf dV

+



σ

(

Ne

)

Ttd S

.

(7)

TheabovevolumeintegralsareapproximatedusingGaussquadraturerulebytransformingitfromglobaltolocalcoordinate system. Here, De

=

di vTNe,theapplied stress ontheboundary alongthe normalist (Neumannbc)which isintegrated

overarea



σ ,thevolumetricforceis f ,integratedoveranelementalvolume



e.Theoperatordi vT isgivenby

di vT

=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

∂ ∂x 0 0 0 ∂ ∂y 0 0 0 _∂∂_z 0 _∂∂_{z ∂}∂_y ∂ ∂z 0 ∂ ∂x ∂ ∂y ∂ ∂x 0

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

.

(8)

Equation (7) can be written based on the stiffness matrix( Ae

=



eDe

T_{C D}

e) and rhs force vector (

˜

fe

=



eNe T_{f dV}

₊



σ

(N)

T_{td S)as} Aeu

˜

e

= ˜

fe

.

(9)

Byassemblingthestiffnessmatrixforeachelement,thegloballinearsystemisformedas

Amechu

˜

= ˜

f

.

(10)

After computingelasticdeformation,stressesandstrainsarecomputedusingconstitutivelawfromEquation (1).Next,the nonlinearelasticdeformation(creep)isrevisited.

3. Creepformulation

Creepisthetendencyofsolidmaterialtodeformpermanentlyundercertainloadthatdependontimeandtemperature. Typicalcreepprocesshasthreephaseswhichareprimarycreep(creepratedecreasingovertime),secondarycreep(constant creep rate) and tertiary creep (increasing creep rateuntil failure) asshown inFig. 1. In thiswork, tertiary creep isnot modelled.Whenthestrainsareinfinitesimal,thetotalstrain

ε

t canbesplitas

ε

t

=

ε

el

+

ε

cr

+

ε

th

.

(11)

Here,

ε

el iselasticstrain,

ε

thisthermalstrainand

ε

cristhecreepstrain.Inthiswork,timeindependentplasticstrainslike

viscoplasticity andthermalstrainsarenot included.Themodified creepequation whichinvolvespowerlawandBingham bodyare[26,24]

ε

cr

=

X



t

Fig. 2. Schematic diagrams of classical creep and relaxation creep methodologies [62].

Dependingonthetype ofmaterialandformulationeitherpowerlaw( X

=

1,Y

=

0)orBurgersmodel( X

=

0,Y

=

1)can bechosen.Here, P

= ˙

ε

cr

=

Ae− Q R T

_σ

n−1 vM stm−1 (13a) B

=

ε

cr

=

s



1 2G1

+

1 2G2

(

1

−

e −G2t η₁

₎

₊

t 2

η

2



(13b)

Materialconstantsare A

,

n

,

m inthepowermodelandintheBurgersmodeltheparametersG1

,

G2 aretheshearmodulus,

η

1

,

η

2 isthe viscositycoefficient andt is the time. Here s is deviatoricstress and

σ

vM isvon Mises stress. TheBurgers

model for3D stress state is obtainedas seriesof elements which representHookean body (elastic, spring), Kelvin body (damperandspringinparallel)andNewtonbody(dashpot/damper).Sothestressesoftheindividualelementsaregivenby

sH

=

2G1

ε

H (14a)

sK

=

2G2

ε

K

+

2

η

1

ε

˙

(14b)

sN

=

2

η

2

ε

˙

N

.

(14c)

Thenon-linearstrainsareincorporatedinthegoverningequationbypluggingEquation (11) inEquation (2),

∇ · (

C

:

ε

el

)

= ∇ · (

C

: (

ε

t

−

ε

cr

))

=

f

.

(15)

Thedeformationiscomputedby

∇ · (

C

: ∇

su

)

= ∇ · (

C

:

ε

cr

)

+

f

.

(16)

Usingtheprincipalofvirtualdisplacements,whichstatesthatsumofexternalandinternalworkis0[21],thefinaldiscrete equationcanbewrittenasmodifiedformofEquation (7),i.e.,



e DeTC Deu

˜

edV

=



e NeTf dV

+



σ

(

Ne

)

Ttd S

+



e DTC

ε

crdV (17)

Equation (17) canberewrittenas

Au

=

Fb

+

Fcr

.

(18)

Where, Fcr isthefictitiouscreepforcesandFbisthebodyforces.

3.1. Numericalmethodology

The creepstrain

ε

cr canbe modelledusingclassical creepandrelaxationcreep methodology[62]. Classicalcreepisa

modellingapproach,whenthe stressisassumedtobe constantandthecreepstrain accumulateswithrespecttotime. If thestressisconstant,theelasticstrainremainsconstantwhichwouldmeantotalstrainisaccumulatedovertime.Whereas, therelaxationcreepisanapproachwhenthetotalstrainisassumedtobeconstant,andascreepstrainaccumulates,elastic strain

ε

el reducesover time.Thiswouldimplythatthestress reduceswithtime.Fig.2showstheschematic diagramsof

classicalandrelaxationcreep.

Thesenumericalmethodologiesaresimilarinnaturetothestressandstraincontrolledcompressionexperimentsonthe rocks,whereeitherstressorstrainiscontrolledontherocksamplewhentheexperimentsareperformed.

ThealgorithmsforclassicalandrelaxationcreeparedepictedinAlgorithm1andAlgorithm2respectively.

Thedevelopedcomputationalframeworkallowstousebothimplicitandexplicitformulation.Theresultsshownhereare obtainedfromimplicitformulationwhich ensures computationalstability even forhighertime stepsizes. Thenon-linear creepformulationisdiscretisedimplicitlyas,

Algorithm1:Classicalcreepalgorithm.

1: Input materialandcreepparameters

2: Calculate u,σandεelattime=0,i=0.Initializeε0cr=0 3: Compute thetotalnumberoftimestepsN=Tsimulation/ t 4: for i<N do

5: ComputedeviatoricandvonMisesstress si j,σvM 6: Calculatecreepstrainattimestepi :εi

cr=X



tP+Y B

7: Computetotalstrainattimestepi :εi t=εe+εicr 8: Seti=i+1

9: return u,σ,εt,εel,εcratlasttimestep

Algorithm2:Relaxationcreepalgorithm.

1: Input materialandcreepparameters

2: Calculate u,σ andεelattime=0,i=0.Initializeε0cr=0 andεt=εel 3: Compute thetotalnumberoftimestepsN=Tsimulation/ t

4: for i<N do

5: ComputedeviatoricandvonMisesstress si j,σvM 6: Calculatecreepstrainattimestepi : εi

cr=X



tP+Y B

7: Calculateelasticstrainεi

el=εt−εicr 8: Calculatestressσi_=C:εi

el 9: Seti=i+1

10: return u,σ,εt,εel,εcratlasttimestep

In the classical creep formulation, the stress remains constant with respectto time. However, in relaxation thestress is dependentontime.Ingeneral, thestressterminthecreepmodelisobtainedbasedonthenewtimestep

(

i

+

1

)

.Inthis case,Newton-Raphsoniterativeapproachisemployedtosolveforthedisplacementiteratively.Thenonlinearresidualreads as



i+1

_{= ∇ · (}

_C

_:

_ε

i+1

cr

)

+

f

− ∇ · (

C

: ∇

sui+1

).

(20)

Tofindui+1_{theaboveequationisfurtherlinearized,i.e.,}



i+1

_{≈ }

ν+1

_{≈ }

ν

₊

∂



∂

u





ν

δ

uν+1

,

(21)

whichisthensolvediterativelyuntilconvergenceisachieved,i.e.



i+1

=

0.TheJacobian J

=

∂_∂_u





ν

iscomputedbasedonly ontheelasticpartoftheresidualequation.Thedisplacementincrementiscomputedby

δ

uν+1

= −(

Jν

)

−1



ν

.

(22)

The iterations

ν

are repeated until the residual norm fallsbelow the prescribed threshold, i.e.,

||

ν

||

2

<



r to achieve

convergence.Algorithm3describestheimplicitformulationtomodelcreep.

Algorithm3:Implicitcreepformulation.

1: Initialization materialandcreepparameters,i=0,ε0

cr=0

2: Compute elasticdeformation,totalnumberoftimesteps(N)andsettheJacobian( J )terms 3: for i<N do

4: Set iterationcounterν=0,creepforcesFi+1 cr =Ficr 5: Calculate theresidualn+1_{= ∇ · (C:ε}n+1

cr )+f− ∇ · (C: ∇sun+1).

6: while||ν_||

2<randν<max(ν)do

7: Updatetheresidualν+1

8: Calculatethedisplacementincrementδuν+1_{= −(J}ν₎−1_ν

9: Updatethedisplacementuν+1_=uν_{+ δu}ν+1 10: Calculatestress(σν+1_),strain(εν+1

t ),creepstrainεν+

1

cr andfictitiouscreepforcesFνcr+1 11: Setν=ν+1

12: end 13: Seti=i+1 14: end

4. Multiscalefiniteelementmethodfornonlineardeformationundercyclicloading

Theapproximatefinescaledisplacementu

¯

h _{iscomputedfromcoarsescalesystemu}

_¯

H _{usingProlongationoperator P}h H,

Fig. 3. Multiscalemeshwithcoarseelements.Theblackdotsatthecornersarethecoarseelementsnodesandthecoloured blocksarecoarseelements. Thesmallercubesinsidethecolouredelementsarethefinescaleelements.(Forinterpretationofthecoloursinthefigure(s),thereaderisreferredtothe webversionofthisarticle.)

¯

uh

=

Ph_Hu

¯

H

.

(23)

Fig.3showsthecomputationaldomain.Theblackdotsatthecornersarethenodesofthecoarseelementsandthecoloured blocksarecoarseelements.TheprolongationoperatorisformedfromdisplacementbasisfunctionsNH_{.Thebasisfunctions}

arecomputedbylocalmomentumbalanceequationswithineachcoarsecellsas

∇ · (

C

: ∇

sNH_i

)

=

0 in



u (24a)

∇ · (

C

: ∇

_||sNH_i

)

=

0 on

∂



u (24b)

N_iH

(

xj

)

= δ

i j

∀

j

∈ {

1

...

nuH

}

(24c)

The localconditions aresolved in eachcoarse element



u_i forbasis functions NH_i ofi coarse block. Here

δ

i j isthe

Kro-neckerdeltaandxj isthepositionvectorofcoarsemeshnodeassociatedtoj-thunknown.

∇

sand

∇

_||s denoteapproximate

divergenceandsymmetricgradient operatorsactingintangential directionofcoarseelements boundary(local coordinate system)tosolvelocalboundaryconditionsat

∂

u forfacesandedges,respectively.

Basedonthewirebasketdecompositionoffinescalesystem[63],apermutationmatrixW isemployedtorearrangeall thenodesintheorderofinterior,faces,edgesandvertex.Basedonthatthefinescalesystemcanbewrittenas

ˆ

Ahu

ˆ

= ˆ

fh

.

(25)

DenotingI,F,EandVasinterior,face,edgeandvertexnode,relativetothecoarsemeshconfiguration,allowsforre-stating

thepermutedsystemas

ˆ

Ah

=

WTAhW

=

⎡

⎢

⎢

⎣

ˆ

AI I A

ˆ

I F A

ˆ

I E A

ˆ

I V

ˆ

AF I A

ˆ

F F A

ˆ

F E A

ˆ

F V

ˆ

AE I A

ˆ

E F A

ˆ

E E A

ˆ

E V

ˆ

AV I A

ˆ

V F A

ˆ

V E A

ˆ

V V

⎤

⎥

⎥

⎦

and

ˆ

uh

₌

_WT_u

_¯

h

₌

⎡

⎢

⎢

⎣

ˆ

uI

ˆ

uF

ˆ

uE

ˆ

uV

⎤

⎥

⎥

⎦ ,

f

ˆ

h

₌

_WT_fh

₌

⎡

⎢

⎢

⎢

⎣

ˆ

fI

ˆ

fF

ˆ

fE

ˆ

fV

⎤

⎥

⎥

⎥

⎦

.

TheabovematrixisreducedbyperformingGaussianeliminationtwicewhichgives

⎡

⎢

⎢

⎣

ˆ

AI I A

ˆ

I F A

ˆ

I E A

ˆ

I V 0

ˆ

SF F

ˆ

SF E S

ˆ

F V 0 0 S S

ˆ

E E S S

ˆ

E V 0 0 S S

ˆ

V E S S

ˆ

V V

⎤

⎥

⎥

⎦

⎡

⎢

⎢

⎣

ˆ

uI

ˆ

uF

ˆ

uE

ˆ

uV

⎤

⎥

⎥

⎦ =

⎡

⎢

⎢

⎢

⎣

I 0 0 0

− ˆ

AF IA

ˆ

−_{I I}1 I 0 0

− ˆ

AE IA

ˆ

−I I1

−ˆ

SE FS

ˆ

−F F1 I 0

− ˆ

AV IA

ˆ

−_{I I}1

−ˆ

SV FS

ˆ

−_{F F}1 0 I

⎤

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎢

⎣

ˆ

fI

ˆ

fF

ˆ

fE

ˆ

fV

⎤

⎥

⎥

⎥

⎥

⎦

.

(26) Here,

ˆ

Si j

= ˆ

Ai j

− ˆ

Ai IA

ˆ

−I I1A

ˆ

I j

∀(

i

,

j

)

∈

F

,

E

,

V

×

F

,

E

,

V (27)

ˆ

S Si j

= ˆ

Si j

− ˆ

Si FS

ˆ

−F F1S

ˆ

F j

∀(

i

,

j

)

∈

E

,

V

×

E

,

V (28)

Byapplyingthereducedorderboundaryconditionsitcanbestatedthat



AF Fu

ˆ

F

+ 

AF Eu

ˆ

E

+ 

AF Vu

ˆ

V

=

0 (29)



AE Eu

ˆ

E

+ 

AE Vu

ˆ

V

=

0

.

(30)

The matrices



AF F

,



AF E

,



AF V

,



AE E

,



AE V are computed by modifying the divergence andgradient operators for faces and

edgesrespectivelyandsolvingEquation (24).Theinformationfromthefacesisinterpolatedtofaces,edgesandverticesby reducedboundaryconditions,andinformationfromedgesisinterpolatedtoedgesandvertices.Thismeans

ˆ

SF F

≈ 

AF F

, ˆ

SF E

≈ 

AF E

, ˆ

SF V

≈ 

AF V

, ˆ

S SE E

≈ 

AE E

,

S S

ˆ

E V

≈ 

AE V

.

(31)

ByimplementingEquation (24),themodifiedsystemis

⎡

⎢

⎢

⎣

ˆ

AI I A

ˆ

I F A

ˆ

I E A

ˆ

I V 0



AF F



AF E



AF V 0 0



AE E



AE V 0 0 0 A

ˆ

V V

⎤

⎥

⎥

⎦

⎡

⎢

⎢

⎣

ˆ

uI

ˆ

uF

ˆ

uE

ˆ

uV

⎤

⎥

⎥

⎦ =

⎡

⎢

⎢

⎢

⎣

ˆ

fI

ˆ

fF

ˆ

fE

ˆ

fV

⎤

⎥

⎥

⎥

⎦

.

(32)

Here A

ˆ

V V isthecoarsescalematrixA

ˆ

H.Thecoarsescalesystemissolvedforu

¯

H andallthedeformationscanbecomputed

byexpressingcoarsescaledeformationasu

¯

h

=

Ph_Hu

¯

H with

Ph_H

=

W

⎡

⎢

⎢

⎢

⎢

⎣

− ˆ

A−_{I I}1

( ˆ

AI F



A−F F1

(

AF E



A−E E1



AE V

))

− ˆ

A−I I1A

ˆ

− 1 I E

(

−

A− 1 E E



AE V

)

− ˆ

A−I I1A

ˆ

I V



A−_{F F}1

(

AF E



A−_{E E}1A

ˆ

E V

)

− 

A−_{F F}1



AF V

−

A−_{E E}1



AE V IV V

⎤

⎥

⎥

⎥

⎥

⎦

.

(33)

Therestrictionoperatoristhetransposeofprolongationmatrix,i.e.,

R_hH

= [

Ph_H

]

T

.

(34)

Thecoarsescalesystemisformulatedtosolveforu

¯

H

ˆ

AH

=

RhHA

ˆ

hPhH

,

ˆ

fH

=

R

ˆ

fh

,

and A

ˆ

Hu

¯

H

= ˆ

fH

.

(35)

UsingthecoarsescaledeformationandemployingEquation (23),finescaledisplacementsarecomputed.

Algorithm4:AlgebraicMSFE.

1: Formulate PermutationoperatorW

2: Solve forreducedbcsAF F,AF E,AF V,AE E,AE V

3: Formulate prolongationoperatorPhH

4: Compute restrictionoperatorRH h= (P

h H)

5: Formulate coarsescalesystemAˆH=R_hHAˆhPhHand fH=Rˆfh

6: Solve forcoarsesystemu¯H

7: Map ittofinescalesystemu¯h_=Ph Hu¯

H

Fig.4 showsthe basisfunctionsatthevertex, edge,faceandinteriorcoarse element.Thecoarseningratio alongeach axisisgivenby

Cr

=

Number of fine elements along each axis

Number of coarse elements along each axis

.

(36)

TheMSFEprocedureforinelasticcreepdeformationisshowninAlgorithm4.Thenormalizedaverageerrorbetweenthe fine-scaleandmultiscalesolutionisgivenby



=

||

uFS

−

uMS

||

2

||

uFS

||

2

.

(37)

ThesolutionsobtainedfromMSstrategyemployfewerdegreesoffreedomwhichleadstothereductionofcomputational effortcomparedtofinescalestrategy.

Fig. 4. Examplesofbasisfunctionsobtainedfromalgebraicmultiscaleinuni-axialcompressiontestcasewithCr=5×5×5 alongeachaxisandfinescale domainof10×10×10.

4.1. Iterativemultiscalesolver

When creep strains are included, the solutions obtained frommultiscale formulation are approximate.To reduce the error to a desiredtolerance, an iterativeprocedure is used,employing a fine scale smootherlike ILU or GMRES [64,65]. In thiswork, ILU(0)andGMRES isemployed forfew test casesto compare thehighercomputational efficiencyandalso improveMSsolutionstothedesiredtolerance.Thestoppingcriteriausedhereer isgivenby

Algorithm5:IterativeAlgebraicMSFE.

1: Formulate prolongationandrestrictionoperatorsandinputtolerance(tol). 2: Iterate untiler<tol

3: Multiscalestage1: uν+1 2=Ph H(R H hA h_Ph H)− 1_RH hrν 4: Compute residualrν+1 2= ||f−Ahuν+ 1 2|| 5: Smoothingstage2: uν+1_{= (M} sm)−nsrν+ 1 2 er

=

||

uFS

−

uMS

||

2

Total number of nodes

.

(38)

Whereeristheaverageerrorbetweenfine-scaleandmultiscalesolutions.Inrealisticfieldcases,finescalesolutionwillnot

beavailable.Inthosecases,thestoppingcriterioncanbetherelativeresidual,i.e.,

Stopping criteria =

=

||

f

−

Ahu

||

2

||

f

||

2

.

(39)

The procedure foriterative formulation is described in Algorithm 5. The smoothing stage (stage- 2) can be applied ns

numberoftimes,dependingonthecomplexityoftheproblemandthecomputationalefficiency. 5. Results

5.1. Convergencetest

Tochecktheconsistencyofthefine-scale3DFEMformulation,asynthetictestcaseformechanicalequilibriumis devel-opedwhichisshowninEquation (40)

Fig. 5. Order of convergence for the displacements computed using Energy norm.

Fig. 6. Synthetictestcase: Schematicofthesynthetictestcasewhenahorizontalshearingloadisappliedonthetopface.Restofthefaceshavestressfree boundaryconditions.TheparametersusedinthistestcasearepresentedinTable2.

Table 2

Parametersusedtomodelcreepinsynthetictestcase.

# Value # Value # Value

dt 1.65 years [day] Q 80000 [66] [J/mol] n 4 [-]

a 5e-24 [-] T 290 [K] E 5 [GPa] [67] m 0.21 ν 0.25 [-] Traction 190 [MPa] Cells 25×25×25 Cr 5×5×5 Dimensions 10×10×10 [m] ux

=

sin



_π

_x L



uy

=

cos



_π

_y W



uz

=

sin



_π

_z H



(40)

Here L

,

W

,

H arethedimensionsoftherockinx

,

y

,

z directionsrespectively.Usingtheforcebalance,thevolumetricforces are computed. Accordingly, these forcesact as righthand side vector forsolving the globalsystemEquation (10). Roller boundary conditions are applied to all the side faces. The bottom face is fixed with no displacement. Elasticmoduli of

E

=

105 _{and the Poisson ratio of}

_ν

₌

₀

_.

_{1 are considered. The error is computed using the 2nd norm where dV is the}

volumeofanelement.Fig.5showsthe2nd-orderaccuracyinspace.Theerroriscomputedby



L2

=











e

(

uanalytical

−

unumerical

)

2dV

.

(41)

5.2. Syntheticsheartestcase

Inthissubsection,asynthetictestcaseismodelledwithcreeptoshowthecapabilitiesofthecomputationalframework andunderstandtheinfluenceofcreeponthebodydeformation.Rocksaltmaterialpropertiesarechosenformodellingthis testcase.Fine-scaleandmultiscalestrategiesareemployedtothismodel.Theparametersusedinthistestcaseareshown inTable2,seealsoFig.6fortheschematicdiagram.Theparametersarechosenartificially,todemonstratethedeformation causedduetocreep.

Ahorizontal shearingforce isapplied onthe topface, thebottomface isfixed andrestofthefaces are tractionfree. Fig. 7b showsthe snapshot of the deformedbody attime

=

0, obtainedfrom both fine-scaleand multiscale strategies. Fig.7bshowsthefinescaleandmultiscalesolutionsobtainedaftertime

=

1.5years.Thisdeformationoccursduetocreep. The error between the FS andMS results ascalculated fromEquation (37) is lessthan 10%. Iterative MS is not further employedinthistestcase.ThesolutionobtainedfromMSstrategyisnotaccuratebecausethereducedboundaryconditions havenot providedanaccuratelocalisationcondition.

Fig. 7. Synthetictestcase: Theabovesnapshotsshowthedeformedbodywhenaforce(x-direction)isappliedonthetopfaceandthebottomfaceis fixed.Fine-scalegridhas25×25×25 cells,andthecoarseningratiois5×5×5.Fig.7ashowstheresultsattime=0forfine-scale(Surface)andMS (Wire-frame).Fig.7bshowstheresultsattime=dt forfine-scale(Surface)andMS(Wire-frame)andwhiteoutlineisthefine-scaleresultattime=0. TheparametersusedinthistestcasearepresentedinTable2.

Table 3

Parametersusedtomodelrocksaltcreepforuni-axialcompressiontestcase.Parametershighlightedwith∗arefoundbyleast-squarefittingofexperimental dataforstrainrate.

# Value # Value # Value

dt 1 [day] Q 80000 [66] [J/mol] n∗ 4 [-]

a∗ 3e-13 [-] T 290 [K] E 5 [GPa] [67]

m∗ 0.21 ν 0.25 [-] [67] Traction 0.1-0.3 [MPa] [68]

Cells 30×30×30 Cr 10×10×10 Dimensions 12×7×7 [cm3]

5.3. Uni-axialcompressionofrocksalt

Uni-axialcompressionofrocksaltisinvestigatedinthissubsection.Adistributedloadisimposedonthetopfaceandall thefourfaces (north,east,westandsouth)ofthedomainaresubjectedtorollerconstraints.Thebottomfaceisfixed.The schematicisshowninFig.9a.AdistributedloadP of0.1-0.3MPaisimposedforaperiodof8months.TheresultsforP

=

0.2MPaarecompared withexperimental results[68].Leastsquare fitisusedonthestrain ratetoobtaintheparameters ofthecreep formulation.Equation (13a) isusedhereasa constitutiveequation.The parametersusedto modelcreepare presentedinTable3.Theexperimentsareconductedoncylindricalspecimens,howevertheeffectofcurvatureisignoredin thiscaseandonlythedimensionsaretaken.

Firstly, classical andrelaxationcreep methodologyis employed tounderstand uni-axialcompression.Fig. 8 showsthe resultsofbothclassicalandrelaxationcreep.Fig.8aandFig.8bshowthevariationofabsolutestress(

σ

zz

)

,absolutestrains

(

ε

t,zz,

ε

el,zzand

ε

cr,zz)withtimeforclassicalcreep.Here,wecanseethatthestressisconstantwithtimewhichmeansthat

elasticstrainisconstantwithtime.Theabsoluteofthetotalstrainisincreasingduetoincreasing creepstrain.Fig.8cand Fig.8d showthe variationofabsolutestress (

σ

zz

)

,absolutestrains(

ε

t,zz,

ε

el,zz and

ε

cr,zz) withtime forrelaxationcreep.

Inthiscase, thetotalstrain isconstant withtime, sincetheabsoluteofthecreep strainis increasingtheelastic strainis reducing withtime. Since the elasticstrain is reducing withtime, theabsolutestress is reducing withtime. Theseplots depictthesamemethodologyasshowninFig.2.

Theresultsare showninFig.9are obtainedfromemployingclassical creep.Alltheresultsareplottedfora nodethat islocatedinthecentre ofthetop face.ThevariationofstrainratewithtimeisshowninFig.9b.Itcanbe seenthatthe strain rateincreaseswhentheimposedloadincreases.SimilarobservationisseeninFig.9c.Least squarefitisconducted on strainrateandaccordinglythe strainis computed.There issome difference betweentheexperimental andnumerical resultsforvariationofstrainwithtimeFig.9c.Thisisbecauseitisassumedthatthepropertiesoftherocksaltareassumed tobehomogeneous.However,Young’smodulusandPoissonratiovarydepending onthesaltcomposition.Thepresenceof saltofvaryingcompositionswouldaffectthefittingparametersofthecreepconstitutiveequation.Investigatingtheeffects ofthecompositionofrocksaltontheoverallbehaviour of thenumericalmodelisbeyondthescopeofthiswork. Fig.9d showsthevariationofdeformationwithtimefordifferentloadingconditionsshownforbothfine-scaleandmultiscale.It canbeseenthattheresultsoffine-scaleandmultiscaleoverlap.Astheimposedloadincreases,theelasticdeformationand deformationaccumulatedovertimeduetocreepincreases.Thisisasituationinwhichtheboundaryconditionimposedon localbasisfunctionsisvery accurate,leading toaccuratemulti-scalesolutions.Notethat,the multiscalesolutionaccuracy dependshighlyontheaccuracyofthelocalbasisfunctions.Whenthereducedboundaryconditionhasprovidedanaccurate localisationcondition,theMSsolutionsareaccurate.

Fig. 8. Theabovefiguresshowthevariationofstrainandstresswithtimewhenclassicalcreep(Figs.8a,8b)andrelaxationcreep(Figs.8c,8d) methodolo-giesareemployedforuniaxialcompressiontestcase.Thefigures(8band8d)showthe variationofelastic,creepandtotalstrainwithtime.Thealgorithms tomodeltheseareshowninAlgorithm1andAlgorithm2respectively.

Fig. 9. Uni-axialcompression: Theaboveplotsshowthevariationofcreepstrainrate(Fig.9b),totalstrain(Fig.9c)anddeformation(Fig.9d)withtimefor rocksaltwhenitundergoesinelasticcreep.TheexperimentalresultsforverticalloadingatthetopofP=0.2MPacase[68] arecompared withnumerical results.ThenumericalresultsareshownfortractionP=0.15to0.25MPa.Fig.9dcomparesthefine-scaleandmultiscaleresultsforcoarseningratioof 1000.ThedetailsofthenumericalmodelarepresentedinTable3.

Fig. 10. Tri-axialcompression: Thefiguresshowthevariationofstrain,strainrate(Fig.10b),deformation(Fig.10c)withtimeforcoalspecimenwhen itundergoescreep.Theexperimentalresults[26] arecomparedwithnumericalresultsusingparametersshowninTable4.Thedeformationresultsare comparedforFS,MSanditerativeMS.TheconvergenceplotobtainedbyusingiterativeMSisshowninFig.10dfordeformation(uz).Theconvergenceplot

isshownforstage2oftheiterativemultiscalesolverfordifferentnumberofsmoothingstagesusingonlyILU(0)ns=1,2,3 orusingonlyGMRESwithout

anypreconditioner.

Table 4

Parametersforcreepformulationofcoalusedtosimulatetri-axialcompressionwithcreepusingBinghammodel.

# Value # Value # Value

dt 60 [s] scr _{44.99 [MPa] [26] Confining stress 30 [MPa] [26]}

G1 1.91 [GPa] [26] G2 79.01 [GPa] [26] E 4 [GPa] [26]

η1 2.6 [GPa h] [26] η2 40.18 [GPa h] [26] ν 0.25 [-] [26]

Cells 20×20×20 Cr 10×10×10 Load 142 [MPa]

5.4. Tri-axialcompression

In thissubsection,influence ofcreep on tri-axialloading is studied.Fig. 10a showstheschematic of tri-axialloading wherea confiningstress isimposed onallthefourfaces (north,east,southandwest),thebottomis fixedandavertical loadisimposedonthetopface.Binghambodyformulationisemployedhere,asinEquation (13b),tostudytheinfluenceof creeponcoalspecimen.TheparameterschosentostudythisformulationarepresentedinTable4.Theexperimentalresults [26] are compared with the numerical results.The experimentis conducted for a period of 45 minutes. The numerical resultsarecomparedforfine-scale,multiscaleanditerativemultiscalesimulations.

Fig.10bshowsthevariationofcreepstrain andcreepstrain ratewithtime.Hereitcanbe seenthattheexperimental datafitswell withthenumericaldatafortheparameters chosen.However, theexperimentaldataalsoshowsthetertiary creep phasewhichisnotcapturedby thenumericalformulation.Tomodeltertiary creep,adifferentdamageconstitutive formulationhastobeemployed thatcanexplain thecouplingofhighlynon-lineareffectswhichleads tomaterialfailure. Thoughitisofcriticalimportance,itisbeyondthescopeofthiswork.Thevariationofdeformation(uz)withtimeisshown

inFig.10cforanodethatislocatedinthecentreofthetopface.Itcanbeseenthatthedeformationincreaseswithtime fortheimposedload.Theprimary andsecondarystageofcreepcanalsobeseeninthisformulation.Thenormalizederror computedfromEquation (37) islessthan2% forthechosencoarseningratio.

Theresultsoffine-scaleandmultiscaleareslightlydifferent.Multiscalesolutionisfurtherimprovedbyemploying itera-tivemultiscalestrategy. ILU(0)andGMRESarecomparedfortheirefficiencyofsmoothingini-MSstrategy. Fig.10dshows thevariationoferrorfordeformation(uz)withnumberofiterationstoachieveconvergence.Thetolerancewassetto1e-10.

Table 5

Parametersforcreepformulationofsandstoneusedtosimulatetri-axialcycliccompression.Parametershighlightedwith∗arefoundbyleast-squarefitting ofexperimentaldataforstrain.

# Value # Value # Value

dt 24 [mins] Q 5000 [24] [cal/mol] Confining stress 25 [MPa] [69]

n∗ 1.6 [-] Dimensions 5×5×10 [cm] E 40 [GPa] [69]

m∗ 1.005 ν 0.32 [-] [69] a∗ 1.5e-19

Cells 30×30×30 Cr 10×10×10 F 265 [MPa] [69]

Fig. 11. Tri-axialcyclicloading: Fig.11ashowsthevariationoftotalstrainwithtimecomparedwiththeexperimentalresults[69] forsandstonerock, andFig.11bshowsthevariationofdeformation(uz)withtimeforbothfinescaleandmultiscaleresults.Theparametersusedforthisformulationare

presentedinTable5.

TheconvergenceplotiscomparedforILU(0)withvaryingsmoothingstagesns

=

1

,

2

,

3 andGMRESwithoutany

precondi-tioner.ItcanbeseenthatGMRESwithoutanypreconditionerconvergestothesolutioninlessthan5iterations,atafaster ratecomparedtoILU(0).Thedeformation obtainedfromMSisnot verydifferentfromFS.However, whentheproperties ofrockareheterogeneous anddifferentloadingconditiontheconvergencerateoftheiterativesolverswillvary.

5.5. Tri-axialcyclicloadingofsandstone

Inthissubsection,a cyclicloadisimposed whenthe rockisconfinedby a stressinall the 4sides. Theexperimental results[69] arecomparedwithnumericalresultsforsandstone.Theparameters chosenforthisformulationare presented inTable5.Investigationoftheeffectofvaryingmaterialpropertiesonthefittingparametersofcreepisbeyondthescope ofthiswork.Thesimulationisconductedforaperiodof210minuteswith4cyclesforavertical loadof265MPa. Equa-tion (13a) isusedbyemployingclassicalcreepinthisstudy.

Fig.11showsthevariationoftotallongitudinalstrainandlongitudinaldeformationforfine-scaleandmultiscale strate-gies. Fig.11acomparestheexperimentalandnumericaltotalstraindata.Leastsquare fitisemployedontheexperimental dataonlyonthefirstcycle.Thefittingparameterm

=

1

.

005 isverycloseto1.Thiswouldmeanthatthestrain rateis al-mostnotdependentontime.Inaddition,itcanbeseenthattheexperimentalandnumericaldataofstrainmatchformore than2cycles.Theobservederrorbetweenthenumericalandexperimentaldataisslightlyhigherinthelastcycleafter300 hours. We thinkthat thismightbe becausedueto cyclicloadingtheremight beopening andclosingofmicro-cracksor expansion ofvoids.Toincorporatethis,the fittingparameters ofthecreep constitutiveformulationwoulddependonthe type ofloadinglikefrequencyandamplitudeofloads.Itisofcriticalimportance tofurtherinvestigatethiseffect.Dueto lackofliterature,thisisbeyondthescopeofthiswork.

Fig. 11bcompares the deformation obtained fromfine-scale andmultiscale strategy. The multiscale solutions are ob-tainedusingmuchfewerdegreesoffreedom.Assuchitclearlyoutperformsthefine-scalesimulationstrategy.Thevariation ofthedeformationofthenodethatislocatedinthecentreofthetopfaceisshownhere.Itcanbeseenthat,the approx-imate solutionsobtainedfrommultiscalesolutions areveryclosetofine-scalesolutions.Hence forthistest case,iterative multiscalestrategyisnotfurtheremployed.

5.6. Influenceofcreeponplainstrainsubsidence

Subsurfaceenergystoragetechnologywillinvolvecyclicinjectionandproductionofthereservoir.Duetothis,subsidence andupliftoftheearths toplayermightoccur.Tocomprehendthis,theinfluenceofcreeponsubsidenceanduplifthasto bestudied.Inthiswork, theinfluenceofnon-lineardeformationonlandsubsidencewhichiscausedduetoinjectionand productionofreservoirismodelled inheterogeneousmediaby assumingplainstrain inthe domain(x

−

z plane). [70,60] gives you the propertiesof the geological formation inside thedomain. The domain mainly consistsof clayedrock. The

Fig. 12. Plainstrainsubsidence: Fig.12ashowstheschematicoftheplainstrainsubsidencetestcasewherethepressurizedH2gasisstoredinthereservoir (smallbox)insidethegeologicaldomain.Thedomainisconstrainedbyrollerboundaryconditionatnorth,east,southandbottom.Fig.12bdepictsthe distributionofYoung’smodulusinthex-zplanewithinthegeologicaldomainasdescribedbyEquation (44).Thegreenboxshowsthelocationofthe reservoir.

reservoir is modelled in the span of 10 km andfor a height of 3 km. The elasticity modulus is obtained fromvertical uni-axialcompressibility(cM)whichisgivenby[71],i.e.,

cM

=

0

.

01241

|

σ

Y

|

−1.1342

.

(42)

HerecM and

σ

_Y areexpressedin[bar−1]and[bar]respectively.Here

σ

_Y istheverticaleffectivestress.Theverticaleffective

stressisobtainedassuperpositionoftotalverticalstress

σ

Y andhydrostaticpressure p,i.e.,

σ

_Y

=

σ

_Y

+

p

=

0

.

12218

|

y

|

1.0766

+

0

.

1

|

y

|.

(43)

TheYoung’smoduluscanbeexpressedasfrom[71],i.e.,

E

=

(

1

−

2

ν

)(

1

+

ν

)

(

1

−

ν

)

cM

.

(44)

TheschematicofthistestcaseisshowninFig.12a.Here,rollersaresubjectedtothebottom,westandeastboundaryand thenorthboundaryistractionfree.Thehydrogengashastobestoredatrighttemperatureandpressureforfeasibility.In thenaturalstateoftheearth,lithostaticpressurewouldbepresentinsidethegeometricaldomain.However,itisassumed thatthelithostaticpressurewouldnotcausenonlineartimedependentdeformation.Adeviationfromtheequilibriumstate, which involves injectionandproduction ofthe reservoir wouldundergononlinear time dependent deformation. Assume that thenaturalstateofthereservoirhasapressureof P0

=

150 barwithoutanysubsidence.Theeffectivepressure P is

givenby

P

=

P

−

P0

.

(45)

Inthiswork,thepressure P isvariedbetween50and110barintheformofcyclicloadasshowninFig.12b.If P

=

0 bar, then P

=

P0whichmeansthattheearthisinequilibriumandthereisnosubsidence.ForP

=

50 bar( P

= −

100 bar)higher

subsidenceisexpectedcomparedto P

=

110 bar( P

= −

40 bar)becauseofthehigherpressuredrop.Thetemperatureis chosen tobe373K[72],withminimummicrobialandbio-chemicalreactions [73].Elasticitymodulus degradationdueto possiblereservoirrockin-elasticityisneglected.Therockmaterialisassumedtobeisotropic.Classicalcreepmethodology isemployedhere.

The fine-scalegrid contains 110

×

110 cells,while themultiscale solveremploys 22

×

22 cells,havingthe coarsening ratioof5

×

5.TheparametersforcreepformulationareobtainedfromprevioustestcaseasshowninTable5.Thetime-step (dt)sizeischosentobe4monthsandtheentiresimulationisconductedforaperiodof10years(30

×

dt).Thetime-step sizeisassumedtobeofsimilarorderasMaxwelltimewhichisgivenbyratioofshearviscositytoitsshearmodulus.The uncertaintyinviscosityduetodifferentgrainsizesisnotconsidered here.Poissonratiois

ν

=

0

.

3 andLame’sparameters areobtainedfromYoungsmodulusinEquation (44).

Fig.14showstheresultsobtainedforplainstrain subsidencetestcasewhennonlineartime-dependentdeformationis incorporated.Fig.13ashowsthevariationofsubsidenceinsidethedomainfordifferenteffectivepressures P

= −

40

,

−

100 barwhentime

=

0.Themultiscalesolutionsareobtainedusingmuchfewerdegreesoffreedom.Thisreducesthe compu-tationaleffortandoutperformsthefine-scalesimulationstrategy.Itcanbeseenthatthemultiscalesolutionandfine-scale solutionhavesimilarmagnitudeexceptforthesmoothkinkwhichisnotobtainedinthemultiscalesolution.Thisisbecause thenumberofcellsoccupiedbythereservoirisclosetothenumberoffine-scaleelementsinacoarseelement.Ifrequired, aniterativeMSstrategycanbeemployedtoimprovethesolution.Overall,thequalityofthemultiscalesolutionisgood.It canbeseenthat thelower istheeffectivepressure,thehigheristhemagnitudeofthesubsidenceduetohigherpressure dropcomparedtothenaturalstate.Theinfluenceofcreepisstudiedbyinjectionandproductioninsidethereservoirwhich variesthepressuresasshowninFig.13bintheformofconstantloadandcyclicloadfor10years.

Fig. 13. Plainstrainsubsidence: TheabovefiguresshowtheresultsforplainstrainsubsidenceusingtheparametersforcreepfromTable5.Thefinescale gridis110×110.Themultiscalesolutionsareshownforacoarseningratioof5×5.Fig.13ashowsthesubsidenceplottedforfine-scaleandmultiscale resultsforP = −40 and −100bar at time=0.

Fig. 14. Plainstrainsubsidence: TheabovefiguresshowtheresultsforplainstrainsubsidenceusingtheparametersforcreepfromTable5.Thefinescale gridis110×110.Themultiscalesolutionsareshownforacoarseningratioof5×5.Constantpressure( P = −100 bar)andcyclicpressure( P = −100 to P = −40 bar)isimposedinsidethereservoirasshownin(Fig.13b)andtheinfluenceofcreeponsubsidenceisshowninFig.14aandFig.14brespectively. Thevariationoftheaverageerror(Equation (37))betweenfine-scaleandmulti-scalesolutionsforux,uzateachtimestepisshown inFig.14candFig.14d.

Fig.14ashowsthesubsidenceinsidethedomainatt

=

0(without creep)andt

=

10years(withcreep)ataconstant load foreffective pressure of P

= −

100 bar. It can be seen that the magnitude of subsidencehas increased att

=

10 years comparedto t

=

0. Duetocreep, thedeformation increaseswhichalso increasesthemagnitudeofthe subsidence caused.The deviationbetweenfine-scaleandmultiscalesolutions hasincreasedwhich isbecausethebasis functionsare constructed based on linear elasticityand the magnitudeof the fictitious creep forces Fcr increaseswith time. Fig. 14b

showsthevariationofsubsidenceinsidethedomainwhenthereservoirisinjectedandproducedcyclically.Themagnitude ofsubsidencecausedonlyduetocreepincyclicloadwillbelesserthanconstantload(ofhighermagnitude).Thisisbecause aconstantloadofhighermagnitudewillresultinhighermagnitudesofstresswhichfurtherresultsintheaccumulationof highermagnitudeofcreepstrainanddeformation.ItcanbeseenthatthemagnitudeoferrorbetweenFSandMSsolutions for P

= −

40 baratt

=

0(asseeninFig.13a)andatt

=

10yearsincyclicloading(asseeninFig.14a)islower.Theerror

6. Conclusion

Inthiswork, aniterativemultiscalenonlinearmodelling strategywasproposedwhichcanmodeltime-dependent non-lineargeomechanicaldeformationusingfiniteelementmethod(FEM).Thedevelopedcomputationalframeworkwouldbeof criticalimportance tomodelsubsurfaceenergystorage whichinvolvesnonlineardeformation. Constitutiveformulaewere chosentomodelthecreepstrainsofdifferentrockmaterialsfromtheavailableliterature.Employingtheleastsquarefiton theexperimentaldata,thefittingparametersofcreepstrainwerecomputed.Classicalcreepandrelaxationcreepwereboth modelledincorporated.Duetothecomputationalstability,an implicittimediscretizationstrategywas followed.Multiscale finiteelement (MSFE)method,formulatedandimplementedalgebraically,was thendevelopedbased onlocallycomputed elastic basis functions. Tofurther improvethe approximatemultiscale solution,2 stage iterativestrategy is developedto combinethemultiscalesystemwithasmootheratfinescale.

Severalnumericalexperimentswereconductedwithsyntheticandrealisticmaterialparameters,tounderstandthe non-lineartime-dependentdeformation.First,creepdeformationofrocksaltunderuni-axialdeformationisstudiedunderboth classical and relaxationcreep. It was found that the fine-scale andmultiscale solutions matchwell. Employing Bingham creepformulation,tri-axialcompressionofcoalisstudied.Thedifferencebetweenfine-scalemultiscalesolutionwasfound tobelessthan2%.Differentiterativemultiscalesmoothingstrategieswerecomparedtoobtainfine-scalesolutionofdesired accuracy.Sincesubsurfacestoragetechnologyinreservoirs involvescyclicloading,assuch tri-axialcyclicloadingon sand-stone was alsostudied.In agreementwiththe previous cases,thedifference betweenfine-scale andmultiscalesolutions was lessthan1%. Theinfluenceofnon-lineardeformation wasfurtherstudied onrealisticreservoirs whichinvolvedland subsidence inheterogeneous fields.Without creep,the multiscalesolutions show a similar magnitudeasfine-scale solu-tions.Forthestudiedrealisticcase,thequalityofmultiscalesolution(withnoiterations)wassatisfactoryevenwhencreep isincorporated.

Thisstudyallows forquantificationoftheinelastictime-dependent creepdeformationforgeoscienceapplications, spe-ciallythoserelevanttosubsurfacestorage.Italsomakesitpossibletorepresentthecomplexnonlinearandheterogeneous systems with much fewer degrees of freedom than the fine-scale classical approaches. As such, it sheds new lights in advancingtheaccuracyandspeedupofthequantitativeanalysesofcomplexnaturalporousmediabehaviour.

Futureworkincludesporo-inelasticsystemstostudytheinfluenceoffluctuatingporepressureonthecyclicbehaviourof therockmaterials.Itwas alsofoundinsubsection5.5thatasthenumberofcyclesincreases, theexperimentaldatastarts todeviatefromthefittingdata.Thisneedstobeinvestigatedfurtherwhenthetime-scaleisinyears.Propagationoffaults andfractureswillbeincorporatedinsidethedomaintostudytheirinfluenceoninelasticdeformationwhenacyclicloadis imposed.

CRediTauthorshipcontributionstatement

KishanRameshKumar: Conceptualization,Methodology,Software,Writing–originaldraft. HadiHajibeygi: Conceptual-ization,Fundingacquisition,Methodology,Supervision,Writing–review&editing.

Declarationofcompetinginterest

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

Acknowledgements

ThefinancialsupportofNWO-TTWViDigrant(projectADMIRE17509)isacknowledged.Theauthorswouldliketothank allthemembersoftheDARSim(DelftAdvancedReservoirSimulation)researchgroupandADMIREprojectmembersatTU Delft,forthefruitfuldiscussions.

References

[1]C.Coutanceau, S. Baranton, T. Audichon,Hydrogen production from water electrolysis,in: HydrogenElectrochemical Production, Elsevier,2018, pp. 17–62.

[2]S.Bauer,A.Dahmke,O.Kolditz,Subsurfaceenergystorage:geologicalstorageofrenewableenergy—capacities,inducedeffectsandimplications, Envi-ron.EarthSci.76(2017)1–4.

[3]M.Budt,D.Wolf,R.Span,J.Yan,Compressedairenergystorage–anoptionformediumtolargescaleelectrical-energystorage,in:CUE2015- Applied EnergySymposiumandSummit2015:LowCarbonCitiesandUrbanEnergySystems,EnergyProc.88(2016)698–702.