# The Norton-Based Double Power Creep Model for Rock Salt

## Full text

(1)
(2)

PLAXISbv, The Netherlands

R.B.J. Brinkgreve

Delft University of Technology & PLAXISbv, The Netherlands

L. Zampich

PLAXISbv, The Netherlands

Windows®is a registered trademark of the Microsoft Corporation.

PLAXIS is a registered trademark of the PLAXIS company (Plaxis bv).

Plaxis bv P.O. Box 572, 2600 AN DELFT, Netherlands

Fax: +31 (0)15 257 3107; E-mail: info@plaxis.nl; Internet site: www.plaxis.nl

These manuals may not be reproduced, in whole or in part, by photo-copy or print or any other means, without written permission from Plaxis bv.

(3)

1 Introduction 3

2 Theory and implementation of the N2PC Rock Salt model 4

2.1 Norton’s creep law 4

2.2 General non-isothermal framework of the N2PC Rock Salt model law 4

2.3 Model Implementation 6

2.4 Summary of the model 7

2.5 Comments on the model and its PLAXIS implementation 8

3 Model Verification 10

3.1 Unconfined uniaxial relaxation test ignoring temperature calculation 11 3.2 Unconfined uniaxial relaxation test at different temperatures 15

3.3 Constant Strain Rate Unconfined 19

4 Case study: a salt cavern for hydrocarbon storage 25

4.1 Model definition 25

4.2 Model calculations 32

4.3 Numerical results 37

5 Conclusions 43

(4)

1 INTRODUCTION

Rock, especially rock salt, exhibits time-dependent behaviour. Creep of rock and rock salt is related to various applications in civil engineering, the mining industry or petroleum engineering, such as energy storages, wellbore stability in oil and gas extraction, radioactive waste facilities design, etc. (Cornet, Dabrowski & Schmida, 2018). For instance, drilling a deep borehole (around 5000m of depth) for oil extraction could face a serious instability problem when going through a creeping salt layer (see Dusseault, Maury, Sanfilippo & Santarelli (2004)).

Most of the underground applications are concerned with large time scales (tens to thousands of years or even longer) and thermal effects at great depth (hundreds to thousands of meters of depth). The Norton power creep law (Norton, 1929) is widely used for describing the steady state creep behaviour, which is relevant for those large time scales. In this context, an enriched version of the Norton law, the Norton-based Double Power Creep (N2PC) model for Salt Rock, including temperature effects, was implemented in PLAXIS software. This model will be presented hereafter. Although the model is particularly adapted for creep of rock salt, it can be used also for other types of rock which exhibit similar creep behaviour.

(5)

2 THEORY AND IMPLEMENTATION OF THE N2PC ROCK SALT MODEL

2.1 NORTON’S CREEP LAW

The Norton creep law (Norton, 1929) describes the creep strain as a function of the stress level as follows:

˙ ε˙νp ε˙νp ενp= A qN 3sss 2q (2.1) whereq = r 3

2s : ss : ss : s the Von-Mises shear stress. A and N are two model parameters. This law is widely used in the literature, especially to describe the creep behaviour of rock salt. The parametersA and N can be easily calibrated by plotting in the logarithm scale the deviatoric stress versus the creep strain rate.

Note that to be consistent, the unit ofA has to be: 1 / time unit / stress unitN

For instance, if N = 3.5 then the unit ofA is: 1 / time unit / MPa3.5

2.2 GENERAL NON-ISOTHERMAL FRAMEWORK OF THE N2PC ROCK SALT MODEL LAW

The N2PC Rock Salt model law is formulated based on the original Norton’s law, see Eq. 2.1. The whole constitutive model is written in terms of the generalised effective stressσσσ000

(accounting for both saturation degree and Biot’s pore pressure coefficient) as defined in Plaxis (2017) (Eq 2.39). This amounts to saying that the constitutive law described below models the relationship between the strains and effective stresses. Note that in the particular case of Salt, the material is often considered non-porous due to its very low porosity (only around 1%) (Fokker, 1994). In the following, the prime symbol is dropped for convenience and the effective stress tensor will be simply denotedσσσ. This effective stress tensor can be decomposed into an isotropic part and a deviatoric part:

σ σ

σ= σmδδδ+ sss (2.2)

whereσmis the mean effective stress,sss the deviatoric stress tensor and δδδ the

second-order identity tensor.

The rock is modelled as an isotropic thermo-elasto-viscoplastic material of which the total strain tensor is decomposed into an elastic part, a thermo-elastic part and a viscoplastic creep part:

ε ε

ε= εεεeee+ εεεTTT+ εεενpνpνp (2.3)

(6)

˙ εeν = 1 K σ˙m (2.5) ˙ ee d ˙ ee d ˙ ee d = 1 2Gs˙s˙s˙ (2.6)

whereε˙ev is the volumetric elastic strain;eedis the deviatoric elastic strain tensor;K and G are the elastic drained bulk and shear moduli, respectively. The overdot notation denotes time derivative of an arbitrary quantity˙r = dr

dt. Remember that the bulk and shear moduli are related by:

K = E

3(1 − 2ν) =

2(1 + ν)

3(1 − 2ν) G (2.7)

The thermo-elastic strains (due to thermal expansion) are already taken into account in PLAXIS 2D together with the Thermal module. They are thus not considered in the constitutive model. For completeness, they are recalled here:

