Rock mass response for lined rock caverns subjected to high internal gas pressure

Davi Rodrigues Damasceno,Johan Spross,Fredrik Johansson

Division of Soil and Rock Mechanics,KTH Royal Institute of Technology,SE-100 44,Stockholm,Sweden

Keywords: Lined rock caverns (LRCs)High pressure Rock mass response In situ stress condition Cavern shape Gas storage

ABSTRACT The storage of hydrogen gas in underground lined rock caverns(LRCs)enables the implementation of the first fossil-free steelmaking process to meet the large demand for crude steel.Predicting the response of rock mass is important to ensure that gas leakage due to rupture of the steel lining does not occur.Analytical and numerical models can be used to estimate the rock mass response to high internal pressure;however,the fitness of these models under different in situ stress conditions and cavern shapes has not been studied.In this paper,the suitability of analytical and numerical models to estimate the maximum cavern wall tangential strain under high internal pressure is studied.The analytical model is derived in detail and finite element (FE) models considering both two-dimensional (2D) and threedimensional (3D) geometries are presented.These models are verified with field measurements from the LRC in Skallen,southwestern Sweden.The analytical model is inexpensive to implement and gives good results for isotropic in situ stress conditions and large cavern heights.For the case of anisotropic horizontal in situ stresses,as the conditions in Skallen,the 3D FE model is the best approach.

1.Introduction

The development of modern society demands an increase in the global production of crude steel;however,conventional blast furnace processes cannot meet this production demand within the international climate agreements.To bring the steelmaking process in line with these agreements,the Hydrogen Breakthrough Ironmaking Technology (HYBRIT) project,i.e.a joint venture between major mining,steelmaking,and energy companies in Sweden,aims to implement the first steelmaking process that is completely fossil-free.With this new technology,renewably produced hydrogen gas will replace fossil fuels used in the blast furnace and other parts of the process.Large quantities of hydrogen gas must be readily available in a storage,which raises concerns regarding safety and cost.As discussed in a technical report by Sofregaz US Inc.(1999) and in a prestudy report by Johansson et al.(2018),underground storages using the lined rock cavern (LRC) concept should be a viable option.Alternatives for the storage of compressed hydrogen gas are gas holders,spherical pressure vessels,and pipe storages,as presented by Andersson and Gr?nkvist(2019);however,the maximum allowed pressures for such alternatives are much lower than that of LRCs.

In Skallen,southwestern Sweden,an LRC for the storage of pressurized natural gas was successfully constructed and put into commercial operation in 2002 (Glamheden and Curtis,2006;Tengborg et al.,2014).Small-scale laboratory experiments have also been conducted to predict the response of such rock caverns(Tunsakul et al.,2013,2017;Jongpradist et al.,2015).Stille et al.(1994) and Johansson (2003) proposed design principles for the LRC concept.The integrity of the storage is sustained by the confining pressure of the surrounding rock mass,while the steel lining maintains gas-tightness.Fig.1 shows the main lining components for an LRC.The LRC for the HYBRIT project will be constructed in Sweden,where high in situ stresses and hard jointed rock masses are typical geological conditions.

Fig.1.Sketch of the LRC concept,highlighting the main lining components.

For the design principles of the LRC,Johansson (2003) showed that excessive straining of the steel lining is avoided by placing the cavern at sufficient depth,which leaves fatigue of the steel lining to be in the critical failure mode.The steel-concrete lining is thin compared to the size of LRC;Johansson (2003) showed that the rock mass,as a consequence,takes essentially all deformations of a large LRC,while the lining’s displacement only follows that of the rock cavern wall.The strain range in the rock cavern wall may therefore be used for the fatigue design of the steel lining.Mansson et al.(2006)’s observations of field tests from Skallen showed that the gas temperature variation of the natural gas inside the LRC would reduce the deformation range at the cavern wall within pressure cycles.This is due to the resulting deformations from thermal expansion and contraction at the cavern wall which are in the opposite directions compared to that of the deformations induced from the variation in the storage internal pressure.Notably,such temperature variation reduces the amount of injectable and withdrawable gas volume,thus using a gas circulating cooling/heating system can be favorable to decrease the temperature variation.

An analytical model to estimate the maximum strain in the rock cavern wall caused by gas pressurization was presented in a Swedish report (Johansson et al.,1995),without detailed derivations.The main simplifications of this analytical model were that the in situ stresses are isotropic and the cavern is long with spherical cupolas,which results in the maximum strain being the tangential strain located at the cavern mid-height (horizontal cross-section).This analytical model is computationally inexpensive and relatively easy to implement.Despite the inherent simplifications associated with analytical models,these models are of importance for the conceptual understanding of the rock mass behavior.For pressurized rock caverns,however,the present literature lacks such analytical models.

