Journal of Engineering
Volume 2017 (2017), Article ID 4723017, 16 pages
https://doi.org/10.1155/2017/4723017

Research Article

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

Lauri Kalle Tapio Uotinen and Topias Kalle Aleksi Siren

Aalto University, Rakentajanaukio 4 A, 02150 Espoo, Finland

Correspondence should be addressed to Lauri Kalle Tapio Uotinen; lauri.uotinen@aalto.fi

Received 3 November 2016; Revised 26 December 2016; Accepted 27 December 2016; Published 31 January 2017

Academic Editor: İlker Bekir Topçu

Copyright © 2017 Lauri Kalle Tapio Uotinen and Topias Kalle Aleksi Siren. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 (Figure 1). 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 (Figure 1) 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].

(a)

(b)

Figure 2: The end of the POSE niche (a) and the main dimensions of the experiment (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 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 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 criterion or von Mises criterion . The criterion is given by

where

, , , and 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 , a failure state on the compressive meridian , uniaxial tensile strength , and biaxial compressive strength . 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 Figure 2(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
[MPa]
Dip direction [°]
Dip [°]
[MPa]
Dip direction [°]
Dip [°]
[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
85

Figure 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 -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
[MPa]
Dip direction [°]
[MPa]
Dip direction [°]
[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 () 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 Table 3 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 conductivity [W/m·K]
Density [kg/m³]
Heat capacity [J/kg·K]
Poisson’s ratio
Young’s modulus [GPa]
Coefficient of thermal expansion [1/K]

Rock mass (3D)

Pegmatite (2D)

Concrete C35/45

Insulation

[26]. [27]. [32]. [28]. [29]. Based on quality assurance tests for sprayed concrete. [30].

Manufacturer quote. Arbitrary value. Insufficient amount of test results for PGR. Property for rock mass. Property for veined gneiss. Property 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

Concrete C35/45

[31]. [34]. Calculated value [33]. Property for pegmatitic granite. Property for rock mass.

Table 5: Ottosen plasticity criterion parameters for C35/45 concrete.

Material

Concrete C35/45

Calculated value [33]. Calculated 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

MPa

Friction angle

°

Tensile strength

MPa

Mode I fracture toughness

MPa

Mode II fracture toughness

MPa

Cohesion, tensile & shear

MPa

Friction angle, tensile & shear
,

°, °

Dilatation angle, tensile & shear
,

,
°, °

Normal stiffness, tensile & shear

GPa/m

Shear stiffness, tensile & shear

GPa/m

[34]. Modified after [35]. [32]. [33]. Derived by [5] after [36] and [37]. Calculated after relation by [38].

2.4. Initial Values

The initial values for convection and for thermal calculation initial values are shown in Table 7. 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

(13.5)
K [°C]

Air pressure

(1.05)
kPa [bar]

Rock mass near field temperature

(18)
K [°C]

Rock mass far field temperature

(11.5)
K [°C]

Concrete temperature

(14)
K [°C]

[26]. Air pressure at −345 m. Estimation 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 Figure 3 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 Figure 6. 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. Figure 7(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.

(a)

(b)

(c)

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).

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 (Figure 7(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 (Figure 8(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.

(a)

(b)

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).

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 Figure 6.

The heat initially increases at a relatively steady rate of 12°C per week, with the rate slowing down towards week 9 (Figure 9). 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 (Figure 9).

No damage in the concrete is expected during the first three weeks (Figure 10). Eurocode 2 [33] defines the characteristic strength ( 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 ( 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 (Figure 6). 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 (Figure 11). 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 (Figure 12), 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 (Figure 13). 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.

(a)

(b)

(c)

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.

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 (Figure 15(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.

(a)

(b)

(c)

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.

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 (Figure 15(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.

(a)

(b)

(c)

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.

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 (Figure 17(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.

(a) Without concrete support

(b) With concrete support

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.

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].

Appendix

A. Ottosen Four-Parameter Failure Criterion for Concrete

A.1. Failure States

To calibrate the four material parameters , , , and , four arbitrary failure stress states must be obtained from multiaxial testing. These can be chosen for convenience compressive strength , a failure state on the compressive meridian , uniaxial tensile strength , and biaxial compressive strength . The two first states are located on the compressive meridian () and the two latter states are located on the tensile meridian ().

A.2. Equations

As shown in [14], using the meridian equations and the four failure states, a set of equations can be reached where

The failure criterion is then

where the function is defined by

A.3. Approximated Empirical Values for Concrete

Sometimes not enough parameters are known. When the concrete grade is within C20–C50 ( MPa) the following approximation is valid

where

Based on [22, 23] the biaxial compressive strength can be approximated as

For the failure state on the compressive meridian, based on [24, 25], it can be assumed that

A.4. Eurocode Based Calculated Values

Using the compressive and tensile strength in Eurocode 2 [33] and the approximations for biaxial strength and compressive meridian above, values for Eurocode concrete strengths can be calculated and are shown in Tables 8 and 9.

Table 8: Ottosen parameters for Eurocode concrete strength classes C12–C40.

MPa
12
16
20
25
30
35
40

MPa
15
20
25
30
37
45
50

MPa
20
24
28
33
38
43
48

MPa
1.57
1.90
2.21
2.56
2.90
3.21
3.51

1.8654
1.8398
1.8544
1.8967
1.9508
2.0098
2.0704

4.1943
4.1510
4.1756
4.2473
4.3390
4.4388
4.5414

14.6958
14.5640
14.6389
14.8571
15.1362
15.4402
15.7528

0.9907
0.9903
0.9905
0.9911
0.9918
0.9925
0.9932

Table 9: Ottosen parameters for Eurocode concrete strength classes C45–C90.

MPa
45
50
55
60
70
80
90

MPa
55
60
67
75
85
95
105

MPa
53
58
63
68
78
88
98

MPa
3.80
4.07
3.97
4.13
4.41
4.66
4.88

2.1311
2.1911
2.5393
2.6690
2.9273
3.1833
3.4362

4.6441
4.7456
5.3351
5.5546
5.9920
6.4252
6.8533

16.0657
16.3749
18.1719
18.8415
20.1757
21.4978
22.8047

0.9938
0.9943
0.9967
0.9974
0.9984
0.9990
0.9995

Disclosure

The research was carried out without involvement of the funding sources.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was partially funded by Aalto Doctoral Programme in Engineering, Posiva Oy, and YTERA Doctoral Programme for Nuclear Engineering and Radiochemistry and the Academy of Finland (no. 297770). The authors thank Matti Hakala (Rock Mechanics Consulting Finland) for advice concerning the thermomechanical modelling, Ari Hartikainen (Aalto University) for technical help in running the preliminary calculations, Johanna Tikkanen (TUKES) for advice regarding concrete thermomechanical properties, and Pasi Marttila (COMSOL) for support with plastic thermomechanically fully coupled direct transient modelling.

References

SKBF/KBS. KBS 3—Final storage of spent nuclear fuel—KBS-3, I general, Art716-1, 1983, http://www.skb.com/publication/4397/

J. C. Andersson, “Äspö Hard Rock Laboratory. Äspö pillar stability experiment, final report. Rock mass response to coupled mechanical thermal loading,” SKB Technical Report TR-07-01, Swedish Nuclear Fuel and Waste Management, Stockholm, Sweden, 2007, http://www.skb.com/publication/1353659/TR-07-01.pdf.

J. Valli, M. Hakala, T. Wanne, P. Kantia, and T. Siren, “ONKALO POSE experiment—phase 3: execution and monitoring,” Working Report 2013-41, Posiva, Eurajoki, Finland, 2014, http://www.posiva.fi/files/3646/WR_2013-41.2.pdf.

L. Uotinen, T. Siren, D. Martinelli, and M. Hakala, “In-situ experiment concerning thermally induced spalling of circular shotcreted shafts in deep crystalline rock,” in Proceedings of the World Tunnel Congress (WTC '13), Geneva, Switzerland, May 2013, https://research.aalto.fi/files/10378012/WTC2013_full_paper_v15.pdf.

T. Siren, L. Uotinen, M. Rinne, and B. Shen, “Fracture mechanics modelling of an in situ concrete spalling experiment,” Rock Mechanics and Rock Engineering, vol. 48, no. 4, pp. 1423–1438, 2015.

R. Glamheden, B. Fälth, L. Jacobsson, J. Harrström, J. Berglund, and L. Bergkvist, “Counterforce applied to prevent spalling,” SKB Technical Report TR-10-37, Swedish Nuclear Fuel and Waste Management, Stockholm, Sweden, 2010, http://www.skb.se/upload/publications/pdf/TR-10-37.pdf.

E. Johansson, T. Siren, M. Hakala, and P. Kantia, “ONKALO POSE experiment—phase 1 & 2: execution and monitoring,” Working Report 2012-60, Posiva, Eurajoki, Finland, 2014, http://www.posiva.fi/files/3641/WR_2012-60.pdf.

B. Shen, O. Stephansson, and M. Rinne, Modelling Rock Fracturing Processes: A Fracture Mechanics Approach Using FRACOD, Springer, Berlin, Germany, 2014.

B. Shen and O. Stephansson, “Modification of the G-criterion for crack propagation subjected to compression,” Engineering Fracture Mechanics, vol. 47, no. 2, pp. 177–189, 1994.

O. Mohr, “Welche umstände bedingen die elastizitätsgrenze und den bruch eines materials,” Zeitschrift des Vereins Deutscher Ingenieure, vol. 46, no. 1524–1530, pp. 1572–1577, 1900.

C. A. Coulomb, Essai sur une Application des Règles de Maximis & Minimis à Quelques Problèmes de Statique, Relatifs à l'Architecture, De l'Imprimerie Royale, 1776.

D. C. Drucker and W. Prager, “Soil mechanics and plastic analysis or limit design,” Quarterly of Applied Mathematics, vol. 10, no. 2, pp. 157–165, 1952.

N. S. Ottosen, “A failure criterion for concrete,” Journal of the Engineering Mechanics Division, ASCE, vol. 103, no. 4, pp. 527–535, 1977.

N. S. Ottosen and M. Ristinmaa, The Mechanics of Constitutive Modeling, Elsevier, Oxford, UK, 2005.

M. Zyczkowski, Combined Loadings in the Theory of Plasticity, Springer Science & Business Media, 1981.

E. Hoek and E. T. Brown, “Empirical strength criterion for rock masses,” Journal of the Geotechnical Engineering Division, ASCE, vol. 106, no. 15715, pp. 1013–1035, 1980.

E. Hoek and E. T. Brown, Underground Excavations in Rock, Instn Min. Metall, London, UK, 1980.

P. Menetrey and K. J. Willam, “Triaxial failure criterion for concrete and its generalization,” ACI Structural Journal, vol. 92, no. 3, pp. 311–318, 1995.

K. J. Willam and E. P. Warnke, “Constitutive models for the triaxial behavior of concrete,” in Proceedings of the International Association for Bridge and Structural Engineer Seminar on Structures Subjected to Triaxial Stresses, vol. 19, pp. 1–31, Bergamo, Italy, 1974.

B. Bresler and K. S. Pister, “Strength of concrete under combined stresses,” ACI Journal Proceedings, vol. 55, no. 9, pp. 321–345, 1958.

W. F. Chen and D. J. Han, Plasticity for Structural Engineers, Springer Science & Business Media, New York, NY, USA, 2012.

H. Kupfer, H. K. Hilsdorf, and H. Rusch, “Behavior of concrete under biaxial stresses,” Journal of American Concrete Institute, vol. 66, no. 8, pp. 656–666, 1969.

G. Shickert and H. Winkler, Results of Tests Concerning Strength and Strain of Concrete Subjected to Multiaxial Compressive Stresses, Deutscher Ausschuss für Stahlbeton, Berlin, Germany, 1977.

G. G. Balmer, Shearing Strength of Concrete Under High Triaxial Stress-Computation of Mohr's Envelope as a Curve, Branch of Design and Construction, US Bureau of Reclamation, 1949.

F. E. Richart, A. Brandtzaeg, and R. L. Brown, “A study of the failure of concrete under combined compressive stresses,” Engineering Experiment Station 185, Bulletin, 1928.

M. Hakala and J. Valli, “ONKALO POSE experiment—phase 3: 3DEC prediction,” Working Report 2012-58, Posiva, Eurajoki, Finland, 2013, http://www.posiva.fi/files/3160/WR_2012-58.pdf.

I. Kukkonen, L. Kivekäs, S. Vuoriainen, and M. Kääriä, “Thermal properties of rocks in Olkiluoto: results of laboratory measurements 1994–2010,” Working Report 2011-17, Posiva, Eurajoki, Finland, 2011, http://www.posiva.fi/files/1504/WR_2011-17_web.pdf.

U. Åkesson, “Laboratory measurements of the coefficient of thermal expansion of olkiluoto drill core samples,” Working Report 2012-14, Posiva, Eurajoki, Finland, 2012, http://www.posiva.fi/files/2825/WR_2012-14web.pdf.

SRMK C4 National Building Code of Finland: Thermal Insulation: Guidelines 2003, Ministry of the Environment, Helsinki, Finland, 2003, http://www.ym.fi/download/noname/%7B09D023C7-7562-449E-86A3-859BD06C7BB9%7D/31549.

A. M. Neville, Properties of Concrete, Longman Group, London, UK, 1995.

Posiva, “Olkiluoto site description 2011,” Posiva Report 2011-02, Posiva, Eurajoki, Finland, 2012, http://posiva.fi/files/3002/POSIVA_2011-02.pdf.

Posiva, “Olkiluoto site description 2008,” Posiva Report 2009-01, Posiva, Eurajoki, Finland, 2009, http://www.posiva.fi/files/1035/Posiva2009-01_web_part1.1rev.pdf.

EN, “Eurocode 2: design of concrete structures. Part 1-1: general rules and rules for buildings,” EN 1992-1-1, European Committee for Standardization, CEN, Brussels, Belgium, 2004.

T. Siren, “Fracture mechanics prediction for Posiva's Olkiluoto spalling experiment (POSE),” Working Report 2011-23, Posiva, Eurajoki, Finland, 2011, http://www.posiva.fi/files/1703/WR_2011-23_web.pdf.

T. Siren, “Fracture toughness properties of rocks in Olkiluoto: laboratory measurements 2008-2009,” Working Report 2012-25, Posiva, Eurajoki, Finland, 2012, http://www.posiva.fi/files/2854/WR_2012-25web.pdf.

J. Davies, “Numerical study of punch-through shear specimen in mode II testing for cementitious materials,” International Journal of Cement Composites and Lightweight Concrete, vol. 10, no. 1, pp. 3–14, 1988.

R. N. Swamy, “Fracture mechanics applied to concrete,” in Developments in Concrete Technology, Part I, F. D. Lydon, Ed., pp. 221–281, Applied Science, London, UK, 1979.

H. W. Reinhardt, J. Ošbolt, X. Shilang, and A. Dinku, “Shear of structural concrete members and pure mode II testing,” Advanced Cement Based Materials, vol. 5, no. 3-4, pp. 75–85, 1997.

S. L. Crouch, “Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution,” International Journal for Numerical Methods in Engineering, vol. 10, no. 2, pp. 301–343, 1976.

Y. J. Liu and Y. X. Li, “Revisit of the equivalence of the displacement discontinuity method and boundary element method for solving crack problems,” Engineering Analysis with Boundary Elements, vol. 47, pp. 64–67, 2014.