˙ εεε˙˙TTT =      αTx · · · · · · · αTy · · · · · · αTz      ˙ T˙T˙T (2.8)

whereT denotes the temperature, and αTx,αTy andαTzare the anisotropic thermal

expansion coefficients.

The viscoplastic strains are modelled as follows: ˙ εεε˙˙νpνpνp=| ˙ενp|3sss 2q (2.9) whereq = r 3

2s : ss : ss : s the Von-Mises shear stress and | ˙ε

νp| is the intensity of the creep

strain rate, defined by the following double-power law:

| ˙ενp| = A∗1 q q0 N1 + A∗2 q q0 N2 (2.10) whereN1andN2are two stress exponent parameters, whileA∗1andA∗2are two

viscosity-like parameters, which are generally temperature-dependent:

A∗1= A1e −Q1 R 1 T ! (2.11) A∗2= A2e −Q2 R 1 T ! (2.12) In whichA1andA2are two reference viscosity parameters,R = 8.314 J.mol-1K-1is the

universal gas constant andQ1andQ2are two creep activation energy parameters

(J.mol-1). SinceR is a constant independent of material, in practice it is only important to know the values of Q1

R and Q2

(7)

NOTE: Note that in Eq. 2.11 and Eq. 2.12 the temperature T must be the absolute temperature (in the unit Kelvin).

Note that in the standard form of the empirical Norton’s law (Eq. 2.1) the unit of the viscosity-like parameterA is “material dependent” and a bit confusing. In the current formulation Eq. 2.10 the reference stressq0is used to render the stress term under the

exponents q q0

dimensionless, and thereby to get a proper and clearer unit for the viscosity-like parametersA1andA2(in the unit 1/time unit). Depending on the context

(the range of stress value, the calibration process of model parameters in laboratory, etc.),q0may take for instance 1 MPa or 100 kPa (atmospheric pressure). Certainly, this

reference stressq0can take any other positive value, but the model parametersA1and

A2have to be modified correspondingly.

2.3 MODEL IMPLEMENTATION

The general (non-isothermal) version of the N2PC Rock Salt model was implemented in PLAXIS as a User-Defined Soil Model (UDSM). The model can be used for analysing both the general non-isothermal and isothermal cases (without thermal analysis). Some main features of the stress integration scheme used for the N2PC Rock Salt model are described hereafter.

Assuming that all the values at time “n” (denoted by a subscript “n”), as well as the time increment∆tn+1, the total strain increment∆εεεn+1n+1n+1, and the temperature increment∆Tn+1

are known. Those incremental values emanate from the general Finite Element

calculation of the PLAXIS kernel. The role of the constitutive model is to update correctly the effective stressσσσn+1n+1n+1, which can be classically computed as follows:

σn+1 σn+1 σn+1= σσσnnn+ D : h ∆εεεn+1n+1n+1− ∆εεTn+1 T n+1 εTn+1− ∆εεενpn+1νpn+1νpn+1i (2.13) where D is the elastic stiffness tensor, ∆εT

n+1

εTn+1 εT

n+1is the thermal expansion strain increment

and∆εεενpn+1νpn+1νpn+1is the creep strain increment. As aforementioned, the thermal strain∆εT

n+1

εT n+1

εT

n+1is automatically calculated inside the

PLAXIS programs and therefore not needed to be considered at the constitutive level. The N2PC Rock Salt model only aims at evaluating the component∆εεενpn+1νpn+1νpn+1, which can be computed based on a generalised trapezoidal rule:

∆εεεn+1νpνpn+1n+1νp = ∆tn+1 (1 − θ) ˙εε˙νpn νp n ˙ ενpn + θ ˙εεε˙˙νpn+1νpn+1νpn+1 (2.14) where0 ≤ θ ≤ 1 is a stress integration scheme parameter, while ˙εεε˙˙νpn+1νpn+1νpn+1is the creep strain rate evaluated at the unknown (targeted) stress levelσσσn+1n+1n+1. As the creep strain rate is a

non-linear function of stress (Eq. 2.10), the integration Eq. 2.13 is a non-linear problem. The choice of integration scheme depends on two important criteria: stability condition and accuracy condition.

(8)

therefore chosen in the implementation of the N2PC Rock Salt model: θ= 1; ∆εεενpn+1νpn+1νpn+1= ∆tn+1εεε˙˙˙νpn+1νpn+1

νp

n+1 (2.15)

In this case of implicit scheme, asεεε˙˙˙νpn+1νpn+1νpn+1is a function of the unknown stress variable ˙

εεε˙˙n+1νpνpn+1n+1νp (Tn+1, σσσn+1n+1n+1), a Newton Raphson iteration process is used to solve the non-linear Eq.

2.13 to get the solution ofσσσn+1n+1n+1.

For the accuracy condition, it is well-known that for time-dependent issues, the accuracy is time-step dependent. Both implicit and explicit schemes are first-order accurate, and only Crank-Nicholson schemes are second-order accurate. To ensure the accuracy condition, an automatic time-stepping is used, based on the evaluation of the local truncation error inspired by Mahnken & Stein (1989). The global time step (which can be very high) is subdivided into smaller time steps to ensure a small enough local error of creep strain within each step integration.

Note that in general, the time step criterion for accuracy is less restrictive than that for the stability criterion. Therefore even with sub-stepping, an implicit scheme allows using much larger time steps than an explicit scheme. Moreover, the stability condition is generally a property of the system of equation and should be imposed at all the stress points, whereas accuracy based time-stepping is only used for the stress points which exhibit high creep responses.