In order to study tunnels and caverns with complex geometries and boundary conditions,numerical models are often used (Jing and Hudson,2002;Jing,2003) and such models have also been used to estimate the response of LRCs.Glamheden and Curtis(2006) used finite difference models with two-dimensional (2D)simplifications of plane strain and axisymmetry to simulate the excavation convergence of the Skallen LRC horizontal(mid-height)and vertical cross-sections,respectively.Furthermore,Lu (1998)used a three-dimensional (3D) finite element (FE) model to simulate field pressure tests in the pilot LRCs in Gr?ngesberg,Sweden,and obtained good agreement with the expected response.

Rock caverns can also be used for compressed air energy storage(CAES).Perazzelli and Anagnostou (2016) investigated different failure modes for such storages under similar operational conditions as that of the LRC.CAES is,however,often constructed in nonsilo shape and operated using lower cavern gas pressures(e.g.Kim et al.,2012,2013;Rutqvist et al.,2012;Park et al.,2013;Carranza-Torres et al.,2017;Zhou et al.,2017,2020).Nevertheless,to ensure the structural safety of LRC designs,there is a need to further study the rock cavern response under different in situ stress conditions and cavern shapes.

In this paper,an analytical model for the tangential strain in the rock cavern wall due to the internal gas pressure is derived in detail,based on the work by Johansson et al.(1995).The solutions from analytical and FE models are verified against field measurements from the LRC in Skallen.The maximum tangential strain in the cavern wall is studied for the presented models,considering different conditions of in situ stress anisotropy and cavern shapes.Geological conditions of good-quality rock mass and high cavern gas pressure are considered.

2.Analytical model for the rock mass response

In this section,the rock mass response caused by cavern gas pressurization is derived and the equations from Johansson et al.(1995) are updated to a more general and concise formulation.(Another analytical solution specifically for shallow cavities excavated in cohesive ground was presented by Carranza-Torres et al.(2017)).All these equations work similarly to the response curve of the rock mass in the convergence-confinement method (Hoek and Brown,1980;Brown et al.,1983).However,due to the pressurization at the cavern wall,the rock mass is displaced outward from the cavern center.

2.1.Governing equations

Positive sign convention is used for the direction of compression,where the cavern gas pressure,pi,is applied.The gas temperature increase during gas pressurization decreases the deformations at the rock cavern wall;therefore,for a more conservative design,this effect is not considered in this work.It is further assumed that the rock mass properties are homogenous,and that the horizontal in situ stresses (σo) are isotropic.If the cavern geometry is symmetric and long in the axial direction and the in situ stress does not change with depth,a plane-strain assumption can be used for 2D simplification and cylindrical coordinates,as shown in Fig.2.The maximum magnitude of the cavern wall tangential strain,εθi,is expected at the cavern midheight under the condition that the cupolas are spherical,the cavern is positioned at a sufficient depth,and the lining interaction is efficient.These assumptions avoid a concentration of strains caused by sharp edges in the geometry,uplift of the overlying rock mass,and concrete cracking.The paragraphs below use a small wedge element located at a radial distance (r) from the cavern center for the mathematical description,as shown in Fig.2c.This wedge element is subjected to radial and tangential stresses σrand σθ,respectively,which results in a radial displacementu.

To derive the analytical solution for εθi,the governing equations for stress equilibrium,strain-displacement compatibility,and constitutive relationship,are needed.Force equilibrium in the radial direction for the wedge element in Fig.2c,results in

Considering both drand dθ to be small,Eq.(1)can be simplified into the stress equilibrium:

Based on the undeformed and deformed shapes of the wedge element in Fig.2c,strain-displacement compatibility is set by

Fig.2.Location of internal pressure,pi,and internal radius,ri,at the cavern for both(a)the vertical cross-section and(b)the mid-height horizontal cross-section.Internal pressure and stresses are represented by red arrows.(c)A small wedge element located in the rock mass at a radial distance r from the cavern center is subjected to stresses and deformation.The radial and tangential stresses on the wedge element are shown in the undeformed shape and result in the deformed shape (in blue).

where εrand εθare the total radial and tangential strains,respectively.

An elastoplastic constitutive relationship is used to describe the rock mass response.Assuming small strains,the total strains εrand εθcan be decomposed into their elastic and plastic strain components as

