# Elastoplastic Modelling of an In Situ Concrete Spalling Experiment Using the Ottosen Failure Criterion

Authors: Lauri Kalle Tapio Uotinen; Topias Kalle Aleksi Siren
Journal: Journal of Engineering
 (2017)

---

## Abstract

An in situ concrete spalling experiment will be carried out in the ONKALO rock characterization facility. The purpose is to establish the failure strength of a thin concrete liner on prestressed rock surface, when the stress states in both rock and concrete are increased by heating. A cylindrical hole 1.5 m in diameter and 7.2 m in depth is reinforced with a 40 mm thin concrete liner from level −3 m down. Eight 6 m long 4 kW electrical heaters are installed around the hole 1 m away. The experiment setup is described and results from predictive numerical modelling are shown. Elastoplastic modelling using the Ottosen failure criterion predicts damage initiation on week 5 and the concrete ultimate strain limit of 0.0035 is exceeded on week 10. The support pressure generated by the liner is 3.2 MPa and the tangential stress of rock is reduced by −33%. In 2D fracture mechanical simulations, the support pressure is 3 MPa and small localized damage occurs after week 3 and damage process slowly continues during week 9 of the heating period. In conclusion, external heating is a potent way of inducing damage and thin concrete liner significantly reduces the amount of damage.


---

## Body

## 1. Introduction

In the spent nuclear fuel disposal using the KBS-3 concept [1] the disposal holes and tunnels will be filled with swelling bentonite clay to provide support pressure to the deposition hole and tunnel surface. The suppressed rock damage directly decreases the number of transport routes for radionuclides and has a positive effect when assessing the safety of nuclear waste disposal. The support pressure also has significant effect when assessing rock reinforcement capacity. Another use case for the support pressure of thin concrete liners is the active prevention of stress induced damage such as spalling, rockburst, or strainburst. In such cases, the support pressure mobilizes the internal friction and inhibits dilatation associated with the rock mass damage onset.Shotcrete (sprayed concrete) or concrete lining with steel wire mesh are the reinforcement options selected to protect the shafts leading to the Olkiluoto spent nuclear fuel repository in Finland. Sprayed concrete is considered for the main tunnels in the deposition panels. Decay heat from the canister panels will first reach the main tunnels in deposition level and then the access tunnel and finally the shafts located in the technical area (Figure1). The three shafts are ventilated and temperature-controlled. The shafts are over 400 meters deep and will house important technical connections and will be used to transport personnel, materials, and possibly spent nuclear fuel canisters. They are difficult to maintenance without interrupting operation. The rock mass surrounding the shafts is loaded by the excavation shape concentrated regional in situ stress and changes in the stress state (e.g., thermally induced or excavation induced) may initiate spalling. In such events, support pressure generated by the concrete liner may reduce or prevent the damage or contain the spalled rock pieces [2].Figure 1
Layout of ONKALO with POSE niche, access tunnel, the shafts, and the technical area.The failure strength of sprayed concrete or concrete liner when subjected to the pressure of the surrounding heated rock is unknown. The liner forms a thin cylinder shell, which will be subject to almost hydrostatic stress. Together with the external rock mass, the liner forms a composite structure. The elastic moduli of rock mass and concrete are significantly different. When both materials are subjected to the same external stress the strains will be different and debonding or slipping may occur. This may damage the concrete-rock adhesion that holds the liner to the rock wall.To study the damage process of the concrete liner an in situ concrete spalling experiment (ICSE) is planned. Preliminary numerical modelling is needed to plan the experiment and the required instrumentation. Further numerical modelling is required to carry out a prediction-outcome study. The intent of the in situ experiment is to define the failure strength of concrete on a prestressed rock surface, when the stress state of both rock and concrete is increased by heating.The experiment will be carried out in the ONKALO rock characterization facility in Finland. The location is the POSE niche at a depth of −345 m at chainage 3650 (Figure1) where Posiva’s Olkiluoto Spalling Experiment was previously conducted. The third experiment hole (ONK-EH3) located at the back of the niche will be used (Figure 2). In the preceding third phase of POSE [3], the ONK-EH3 was heated from inside to study thermally induced damage caused by the excess heat, as a simulation of the heat from a spent nuclear fuel canister inside a canister hole. No significant damage was observed [3] and the hole is to be reused in the ICSE and heated to a higher temperature until failure occurs. Early scoping results for the ICSE predicted that eight 4 kW heaters could produce damage both in sprayed concrete and in surrounding rock mass within 9 weeks of heating period [4].Figure 2
The end of the POSE niche (a) and the main dimensions of the experiment (b).
(a) (b)In ICSE, the ONK-EH3 is used to model the round shafts. Vertical heaters will be placed outside the experiment hole to simulate the external decay heat from the canister panels. The acoustic emission sensor array located in four holes around the ONK-EH3 hole from the preceding experiment will be reused. New arrays of strain gages and thermal sensors are installed to monitor the behavior of both rock surface and the concrete liner. Numerical modelling was used extensively:(1)
Preliminary  3D FEM modelling was used on a 1/8th subgeometry to determine amount of heater elements and location of heater holes and to plan the heating power pattern [4].(2)
Full-scale 3D FEM modelling was used to include the effect of the above tunnel, to add the in situ stress effect, and to plan the monitor locations and probe the expected responses [5].(3)
2D fracture mechanics modelling was used to study the failure onset and propagation and to produce a quantified effect of the reduction provided by the concrete liner [5].Thermally induced spalling can occur in highly stressed rock, for example, spent nuclear fuel repositories in deep crystalline rock. Concrete support may be used to provide support pressure. The interaction of the concrete liner and surrounding rock mass during external heating is a thermomechanical transient problem. The problem must be studied in detail to confirm the validity of the chosen design methods and to establish a factor of safety.The spalling reducing effect of support pressure has been studied in ÄSPÖ Hard Rock Laboratory in Sweden. Glamheden et al. [6] conducted the Counterforce Applied to Prevent Spalling (CAPS) in situ experiment in which 4 holes were filled with LECA pellets and 4 holes were left unconfined. It was observed that the support pressure created will decrease the thermal-induced spalling, even at very small support pressures of 10–30 kPa with minor deformation. In Äspö Pillar Stability Experiment (APSE) first experiment hole was left unconfined and second hole was confined with a rubber bladder. APSE indicated that spalling could be totally controlled with approximately 150 kPa confining pressure [2]. A similar study was conducted in ONKALO in the second phase of POSE experiment, where the confinement was created using sand filling. The second phase of POSE experiment concluded that confined hole actually had more damage compared to unconfined hole [7]. However, this might result from the complicated geology of ONKALO. In the third phase of POSE the ONK-EH3 hole was totally filled with aluminium oxide providing similar support pressure compared to sand. The experiment resulted in minor damage only [3] enabling the ICSE, where the spalling suppressing effect of support pressure can be further investigated by applying concrete liner to only the bottom part of the hole.Confinement increases the compressive strength of concrete. The tubular thin concrete liner is a self-confining continuously arched shell structure. Outwards expansion is prevented by the surrounding rock mass and expansion downwards is prevented by the rock mass below. The top expansion is only partially inhibited and the floor of the tunnel may rise. This creates geometric asymmetry in the experiment but also enables the study of the effect of free and restricted movement.The 2D fracture mechanics model uses F-criterion [8] modified from G-criterion [9] for fracture propagation and Mohr-Coulomb criterion [10, 11] to determine fracture initiation. However, the Mohr-Coulomb criterion lacks the effect of the intermediate stress σ2 and it has sharp corners in the deviatoric plane, which can cause numerical problems and make the flow rule indeterminate. The Drucker-Prager criterion [12] is circular in the deviatoric plane and lacks the cos⁡3θ parameter to allow for effect of the deviatoric stress. The Ottosen [13] four-parameter failure criterion contains all the three stress invariants explicitly. At high hydrostatic stress it is semicircular and at low hydrostatic stress it has a more triangular shape.In the planned experiment a nonzero intermediate stress is expected and while the thermally induced stress is hydrostatic, a secondary deviatoric effect can arise from the above tunnel geometry. The hydrostatic stress of the concrete liner will begin at zero and increase continuously during the experiment until failure occurs. Therefore, the Ottosen criterion was deemed necessary to be used with the finite element models. The Ottosen criterion may be reduced back to Drucker-Prager criterionA=0,λ=constant or von Mises criterion A=0,B=0,λ=constant. The criterion is given by(1)AJ2σc2+λJ2σc+BI1σc-1=0,where(2)λ=K1cos⁡13arccos⁡K2cos⁡3θifcos⁡3θ≥0K1cos⁡π3-13arccos⁡-K2cos⁡3θifcos⁡3θ≤0.A, B, K1, and K2 are dimensionless material parameters, which can be obtained from laboratory testing data following [14] using four failure points on the failure surface, for example, compressive strength fcm, a failure state on the compressive meridian I1,J2, uniaxial tensile strength fctm, and biaxial compressive strength fcbc. Appendix A provides the formulae and calculated parameter values for Eurocode concrete strength classes.Ottosen [13] is not the only multiparameter octahedral shear failure criterion and [15] has identified more than 30 similar criteria. Models with sharp corners are ignored due to lack of explicit derivate at sharp corners, which causes computational problems. Notable alternatives suitable for numerical calculations include Hoek-Brown [16, 17], Menetrey-Willam [18], Willam-Warnke [19], and Bresler-Pister [20]. Hoek-Brown is a purely empirical fit and it ignores the effect of the intermediate stress. Menetrey-Willam is a three-parameter failure criterion based on Hoek-Brown with one shape parameter for the deviatoric surface shape. Menetrey-Willam improves on Hoek-Brown by accounting for the intermediate stress and allowing the deviatoric section to gradually change from triangular to circular. Willam-Warnke and Bresler-Pister are covered in detail by [21]. Bresler-Pister implements curved envelopes but has only one failure meridian and no deviatoric effects. The three-parameter Willam-Warnke produces a smooth solutions surface but retains the linear envelopes in the meridian planes. The five-parameter Willam-Warnke has all the desired mechanical properties and it produces a slightly better fit than the Ottosen, but it uses two empirically determined parameters. The five-parameter Willam-Warnke is the recommended model, when enough true triaxial laboratory data is available. For this research, Ottosen is the simplest suitable constitutive model. The four parameters are derived based on the uniaxial compressive strength and selected literature sources on concrete material testing [22–25] as shown in Appendix A.