2.4 SUMMARY OF THE MODEL

2.4.1 GENERAL NON-ISOTHERMAL VERSION

The general version of the N2PC Rock Salt model allows modelling creep responses together with thermal analyses. In this case, it consists of ten parameters (Figure 2.1), including:

1. Elastic shear modulusG (stress unit) and Poisson’s ratio ν (see Eq. 2.5 and Eq. 2.7). 2. Creep stress exponentsN1,N2(see Eq. 2.10). Typical values range from 1 to 8.

3. The viscosity like parametersA1,A2(1/time unit, see Eq. 2.11 and Eq. 2.12).

4. Thermal-related parameters Q1

R , Q2

R (Kelvin unit, see Eq. 2.11 and Eq. 2.12). NOTE: Since Eq. 2.11 and Eq. 2.12 are valid for an absolute pressure only,

Q1

R and Q2

R must always be given in Kelvin. In the PLAXIS project properties, the user may use any unit of Temperature; however, the properties Q/R must always be specified in Kelvin.

5. The reference stress parameterq0(stress unit). As explained earlier, this is not a

material constant, but only a parameter to get a proper unit for the parametersA1,A2.

It is recommended, for simplicity, to chooseq0= 1 unit of stress. However, any other

positive value can be chosen forq0butA1,A2have to be changed correspondingly.

Example: if stresses are expressed in MPa, we may chooseq0= 1 MPa, i.e. the

material propertiesA1,A2are fitted using creep test results for strain rate and the

dimensionless stress q

(9)

Figure 2.1 Model Parameters

but the parameterA1andA2have to be changed toA1/(0.1)N1andA2/(0.1)N2.

6. All the parameters should be positive to ensure mathematically and physically realistic results. If negative values are used, an abort message will be displayed, and the calculation will be stopped.

2.4.2 IGNORING TEMPERATURE VERSION

The model can also be used without thermal analysis. This is where the temperature calculation is ignored (for instance when the temperature change is not significant). In this case, the following points should be considered:

1. The temperature-related parameterQ1

R, Q2

R are not needed (see Figure 2.1). Their values will not affect the calculation results.

2. Only Eq. 2.10 is considered, i.e. the parametersA∗

andA2.

NOTE: Although Eq. 2.11 and Eq. 2.12 are not explicitly considered in this case, it is implicitly understood that A∗

1 andA∗2 are the values evaluated at

an assumed constant temperature T0. The user has to provide the correct

values forA∗1andA∗2at that constant temperature.

(10)

1. Only steady-state creep is considered. This is more or less relevant for the time scale considered in underground engineering applications. Although only steady-state creep is considered at the material level, the model may also describe the “geometrical transient creep” due to stress redistribution around the cavity as reported in previous works (Bérest, Karimi-Jafari & Brouard, 2017). This is demonstrated by the numerical example presented in Section 4.

2. No creep threshold is considered, i.e. creep occurs once deviatoric stresses are non-null. This is more or less consistent with rock salt but also with other rock materials since detecting a stress threshold for creep is not always straightforward in a laboratory.

3. No volumetric creep strain is considered, i.e. creep is isochoric. This is often observed for rock salt.

4. The creep rate is increased with the temperature.

5. The creep law is an extension of the original Norton’s law with two power branches. This provides more flexibilities allowing to better reproduce the rock behaviour (Gunther, Salzer, Popp & Ludeling, 2015). In particular, rock salt exhibits different responses at low and high stress levels. At low stresses, the pressure-solution mechanism is dominant, which can be described by a “linear” power law with the exponent equal to unity (Fokker, 1994). At higher stresses, for instance, higher than 5 MPa, a dislocation mechanism governs the creep behaviour and a higher value of the stress exponent should be used (Bérest, Ghoreychi, Hadj-Hassen & Tijani, 2012) (Gunther, Salzer, Popp & Ludeling, 2015).

6. A robust and efficient implementation was performed with a fully implicit scheme with an automatic local sub-incrementing procedure. This allows using a relatively large time step, and only sub-steps at the points of the structure that are subjected to significant creep loading to satisfy the accuracy condition.

7. This UDSM is developed in and compatible with PLAXIS 3D 2018 and PLAXIS 2D 2019. When using older versions, please set the unit of Temperature to Kelvin exclusively and set Q1

R = Q

2

(11)

3 MODEL VERIFICATION

Some basic testing is presented in this section to verify the model implementation, for different PLAXIS products, including: PLAXIS 2D (2019), PLAXIS 3D (2018) and the SoilTest facility. They are based on the simulation of typical homogenous laboratory tests (with simple geometry).

The unit of time is set to hours, that of temperature is set to◦C. The stress unit is chosen as kPa. The creep parameters were calibrated based on stress results in MPa, therefore the unit reference stress parameter is chosen asq0= 1 Mpa = 1000 kPa.

To run the model, the N2PC64.dll containing the N2PC Rock Salt model is placed in the UDSM sub-folder of the PLAXIS program directory before the program is started. Let’s consider two types of rock salt with the same material properties except for the creep parameters. The other parameters are presented in Table 3.1.

Table 3.1 Parameters for verification examples

Parameters Salt 1A (general non-isothermal analysis) Salt 1B (ignoring temperature analysis) Salt 2 (constant strain rate test)