where εreand εθe are the radial and tangential elastic strains,respectively;and εrpand εθpare the radial and tangential plastic strains,respectively.These elastic and plastic strains depend on the principal stresses.For the case of cavern wall pressurization (pi>σo),i.e.the opposite stress path direction compared to the rock mass response curve in the convergence-confinement method,the maximum principal stress,σ1,and minimum principal stress,σ3,are equal to σrand σθ,respectively.

The elastic behavior is hypothesized to be linear.With the plane-strain assumption(i.e.zero strain in the axial direction),the stress-strain relationship is defined as

whereMmis the constrained modulus of the rock mass;and λmis Lamé’s first parameter for the rock mass,which can be written as functions of the Young’s modulus of the rock mass,Em,and Poisson’s ratio of the rock mass,vm,such thatMm=Em(1-vm)/[(1+vm)(1-2vm)] and λm=Emvm/[(1+vm)(1 -2vm)].

In terms of plastic behavior,εrpand εθpare defined by the flow rule,which is set by

where λpis the plastic multiplier andqpis the plastic potential for non-associated flow.The expression forqpis derived based on criteria for tensile yield,Fyt,and shear yield,Fys.In this paper,Fytis represented by the Rankine failure criterion andFysby the Mohr-Coulomb failure criterion,given by

whereftmandfcmare the tensile strength of the rock mass and unconfined compressive strength of the rock mass,respectively.The parameterfcmis estimated asfcm=am(km-1),whereamandkmare the attraction of the rock mass and passive earth pressure coefficient of the rock mass,respectively.The parametersamandkmare functions of the cohesion of rock mass,cm,and its friction angle,φm,i.e.am=cm/tan φmandkm=tan2(45+φm/2).

Based on the governing equations in Eqs.2-12,three response types were derived for the rock mass:compressive elastic behavior,tensile yield (crack opening) behavior and shear yield behavior.Fig.3 presents the conceptual model.

Fig.3 shows yielding paths I and II with Mohr circles at yielding initiations and possible rock mass responses.For yielding path I,the rock mass deforms elastically at lowpiand yields by tension and shear with increasingpi.If the rock mass is weak,no tensile yield takes place following yielding path II.The parameters σrtand σrsare the radial stresses at the tensile and shear yield boundaries,respectively,located at radiirtandrs,from the cavern center.The radial displacements in each deformation zone are derived in the following Sections 2.2-2.4 and a general approach to estimate εθiis summarized in Section 2.5.

2.2.Compressive elastic behavior

Fig.3.Conceptual model showing(a)two possible stress paths and corresponding Mohr circles at yielding initiation,and(b)development of rock mass responses for yielding path I.For yielding path II,tensile yield does not take place.

As Fig.3b shows,compressive elastic behavior takes place in rock volumes located betweenrtandr→+∞,wherertcan be larger thanrsfor yielding path I and equal torsfor yielding path II.To describe the rock mass under compressive elastic behavior,the strain-displacement compatibility Eqs.(3) and (4) are applied in the stress-strain relationship of Eqs.(7) and (8),resulting in

whereucis the radial displacement in the compressive elastic zone from all deformations,i.e.from virgin stress,excavation,and cavern gas pressurization.Inserting σrand σθinto the stress equilibrium of Eq.(2) results in a homogeneous second-order differential equation:

The general solution is given by

whereC1andC2are the constants of integration.Insertingucinto Eq.(13) and applying the boundary conditionsandresult in the stresses in the compressive elastic zone:

The specific solution for the differential equation is then obtained by

Only radial displacement from cavern gas pressurization is of interest,whileucalso accounts for deformations from both virgin stress and excavation.Conditions of high-quality rock mass are generally recommended for the LRC concept,thus the deformation caused by excavation is typically elastic.Therefore,dilatancy of the rock mass volume may occur only during cavern gas pressurization.Using Eq.(19)withgives the deformations from virgin stress and excavation and allows the radial displacement prior to cavern pressurization,u0,to be determined:

The radial displacement caused by cavern gas pressurization,,can then be obtained by subtracting=uc-u0,resulting in

2.3.Tensile yield behavior

Tensile yield behavior is caused by the opening of rock cracks in the zone betweenrsandrt(Fig.3b).The radial strength against tensile yieldingfbt,as Fig.3a shows,is calculated at the boundary between the compressive elastic and tensile yield zones(i.e.atr=rt) by inserting the tangential stress σθ(Eq.(18)) into the yield criterionFyt(Eq.(11)),resulting in

