Numerical Parametric Studies on the Stress Distribution in Rocks around Underground Silo

Abstract

The underground silo was constructed as a facility for the disposal of low- and intermediate-level radioactive waste. It is divided into three parts: the upper-dome core, the lower-dome core, and the cylindrical-space core. Numerical parametric studies on the stress distribution occurring in the surrounding rocks around the underground silo are presented in this paper. It is assumed that the soil layer was distributed to a depth of −4.3 m from the ground level, the weathered rocks were distributed to a depth of −9.5 m from the bottom of the soil layer, and the rocks were distributed in the lower part of the weathered rocks. A 2D axial symmetric finite element model was considered for the numerical analysis of the underground silo. A 3D finite element model was used to verify the reliability of the 2D axial symmetric model. Finite element analysis was carried out under various ratios of in situ horizontal stress to vertical stress (Ko). The numerical results obtained through these analyses include detailed stress states in the p–q and octahedral planes at key locations of finite element models around an underground silo. Contours of safety factor distributions are also presented to evaluate the overall structural safety of the surrounding rock mass, which is the main supporting body of the underground silo.

1. Introduction

In Korea, the Wolsong Low and Intermediate Level Radioactive Waste Disposal Center (WLDC) at Gyeongju in North Gyeongsang province has a plan to construct a facility with a total capacity of 800,000 drums. The WLDC site is located in the southeast of the Korean Peninsula (see Figure 1) [1]. The first phase of the facility, an underground silo for the disposal of low- and intermediate-level radioactive waste (LILW) with a capacity of 100,000 drums, was completed in 2014 (see Figure 2) [2,3,4]. Low-level radioactive waste mainly consists of items such as filters, work clothes, gloves, and replacement parts for devices used routinely in nuclear power plants. Intermediate-level radioactive waste includes radioisotope (RI) waste generated by hospitals, industries, and research institutes. Low- and intermediate-level radioactive waste must be managed safely and according to strict protocols for a certain period of time in accordance with government regulations and guidelines. Such waste is permanently disposed of in the WLDC.

It is essential to conduct a safety evaluation that comprehensively considers various types of impact on sites under construction or on the operation of nuclear facilities such as nuclear power plants, low- and intermediate-level radioactive waste disposal facilities, and spent fuel storage [5,6,7,8]. The Wolsong LILW facility’s first phase was constructed 130 m below sea level, and consists of six silos, 50 m in height and 23.6 m in diameter (see Figure 3) [3,9]. Many studies have been carried out and published to verify its safety assessment since the facility was completed [9,10].

Large-scale earthquakes occurred in 2016 and 2017 in Gyeongju and Pohang, respectively. Both arears are near the WLDC, and interest in the stability of the underground silo increased significantly [11,12]. However, unfortunately, there has been little research on the stress distribution in the rocks around the underground silos used as LILW disposal facilities [13].

Therefore, this paper presents numerical parametric studies on the stress distribution occurring in the surrounding rock mass around an underground silo which was constructed as a facility for the disposal of low- and intermediate-level radioactive waste. The 2D axial symmetric finite element model was used. The 3D finite element model was also used to verify the reliability of the 2D axial symmetric model. Finite element analyses of the underground silo were performed under various changes in the value of Ko, and the numerical results obtained through these analyses were also examined.

2. The Underground Silo

2.1. Layout of the Underground Silo

A conceptual drawing of the underground silo is presented in Figure 4 [1,14]. Its engineering barrier system is composed of a concrete silo, waste packages, disposal containers, and backfill. The walls of the silo are circular and the roof is in the form of a dome. The inside diameter is 23.6 m and the height of the wall is 35 m. The silo dome’s diameter is 30 m and its height is 17.4 m. The upper part of the underground silos is located at −80 m below sea level, and the lower part is located at −130 m below sea level.

2.2. In-Situ Stress State and Material Properties

According to the results of a detailed geological survey of the site area with a radius of 1 km centered on the WLDC site in Wolsong, the geology of the WLDC site is mainly composed of tertiary plutonic rocks, cretaceous sedimentary rocks, and intrusive rocks (see Figure 5) [15,16]. More detailed geological survey data are described in [16].