## 2. Initial Data

### 2.1. ICSE Dimensions

The general layout of the experiment is shown in Figure2(a). The main dimensions of the experiment are shown in Figure 2(b). The plan view of the area with hole labels is shown in Figure 3. The ICSE was modelled using one 2D model and two 3D models. The 2D simulation was carried out in an infinite horizontal plane crossing the hole at −3 m depth. The 3D preliminary simulation was carried out as a subgeometry of a finite rock mass cylinder with height of 10.5 m and radius of 20 m and the final 3D model was contained in a rock mass cube with dimensions 100 m × 100 m × 100 m. The pilot hole, which acted as a drilling guide for the drilling the larger experiment hole, is 9.27 m in length and 300 mm in diameter. The main hole (ONK-EH3) is 7.2 m deep and 762 mm in radius. The concrete liner is 40 mm thick and it covers the bottom 4.2 m of the hole surface (Figure 2(b)). The top 3.0 m is unreinforced. To provide a level plane for the borehole drilling, a 30 cm thick concrete slab was cast in a 7 m × 7 m square centered to the hole center. Eight heater holes are located 1.762 m away from ONK-EH3 in a symmetrical radial array and they are 8 m deep and 76 mm in radius. The four acoustic emission holes are 7.5 m deep and 76 mm in radius. Four of the temperature measurement holes are 10 m deep each and the fifth hole T12 is 7 m deep. Finally, the floor of the experiment area is covered with two perpendicular layers of Paroc eXtra XS 10 cm thick rock wool insulation for a total thickness of 200 mm in an 8 m × 8 m area.Figure 3
The resulting heater pattern (HH) and preexisting temperature monitoring holes (T) and acoustic emission monitoring holes (AE). The inner 7 m × 7 m square shown is the levelling concrete slab and the outer 8 m × 8 m square is the insulation layer.


### 2.2. In Situ Stress State

There are two applicable in situ stress state interpretations, which will be referred to as “EDZ and VT1” and “ONK-EH3” hereafter. EDZ and VT1 interpretation originates from the original smaller research niche 30 m South from the ICSE. The newer ONK-EH3 interpretation is based on measurements from the ICSE area. The stress states are reported in more detail in [26], shown in Table 1, and illustrated in Figures 4 and 5. The ONK-EH3 is mostly located in pegmatite for the top four meters and then gradually introducing veined gneiss in the hole bottom.Table 1
EDZ and VT1 and ONK-EH3 in situ stress state interpretations at depth −345 m.
Stress interpretation σ 1 [MPa] Dip direction [°] Dip [°] σ 2 [MPa] Dip direction [°] Dip [°] σ 3 [MPa] Dip direction [°] Dip [°] EDZ and VT1 25.1 166 1 17.1 256 10 12.3 68 80 ONK-EH3 18.2 120 2 15.6 210 4 8.9 3 85Figure 4
EDZ and VT1 stress state interpretation [26].Figure 5
ONK-EH3 stress state interpretation [26].The 2D in-plane stresses used in 2D modelling are acquired from previous 3D modelling work [26] including the effect of the tunnel above. The used in-plane stresses are shown in Table 2. The 2D modelling was carried out in the plane of horizontal principal stresses with y-axis being the larger principal stress. The 2D modelling was carried out in plane strain condition.Table 2
The in-plane stresses below the tunnel at experiment hole location.
Stress interpretation σ y y [MPa] Dip direction [°] σ x x [MPa] Dip direction [°] σ z z [MPa] EDZ and VT1 −3 m 25.1 166 21.5 256 4.0 ONK-EH3 −3 m 23.0 120 15.0 210 3.0


### 2.3. Material Properties

The parameters used are the same as those in thermomechanics prediction for the third phase of POSE experiment by [26] complemented with data from [27–31]. The rock mass is assumed to be homogeneous, isotropic, linearly elastic, and mostly pegmatitic granite. For the 3D FEM modelling, more generalized rock mass properties are used, as the rock material is migmatitic gneiss containing both gneissic and granitic sections. Young’s modulus is reduced from 60 GPa [31] to 53 GPa [26] to account for rock discontinuities. The difference in Young’s moduli between rock and concrete will cause a jump in the strain at the material boundary.For the fracture mechanics prediction, the thermal rock properties [27, 28] of pegmatitic granite are used; however the elastic properties are the same for both modelling approaches due to insufficient amount (n=13) of test data for pegmatitic granite. The thermal expansion coefficients of rock and concrete are close to each other and there is no significant difference in heat capacity. The insulation was given arbitrary low values for mechanical parameters to prevent it from contributing to the mechanical solution.The material properties for rock, concrete, and insulation are shown in Table3 and the material strengths for rock and concrete are shown in Table 4. The concrete was modelled with the Ottosen [13] plasticity (Appendix A) and the used parameters are show in Table 5. A bilinear stress-strain relation as described in [33] was used and the concrete is assumed damaged when the effective plastic strain exceeds 0.00175, which corresponds to total strain exceeding 0.0035.Table 3
Material properties used in simulations.
Material Thermal conductivityk [W/m·K] Densityρ [kg/m3] Heat capacityCp [J/kg·K] Poisson’s ratioν Young’s modulusE [GPa] Coefficient of thermal expansionα [1/K] Rock mass (3D) 3.33 a  R M 2635 b  P G R 716 a  R M 0.25 c  V G N 53 a  R M 9.76 e - 6 d  V G N Pegmatite (2D) 3.20 b  P G R 2635 b  P G R 689 a  P G R 0.25 c , ∗  V G N 53 a  R M 7.20 e - 6 d  P G R Concrete C35/45 1.70 e 2200 f 840 g 0.20 b 34 b 10 e - 6 e Insulation 0.0353 h 28.5 h 840 h 0.10 i 0.1 i 10 e - 10 i a[26]. b[27]. c[32]. d[28]. e[29]. fBased on quality assurance tests for sprayed concrete. g[30]. hManufacturer quote. iArbitrary value. ∗Insufficient amount of test results for PGR. RMProperty for rock mass. VGNProperty for veined gneiss. PGRProperty for pegmatitic granite.Table 4
Material strengths based on measured values and material properties.
Material Uniaxial compressive strength, mean [MPa] Uniaxial compressive strength, characteristic [MPa] Uniaxial tensile strength, mean [MPa] Uniaxial tensile strength, characteristic [MPa] Rock mass U C S m = 102 a  P G R U C S k = 65.5 b  R M U T S m = 12 b  R M Concrete C35/45 f c m = 43 c f c k = 35 f c t m = 3.2 c f c t k , 0.05 = 2.2 c a[31]. b[34]. cCalculated value [33]. PGRProperty for pegmatitic granite. RMProperty for rock mass.Table 5
Ottosen plasticity criterion parameters for C35/45 concrete.
Material f c m A B K 1 K 2 ε c 3 ε c u 3 Concrete C35/45 43  M P a a 2.0098 b 4.4388 b 15.4402 b 0.9932 b 0.00175 a 0.0035 a aCalculated value [33]. bCalculated as shown in Appendix A.The fracture mechanics parameters from [34, 35] are shown in Table 6 complemented with data from [32, 33]. There was no laboratory testing data available for the concrete fracture mechanics parameters and values were derived by [5] from literature relationships by [36–38]. For concrete, the lowest principal stress is close to zero and cohesion will determine the strength. A nonzero friction angle was still used for better numerical stability [5].Table 6
Fracture mechanics modelling parameters.
Property Symbol Pegmatitic granite Concrete C35/45 Unit Cohesion c 12.9 a 17.5 e MPa Friction angle ϕ 47 a 1 e ° Tensile strength σ t 12 a 2.2 d MPa Mode I fracture toughness K I C 1.96 b 0.31 e MPam Mode II fracture toughness K I I C 3.30 b 1.84 f MPam Cohesion, tensile & shear c 10 a 10 a MPa Friction angle, tensile & shear ϕ t, ϕs 35 c , 35 c 35 c ,  35 c °, ° Dilatation angle, tensile & shear ψ t, ψs 2.5 c ,  2.5 c 2.5 c, 2.5c °, ° Normal stiffness, tensile & shear k n 20,000 a 20,000 a GPa/m Shear stiffness, tensile & shear k s 2,000 a 2,000 a GPa/m a[34]. bModified after [35]. c[32]. d[33]. eDerived by [5] after [36] and [37]. fCalculated after relation by [38].