Based on Eq.(22) and assuming that the residual tensile strength of the rock mass is negligible,the magnitude of the residual radial stress atrtis given by σrt=2σo.

The tangential stress in the tensile yield zone is calculated using a negligible (zero) residual tensile strength inFyt(Eq.(11)),resulting in

Setting up a force equilibrium in the radial direction for the tensile yield zone,where σθ=0 and,yields a radial stress of

The radius of the tensile yield boundary is found by reformulating Eq.(24)with,giving

To derive the radial displacement,the plastic potential based on Eq.(23),i.e.qp=σθ,is inserted in the flow rule in Eqs.(9)and(10).Note that plastic strain takes place only in the tangential direction;therefore,the total radial strain in the tensile yield zone is equal to the radial elastic strain alone,which is found by rearranging Eqs.(7)and (8) as

Applying the strain-displacement compatibility Eq.(3) and inserting the stresses σr(Eq.(24))and(Eq.(18))result in the following expression:

whereutis the radial displacement in the tensile yield zone from all deformations,i.e.from virgin stress,excavation,and cavern gas pressurization.The formulation of Eq.(27) using σθ=2σo-σrtinstead of inserting σθ=0 ensures that the differentiation of the resulting radial deformation,i.e.for the calculation of the radial strain Eq.(3),has continuity between the tensile yield and compressive elastic behaviors.Note that the condition σθ=0 from Eq.(23) holds when σrt=2σoin the tensile yield zone.This is important in order to concisely formulate the final analytical solution in Section 2.5.Integrating Eq.(27)while setting the condition that(from Eq.(19))at the zone boundaryrtgives

In a manner similar to(Eq.(21)) in the compressive elastic zone,radial displacement in the tensile yield zone from cavern gas pressure is obtained by subtracting=ut-u0:

2.4.Shear yield behavior

The shear yield zone may develop through two different yielding paths(Fig.3a).From the compressive elastic behavior,the lowest magnitude out offbt(Eq.(22))and the corresponding radial strength against shear yieldingfbsdecide which yield behavior occurs first,i.e.yielding by tension followed by shear(yielding path I) or only yielding by shear (yielding path II).Mathematically,yielding path I can be formulated such that iffbt≤fbs,then tensile yielding takes place whenfbtfcm.Meanwhile,yielding path II occurs iffbsfbs.

If shear yielding takes place after tensile yielding(yielding path I),the confining stress at the boundary between the tensile yield and shear yield zones is σθ=0.From the yield criterionFys(Eq.(12)),it can then be noted that the radial stress at the shear yield boundary is equal tofcm.The residual radial stress atrsis defined by the residual parameters,i.e.σrs=fcmr=amr(kmr-1),whereamr=cmr/tan φmrandkmr=tan2(45+φmr/2),in whichfcmris the residual unconfined compressive strength of the rock mass,amris the residual attraction of the rock mass,kmris the residual passive earth pressure coefficient of the rock mass,cmris the residual cohesion of the rock mass,and φmris the residual friction angle of the rock mass.

Alternatively,shear yielding may initiate without any previous tensile yielding (yielding path II),which is modeled by setting the tensile yield zone to zero thickness,i.e.rt=rs.The upper limit for compressive elastic behavior is then given by the radial strength at the shear yield boundary,fbs.The expression forfbsis derived by inserting the stresses σr(Eq.(17)) and σθ(Eq.(18)) into the yield criterionFys(Eq.(12)) at the boundary between the compressive elastic and shear yield zonesrs,which gives

If the shear yield zone develops,the residual radial stress at the boundaryrsis given by Eq.(30),using the residual parameter,i.e.σrs=(fcmr+2kmrσo)/(kmr+1).

The tangential stress in the shear yield zone is determined by isolating σθinFys(Eq.(12))with the residual parameters,resulting in

To derive the radial stress,σθis then inserted into the stress equilibrium Eq.(2),resulting in the first-order differential equation:

which has the general solution:

The constant of integrationC1is determined using=pias the boundary condition at the cavern wall,giving

The radial distance to the outer boundary of the shear-yielded rock mass is calculated using Eq.(34) with=σrs,giving

In order to derive the radial displacement in the shear yield zone,the plastic potentialqpmust be set usingFys(Eq.(12)) with the dilatancy parameter,i.e.qp=σr-aψ(kψ -1)-kψσθ=0,in whichaψandkψare the attraction and passive earth pressure coefficient expressed using the dilatancy angle of the rock mass(ψm)in place of the friction angle,respectively.Usingqpin the flow rule(Eq.(9)and(10))results in the following plastic strain relationship:

The parameterkψin Eq.(36) expresses the incremental ratio between plastic strains as

Following Eqs.(5) and (6),the plastic strains in Eq.(36) can be reformulated in terms of total strains and elastic strains,where the elastic strains are calculated at the shear yield boundaryrs.Based on this and applying the strain-displacement compatibility in Eqs.(3)and(4)for the total strains at the shear yield zone,a first-order differential equation is obtained:

whereusis the radial displacement in the shear yield zone;and εrersand εθers are the radial and tangential elastic strains at the shear yield boundary,respectively.The total radial strain atr=rsis purely elastic for both yielding paths;therefore,εrersis calculated using the strain-displacement compatibility of Eq.(3) in the expression for(Eq.(29)) giving

On the other hand,the total tangential strain atr=rsis elastoplastic for yielding path I since tensile yielding takes place.To obtain an analytical solution for Eq.(38),only the elastic part of the tangential strain atr=rsis used for εθers.The tangential elastic strain εθersis approximated by first subtractingu｛r=rt｝c-u0(Eq.(19)and(20))and then applying the strain-displacement compatibility Eq.(4),withr=rs,resulting in

The elastic strains εrersand εθers are independent ofr,and therefore the general solution for Eq.(38) is

Using the condition=-εθrsrs,where εθrsis the total tangential strain at the shear yield boundary,in Eq.(41) results in the specific solution:

The εθrsis expressed by inserting(Eq.(29)) in the straindisplacement compatibility (Eq.(4)),usingr=rs,which gives

Since Eq.(42) is a function of strains which have already been corrected for virgin stress and excavation,the calculatedusis the displacement caused by cavern gas pressurization alone.

2.5.Calculation of cavern wall tangential strain

Note specifically that the derivations presented here have been organized carefully so that Eq.(42)can,in fact,be used to calculate the cavern wall radial displacement including both yielding paths(Fig.3a) and allpi.The expression for the sought εθiis obtained using the strain-displacement compatibility,Eq.(4),inus(Eq.(42))withr=ri,which results in

In summary,to calculate εθiusing Eq.(44),first the radial stresses σrtand σrsmust be defined for all possible responses,as shown in Table 1,following the conceptual model in Fig.3.The yielding sequence depends on the magnitudes offbt(Eq.(22)) andfbs(Eq.(30)),and the responses depend,additionally,onpi,fcmandfcmr.Then,the radiirt(Eq.(25))andrs(Eq.(35))can be calculated,followed bykψ (Eq.(37)),εrers(Eq.(39)),εθers (Eq.(40))and εθrs (Eq.(43)) which are all inputs to obtain εθi.For conditions of hard rock masses,the analytical model(Eq.(44))has a good agreement with the magnitudes of εθi calculated from a 2D FE model.

Table 1 Definition of radial stresses σrt and σrs for possible rock mass responses.

3.FE models for the rock mass response

To investigate the effect of anisotropic in situ stresses and cavern shapes,numerical models were constructed using FE analysis,via the commercial software package Abaqus?,which allowed εθi to be estimated.Both 2D and 3D FE models were used in this paper,as shown in Fig.4.The lateral dimensions of the models were set to 1000 m to avoid boundary effects.For the 2D FE model,plane strain assumption was used at the mid-height horizontal cross-section of the LRC with anisotropic horizontal in situ stresses.The 3D FE model included the height of the LRC and assumed perfectly spherical cupolas.The vertical in situ stress,σv,increases with depth and both maximum and minimum horizontal in situ stresses are defined by multiplying stress ratios to σv.Roller and pinned boundary conditions were used in the FE models as indicated in Fig.4.The resulting state of stress after the excavation of the cavity was imported into the FE models before the pressurization of the cavern wall.The material model for the rock mass may fail in tension and/or shear following the Rankine and Mohr-Coulomb failure criteria,respectively.

4.Verification of models using field measurements

For the development of the LRC concept,a large-scale demonstration storage was commissioned in Skallen,southwestern Sweden.The rock mass in the area is dominated by good-quality gneiss,i.e.GSI≈75,as presented by Glamheden and Curtis (2006).The demonstration LRC was constructed between 1998 and 2002.Extensive testing was performed to validate the LRC concept for the storage of natural gas.This storage was then used for commercial operation and merged to the Swedish gas grid.

Fig.4.Meshed geometry with boundary conditions and loads used in analysis for (a) 2D FE and (b) 3D FE models.