Shear modulus G 10 GPa

Poisson’s ratioν 0.25 1st creep exponentN1 3.1 2 1st creep parameterA1 2.965 ·10-41/h A∗1= 0.25 · 10−91/h 11.86 ·10-41/h 2nd creep exponentN2 1 -2nd creep parameterA2 1.089 ·10-31/h A∗2= 1 · 10−81/h 0 1st temperature-related term Q1 R 4100 K - 4100 K 2nd temperature-related term Q2 R 3400 K - 3400 K

The parameters for Salt 1 are chosen in agreement with rock salt experimental data. The creep behaviour is typically described by two branches. The first branch represents dislocation creep and is modelled by a power law with a typical value of the exponent N1= 3.1. The second represents pressure-solution creep (dominant at low stress-level)

and can be modelled by a linear law, i.e.N2= 1. To illustrate different applications of the

model implemented in PLAXIS, two material sets are considered based on the same material properties of Salt 1. The Salt 1A includes all the material properties of Salt 1, which can be used for general non-isothermal analysis. The second set for Salt 1B corresponds to the creep parameters of Salt 1 at20◦C which is used when temperature

calculation is ignored. The temperature related parameters Q1

R and Q2

R are not needed. Moreover, as explained before,A∗1andA∗2are needed instead ofA1andA2. In this case,

they may be estimated at20◦C based on Eq. 2.11 and Eq. 2.12 using the data ofA1,A2,

Q1

R and Q2

(12)

facilitate the comparison against the analytical solution for constant strain rate testing, where the closed-form solution is only obtained for the single power law withN1= 2. This

is explained in detail in Section 3.3.

Note that, to facilitate comparisons, the gravity effect is ignored by setting the unit weight equal to zero. Furthermore, all the materials are drained and no water pressure effect is taken into account in these tests.

3.1 UNCONFINED UNIAXIAL RELAXATION TEST IGNORING TEMPERATURE CALCULATION

The uniaxial relaxation test is simulated without considering temperature analysis, i.e. for Salt 1B. This option can be used for both PLAXIS 2D, PLAXIS 3D and SoilTest facility.

3.1.1 GENERAL DESCRIPTION AND CLOSED FORM SOLUTION The sample is subject to an initial elastic loadingεini

1 corresponding to an initial stress

valueσ1ini = −qini where “1” indicates the loading axis. The configuration of the uniaxial test allows to write the relationship between the applied axial stressq = −σ1and the

axial strain as follows: ˙

σ1= E ˙ε1− ˙ενp1



(3.1) whereE = 2G(1 + ν) is the Young’s modulus.

1 = −0.04% 0 hours

Relaxation phase ε˙1= 0 30.000 hours

The relaxation conditionε˙1= 0 and the creep constitutive equations Eq. 2.9 and Eq. 2.10

allow to rewrite Eq. 3.1:

˙ q + E " A∗1 qN1 0 ! qN1+ A ∗ 2 qN2 0 ! qN2 # = 0 (3.2)

In the case whereN2= 1, some mathematical techniques (see for instance Polyanin &

Zaitsev (2017)) allow solving the non-linear ordinary differential Eq. 3.2 analytically to get the evolution of the axial deviatoric stress:

q =    qini1−N1+A ∗ 1 A∗2q 1−N1 0  e −(1 − N1)EA ∗ 2t q0 −A ∗ 1 A∗2q 1−N1 0   1 1 − N1 (3.3)

(13)

3.1.2 SOILTEST SIMULATION

Figure 3.1 Modelling of unconfined relaxation test by SoilTest facility of PLAXIS

The SoilTest simulation is simply based on stress calculation without any geometry effect. The relaxation test configuration can be defined in the General tab sheet. Two Steps were used for the instantaneous loading phase, while only ten time steps to simulate the relaxation phase, which lasts for a very long time (30.000 hours, namely around years). The typical result of the relaxation test is often expressed by plotting a stress - time graph. This can be done by Right click→ Add custom chart in the SoilTest user interface (Figure 3.1).

3.1.3 PLAXIS 2D SIMULATION

The same test was simulated using PLAXIS 2D. An axisymmetric 2D model was built to represent the test condition. The model is 1 m by 1 m in the XY plane. The boundary conditions are: fixities at the left side (symmetry axis of the sample) and at the base of the model, and an imposed displacement equal to0.4 · 10−3m (corresponding to εini

1 = 0.04%) at the top to represent the loading condition (Table 3.2). Only two elements

(14)

Figure 3.2 Modelling of unconfined relaxation test by PLAXIS 2D: Axisymetric model mesh and boundary conditions

Figure 3.3 Modelling of unconfined relaxation test by PLAXIS 2D: Calculation Type Ignore temperature Thermal calculation type (Figure 3.3). The other options keep their default values.

3.1.4 PLAXIS 3D SIMULATION

PLAXIS 3D can also be used for simulating the triaxial stress state of the test. Considering the homogenous stress state characteristic of the test, a cubic sample is modelled. The boundary conditions and mesh are shown in Figure 3.4. A surface displacement equal to0.4 · 10−3m is applied at the upper boundary of the model. The base of the model is vertically fixed. Due to symmetry, the left and back sides are normally fixed (i.e. fixed in the direction perpendicular to their plane), while the other boundaries are free due to the unconfined condition of the test.

(15)

Figure 3.4 Modelling of unconfined relaxation test by PLAXIS 3D: Model mesh and boundary conditions