### 2.4. Initial Values

The initial values for convection and for thermal calculation initial values are shown in Table7. The initial temperature of the rock is expected to be 18°C after the third phase of POSE experiment. The in situ temperature was set to 13.5°C based on temperature measurements in the niche. The air temperature and air pressure were only used in the preliminary calculations to study the effect of convection to the overall solution. The starting temperatures affect how long the heating period lasts but have only a small impact on the induced stresses.Table 7
Initial values.
Property Symbol Value Unit Air temperature T a i r 286.65 a (13.5) K [°C] Air pressure ρ a i r 105 b (1.05) kPa [bar] Rock mass near field temperature T r m 291.15 a (18) K [°C] Rock mass far field temperature T r m , f f 284.65 a (11.5) K [°C] Concrete temperature T s c 287.15 c (14) K [°C] a[26]. bAir pressure at −345 m. cEstimation based on ambient air temperature.


### 2.5. Heating Power Plan

The heater amount and positions were determined iteratively using the preliminary model and requiring an evenly distributed temperature in the main hole surface reaching the target temperature of 120°C within a time window of two months. The resulting heater arrangement is shown in Figure3 with heater holes denoted by HH1-8. Eight heaters were placed symmetrically around the main hole one meter away from the experiment hole. The arrangement was slightly rotated clockwise to provide maximum clearance from existing measurement holes. The heater holes are 8 meters deep and the heater elements are 6 meters long.The heating patterns used are shown in Figure6. In the 2D fracture mechanics analysis the heat power (on right axis) is increased in three equal size steps over 9 weeks. For the 3D preliminary analysis, the output (on left axis) is higher and the heaters are turned off on week 9. For the full model, after 9 weeks, the heaters are turned on with maximum power output of 4000 W per heater for the last 7 weeks. During the in situ experiment, if there is sufficient damage by visual observation before end of week 9, the heaters will be shut down. If the hole has not suffered observable damage, the heaters are turned on with maximum power for a maximum of 7 weeks or until observable damage has occurred.Figure 6
The heating power plan per heater unit.


## 2.1. ICSE Dimensions

The general layout of the experiment is shown in Figure2(a). The main dimensions of the experiment are shown in Figure 2(b). The plan view of the area with hole labels is shown in Figure 3. The ICSE was modelled using one 2D model and two 3D models. The 2D simulation was carried out in an infinite horizontal plane crossing the hole at −3 m depth. The 3D preliminary simulation was carried out as a subgeometry of a finite rock mass cylinder with height of 10.5 m and radius of 20 m and the final 3D model was contained in a rock mass cube with dimensions 100 m × 100 m × 100 m. The pilot hole, which acted as a drilling guide for the drilling the larger experiment hole, is 9.27 m in length and 300 mm in diameter. The main hole (ONK-EH3) is 7.2 m deep and 762 mm in radius. The concrete liner is 40 mm thick and it covers the bottom 4.2 m of the hole surface (Figure 2(b)). The top 3.0 m is unreinforced. To provide a level plane for the borehole drilling, a 30 cm thick concrete slab was cast in a 7 m × 7 m square centered to the hole center. Eight heater holes are located 1.762 m away from ONK-EH3 in a symmetrical radial array and they are 8 m deep and 76 mm in radius. The four acoustic emission holes are 7.5 m deep and 76 mm in radius. Four of the temperature measurement holes are 10 m deep each and the fifth hole T12 is 7 m deep. Finally, the floor of the experiment area is covered with two perpendicular layers of Paroc eXtra XS 10 cm thick rock wool insulation for a total thickness of 200 mm in an 8 m × 8 m area.Figure 3
The resulting heater pattern (HH) and preexisting temperature monitoring holes (T) and acoustic emission monitoring holes (AE). The inner 7 m × 7 m square shown is the levelling concrete slab and the outer 8 m × 8 m square is the insulation layer.


## 2.2. In Situ Stress State

There are two applicable in situ stress state interpretations, which will be referred to as “EDZ and VT1” and “ONK-EH3” hereafter. EDZ and VT1 interpretation originates from the original smaller research niche 30 m South from the ICSE. The newer ONK-EH3 interpretation is based on measurements from the ICSE area. The stress states are reported in more detail in [26], shown in Table 1, and illustrated in Figures 4 and 5. The ONK-EH3 is mostly located in pegmatite for the top four meters and then gradually introducing veined gneiss in the hole bottom.Table 1
EDZ and VT1 and ONK-EH3 in situ stress state interpretations at depth −345 m.
Stress interpretation σ 1 [MPa] Dip direction [°] Dip [°] σ 2 [MPa] Dip direction [°] Dip [°] σ 3 [MPa] Dip direction [°] Dip [°] EDZ and VT1 25.1 166 1 17.1 256 10 12.3 68 80 ONK-EH3 18.2 120 2 15.6 210 4 8.9 3 85Figure 4
EDZ and VT1 stress state interpretation [26].Figure 5
ONK-EH3 stress state interpretation [26].The 2D in-plane stresses used in 2D modelling are acquired from previous 3D modelling work [26] including the effect of the tunnel above. The used in-plane stresses are shown in Table 2. The 2D modelling was carried out in the plane of horizontal principal stresses with y-axis being the larger principal stress. The 2D modelling was carried out in plane strain condition.Table 2
The in-plane stresses below the tunnel at experiment hole location.
Stress interpretation σ y y [MPa] Dip direction [°] σ x x [MPa] Dip direction [°] σ z z [MPa] EDZ and VT1 −3 m 25.1 166 21.5 256 4.0 ONK-EH3 −3 m 23.0 120 15.0 210 3.0


## 2.3. Material Properties