Mansson et al.(2006)presented the measured average εθifrom mini-extensometers(M1-M36),which were placed at mid-height in the concrete layer around the cavern wall of the LRC in Skallen.These data were used to verify the analytical and FE models presented in Sections 2 and 3.

Geometrical specifications for the models were obtained from Mansson et al.(2006).The characterization of rock mass properties for the geological conditions of this LRC was provided by Glamheden and Curtis (2006) and used here as input parameters for simulation.The shear yield limit,i.e.fbt≤fbsandfcm=50?3 MPa,was not exceeded for the consideredpi=25 MPa and the rock mass yielded only by tension.(In case shear yield would take place,residual strength properties could for example be estimated from Hoek and Brown(1997)or Cai et al.(2007)and the dilatancy angle from Alejano and Alonso(2005)).The expressions for σHand σh(Table 2) were presented by Glamheden and Curtis (2006) and these horizontal in situ stresses are estimated at the cavern midheight section,i.e.z=140?5 m corresponding to the location of the mini-extensometers.For the analytical solution the average horizontal in situ stress σo=(σH+σh)/2 atz=140?5 m is used.A summary of the input parameters for the analytical and FE models is shown in Table 2.

Table 2 Summary of input parameter used for the analytical and FE models.

The deformation behavior of the rock cavern depends on the stresses that develop around the storage during gas pressurization.In Fig.5,the distributions of σrand σθat the directions of principal horizontal in situ stresses are shown for the analytical,2D FE and 3D FE models.Note that the analytical model can only account for isotropic in situ stresses while anisotropic in situ stresses are used in the FE models;the results are therefore expected to deviate.The magnitudes of the radial stresses at the cavern wallr/ri=1 are equal to the applied gas pressure(σr=pi)and decrease rapidly in magnitude toward the respective in situ stresses.The tangential stress σθdecreases in magnitude toward the cavern wall (r/ri=1),and σθ=0 if tensile yield takes place.

Fig.6 shows,for increasingpi,a comparison between the average εθi from mini-extensometers M1-M36 and the corresponding magnitudes calculated from the analytical and FE models.The shaded areas represent the strain ranges between the minimum and maximum εθi for the FE models.As the largest magnitudes of εθitend to concentrate at only a small area of the cavern wall(in the σHdirection),the average εθihas a magnitude closer to the minimum εθi.Good agreement was observed between the field measurements and the 3D FE model.The field data fall between the maximum and average εθiestimated from the 3D FE model (blueshaded area).The analytical εθi has a good agreement with the average εθiestimations from the 2D FE model.

The 3D FE model was capable of estimating more realistic εθi magnitudes than the 2D model;however,the average tangential strain is slightly underestimated by the 3D model compared to the measured values.In that case,a more conservative approach would be investigating the maximum magnitude of εθi estimated from the 3D model.The effects of in situ stress anisotropy and geometry of the cavern were the main simplifications used in the 2D model and the influence of these assumptions is further discussed in Section 5.

5.Influences of in situ stress condition and cavern shape

In addition to the in situ stress condition and cavern shape,the location of the maximum strain concentration at the cavern wall depends on whether uplift of the overburden occurs.Therefore,to allow model comparison of tangential strains at the cavern midheight,a minimum required rock cover was defined for the 3D FE model which makes the influence from uplift negligible.The value of the required rock cover was defined to be the depth where the maximum cavern wall tangential strain,max|εθi|,is located at the mid-height cross-section for the relevant average horizontal in situ stress,σo.The influences of in situ stress condition and cavern shape was then studied for the estimation of max|εθi| using the analytical and FE models.

Fig.5.Calculated magnitudes of σr and σθ in the directions of (a) σH and (b) σh at distance from the gas storage for analytical,2D FE and 3D FE models.

Fig.6.Comparison between field data(mini-extensometers M1-M36) and analytical,2D FE and 3D FE models for tangential strains at the cavern wall.The field data plot between the minimum and maximum tangential strains of the 3D FE model,which are indicated by the blue shaded area.

5.1.Definition of minimum rock cover

In order to investigate how the rock cover,hrc,influences the location of strain concentrations around the 3D cavern wall,the max|εθi| was studied under different in situ stress conditions.For this purpose,σovaries between σv,2σvand 3σv,assuming isotropic horizontal in situ stresses(σH=σh).Sincepi,riand the rock mass properties are fixed in the analysis,hrcand σoare the only parameters that influence the resistance against uplift.Therefore,max|εθi| was analyzed with respect to these two parameters as shown in Fig.7,to determine the minimum requiredhrcfor the 3D FE model.At shallow depth,strain concentrates at the upper cupola for σo=2σvand σo=3σv.For σo=σv,the max|εθi|takes place at the mid-height of the cavern,even at shallow depth.The robust depth for the LRC,following the guidelines by Johansson(2003),for the three curves was found athrc=216 m;therefore,this value is used for the minimum rock cover in the following analysis.In the real case,shallower depths could be allowed depending on the design requirements and prevailing geological conditions.