It was reported in an investigation of the site’s characteristics that the ratio of in situ horizontal stress to vertical stress (Ko) was approximately between 1.17 and 1.92 [17]. The numerical analysis presented here was therefore performed to predict the stress distribution occurring in rocks around the silo for three cases (Ko = 0.5, 1.0 and 2.0). The soil layer was assumed to be −4.3 m deep from ground level, and the weathering rock was assumed to be −9.5 m deep, starting from the bottom of the soil layer. It was assumed that the rock was distributed from the lower part of the weathering rock to the deep layer. Typical values for the material properties of geomaterials used for this finite element modeling were assumed, as shown in

Table 1. Typical material properties of geomaterials.

Ground Layer	Unit Weight (kN/m3)	Young’s Modulus (MPa)	Poisson’s Ratio	Internal Friction Angle (degree)	Cohesion (MPa)
Soil layer	18.56	0.124 × 104	0.33	30	0.35
Weathering rock	20.52	0.342 × 104	0.30	38	2.08
Rock	26.28	8.260 × 104	0.27	43	2.26

[13]. All geomaterials were assumed to be homogeneous and isotropic.

3. Theoretical Background

3.1. Material Model Description

The Generalized Hoek and Brown Model [18,19,20] was applied to this numerical analysis. This model represents the constitutive relationship of rocks or soils. The failure surface is described in its generalized form by the following equation:

$$\mathrm{F}\left(\mathrm{p}, \mathrm{q}, \mathsf{\theta }\right)=\mathrm{q}-\{{(\mathsf{\alpha }+\mathsf{\beta }\mathrm{p})}^{\mathrm{n}}+\mathsf{\kappa }\} \mathrm{R}\left(\mathsf{\theta }\right)=0$$

(1)

where the stress invariants (p, q, and $\mathsf{\theta }$) are defined in Appendix A along with the expression for $\mathrm{R}\left(\mathsf{\theta }\right)$.

The function $\mathrm{R}\left(\mathsf{\theta }\right)$ denotes the shape of the yield surface projected onto the octahedral plane. The effect of the change of parameter k on the shape of the yield surface is shown in Figure 6 and Figure 7. The value k denotes the ratio of the shear strength in the triaxial extension state to the shear strength in the triaxial compression state at the same mean pressure. It can vary from 0.5 to 1.0 and is a measure of the influence of the intermediate principal stress on the yield surface. When it is unity, $\mathrm{R}\left(\mathsf{\theta }\right)$ is circular, representing a von Mises or Drucker–Prager failure model. When it is less than unity, $\mathrm{R}\left(\mathsf{\theta }\right)$ is a smooth cornered approximation to the Mohr–Coulomb failure envelope [21,22,23]. The parameter n in Equation (1) determines the shape of the yield surface in the p–q plane.

When n = 0, the strength envelope decreases to the von Mises or Tresca yield surface, and the shear strength shows a constant value for mean pressure. When n = 1/2, the strength envelope represents a Hoek and Brown failure surface [22]. This nonlinear failure model means a multidimensional generalization of the original one-dimensional axial symmetric Hoek and Brown model, based on extensive laboratory experiments and field data [24]. When n = 1, the strength envelope in the p–q plane is representative of the Drucker–Prager or Mohr–Coulomb failure surface and shear strength is linearly proportional to the mean pressure.

The parameters $\mathsf{\alpha }$, $\mathsf{\beta },$ and $\kappa$ of Equation (1) define the failure envelope in the p–q plane as values determined in laboratory tests. Recommended relationships for determining these parameters for various types of materials are listed in

Table 2. Material constants in the Generalized Hoek and Brown Model.

	n = 0 von Mises or Tresca	n = 1/2 Hoek and Brown	n = 1 Mohr–Coulomb or Drucker–Prager
$\mathsf{\alpha }$	N/A	$\left(m^2/36+s \right)·{\sigma }_c^2$	1000
$\mathsf{\beta }$	N/A	$\mathrm{m}· {\sigma }_c$	$6\sin \varnothing /(3-\sin \varnothing )$
κ	${\mathrm{q}}^{\prime }-1$	$\left(1/6\right)·\left(m·{\sigma }_c\right)$	$[3(1-\sin \varnothing )/(3-\sin \varnothing )]·{\sigma }_c-1000$