(16)

3.1.5 RESULTS COMPARISON

The typical results for the stress-time relationship in the relaxation test obtained from different simulations are compared (Figure 3.6). It is readily seen that all the simulations provide accurate results compared with the analytical solution (Eq. 3.3).

Figure 3.6 Stress relaxation curve in Unconfined relaxation test: Comparisons of simulation results by SoilTest, PLAXIS 2D and PLAXIS 3D versus Analytical solution

Moreover, only a small number of time steps (around ten steps) are needed to get this accurate result. All the time-stepping processes are automatic and very practical. This illustrates the robustness of the implementation.

3.2 UNCONFINED UNIAXIAL RELAXATION TEST AT DIFFERENT TEMPERATURES Now to illustrate the thermal effect included in the model, the same test configuration is repeated accounting for different (constant) temperatures. Simulations are made by means of PLAXIS 2D, which is equipped with the Thermal module.

In this case, the closed form solution Eq. 3.3 is still valid, but Eq. 2.11 and Eq. 2.12 are also involved in enabling the temperature dependency.

3.2.1 PLAXIS 2D SIMULATIONS WITH THERMAL EFFECT

The parameter set Salt 1A (Table 3.1) is used to simulate the relaxation tests at20◦C

(first test) and80◦C (second test). The full data ofA

1,A2,Q1

R and Q2

R are needed. A constant temperature is generated during the initial phase via Earth gradient mode imposed in Model condition (Figure 3.7).

(17)

Figure 3.7 Earth gradient imposed in the initial phase of the test

(18)

Figure 3.9 Temperature boundary condition of80◦C applied during Phase 4 (2nd test).

Figure 3.10 Calculation type for the 2nd relaxation test at80◦C

(19)
(20)

3.2.2 RESULTS COMPARISON

Again, the typical results for the stress-time relationship in the relaxation test obtained from different simulations are compared (Figure 3.12). It is readily seen that all the simulations provide accurate results compared with the analytical solution Eq. 3.3 together with condition Eq. 2.11 and Eq. 2.12. An increase of temperature leads to an increase of the relaxation rate. Moreover, the simulation result using Salt 1A at20◦C coincides with that obtained previously using Salt 1B as expected (because Salt1B corresponds indeed to the reference temperature of20◦C).

Figure 3.12 Stress relaxation curve in unconfined relaxation tests at20◦C and80◦C

3.3 CONSTANT STRAIN RATE UNCONFINED

3.3.1 GENERAL DESCRIPTION AND CLOSED FORM SOLUTION

The constant strain rate test on rock salt sample (at constant temperature) is now simulated. This test is commonly used to study the long-term time-dependent behaviour of rocks (Hudson, 1993). The test configuration is similar to the above unconfined relaxation test. However, the sample is now subject to a constant strain rate

˙

ε1= −Vdef <0 (negative sign for compression) during the whole phase to reach the

(21)

instead ofε˙1= 0. The governing ordinary differential Eq. 3.2 becomes: ˙ q + E " A∗1 qN1 0 ! qN1+ A ∗ 2 qN2 0 ! qN2 # = EVdef (3.4)

This non-homogenous non-linear equation is not easy to solve for arbitrary values of N1andN2. This normally requires advanced techniques which may provide closed-form

but cumbersome results. The general solving is out of the scope of this report, which aims at verifying our model implementation. Therefore, we only take into account the case where only one power-component is active by settingA∗

2= 0, and the active

component has the power exponentN1= 2. The simplified equation can be rewritten as

follows: dq q2− q0 s Vdef A∗1 !2 = −E  A∗1 q20  dt (3.5)

The integration technique provides the solution, assuming that the initial valueqini = 0:

q = q0 s Vdef A∗1   1 − 2 1 + e 2Et q0 √ VA∗ 1    (3.6)

Note thatA∗1is temperature-dependent via Eq. 2.11 and Eq. 2.12.

3.3.2 PLAXIS 2D SIMULATIONS WITH THERMAL EFFECT

The test can be simulated without considering temperature dependency using both PLAXIS 2D, PLAXIS 3D and SoilTest facility. Nonetheless, to clarify the temperature effect, PLAXIS 2D simulations with temperature option were conducted. The Salt 2 material parameters (Table 3.1) are used (see Figure 3.13). Four tests are simulated for two different strain rates at two different temperatures of20◦C and60C (Table 3.3). The

desired total strain level to be reached during the tests is 1 %. To this end, a

displacement of -0.01 m is imposed in the upper boundary of the sample (Figure 3.14). Table 3.3 Testing condition for four constant strain rate compression tests

Strain rate Temperature

Test 1 Vdef= 10−3hour−1= 2.78 · 10−7s−1 20◦C

Test 2 Vdef= 10−4hour−1= 2.78 · 10−8s−1 20◦C

Test 3 Vdef= 10−3hour−1= 2.78 · 10−7s−1 60◦C

Test 4 Vdef= 10−4hour−1= 2.78 · 10−8s−1 60◦C

Apart from the steady-state thermal flow option presented in Section 3.2, other methods may be used to generate a constant temperature in a cluster. In this example, the temperature of20◦C is generated on Test 2 by assigning a uniform temperature profile in Thermal Condition then using earth gradient type together with Plastic calculation (Figure 3.15).

(22)

Figure 3.13 Model Parameters of Salt 2 for constant strain rate tests