5.2.Comparisons between analytical and FE models

As a reference for this analysis,the steel lining resistance to fatigue is estimated to be ΔεR≈6‰,based on the Eurocode 3(CEN,2005),assuming 600 filling-emptying cycles and a welding quality category of 80 MPa.The analyzed in situ stress conditions were defined by assuming different degrees of anisotropy between average horizontal and vertical in situ stresses,i.e.σo=σv;2σv;3σv.Furthermore,anisotropy between maximum and minimum horizontal in situ stresses was also considered,letting σH=σh;2σh;3σh.This anisotropic in situ stress condition was obtained by rearranging the expression σo=(σH+σh)/2 into

The influence of cavern shape is studied forhc/36=1;2;3,wherehcis the cavern height and 36(m)is the fixed diameter of the cavern,thushc/36=1 represents a perfectly spherical shape.Consequently,the maximum cavern height for the simulations is max(hc)=3×36=108 m.To avoid concentration of strains in the cupolas of the 3D FE model,the cavern mid-height depth,zm,was fixed for all simulations atzm=hrc+max(hc)/2=270 m,wherehrc=216 m from Section 5.1.The vertical in situ stress at this depth had a fixed value of=7 MPa.

The comparison between analytical and FE models with respect to in situ stress condition and cavern shape is shown in Fig.8.For these simulations,no shear yielding of the rock mass was observed.To allow comparison with the FE models,the analytical solutions which considers only the case σh=σH,are also shown for the cases σH=2σhand 3σh.

From Fig.8,the following observations can be made:

(1) The magnitude of max|εθi| decreases with increasing σo(from left to right columns of Fig.8) and increases with increasing in situ stress anisotropy in the horizontal direction(from top to bottom rows of Fig.8).

(2) For the cases of σH=σh(first row of Fig.8),the results from the analytical and 2D FE simulations match well.

(3) Ifhc/36 approaches three and σH=σh(first row of Fig.8),the results from the 3D FE simulation match well with that from the analytical and 2D FE simulations,which consider the plane strain assumption.For σH≥2σh,the 2D FE simulation overestimates the max|εθi| compared to the 3D FE simulation,even if the cavern’s height in the FE model is large(hc/36=3).

Fig.7.Study of required rock cover and average in situ stress for the 3D FE model to find the hrc that avoids uplift of the overlying rock mass.The locations of max|εθi|around the rock cavern wall are indicated by the crosses for hrc=72 m and 216 m.

(4) With increasing horizontal in situ stress anisotropy,the result difference between 2D and 3D FE simulations becomes greater.

(5) At higher horizontal in situ stress anisotropies (σH≥ 2σh),the values of max|εθi| from the 3D FE simulation are mostly larger than that from the analytical simulation (except for those close to the spherical shape),which is based on the plane strain assumption.

The results from this analysis show that the estimations of max|εθi|for the analytical and FE simulations are sensitive to the in situ stress condition and cavern shape.The difference between analytical,2D and 3D FE simulations is more pronounced for anisotropic horizontal in situ stress conditions.For the cases of isotropic horizontal in situ stresses,the strain magnitudes are relatively small and the max|εθi| can be well estimated using the analytical solution under the assumed geological conditions.As indicated by the top row of Fig.8,this is the most favorable stress condition for an LRC.

6.Conclusions

This paper compared analytical and FE models for the rock mass response to high internal gas pressure in LRCs.The analytical model was derived and presented in detail to aid in the conceptual understanding of the rock mass behavior for pressurized caverns.Analytical and FE simulation results were verified using field measurements from the LRC in Skallen,southwestern Sweden.The suitability of different models was studied,with respect to the influences of in situ stress condition and cavern shape on the estimation of maximum cavern wall tangential strain.The following conclusions can be drawn:

Fig.8.Estimations of max|εθi|for analytical and FE models at different in situ stress conditions and cavern shapes.The magnitude of σo is shown at each column and the magnitudes of σh and σH are shown inside each individual graph.*The analytical model results from σH=σh are displayed also in the plots for σH=2σh and 3σh as a reference,even though this model cannot account for anisotropic horizontal in situ stresses.