${\mathrm{q}}^{\prime }={\sigma }_1-{\sigma }_3$ where ${\sigma }_1$ and ${\sigma }_3$ are major and minor principal stresses at failure. ${\sigma }_c$ = unconfined compressive strength; $\varnothing$ = internal friction angle. m, s = Hoek and Brown’s material constants as tabulated in Table A1 in Appendix A.

. The empirical material parameters for n = 1/2 are tabulated in

Table A1. Hoek and Brown Material Parameters (m, s).

	Rock Type					Amphibolite
		Dolomite	Mudstone		Andesite	Gabbro
		Limestone	Siltstone	Sandstone	Dolerite	Gneiss
		Marble	Shale	Quartzite	Rhyolite	Norite
Rock Quality			Slate			Quartz-Diorite
Intact
CSIR rating = 100		m = 7	10.0	15.0	17.0	25.0
NGI rating = 150		s = 1	1.0	1.0	1.0	1.0
Very Good Quality
CSIR rating = 85		3.5	5.0	7.5	8.5	12.5
NGI rating = 100		0.1	0.1	0.1	0.1	0.1
Good Quality
CSIR rating = 65		0.7	1.0	1.5	1.7	2.5
NGI rating = 10		0.004	0.004	0.004	0.004	0.004
Fair Quality
CSIR rating = 44		0.14	0.20	0.3	0.34	0.5
NGI rating = 1		0.001	0.0001	0.0001	0.0001	0.0001
Poor Quality
CSIR rating = 23		0.04	0.05	0.08	0.09	0.13
NGI rating = 0.1		0.00001	0.00001	0.00001	0.00001	0.00001
Very Poor Quality
CSIR rating = 3		0.007	0.01	0.015	0.017	0.025
NGI rating = 0.01		0.0	0.0	0.1	0.0	0.0

in Appendix A for several different rock types, as descriptors of rock quality. A detailed description of rock quality is presented in

Table A2. Description of Rock Quality in Table A1.

Intact Rock Samples	Laboratory size specimens free from joints
Very Good Quality Rock Mass	Tightly interlocking undisturbed rock with unweathered joints at 1 to 3 m
Good Quality Rock Mass	Fresh to slightly weathered rock, slightly disturbed with joints at 1 to 3 m
Fair Quality Rock Mass	Several sets of moderately weathered joints spaced at 0.3 to 1 m
Poor Quality Rock Mass	Numerous weathered joints at 30 to 500 mm with sane gouge. Clean compacted waste rock
Very Poor Quality Rock Mass	Numerous heavily weathered joints spaced <50 m with gouge. Waste rock with fines

in Appendix A.

The model includes the classical von Mises, Drucker–Prager, and Mohr–Coulomb failure equations, as well as the empirically based Hoek and Brown failure equation [25,26]. One of its useful features is that the model can use empirical data as a basis for modelling the strength of an in situ rock mass when the in situ strength data are not available. The three-dimensional elasto-plastic matrix for the Generalized Hoek and Brown Model is presented in Appendix A. The model assumes an elastic state below the failure surface, mutually dependent volumetric and deviatoric behavior once it reaches the failure surface, and a perfectly plastic state along the failure surface.

3.2. Safety Factor

The safety factor is defined by the ratio of the ultimate deviatoric stress to the current deviatoric stress at a given mean pressure, as illustrated in Figure 8. It should be noted that the ultimate deviatoric stress depends on not only the mean pressure but also on the Lode angle in the octahedral plane. The recommended content in the state range according to the safety factor is summarized in

Table 3. State range of the surrounding rock according to the safety factor.

Safety Factor (FS)	State of Surrouding Rock
FS = 1	Possibility of failure
1 < FS ≤ 2	Safe but observation required
2 < FS ≤ 10	Safe

. The contour of the safety factor distributions provides valuable information for evaluating the structural safety of the surrounding rock, which is the main supporting body of the underground silo.

4. Finite Element Modeling

There are many kinds of numerical method available to solve geomechanical problems that take into consideration the effect of nonlinearity [27,28]. In this study, finite element analysis [29,30,31,32,33] was used in the numerical modeling of the underground silo and the rock surrounding the underground silo.