linear variation. The time interval of each test is thus deduced corresponding to each constant strain rate. For example, the test with a strain rate of10−4hour−1lasts for 100 hours to reach the desired total strain level (Figure 3.15).

(23)

Figure 3.15 Calculation option setting for Tests 1 and 2

Figure 3.16 Thermal condition to generate a uniform temperature of60◦C in the sample Another way to apply a temperature to a cluster is to set directly the temperature

condition of a soil entity. This way is used to impose the temperature of60◦C on our tests

(24)
(25)

3.3.3 RESULTS COMPARISON

The comparison between the PLAXIS results and analytical solutions Eq. 3.6 together with Eq. 2.11 and Eq. 2.12 are carried out and shown in Figure 3.18. A perfect agreement between the analytical and numerical results is obtained. Both results illustrate the physical processes considered in the N2PC Rock Salt model creep effect is small for fast loading, while non-linear deformation may develop during a slow loading, and these non-linear strains are enhanced with a higher temperature.

(26)

4 CASE STUDY: A SALT CAVERN FOR HYDROCARBON STORAGE

In the following, a case study on a typical deep Compressed Air Energy Storage (CAES) will be examined. An underground cavern for gas storage, excavated in a deep salt layer, is simulated during its operation state.

In general, CAES caverns are created by solution mining. Soft water is injected through a cemented well into a salt layer. Rock salt dissolves in water to create a cavern filled with brine, which will be replaced by compressed air used in energy production. Seasonal exploitations of the stored gas consist of gas withdrawal and injection, relating to depressurisation/repressurisation as well as heating/cooling cycles (Bottcher, Gorke, Kolditz & Nagel, 2017).

Modelling the whole life cycle of the cavern involves complex coupled phenomena such as rock dissolution (normally from the bottom to the top of the cavern), heat exchanges between gas, brine and rock mass, thermodynamic effects of fluids and mechanical behaviour of the rock mass (creep) (Bérest, 2011).

In the following, the behaviour of a salt cavern subject to a depressurisation process is simulated. The solution mining, which can be considered relatively short compared to the operation period of the storage, is not modelled. Instead, it is only represented by an initial stress state, assuming that during that short time, stress redistribution due to creep does not occur significantly. A very common assumption that rock salt is non-porous is followed and only thermo-mechanical coupling is taken into account (Zhu, Shen, Arson & Pouya, 2017), (Bérest, Ghoreychi, Hadj-Hassen & Tijani, 2012), (Fokker, 1994) and (Bottcher, Gorke, Kolditz & Nagel, 2017).

4.1 MODEL DEFINITION

The cavern height is 100 m and its diameter is 100 m. It is located at the depth ranging from 1000 m to 1100 m in an Asse salt formation Figure 4.1. Due to the symmetry of the model and the cavern geometry, an axisymmetric model of 500 m by 500 m is used in PLAXIS 2D (corresponding to a depth from 800 to 1300 m), which is large enough compared to the cavern geometry to avoid boundary effects. Instead of modelling complex behaviour of gas stored in the cavern, gas temperature and gas pressure conditions on the cavern wall as simplified boundary conditions are applied directly. In the PLAXIS 2D model, the cavern is located between y = 200 m and 300 m. The left

boundary (symmetry axis), the base and the right far field boundary are fixed mechanically and closed for heat flow. The life cycle of the cavern includes (Serbin,

´

Slizowski, Urba ´nczyk & Nagy, 2015): • Excavation by salt leaching

• Gas injection at high temperature and high pressure (pressurisation) • Storing of gas at maximum pressure

• Gas withdrawing for energy use (depressurisation) • Idling phase: storing of gas at the minimum pressure Except for excavation, the other phases are repeated seasonally.

(27)

Y X Surface Upper layers (Rocks, Soils, ...) 0 m -800 m -1100 m -1300 m -1000 m 500 m 500 m 0y = 16 MPa Pg CavernTg

Figure 4.1 Deep Salt Cavern model geometry and boundary conditions sections:

• Line 1: (0,0) to (50,15) • Line 2: (50,15) to (50,80)

• Arc: start angle45◦, radius 70 m, segment angle35

Two subsections are also used to get a closed geometry of the cavern as in Figure 4.2. At the cavern wall, a temperature boundary condition and a normal force perpendicular to the wall are set in the tunnel properties as in Figure 4.3.

A general fine mesh is used due to complex phenomena involved. However, the mesh is refined around the cavern wall and coarsened at far field by setting a maximum

coarseness factor of 8 (see Figure 4.4). Moreover, the effect of all the upper layers can be simulated by a line load or by using a fictive material layer which represents an equivalent weight. In this example, we use a 1 m thick layer (see Figure 4.4) composed of a fictive non-porous very stiff material with a unit weight of 1600 kN/m3(see Figure 4.5) to generate a vertical stress boundary condition of 16 MPa at the upper boundary of the model (Figure 4.1). In Figure 4.4, some points are selected for later analysis, including nodes A (0,300), B (0,200), C (50,250.2) and stress point N (50.08, 270.3).

(28)

Figure 4.2 Tunnel designer tool to create the cavern geometry

(29)

Figure 4.4 Model mesh and boundary conditions Table 4.1 Parameters for verification examples

Parameters Asse Rock Salt

Shear modulus G 8.846 GPa

Poisson’s ratioν 0.3

Unit weight 20.4 kN/m3

Void ratio 0.02