The parameters used are the same as those in thermomechanics prediction for the third phase of POSE experiment by [26] complemented with data from [27–31]. The rock mass is assumed to be homogeneous, isotropic, linearly elastic, and mostly pegmatitic granite. For the 3D FEM modelling, more generalized rock mass properties are used, as the rock material is migmatitic gneiss containing both gneissic and granitic sections. Young’s modulus is reduced from 60 GPa [31] to 53 GPa [26] to account for rock discontinuities. The difference in Young’s moduli between rock and concrete will cause a jump in the strain at the material boundary.For the fracture mechanics prediction, the thermal rock properties [27, 28] of pegmatitic granite are used; however the elastic properties are the same for both modelling approaches due to insufficient amount (n=13) of test data for pegmatitic granite. The thermal expansion coefficients of rock and concrete are close to each other and there is no significant difference in heat capacity. The insulation was given arbitrary low values for mechanical parameters to prevent it from contributing to the mechanical solution.The material properties for rock, concrete, and insulation are shown in Table3 and the material strengths for rock and concrete are shown in Table 4. The concrete was modelled with the Ottosen [13] plasticity (Appendix A) and the used parameters are show in Table 5. A bilinear stress-strain relation as described in [33] was used and the concrete is assumed damaged when the effective plastic strain exceeds 0.00175, which corresponds to total strain exceeding 0.0035.Table 3
Material properties used in simulations.
Material Thermal conductivityk [W/m·K] Densityρ [kg/m3] Heat capacityCp [J/kg·K] Poisson’s ratioν Young’s modulusE [GPa] Coefficient of thermal expansionα [1/K] Rock mass (3D) 3.33 a  R M 2635 b  P G R 716 a  R M 0.25 c  V G N 53 a  R M 9.76 e - 6 d  V G N Pegmatite (2D) 3.20 b  P G R 2635 b  P G R 689 a  P G R 0.25 c , ∗  V G N 53 a  R M 7.20 e - 6 d  P G R Concrete C35/45 1.70 e 2200 f 840 g 0.20 b 34 b 10 e - 6 e Insulation 0.0353 h 28.5 h 840 h 0.10 i 0.1 i 10 e - 10 i a[26]. b[27]. c[32]. d[28]. e[29]. fBased on quality assurance tests for sprayed concrete. g[30]. hManufacturer quote. iArbitrary value. ∗Insufficient amount of test results for PGR. RMProperty for rock mass. VGNProperty for veined gneiss. PGRProperty for pegmatitic granite.Table 4
Material strengths based on measured values and material properties.
Material Uniaxial compressive strength, mean [MPa] Uniaxial compressive strength, characteristic [MPa] Uniaxial tensile strength, mean [MPa] Uniaxial tensile strength, characteristic [MPa] Rock mass U C S m = 102 a  P G R U C S k = 65.5 b  R M U T S m = 12 b  R M Concrete C35/45 f c m = 43 c f c k = 35 f c t m = 3.2 c f c t k , 0.05 = 2.2 c a[31]. b[34]. cCalculated value [33]. PGRProperty for pegmatitic granite. RMProperty for rock mass.Table 5
Ottosen plasticity criterion parameters for C35/45 concrete.
Material f c m A B K 1 K 2 ε c 3 ε c u 3 Concrete C35/45 43  M P a a 2.0098 b 4.4388 b 15.4402 b 0.9932 b 0.00175 a 0.0035 a aCalculated value [33]. bCalculated as shown in Appendix A.The fracture mechanics parameters from [34, 35] are shown in Table 6 complemented with data from [32, 33]. There was no laboratory testing data available for the concrete fracture mechanics parameters and values were derived by [5] from literature relationships by [36–38]. For concrete, the lowest principal stress is close to zero and cohesion will determine the strength. A nonzero friction angle was still used for better numerical stability [5].Table 6
Fracture mechanics modelling parameters.
Property Symbol Pegmatitic granite Concrete C35/45 Unit Cohesion c 12.9 a 17.5 e MPa Friction angle ϕ 47 a 1 e ° Tensile strength σ t 12 a 2.2 d MPa Mode I fracture toughness K I C 1.96 b 0.31 e MPam Mode II fracture toughness K I I C 3.30 b 1.84 f MPam Cohesion, tensile & shear c 10 a 10 a MPa Friction angle, tensile & shear ϕ t, ϕs 35 c , 35 c 35 c ,  35 c °, ° Dilatation angle, tensile & shear ψ t, ψs 2.5 c ,  2.5 c 2.5 c, 2.5c °, ° Normal stiffness, tensile & shear k n 20,000 a 20,000 a GPa/m Shear stiffness, tensile & shear k s 2,000 a 2,000 a GPa/m a[34]. bModified after [35]. c[32]. d[33]. eDerived by [5] after [36] and [37]. fCalculated after relation by [38].


## 2.4. Initial Values

The initial values for convection and for thermal calculation initial values are shown in Table7. The initial temperature of the rock is expected to be 18°C after the third phase of POSE experiment. The in situ temperature was set to 13.5°C based on temperature measurements in the niche. The air temperature and air pressure were only used in the preliminary calculations to study the effect of convection to the overall solution. The starting temperatures affect how long the heating period lasts but have only a small impact on the induced stresses.Table 7
Initial values.
Property Symbol Value Unit Air temperature T a i r 286.65 a (13.5) K [°C] Air pressure ρ a i r 105 b (1.05) kPa [bar] Rock mass near field temperature T r m 291.15 a (18) K [°C] Rock mass far field temperature T r m , f f 284.65 a (11.5) K [°C] Concrete temperature T s c 287.15 c (14) K [°C] a[26]. bAir pressure at −345 m. cEstimation based on ambient air temperature.


## 2.5. Heating Power Plan

The heater amount and positions were determined iteratively using the preliminary model and requiring an evenly distributed temperature in the main hole surface reaching the target temperature of 120°C within a time window of two months. The resulting heater arrangement is shown in Figure3 with heater holes denoted by HH1-8. Eight heaters were placed symmetrically around the main hole one meter away from the experiment hole. The arrangement was slightly rotated clockwise to provide maximum clearance from existing measurement holes. The heater holes are 8 meters deep and the heater elements are 6 meters long.The heating patterns used are shown in Figure6. In the 2D fracture mechanics analysis the heat power (on right axis) is increased in three equal size steps over 9 weeks. For the 3D preliminary analysis, the output (on left axis) is higher and the heaters are turned off on week 9. For the full model, after 9 weeks, the heaters are turned on with maximum power output of 4000 W per heater for the last 7 weeks. During the in situ experiment, if there is sufficient damage by visual observation before end of week 9, the heaters will be shut down. If the hole has not suffered observable damage, the heaters are turned on with maximum power for a maximum of 7 weeks or until observable damage has occurred.Figure 6
The heating power plan per heater unit.


## 3. Numerical Simulations

The goals for the preliminary numerical simulations were to determine heater arrangement and heater power scheme and to establish expected values for temperature and stress, which clearly exceed the calculated mean uniaxial strength values for concrete and observed uniaxial compressive strength laboratory test values from Olkiluoto rock. The temperature was designed not to exceed the thermal tolerance of the strain gages at 130°C. Full-scale numerical modelling was carried out to take into account the effect of the hole and the tunnel above and in situ stresses and to produce predicted results at planned instrumentation locations. Fracture mechanics predictions were used to predict failure initiation, propagation, and extent in concrete and in rock. Additionally the fracture mechanics predictions may later aid the interpretation of the acoustic emission results.Several different subgeometries were used in the scoping analyses: The preliminary simulations exploited octosymmetry and a 1/16th sector surrounding the hole was modelled in 3D FEM. 2D fracture mechanical analysis was carried out in a horizontal infinite plane at −3 m level. Finally, the full-scale model houses the POSE niche in a 100 m × 100 m × 100 m rock mass cube.


### 3.1. Preliminary Model

Preliminary simulations were carried out using COMSOL Multiphysics 4.3a Thermal Stress module. The preliminary model stage assumed the entire hole is lined with concrete and the concrete slab, pilot hole, and convection were also included to study their contribution to the solution.The slice model assumes the hole is located in a half-space and only the heating induced stress increase was studied. The half-space assumption is somewhat incorrect and will cause error to the results mainly because of lack of restrictions on upwards movement. The insulation, foundation slab, the pilot hole, and convection of air above the tunnel were included. The concrete and the rock mass were meshed together with shared nodes. Figure7(a) shows the modelled area, which includes one-half of a heater hole. Liner thickness, amount, and location of heater holes and heater power plan were optimized iteratively to produce a high stresses without excessive heating.Figure 7
Location of the 1/16th slice taken and conversion of orthogonal components to cylinders (a), the whole modelling domain with mesh (b), and mesh detail showing the insulation, concrete foundation slab, liner, and half a heater hole (c).
(a) (b) (c)The size of the elements for the liner was determined using a parametric 2D sweep and monitoring the change of temperature and tangential stresses at the inner surface of the 40 mm thick concrete layer. For quadratic elements all element sizes from 4 mm up to 40 mm produced an acceptable result with less than 1% difference in the tangential stresses. For linear elements, two layers of 20 mm elements produced a discretization error of −4.7% underestimating the stresses compared to the quadratic solutions. For temperature, all element sizes from 4 mm up to 40 mm produce an acceptable result with less than 1% difference using the linear shape functions.The preliminary model was meshed with general physics driven mesher at extrafine setting (Figure7(b)). The liner mesh was meshed at extremely fine setting with maximum element size constrained to 30 mm (Figure 7(c)) to ensure two layers of elements to the 40 mm thick layer of concrete liner. The total mesh consisted of 1,001,210 elements and 5,488,561 degrees of freedom. For displacements second-order shape functions were used and for temperature linear shape functions were used.On the hole surface, heat flux was set to zero corresponding to the fully insulated condition. The experience from previous experiments in the POSE niche is that modelling with zero thermal flux corresponds well with observations. The effect of convection was small in previous experiments, but it was included on the tunnel floor and on free surfaces of the insulation. The long sides used the symmetry boundary condition to account for the missing geometry.