One underground silo was examined in this finite element modeling procedure, although six underground silos have been constructed and are currently in operation at the Wolsong LILW disposal facility. The underground silo is divided into three parts: the upper-dome core, the lower-dome core linked to the operating tunnel, and the cylindrical-space core linked to the construction tunnel and used for storing radioactive waste packages.

The size of the analysis domain for 2D axial symmetric finite element modeling was set from the ground-surface level to a depth of 250 m from the silo floor in the vertical plane, and up to 250 m—ten times the diameter of the silo wall—to the left and right in the horizontal plane. It was also set to a thickness of 250 m for 3D analysis. The underground silo was excavated in step-by-step stages when it was under construction. However, the effects of this multi-step excavation of the underground silo site were ignored, since it is located in the hardrock zone. Thus, for this numerical study, it is assumed that the underground silo was excavated in a single step.

A 3D analysis model was constructed to verify the 2D axial symmetric analysis model. The 3D finite element model of the analysis domain, which includes the underground silo and the surrounding geomaterials, is depicted in Figure 9. The numerical results of stress and displacement obtained from carrying out a finite element analysis of this model were examined and compared at key locations around the underground silo, as shown in Figure 10.

5. Numerical Results

5.1. Displacements

Exaggerated deformed shapes resulting from the excavation of the underground silo are shown in Figure 11 for three different cases of Ko values. For the case of Ko = 0.5, the maximum displacement (0.85 mm) occurred in the vertical direction at the central point of the silo’s bottom surface (Location F) since the in situ vertical stresses were twice as high as the in situ horizontal stresses. For the case of Ko = 2.0, on the other hand, the maximum displacement (1.32 mm) occurred in the radial direction at mid-point of the cylindrical wall (Location D) since the in situ horizontal stresses were twice as high as the in situ vertical stresses. For the case of Ko = 1.0, where the in-situ stresses were hydrostatic, the maximum vertical displacement (0.786 mm) at central point of the silo’s bottom surface (Location F) were somewhat higher than the radial displacement (0.652 mm) at the mid-point of the cylindrical wall (Location D).

The displacements at key locations obtained through 2D axial symmetric analysis based on change in the ratio of in situ horizontal stress to vertical stress (Ko) are presented in

Table 4. Displacements at key locations obtained through 2D analysis (unit: mm).

Locations	Direction	Ko = 0.5	Ko = 1.0	Ko = 2.0
A	Vertical	−0.485	−0.375	−0.153
	Radial	0.0	0.0	0.0
B	Vertical	−0.318	−0.234	−0.064
	Radial	−0.076	−0.230	−0.538
C	Vertical	0.284	0.284	0.283
	Radial	−0.208	−0.540	−1.205
D	Vertical	0.087	0.070	0.038
	Radial	−0.318	−0.651	−1.316 (max)
E	Vertical	0.281	0.206	0.055
	Radial	−0.067	−0.191	−0.438
F	Vertical	0.848	0.785	0.658
	Radial	0.0	0.0	0.0

. The displacements at key locations obtained through 3D analysis based on change in the ratio of in situ horizontal stress to vertical stress (Ko) are presented in

Table 5. Displacements at key locations obtained through 3D analysis (unit: mm).

Locations	Direction	Ko = 0.5	Ko = 1.0	Ko = 2.0
A	Vertical	−0.489	−0.379	−0.158
	Radial	0.0	0.0	0.0
B	Vertical	−0.321	−0.235	−0.064
	Radial	−0.078	−0.234	−0.547
C	Vertical	0.286	0.285	0.283
	Radial	−0.209	−0.542	−1.208
D	Vertical	0.087	0.703	0.037
	Radial	−0.319	−0.652	−1.320 (Max)
E	Vertical	0.282	0.206	0.054
	Radial	−0.067	−0.191	−0.439
F	Vertical	0.850 (Max)	0.786 (Max)	0.659
	Radial	0.0	0.0	0.0

. It can be seen that numerical results of the 3D analysis were almost identical everywhere to those of the 2D axial symmetric analysis. Therefore, the reliability of the 2D finite element model was confirmed.

5.2. Stresses