(1) The analytical model allows computationally inexpensive evaluations of the cavern wall tangential strain under high internal pressure.This model gives reliable results for conditions of isotropic in situ stresses and large cavern heights.For a LRC in similar conditions to that of the Skallen storage,i.e.good-quality rock mass (GSI≈75),anisotropic in situ stresses and close to spherical cavern,the analytical model overestimates the measured average tangential strains by about 40%.

(2) The field measurements of the average tangential strains in the cavern wall of the Skallen storage fall within the strain range estimated by the 3D FE model,when considering anisotropic in situ stresses.

(3) For isotropic in situ stress conditions,the analytical and 2D FE models should yield the same results,for good-quality rock mass conditions.However,for poor-quality rock mass conditions,the analytical results underestimate the strain compared with the 2D FE simulation.This is caused by a necessary simplification in the derivation of the analytical model.

A limitation of this study was that a large rock cover had to be chosen to ensure maximum strain concentration at the cavern midheight in order to allow comparison of the rock cavern response between 2D and 3D geometries.Spherical shape of the cavern’s lower cupola was also used in the 3D geometry which disregarded strain concentrations due to irregular shapes.If hydrogen gas is to be stored,the effect of hydrogen embrittlement in the steel lining must be evaluated and the resistance decreased accordingly.Furthermore,this paper only studied the large-scale behavior of the rock cavern,and,for the design of LRCs,any additional strain concentration caused by the opening of rock joints due to cavern pressurization must also be evaluated.

To enable the storage of large quantities of hydrogen gas and facilitate fossil-free steel production,further studies related to LRC are critical as the safety of such storages is of the utmost importance.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work has been conducted as part of the HYBRIT research project RP-1.This research was financially supported by the Swedish Energy Agency (Grant No.42684-2).Thanks to Prof.H?kan Stille and Mr.Jan Johansson for their support with the LRC concept.

List of symbols

amPeak attraction of the rock mass

amrResidual attraction of the rock mass

aψDilatancy attraction of the rock mass

cmPeak cohesion of the rock mass

cmrResidual cohesion of the rock mass

C1First constant of integration

C2Second constant of integration

EmYoung’s modulus of the rock mass

fbsRadial strength at the shear yield boundary

fbtRadial strength at the tensile yield boundary

fcmPeak unconfined compressive strength of the rock mass

fcmrResidual unconfined compressive strength of the rock mass

ftmTensile strength of the rock mass

FytTensile yield criterion

FysShear yield criterion

gGravitational acceleration

hcHeight of cavern

hrcHeight of rock cover

kmPeak passive earth pressure coefficient of the rock mass

kmrResidual passive earth pressure coefficient of the rock mass

kψDilatancy passive earth pressure coefficient of the rock mass

MmConstrained modulus of the rock mass

pMean stress for the stress path diagram

piCavern gas pressure

qRadius of Mohr circle for the stress path diagram

qpPlastic potential for non-associated flow

rRadial distance from cavern center

riCavern wall radius

rsRadial distances from cavern center to the shear yield boundary

rtRadial distances from cavern center to the tensile boundary

uRadial displacement

u0Radial displacement prior to cavern gas pressurization

ucRadial displacement in the compressive elastic zone from all deformations

u′cRadial displacement in the compressive elastic zone caused by cavern gas pressure

usRadial displacement in the shear yield zone

utRadial displacement in the tensile yield zone from all deformations

u′tRadial displacement in the tensile yield zone caused by cavern gas pressure

vmPoisson’s ratio of the rock mass

zDepth from the ground surface

zmDepth at cavern mid-height

εrTotal radial strain

εreRadial elastic strain

εrersRadial elastic strain at the shear yield boundary

εrpRadial plastic strain

εθ Total tangential strain

εθeTangential elastic strain

εθersTangential elastic strain at the shear yield boundary

εθiCavern wall tangential strain

εθpTangential plastic strain

εθrsTangential strain at the shear yield boundary

θ Section angle from the cavern center

λmLamé’s first parameter of the rock mass

λpPlastic multiplier

ρmDensity of the rock mass

σ1Maximum principal stress

σ3Minimum principal stress

σhMinimum horizontal in situ stress

σHMaximum horizontal in situ stress

σoAverage horizontal in situ stress

σrRadial stress

σrsRadial stress at the shear yield boundary

σrtRadial stress at the tensile yield boundary

σvVertical in situ stress

σθTangential stress

φmPeak friction angle of the rock mass

φmrResidual friction angle of the rock mass

ψmDilatancy angle of the rock mass