### 3.2. Full-Scale Model

Full-scale simulations were carried out using COMSOL Multiphysics 4.4 Thermal Stress module. The rock mass was modelled as a 100 m × 100 m × 100 m cube, which contains the POSE niche and the experiment area (Figure8(a)). The liner was modelled in high detail and three enhanced precision planes were created in direction towards the temperature sensors T2 (North) and T4 (West) and heater hole HH4 (Southeast) as shown in Figure 8(b). The pilot hole and convection were ignored based on their limited significance from preliminary modelling stage. The concrete foundation slab was converted into rock mass introducing a small material error near the top surface. Based on preliminary modelling, the concrete-rock interface is without tension and without significant shear. The interface remains in compression through the heating period, but cooling may induce tension. Based on these results, interface elements were not used and the concrete and the rock mass were simply meshed together with shared nodes. The potential for interface damage can be evaluated from the stresses at both sides of the interface. The total mesh consists of 1,323,298 elements and 5,313,195 degrees of freedom. For displacements second-order shape functions were used and for temperature linear shape functions. The transient heating with plasticity for concrete only was modelled with linear shape functions for displacements resulting in 857,156 degrees of freedom. For detailed analysis of the concrete plasticity another mesh variant was created. In this variant, all geometrical features are reduced to bare minimum (only tunnel and the experiment hole) and the cylindrical heaters are replaced with lines. This results in a mesh with 809,779 elements and 423,720 degrees of freedom. The concrete mesh has two layers of elements, totaling 230,111 (28%) linear tetrahedrons, which was found to give less than 5% difference in stress with 2D element size scoping analysis.Figure 8
Rock mass block (100 m × 100 m × 100 m) housing the POSE niche and the experiment area (a) and mesh detail showing the three mesh guidance planes (b).
(a) (b)


### 3.3. Fracture Mechanical Simulations

The Fracod is a 2-dimensional BEM/DDM code, which provides an exact solution for the governing partial equations. The fracture mechanics predictions were carried out using Fracod2D 4.11 as described by [8]. DDM, in contrast to FEM, provides exact basis of calculating the stresses especially related to the crack tips, particularly important when modelling fracture propagation. DDM has been used for years for fracture mechanics applications [39]. In [40] the similarities between DDM and BEM have been discussed and shown that they are of identical nature and can be used for solving fracture mechanics problems.Fracod2D uses fictitious heat source method to calculate thermal stresses [8]. Comparing the 2D approach to 3D results, it is noted that the effective heat is slightly higher in 2D due to the absence of the end effect as the 2D model cannot include out-of-plane geometry. The heat flux at the hole inner surface is assumed to be 5% of the total heating power. The thermal input is constructed for each time step (3 week) as an average from the input power at the beginning and end of each time step. Therefore the heating scheme is smoothed in the beginning of last heating step but retains almost the same total power as presented in Figure 6 with the absence of end effect compensating the power loss in fifth week (power loss is 3.6% of total power used in experiment). The final heating power is kept the same as in experiment heating scheme to not exceed the maximum temperate at any time step. The cooling period was not modelled with Fracod2D. The material is assumed to be continuous, isotropic, homogeneous, and linearly elastic with fracture growth. For the two stress interpretations separate models were calculated using quarter symmetry with fracture initiation element size 60 mm and boundary element spacing of 30 mm at hole surface.


## 3.1. Preliminary Model

Preliminary simulations were carried out using COMSOL Multiphysics 4.3a Thermal Stress module. The preliminary model stage assumed the entire hole is lined with concrete and the concrete slab, pilot hole, and convection were also included to study their contribution to the solution.The slice model assumes the hole is located in a half-space and only the heating induced stress increase was studied. The half-space assumption is somewhat incorrect and will cause error to the results mainly because of lack of restrictions on upwards movement. The insulation, foundation slab, the pilot hole, and convection of air above the tunnel were included. The concrete and the rock mass were meshed together with shared nodes. Figure7(a) shows the modelled area, which includes one-half of a heater hole. Liner thickness, amount, and location of heater holes and heater power plan were optimized iteratively to produce a high stresses without excessive heating.Figure 7
Location of the 1/16th slice taken and conversion of orthogonal components to cylinders (a), the whole modelling domain with mesh (b), and mesh detail showing the insulation, concrete foundation slab, liner, and half a heater hole (c).
(a) (b) (c)The size of the elements for the liner was determined using a parametric 2D sweep and monitoring the change of temperature and tangential stresses at the inner surface of the 40 mm thick concrete layer. For quadratic elements all element sizes from 4 mm up to 40 mm produced an acceptable result with less than 1% difference in the tangential stresses. For linear elements, two layers of 20 mm elements produced a discretization error of −4.7% underestimating the stresses compared to the quadratic solutions. For temperature, all element sizes from 4 mm up to 40 mm produce an acceptable result with less than 1% difference using the linear shape functions.The preliminary model was meshed with general physics driven mesher at extrafine setting (Figure7(b)). The liner mesh was meshed at extremely fine setting with maximum element size constrained to 30 mm (Figure 7(c)) to ensure two layers of elements to the 40 mm thick layer of concrete liner. The total mesh consisted of 1,001,210 elements and 5,488,561 degrees of freedom. For displacements second-order shape functions were used and for temperature linear shape functions were used.On the hole surface, heat flux was set to zero corresponding to the fully insulated condition. The experience from previous experiments in the POSE niche is that modelling with zero thermal flux corresponds well with observations. The effect of convection was small in previous experiments, but it was included on the tunnel floor and on free surfaces of the insulation. The long sides used the symmetry boundary condition to account for the missing geometry.


## 3.2. Full-Scale Model

Full-scale simulations were carried out using COMSOL Multiphysics 4.4 Thermal Stress module. The rock mass was modelled as a 100 m × 100 m × 100 m cube, which contains the POSE niche and the experiment area (Figure8(a)). The liner was modelled in high detail and three enhanced precision planes were created in direction towards the temperature sensors T2 (North) and T4 (West) and heater hole HH4 (Southeast) as shown in Figure 8(b). The pilot hole and convection were ignored based on their limited significance from preliminary modelling stage. The concrete foundation slab was converted into rock mass introducing a small material error near the top surface. Based on preliminary modelling, the concrete-rock interface is without tension and without significant shear. The interface remains in compression through the heating period, but cooling may induce tension. Based on these results, interface elements were not used and the concrete and the rock mass were simply meshed together with shared nodes. The potential for interface damage can be evaluated from the stresses at both sides of the interface. The total mesh consists of 1,323,298 elements and 5,313,195 degrees of freedom. For displacements second-order shape functions were used and for temperature linear shape functions. The transient heating with plasticity for concrete only was modelled with linear shape functions for displacements resulting in 857,156 degrees of freedom. For detailed analysis of the concrete plasticity another mesh variant was created. In this variant, all geometrical features are reduced to bare minimum (only tunnel and the experiment hole) and the cylindrical heaters are replaced with lines. This results in a mesh with 809,779 elements and 423,720 degrees of freedom. The concrete mesh has two layers of elements, totaling 230,111 (28%) linear tetrahedrons, which was found to give less than 5% difference in stress with 2D element size scoping analysis.Figure 8
Rock mass block (100 m × 100 m × 100 m) housing the POSE niche and the experiment area (a) and mesh detail showing the three mesh guidance planes (b).
(a) (b)