This section presents the numerical results for stresses obtained by 2D axial symmetric analysis. The results of computed stresses examined at the key locations are presented in Figure 10. Nine different lines were selected around the silo: a vertical upward line, A1, above dome crown; a vertical upward line, B2; a diagonal line, B3; a horizontal line, B4; a horizontal line, C5, at the dome’s bottom; a horizontal line, D6, at the middle of the storage wall; a horizontal line, E7; a vertical downward line, E8, at the silo’s bottom corner; and a vertical downward line, F9, at the bottom center of the silo. Among all these locations, it was found that deviatoric stresses were most concentrated around two locations: C and D. Therefore, in the following analysis, stress results will be interpreted along the two corresponding lines, C5 and D6.

In Figure 12, Figure 13 and Figure 14, the stress state in the p–q plane is summarized for Ko values of 0.5, 1.0, and 2.0, along with the strength envelopes corresponding to the triaxial compression and extension modes. Using the given material properties ($\varnothing =43$ and $\mathrm{c}=2.26 \mathrm{MPa}$, these two upper- and lower-bound strength envelopes were computed based on the following equations:

$${\mathrm{q}}_{\mathrm{TXC}}=[6\sin \varnothing /(3-\sin \varnothing )]· \mathrm{P}+[6\cos \varnothing /\left(3-\sin \varnothing \right)]· \mathrm{c}=1.765 \mathrm{P}+4.27$$

$${\mathrm{q}}_{\mathrm{TXE}}=[\left(3-\sin \varnothing \right)/(3+\sin \varnothing )]·{\mathrm{q}}_{\mathrm{TXC}}=0.63{ \mathrm{q}}_{\mathrm{TXC}}$$

Generally, the deviatoric shear stresses increased with the distance from the storage wall, but the mean pressures remained more-or-less constant, except for the first three points along line C5 for Ko = 0.5.

For Ko = 0.5, all deviatoric stresses were somewhat below the lower-bound strength envelope in the triaxial extension mode. For Ko = 1.0, deviatoric stress in the first element in line C5 reached the lower-bound strength in the triaxial extension mode. For Ko = 2.0, the deviatoric stresses in the first two elements in lines C5 and D6 were above the lower-bound strength envelope and the first element in line C5 almost reached the upper-bound shear strength in the triaxial compression mode.

The plots in the p–q plane shown in Figure 11, Figure 12 and Figure 13 provide a good evaluation of the stress states along lines C5 and D6 by comparing the upper- and lower-bound strength envelopes. However, these plots do not show the ultimate strength corresponding to each stress point since the ultimate strengths are dependent on both mean pressure and Lode angle, as explained in the material model in Section 3 [34,35]. Therefore, in order to evaluate more accurately, an additional plot showing the stress state in the octahedral plane is necessary.

The stress states in the p–q and t–q planes are summarized in Figure 15, Figure 16 and Figure 17 for Ko values of 0.5, 1.0 and 2.0, respectively. The t–q plane means a θ–q plane failure surface. The t–q plot, as introduced in this study, is essentially the simpler form of deviatoric stress plot in the octahedral plane shown in Figure 6 and Figure 7.

For Ko = 0.5, the stress state was somewhat close to the triaxial compression mode at Location C, but was closer still to the triaxial extension mode at Location D. The safety factors, as defined in Figure 17 in the next section, were 2.61 and 1.58 at Locations C and D, respectively.

For Ko = 1.0, the stress state was very close to the triaxial compression mode at Location C, but was very close to the deviatoric pure shear mode at Location D, with a Lode angel of −1.9°. The safety factors were 1.47 and 1.34 at Locations C and D, respectively.

For Ko = 2.0, the stress state at Location C was almost in the triaxial compression mode, reaching ultimate strength with a safety factor of 1.02. The stress state at Location D was also close to the triaxial compression mode, with a Lode angle of −20.9°, reaching ultimate strength with a safety factor of 1.08. Thus, the silo’s storage wall may be subjected to the onset of shear failure at these regions.

5.3. Safety Factors

In this section, the safety factor is presented using the stress value obtained through 2D axial symmetric analysis. The contours of safety factor distributions provide valuable information for evaluating the structural safety of the surrounding rock mass, which is the main supporting body of the underground silo.

The safety factor distributions according to change in the ratio of in situ horizontal stress to vertical stress (Ko) after single-step excavation of an underground silo are presented in Figure 18, Figure 19 and Figure 20 for Ko values of 0.5, 1.0, and 2.0, respectively.

For Ko = 0.5, the zones with lower safety factors were more extended in the horizontal direction along the storage wall, with a minimum safety factor of 1.42. The stress state in these zones, as shown in Figure 14d, were close to the triaxial extension mode where vertical and tangential stresses were contributing as major stresses and the horizontal stress as a minor stress. Overall, the surrounding rock mass was in a safe condition.

For Ko = 1.0, the zones with lower safety factors were more-or-less uniformly concentrated around the silo, with a minimum safety factor of 1.29. The stress states in these zones, as examined in Figure 15d, were close to the deviatoric pure shear mode where the major stress was tangential stress, the intermediate stress was vertical stress, and the minor stress was horizontal stress. Overall, the surrounding rock mass was in a safe condition.

For Ko = 2.0, the zones with lower safety factors were more widely spread around the silo compared to the cases for Ko values of 0.5 and 1.0. The stress states along the storage wall, as shown in Figure 16d, were closer to the triaxial compression mode. Overall, surrounding rock mass was in a safe condition except along the storage wall, where it almost reached ultimate strength with a minimum safety factor of 1.02.

6. Conclusions

This article presented numerical parametric studies on the stress distribution in rocks around the underground silo for low- and intermediate-level radioactive waste disposal facilities in Korea. The 2D axial symmetric finite element model was used in this study, and a 3D finite element model was also used to verify the reliability of the 2D axial symmetric model. A finite element analysis of an analysis domain covering an underground silo and the surrounding rock mass was performed under three cases of differing ratios of in situ horizontal stress to vertical stress (Ko = 0.5, 1.0 and 2.0). Numerical results obtained through these analyses were presented and examined in detail for displacements, stresses, and safety factors associated with the excavation of an underground silo. To evaluate more accurately the stress state of the surrounding rock, a t–q plot showing the stress state in the octahedral plane was created, in addition to a p–q stress plot. The safety factor was also presented, using the stress value obtained through numerical analysis.

The following conclusions arise from the numerical parametric study presented in this paper.

(1)

The numerical results of the 3D analysis were almost identical in every case to those of the 2D axial symmetric analysis, confirming the reliability of the 2D finite element model.

(2)

For the case of Ko = 0.5, the maximum displacement occurred in the vertical direction at the central point of the silo’s bottom surface, since the in situ vertical stresses were twice as high as the in situ horizontal stresses. For the case of Ko = 2.0, on the other hand, the maximum displacement occurred in the radial direction at the mid-point of the cylindrical wall, since the in situ horizontal stresses are twice higher than the in situ vertical stresses. For the case of Ko = 1.0, where the in-situ stresses were hydrostatic, the maximum vertical displacement at the central point of the silo’s bottom surface was somewhat higher than the radial displacement at the mid-point of the cylindrical wall.

(3)

For Ko = 0.5, all deviatoric stresses were somewhat below the lower-bound strength envelope in the triaxial extension mode. For Ko = 1.0, the deviatoric stress in the first element in line C5 reached the lower-bound strength in the triaxial extension mode. For Ko = 2.0, the deviatoric stresses in the first two elements in lines C5 and D6 were above the lower-bound strength envelope and the first element in line C5 almost reached the upper-bound shear strength in the triaxial compression mode.

(4)

For Ko = 0.5, the stress state was somewhat close to the triaxial compression mode at Location C, but was closer to the triaxial extension mode at Location D. The safety factors were 2.61 and 1.58 at Location C and D, respectively. For Ko = 1.0, the stress state was very close to the triaxial compression mode at Location C, but was very close to the deviatoric pure shear mode at Location D, with a Lode angel of −1.9°. The safety factors were 1.47 and 1.34 at Location C and D, respectively. For Ko = 2.0, the stress state at Location C was almost in the triaxial compression mode, reaching ultimate strength with a safety factor of 1.02. The stress state at Location D was also close to the triaxial compression mode, with a Lode angle of −20.9°, reaching ultimate strength with a safety factor of 1.08. Thus, the silo’s storage wall may be subjected to the onset of shear failure at these regions.

(5)

It is better to compare the numerical results with in situ monitoring results. However, unfortunately, in situ monitoring results have not been reported anywhere. Any future data that is acquired should be compared with the results of this study.