Specific heat capacity 858.6 ·103J/t/K

Thermal conductivity 5.4 W/m/K

Linear thermal expansion 2.8 ·10-51/K

1st Creep exponentN1 2.7

1st Creep parameterA1 0.05 1/day

2nd Creep exponentN2 6.8

2nd Creep parameterA2 0.02 1/day

1st Temperature-related term Q1

R 4800

2nd Temperature-related term Q2

(30)

Figure 4.5 Fictive load layer geometry and properties

(31)
(32)
(33)

4.2 MODEL CALCULATIONS

The two-year simulation is carried out in the following phases. During each year, a cycle of pressurisation, storing, depressurisation and idling is fully modelled. The input loading values are based on the realistic data from Serbin, ´Slizowski, Urba ´nczyk & Nagy (2015).

4.2.1 INITIAL PHASE

During the initial phase, the initial stresses and temperatures are generated, using the K 0 procedure and Earth gradient modes.

In rock mechanics, the ratio between horizontal and vertical in situ stresses may depend on the considered depth. Following Brown & Hoek (1978)K 0 = 1.2 is used,

corresponding to the depth range (from 800 to 1300 m) in the present study. This value is set in the Rock Salt material Initial tab. The valueK 0 = 1.2 means that the horizontal tectonic stress is higher than the vertical stress, which is often met in deep rock formations.

(34)

4.2.2 INITIAL OPERATION STATE PHASE

As aforementioned, the mining process is not modelled. Instead, the mechanical and thermal steady state after excavation is modelled using Plastic with Steady-state thermal flow modes. During this phase, a uniform temperature of 310 K and a gas pressure of 8 MPa is applied to the cavern wall during the storage of the gas.

(35)

4.2.3 PRESSURISATION PHASE Y1 (FIRST YEAR)

Starting from the initial state after excavation, gas is then injected intensively into the cavern for 20 days. This phase simulates the gas injection process. Pressurisation means that the pressure is increased from 8 MPa to 12 MPa and the temperature is also increased from 310 K to 325 K during 20 days, using a linear time-dependent function Heating 1 in PLAXIS 2D as shown in Figure 4.9. A transient calculation with Fully Coupled mode is therefore used.

(36)

4.2.4 STORING PHASE

After the gas is injected, it is stored at the constant pressure and temperature during the remainder of the year (205 days). A fully coupled thermo-mechanical transient analysis is performed (see Figure 4.10).

(37)

4.2.5 DEPRESSURISATION PHASE Y1

The gas is withdrawn during 20 days for energy production. The decrease of gas pressure also leads to a temperature decrease due to the thermodynamic relationship of gas. A fully coupled thermo-mechanical transient analysis is performed (see Figure 4.11), where the temperature is reduced from 325 K to 300 K using a Temperature function Cooling (linear) and the pressure is reduced to 5 MPa.

(38)

4.2.6 IDLING PHASE Y1

The gas is then stored at its minimum pressure of 5 MPa. The cavern is kept idle during the remainder of the year (120 days). A transient fully coupled analysis is therefore performed for this time interval.

This cycle is then repeated for year two (phases denoted with Y2). The only difference is at the 2nd pressurisation phase, where the temperature boundary condition at the cavern wall is modelled by a linear function Heating2 other than the temperature function Heating1 used in phase Pressurisation Y1 (Figure 4.12). This is to model the temperature increase from 300 K (previous idling phase) to 325 K, i.e. interval of 25 K, instead of from 310 K (initial state) to 325 K as in the first pressurisation phase.

Figure 4.12 Temperature boundary condition for Pressurisation Y2 Phase

4.3 NUMERICAL RESULTS

(39)

4.3.1 CONVERGENCE ANALYSIS

It is important to estimate the convergence, mainly the volume loss of the storage due to salt creep. Figure 4.13 shows the radial displacement at the right cavern wall, at point C as indicated in Figure 4.4. It relates to the change of the cavern radius. It is seen that the cavern radius is getting smaller due to the creep effect. The convergence rate is highest once the cavern is kept idle (idling phases Y1 and Y2). This is logical because the gas is kept at minimum pressure and the deviatoric stress is smallest during those phases. Figure 4.14 shows the vertical displacements at the roof and floor of the cavern. It is observed that the floor moves upward, while the roof moves downward, which points out the creep convergence effect. Again, the convergence rate is highest at idling phases. It is also noted that the creep convergence exhibits a transient behaviour although the constitutive model only takes the steady-state creep behaviour into account. This is consistent with the “geometrical transient creep” due to stress redistribution around the cavity as reported in previous work Bérest, Karimi-Jafari & Brouard (2017).

After two years, the cavern height lost 30cm and its radius lost around 15cm. Rough estimations lead to a loss of 0.9% of the cavern volume. The numerical results are at the same order of magnitude as reported in (Serbin, ´Slizowski, Urba ´nczyk & Nagy, 2015).

Figure 4.13 Time evolution of the radial displacement at the right side of the cavern (point C)

4.3.2 TENSILE STRESS ANALYSIS

Apart from convergence problems, the tensile stresses are also a critical issue. This may lead to fracturing of the rock mass and thereby a permeability enhancement which induces gas leakage.

The radial (σxx), axial (σyy) and tangential (σzz) stresses at the point N (see Figure 4.4) in

(40)

Figure 4.14 Time evolution of the vertical displacement at the roof and floor of the cavern (point A and B)