## 3.3. Fracture Mechanical Simulations

The Fracod is a 2-dimensional BEM/DDM code, which provides an exact solution for the governing partial equations. The fracture mechanics predictions were carried out using Fracod2D 4.11 as described by [8]. DDM, in contrast to FEM, provides exact basis of calculating the stresses especially related to the crack tips, particularly important when modelling fracture propagation. DDM has been used for years for fracture mechanics applications [39]. In [40] the similarities between DDM and BEM have been discussed and shown that they are of identical nature and can be used for solving fracture mechanics problems.Fracod2D uses fictitious heat source method to calculate thermal stresses [8]. Comparing the 2D approach to 3D results, it is noted that the effective heat is slightly higher in 2D due to the absence of the end effect as the 2D model cannot include out-of-plane geometry. The heat flux at the hole inner surface is assumed to be 5% of the total heating power. The thermal input is constructed for each time step (3 week) as an average from the input power at the beginning and end of each time step. Therefore the heating scheme is smoothed in the beginning of last heating step but retains almost the same total power as presented in Figure 6 with the absence of end effect compensating the power loss in fifth week (power loss is 3.6% of total power used in experiment). The final heating power is kept the same as in experiment heating scheme to not exceed the maximum temperate at any time step. The cooling period was not modelled with Fracod2D. The material is assumed to be continuous, isotropic, homogeneous, and linearly elastic with fracture growth. For the two stress interpretations separate models were calculated using quarter symmetry with fracture initiation element size 60 mm and boundary element spacing of 30 mm at hole surface.


## 4. Results and Discussion

### 4.1. Preliminary Model

The evenness of the heating in the surface of the hole is affected by the amount of heater holes, heater hole distance, and heater power pattern. By using a coarse mesh, it was found that using 8 vertical heaters one meter away from the hole wall produces an even heating pattern in a reasonable experiment duration. The array was rotated slightly counterclockwise to provide maximum clearance from already existing temperature monitoring and acoustic event sensor holes. The heating pattern chosen is shown in Figure6.The heat initially increases at a relatively steady rate of 12°C per week, with the rate slowing down towards week 9 (Figure9). 100°C is reached early on week 7. This leaves enough time for all the water to evaporate before being converted into steam. At the end of week 9, a temperature of 125°C is reached. Heater holes are at a temperature of 160°C at the hole top and 140°C along the top three meters. The concrete liner reaches 100°C at the beginning of 7th week and reaches 125°C (at −1 m level) at the end of 9th week (Figure 9). After the heating is turned off the temperature quickly levels out and dissipates. The seven-week cooling period brings the temperature back down to 35°C. This can further be accelerated by continuously flushing the hole with cool air.Figure 9
Calculated temperature at the center of the concrete liner (weeks 0–9).Each heater element is 6 meters long and the temperature quickly drops when approaching the hole bottom. At −6 m the maximum temperature reached at 9 weeks is 75°C and at 16 weeks is 155°C. Radially at −3 m level the temperature quickly rises towards the heater hole and quickly drops down to zero change at radius of 7.5 meters from the hole center. There is a small temperature drop (−7°C) in top part of the temperature curves which is caused by the convection and the more conductive concrete foundation slab. The pilot hole effects cannot be seen at the hole bottom, but the cooling effect of the rock mass below is clearly seen and strong (Figure9).No damage in the concrete is expected during the first three weeks (Figure10). Eurocode 2 [33] defines the characteristic strength (fck=35 MPa) as 5% failure fraction. Neglecting the effect of the intermediate stress, after five weeks, by definition, there is a 5% chance of damage and after seven weeks the mean uniaxial strength is exceeded. After nine weeks of heating, the stress will peak at 53 MPa, which is 23% higher than the mean uniaxial strength (fcm=43 MPa) of the concrete. The highest loaded region is at level −3.3 m, which is close to the −3 m monitoring level. The second planned monitoring level at −6 will peak at a much lower stress of 35 MPa, which corresponds to the characteristic uniaxial compressive strength of the concrete liner.Figure 10
Calculated tangential stress at the inner surface of the concrete liner (weeks 0–9).While mean uniaxial strengths are exceeded with a margin, the Ottosen criterion predicts that the structure is not damaged. Therefore an extended heating period was added after week 9 with heaters set to maximum (Figure6). This will exceed greatly the maximum temperature of the strain gages and some sensor loss is to be expected in the actual experiment. However, the temperature follows a predictable pattern and the lost information can be extrapolated.The concrete fails when the tangential stress increase reaches 60 MPa. Damage initiation is predicted on week 10 at −3.5 m depth. The damage continues to propagate as the heating is continued. On week 11 damage is predicted between −1 m and −5.5 m. The stress change decreases towards hole bottom and at −6 m level is only +30 MPa after 9 weeks and it reaches failure during week 13. At the hole top 0.5 m there is a region with horizontal tensile stresses up to −5 MPa.The stresses in the rock are 1.56 times higher due to the difference in elastic modulus (Figure11). The jump at 0.3 m is caused by the change of material from rock mass to concrete foundation slab. A tangential stress increase of +80 MPa is reached at the end of week 9 and the elastic increase reaches up to +190 MPa during the extended heating period. The thermally induced stress acts together with the in situ stress (Table 1). The initial tangential stress around the hole is estimated to be 46,…,54 MPa and the combined stress should be 126,…,134 MPa already on week 9. It is possible that the rock wall may sustain damage as the estimated damage strength is 58,...,102 MPa. The liner produces a support pressure, which maximizes at 2.8 MPa at depth of −3.5 m.Figure 11
Preliminary model tangential stress increase calculated at the inner surface of the rock mass (weeks 0–9).The hole moves upwards 5.2 mm and outwards 0.8 mm at top which causes vertical tensile stresses to the topmost 1.5 m of the hole. The calculated stresses exceed the design tensile strength (1.6 MPa in Finland) for concrete after week 5 but peak out at 2.1 MPa and never reach the characteristic tensile strength (2.2 MPa) level. Some tensile cracks may be observed in the top part of the hole. The maximum inward deflection is 0.3 mm at −4 m level and maximum downward displacement is 0.8 mm at −7.2 m level.


### 4.2. Full-Scale Model

The full-scale model is different from the preliminary model as it includes the true geometry with the tunnel above the hole. Additionally, the hole is lined with concrete from −3 m to −7.2 m when in preliminary model it was fully lined. The in situ stresses are included in the model and two stress state interpretations are analyzed separately. Due to the simplifications, the thermally induced stresses are slightly higher (Figure12), but the damage localization and the concentrating effects of excavation shape can now be studied more accurately.Figure 12
Full model tangential stress increase calculated at the inner surface of the rock mass in the North wall (weeks 0–9).For the EDZ and VT1 in situ stress interpretation, the initial maximum stresses are Northeast and Southwest aligned and rotate towards North-South in the heating phase (Figure13). For concrete, the tangential stresses reach 59 MPa after 9 weeks of heating and 84 MPa after 16 weeks of heating. For rock mass, the tangential stresses reach 115 MPa after 9 weeks of heating and 200 MPa after 16 weeks of heating. For the ONK-EH3 interpretation, the initial stresses are almost North-South and do not change direction during heating (Figure 14). For concrete, the tangential stresses reach 59 MPa after 9 weeks of heating and 85 MPa after 16 weeks of heating. For rock mass, the tangential stresses reach 117 MPa after 9 weeks of heating and 205 MPa after 16 weeks of heating.Figure 13
First principal stress [MPa] at −3 m level resulting from (a) 0 weeks in situ stress, (b) 9 weeks of heating, and (c) 16 weeks of heating using the EDZ and VT1 stress interpretation. The left scale is for the concrete and the right scale is for the surrounding rock mass.
(a) (b) (c)Figure 14
First principal stress [MPa] at −3 m level resulting from (a) 0 weeks in situ stress, (b) 9 weeks of heating, and (c) 16 weeks of heating using the ONK-EH3 stress interpretation. The left scale is for the concrete and the right scale is for the surrounding rock mass.
(a) (b) (c)The plastic damage was studied using a refined version of the full model. In this model the concrete mesh was the priority and all other features were modelled using a reduced degree of precision. Both in situ interpretations produce identical results for concrete. The damage initiates in Northeast segment on week 5 (Figure15(a)). The concrete exceeds the ultimate strain (0.0035) and becomes damaged on week 10 (Figure 15(b)) and continues to deteriorate in North and South until the end of the heating period reaching plastic strain of 0.017 on week 16 (Figure 15(c)). The damaged regions are the North and South quadrants with the most damage in the centers.Figure 15
Effective plastic strain of the concrete layer at −3 m level after 5 weeks (a), 10 weeks (b), and 16 weeks (c) of heating.
(a) (b) (c)For the EDZ and VT1 stress interpretation, the 9-week values are 122 MPa and 86 MPa (−30%) in Southwest and Northwest direction. After 16 weeks this becomes 217 MPa and 143 MPa (−34%). In ONK-EH3 stress interpretation, for the rock mass immediately around the hole the maximum tangential stress will be 123 MPa in North and South walls for the unsupported part and 86 MPa (−30%) for the supported part after 9 weeks. After 16 weeks the stresses are 223 MPa and 144 MPa (−35%). The support pressure generated is 3.2 MPa for 9 weeks and 5.7 MPa for 16 weeks for both interpretations. The rock mass uniaxial strength is 108 MPa and intensive damage can be expected in the unsupported part of the hole.Eight heaters produce an even heating to the concrete layer (Figure15(a)). As in the preliminary model (Figure 9) the temperature quickly drops after depth of 3 meters towards the hole bottom. At 9 weeks the concrete reaches a temperature of 129°C at −3 m level and 45°C at hole bottom −7.2 m (Figure 16(b)). The heater holes reach 185°C after 9 weeks (Figure 15(c)) and 347°C after 16 weeks. The subsequent weeks quickly heat the concrete layer until it reaches a temperature of 269°C. The model is perfectly insulated and has no losses and actual observed temperatures are expected to be lower for the concrete and the main experiment hole.Figure 16
Temperature [°C] distribution after 9 weeks of heating shown in (a) top view, (b) cross section plane, and (c) in HH4-HH8 plane.
(a) (b) (c)


### 4.3. Fracture Mechanical Simulations

In the fracture mechanics prediction, the fracturing initiates for both stress interpretations at first heating steps between weeks 0 and 3 at the concrete surface (Figure17(a) on left) on region with highest tangential stress in concrete and propagates in concrete before the first fractures initiate in rock at the third week (Figure 17(a)). In ONK-EH3 stress interpretation the failure the rock failure starts already after 3 weeks of heating and grows slowly until propagating at high rate just at the end of the experiment. The concrete failure propagates gradually during the heating period, finally surrounding the hole perimeter. The EDZ and VT1 model only suffers minor rock damage and only half of the concrete liner fails to the maximum tangential stress direction increasing with constant rate during the experiment. The displacement at the concrete and rock surface after experiment is about 0.5 mm towards to experiment hole (Figure 17(b), right). No tensile stresses exist in concrete liner although the liner has failed. According to the models, the concrete liner generates up to 3 MPa of support pressure.Figure 17
Predicted development of fracture initiation and propagation during nine weeks of heating for EDZ and VT1 stress interpretation (a) without concrete support (upper row) and (b) with concrete support (lower down). The temperature distribution is shown in top right figure and the displacement distribution and displacement vectors after nine weeks are shown on the bottom right. The temperatures and displacements are similar for both stress interpretations. The fractures are color-coded red green for tension (mode I), red for shear (mode II), and blue for boundary fractures.
(a)
Without concrete support (b)
With concrete support


## 4.1. Preliminary Model

The evenness of the heating in the surface of the hole is affected by the amount of heater holes, heater hole distance, and heater power pattern. By using a coarse mesh, it was found that using 8 vertical heaters one meter away from the hole wall produces an even heating pattern in a reasonable experiment duration. The array was rotated slightly counterclockwise to provide maximum clearance from already existing temperature monitoring and acoustic event sensor holes. The heating pattern chosen is shown in Figure6.The heat initially increases at a relatively steady rate of 12°C per week, with the rate slowing down towards week 9 (Figure9). 100°C is reached early on week 7. This leaves enough time for all the water to evaporate before being converted into steam. At the end of week 9, a temperature of 125°C is reached. Heater holes are at a temperature of 160°C at the hole top and 140°C along the top three meters. The concrete liner reaches 100°C at the beginning of 7th week and reaches 125°C (at −1 m level) at the end of 9th week (Figure 9). After the heating is turned off the temperature quickly levels out and dissipates. The seven-week cooling period brings the temperature back down to 35°C. This can further be accelerated by continuously flushing the hole with cool air.Figure 9
Calculated temperature at the center of the concrete liner (weeks 0–9).Each heater element is 6 meters long and the temperature quickly drops when approaching the hole bottom. At −6 m the maximum temperature reached at 9 weeks is 75°C and at 16 weeks is 155°C. Radially at −3 m level the temperature quickly rises towards the heater hole and quickly drops down to zero change at radius of 7.5 meters from the hole center. There is a small temperature drop (−7°C) in top part of the temperature curves which is caused by the convection and the more conductive concrete foundation slab. The pilot hole effects cannot be seen at the hole bottom, but the cooling effect of the rock mass below is clearly seen and strong (Figure9).No damage in the concrete is expected during the first three weeks (Figure10). Eurocode 2 [33] defines the characteristic strength (fck=35 MPa) as 5% failure fraction. Neglecting the effect of the intermediate stress, after five weeks, by definition, there is a 5% chance of damage and after seven weeks the mean uniaxial strength is exceeded. After nine weeks of heating, the stress will peak at 53 MPa, which is 23% higher than the mean uniaxial strength (fcm=43 MPa) of the concrete. The highest loaded region is at level −3.3 m, which is close to the −3 m monitoring level. The second planned monitoring level at −6 will peak at a much lower stress of 35 MPa, which corresponds to the characteristic uniaxial compressive strength of the concrete liner.Figure 10
Calculated tangential stress at the inner surface of the concrete liner (weeks 0–9).While mean uniaxial strengths are exceeded with a margin, the Ottosen criterion predicts that the structure is not damaged. Therefore an extended heating period was added after week 9 with heaters set to maximum (Figure6). This will exceed greatly the maximum temperature of the strain gages and some sensor loss is to be expected in the actual experiment. However, the temperature follows a predictable pattern and the lost information can be extrapolated.The concrete fails when the tangential stress increase reaches 60 MPa. Damage initiation is predicted on week 10 at −3.5 m depth. The damage continues to propagate as the heating is continued. On week 11 damage is predicted between −1 m and −5.5 m. The stress change decreases towards hole bottom and at −6 m level is only +30 MPa after 9 weeks and it reaches failure during week 13. At the hole top 0.5 m there is a region with horizontal tensile stresses up to −5 MPa.The stresses in the rock are 1.56 times higher due to the difference in elastic modulus (Figure11). The jump at 0.3 m is caused by the change of material from rock mass to concrete foundation slab. A tangential stress increase of +80 MPa is reached at the end of week 9 and the elastic increase reaches up to +190 MPa during the extended heating period. The thermally induced stress acts together with the in situ stress (Table 1). The initial tangential stress around the hole is estimated to be 46,…,54 MPa and the combined stress should be 126,…,134 MPa already on week 9. It is possible that the rock wall may sustain damage as the estimated damage strength is 58,...,102 MPa. The liner produces a support pressure, which maximizes at 2.8 MPa at depth of −3.5 m.Figure 11
Preliminary model tangential stress increase calculated at the inner surface of the rock mass (weeks 0–9).The hole moves upwards 5.2 mm and outwards 0.8 mm at top which causes vertical tensile stresses to the topmost 1.5 m of the hole. The calculated stresses exceed the design tensile strength (1.6 MPa in Finland) for concrete after week 5 but peak out at 2.1 MPa and never reach the characteristic tensile strength (2.2 MPa) level. Some tensile cracks may be observed in the top part of the hole. The maximum inward deflection is 0.3 mm at −4 m level and maximum downward displacement is 0.8 mm at −7.2 m level.


## 4.2. Full-Scale Model