follows the same stress path directly induced by boundary conditions at this wall. During pressurisation, the tangential stress is increasing (less compressive). During

depressurisation, the tangential stress is however not decreasing. This may be explained by the cooling effect (see Figure 4.16) where the temperature drop results in a tensile stress near the cavern wall. Note that in Figure 4.16 the highest temperature is

somewhere inside the rock mass but not at far field because the previous storing phase does not reach the steady state yet. During the creep phases (idling or storing phases), the stress state tends to isotropic (equilibrium).

The cooling effect seen in Figure 4.15 is widely reported in the literature (Bérest, Ghoreychi, Hadj-Hassen & Tijani, 2012), (Bottcher, Gorke, Kolditz & Nagel, 2017), (Bérest, 2011) and (Serbin, ´Slizowski, Urba ´nczyk & Nagy, 2015). At this stress point N, the tangential and axial stresses may be tensile and fracturing may develop at this point. For further investigation, the tangential stress profile at the cross section at the same elevation of this point (270.3 m) is plotted in Figure 4.17. A small tensile zone of around 1m is observed in the vicinity of the storage wall at the end of the 2nd depressurisation phase.

(41)

Figure 4.15 Stresses at point N at the vicinity of the cavern wall

(42)
(43)
(44)

5 CONCLUSIONS

A constitutive model for creep of rocks, especially rock salt, has been implemented in PLAXIS as a UDSM.

(45)

6 REFERENCES

 Bérest, P. (2011). Thermomechanical aspects of high frequency cycling in salt storage caverns.

 Bérest, P., Ghoreychi, M., Hadj-Hassen, F., Tijani, M. (2012). Mechanical Behaviour of Salt VII. Taylor & Francis, Paris.

 Bérest, P., Karimi-Jafari, M., Brouard, B. (2017). Geometrical versus rheological transient creep closure in a salt cavern. Comptes Rendus Mecanique, 345, 735–741.  Bottcher, N., Gorke, U., Kolditz, O., Nagel, T. (2017). Thermo-mechanical

investigation of salt caverns for short-term hydrogen storage. Environmental Earth Science, 76–98.

 Brown, E., Hoek, E. (1978). Trends in relationships between measured in-situ stresses and depth. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 15, 211–215.  Cornet, J., Dabrowski, M., Schmida, D. (2018). Long term creep closure of salt

cavities. International Journal of Rock Mechanics and Mining Sciences, 113, 96–106.  Dusseault, M.B., Maury, V., Sanfilippo, F., Santarelli, F.J. (2004). Drilling around salt:

stresses, risks, uncertainties.

 Fokker, P. (1994). The behaviour of salt and salt caverns. Dissertation, Delft University of Technology.

 Gunther, R., Salzer, K., Popp, T., Ludeling, C. (2015). Steady-state creep of rock salt: Improved approaches for lab determination and modelling. Rock Mech Rock Eng.  Hudson, J. (1993). Comprehensive Rock Engineering: Rock Testing and Site

Characterization. Pergamon.

 Mahnken, R., Stein, E. (1989). Adaptive time step control in creep analysis. International journal for numerical methods in Engineering, 28, 1619–1633.  Norton, F.H. (1929). Creep of Steel at High Temperatures. McGraw-Hill Book

Company.

 Plaxis (2017). PLAXIS Scientific Manual.

 Polyanin, A., Zaitsev, V. (2017). Handbook of Ordinary Differential Equations: Exact Solutions, Methods, and Problems. CRC Press/Chapman and Hall.

 Serbin, K., ´Slizowski, J., Urba ´nczyk, K., Nagy, S. (2015). The influence of

thermodynamic effects on gas storage cavern convergence. International Journal of Rock Mechanics & Mining Sciences, 79, 166–171.

 Simo, J., Hughes, T. (1997). Computational Inelasticity. s.n.  Zhu, C., Shen, X., Arson, C., Pouya, A. (2017). Numerical study of

Updating...

## References

Related subjects :
{"ft_lang_label":"__label__en","ft_lang_prob":0.8480448,"math_prob":0.9364587,"size":40875,"snap":"2021-43-2021-49","text_gpt3_token_len":11043,"char_repetition_ratio":0.13897385,"word_repetition_ratio":0.043550834,"special_character_ratio":0.24822019,"punctuation_ratio":0.13544102,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.96897143,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2021-10-28T15:55:14Z\",\"WARC-Record-ID\":\"<urn:uuid:dcd43453-a811-45a8-a15d-c5d0b5f88e23>\",\"Content-Length\":\"291232\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:12bb4e93-bf4b-4761-a095-2beb4a20a8a7>\",\"WARC-Concurrent-To\":\"<urn:uuid:2f810902-677d-44c2-9167-786da27c365e>\",\"WARC-IP-Address\":\"164.90.139.124\",\"WARC-Target-URI\":\"https://1library.net/document/q05r6g0g-norton-based-double-power-creep-model-rock-salt.html\",\"WARC-Payload-Digest\":\"sha1:HHN5QTNP6FTD7LD6VXKUYV4DDLHJKOHU\",\"WARC-Block-Digest\":\"sha1:W72NJWW23MH2CJERYSXKSICIL6OZAJWV\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2021/CC-MAIN-2021-43/CC-MAIN-2021-43_segments_1634323588341.58_warc_CC-MAIN-20211028131628-20211028161628-00245.warc.gz\"}"}