The full-scale model is different from the preliminary model as it includes the true geometry with the tunnel above the hole. Additionally, the hole is lined with concrete from −3 m to −7.2 m when in preliminary model it was fully lined. The in situ stresses are included in the model and two stress state interpretations are analyzed separately. Due to the simplifications, the thermally induced stresses are slightly higher (Figure12), but the damage localization and the concentrating effects of excavation shape can now be studied more accurately.Figure 12
Full model tangential stress increase calculated at the inner surface of the rock mass in the North wall (weeks 0–9).For the EDZ and VT1 in situ stress interpretation, the initial maximum stresses are Northeast and Southwest aligned and rotate towards North-South in the heating phase (Figure13). For concrete, the tangential stresses reach 59 MPa after 9 weeks of heating and 84 MPa after 16 weeks of heating. For rock mass, the tangential stresses reach 115 MPa after 9 weeks of heating and 200 MPa after 16 weeks of heating. For the ONK-EH3 interpretation, the initial stresses are almost North-South and do not change direction during heating (Figure 14). For concrete, the tangential stresses reach 59 MPa after 9 weeks of heating and 85 MPa after 16 weeks of heating. For rock mass, the tangential stresses reach 117 MPa after 9 weeks of heating and 205 MPa after 16 weeks of heating.Figure 13
First principal stress [MPa] at −3 m level resulting from (a) 0 weeks in situ stress, (b) 9 weeks of heating, and (c) 16 weeks of heating using the EDZ and VT1 stress interpretation. The left scale is for the concrete and the right scale is for the surrounding rock mass.
(a) (b) (c)Figure 14
First principal stress [MPa] at −3 m level resulting from (a) 0 weeks in situ stress, (b) 9 weeks of heating, and (c) 16 weeks of heating using the ONK-EH3 stress interpretation. The left scale is for the concrete and the right scale is for the surrounding rock mass.
(a) (b) (c)The plastic damage was studied using a refined version of the full model. In this model the concrete mesh was the priority and all other features were modelled using a reduced degree of precision. Both in situ interpretations produce identical results for concrete. The damage initiates in Northeast segment on week 5 (Figure15(a)). The concrete exceeds the ultimate strain (0.0035) and becomes damaged on week 10 (Figure 15(b)) and continues to deteriorate in North and South until the end of the heating period reaching plastic strain of 0.017 on week 16 (Figure 15(c)). The damaged regions are the North and South quadrants with the most damage in the centers.Figure 15
Effective plastic strain of the concrete layer at −3 m level after 5 weeks (a), 10 weeks (b), and 16 weeks (c) of heating.
(a) (b) (c)For the EDZ and VT1 stress interpretation, the 9-week values are 122 MPa and 86 MPa (−30%) in Southwest and Northwest direction. After 16 weeks this becomes 217 MPa and 143 MPa (−34%). In ONK-EH3 stress interpretation, for the rock mass immediately around the hole the maximum tangential stress will be 123 MPa in North and South walls for the unsupported part and 86 MPa (−30%) for the supported part after 9 weeks. After 16 weeks the stresses are 223 MPa and 144 MPa (−35%). The support pressure generated is 3.2 MPa for 9 weeks and 5.7 MPa for 16 weeks for both interpretations. The rock mass uniaxial strength is 108 MPa and intensive damage can be expected in the unsupported part of the hole.Eight heaters produce an even heating to the concrete layer (Figure15(a)). As in the preliminary model (Figure 9) the temperature quickly drops after depth of 3 meters towards the hole bottom. At 9 weeks the concrete reaches a temperature of 129°C at −3 m level and 45°C at hole bottom −7.2 m (Figure 16(b)). The heater holes reach 185°C after 9 weeks (Figure 15(c)) and 347°C after 16 weeks. The subsequent weeks quickly heat the concrete layer until it reaches a temperature of 269°C. The model is perfectly insulated and has no losses and actual observed temperatures are expected to be lower for the concrete and the main experiment hole.Figure 16
Temperature [°C] distribution after 9 weeks of heating shown in (a) top view, (b) cross section plane, and (c) in HH4-HH8 plane.
(a) (b) (c)


## 4.3. Fracture Mechanical Simulations

In the fracture mechanics prediction, the fracturing initiates for both stress interpretations at first heating steps between weeks 0 and 3 at the concrete surface (Figure17(a) on left) on region with highest tangential stress in concrete and propagates in concrete before the first fractures initiate in rock at the third week (Figure 17(a)). In ONK-EH3 stress interpretation the failure the rock failure starts already after 3 weeks of heating and grows slowly until propagating at high rate just at the end of the experiment. The concrete failure propagates gradually during the heating period, finally surrounding the hole perimeter. The EDZ and VT1 model only suffers minor rock damage and only half of the concrete liner fails to the maximum tangential stress direction increasing with constant rate during the experiment. The displacement at the concrete and rock surface after experiment is about 0.5 mm towards to experiment hole (Figure 17(b), right). No tensile stresses exist in concrete liner although the liner has failed. According to the models, the concrete liner generates up to 3 MPa of support pressure.Figure 17
Predicted development of fracture initiation and propagation during nine weeks of heating for EDZ and VT1 stress interpretation (a) without concrete support (upper row) and (b) with concrete support (lower down). The temperature distribution is shown in top right figure and the displacement distribution and displacement vectors after nine weeks are shown on the bottom right. The temperatures and displacements are similar for both stress interpretations. The fractures are color-coded red green for tension (mode I), red for shear (mode II), and blue for boundary fractures.
(a)
Without concrete support (b)
With concrete support


## 5. Conclusions

The results from the three models, the preliminary subgeometry model, the full-scale 3D model, and the fracture mechanics 2D model, agree well and each predict damage in rock mass without the thin concrete liner and reduced damage with the liner. During the heating stage, the concrete-rock interface remains in compression and no significant shear is expected. Therefore, no adhesion loss at the interface is expected. However, the cooling period does induce tensile stresses. The cooling adhesion loss takes place after plasticity has already occurred and it could not be modelled reliably. It remains unclear what effect the loss of adhesion at the concrete-rock interface will have during cooling. The simulation results helped to design the in situ experiment and more importantly provided blind predictions, which will be compared to the observations from the actual in situ concrete spalling experiment as a part of a prediction-observation study.The thin concrete layer requires extremely small elements to account for plastic effects. For the preliminary model, it was demonstrated that the geometry could be simplified by exploiting symmetry conditions and using a 1/16th subgeometry. The results show that the concrete foundation slab and the convection of air have minimal effect on temperature and can be ignored. Based on the findings, the full model could be constructed more effectively. The full model shows 59 MPa tangential stress for the concrete liner at 9 weeks and 84 MPa at 16 weeks. For the rock mass tangential stress of 117 MPa is reached at week 9 and 205 MPa at 16 weeks. These values are well above the mean uniaxial compressive strength of concrete (43 MPa) and pegmatitic granite (102 MPa) and observable damage is expected.The Ottosen failure criterion was used in the full 3D model to include the effect of the intermediate stress. Compared to uniaxial models or Mohr-Coulomb the Ottosen predicts higher strength and more time until the ultimate strain limit of 0.0035 is reached. The top surface of the thin sprayed concrete liner was unconstrained, so the top part of the liner may be considered biaxial while the bottom part will develop a triaxial stress state. A method is provided to establish the Ottosen failure criterion parameters using only the cylindrical uniaxial compressive strength. The parameters predicted will be compared to material property measurements from the actual experiment.The fracture mechanics code Fracod2D used a modified heating pattern with similar heat load mapped to 2D. The FM modelling predicts the concrete is damaged during the first three weeks. The damage in rock is suppressed until late steps in heating week 9. The damage begins after three weeks of heating and continues as small localized damage until end of the experiment when large cracks occur and interact. There are two stress interpretations for the experiment area. The concrete is thoroughly damaged in the first interpretation while in the second interpretation the damage is localized perpendicular to the highest stress direction. It can be concluded that an anisotropic stress state causes significantly more damage than a more isotropic state.The bottom part of the hole is lined with concrete and the top part is unsupported. The COMSOL FEM full-scale model predicts support pressures up to 3.2 MPa, which is close to the support pressure of 3.0 MPa predicted by the Fracod2D fracture mechanics modelling. The presence of concrete liner reduces the tangential stress by −33%. In Äspö Pillar Stability Experiment there were indications that spalling could be controlled with approximately 150 kPa confining pressure [2] and later CAPS experiment results indicated that even support pressures of 10–30 kPa will decrease the thermal-induced spalling [6]. In the fracture mechanics modelling, the concrete liner was removed to test this hypothesis and it increased the depth of spalling significantly, thus agreeing with observations in [2, 6].



---

*Source: 4723017-2017-01-31.xml*