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

: The underground silo was constructed as a facility for the disposal of low ‐ and interme ‐ diate ‐ 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 occur ‐ ring 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 distrib ‐ uted in the lower part of the weathered rocks. A 2D axial symmetric finite element model was con ‐ sidered 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 loca ‐ tions 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.


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.

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.

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 [13]. All geomaterials were assumed to be homogeneous and isotropic.

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: where the stress invariants (p, q, and θ) are defined in Appendix A along with the expression for R θ . The function R θ 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 Figures 6 and 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, R θ is circular, representing a von Mises or Drucker-Prager failure model. When it is less than unity, R θ 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 α, β, and 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. The empirical material parameters for n = 1/2 are tabulated in Table A1 in Appendix A for several different rock types, as descriptors of rock quality. A detailed description of rock quality is presented in Table A2 in Appendix A.
where and are major and minor principal stresses at failure. = unconfined compressive strength; ∅ = internal friction angle. m, s = Hoek and Brown's material constants as tabulated in Table A1 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.

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

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.

Displacements
Exaggerated deformed shapes resulting from the excavation of the underground silo are shown in Figure 11  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. 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. 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.

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 Figures 12-14 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 lowerbound 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 Figures 11-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 Figures 15-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 Figures 6 and 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.

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 Figures 18, 19 and 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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Generalized Hoek and Brown Model Description
Appendix A. 1

. Failure Surface
The failure surface in Equation (1) in Section 3 is written in terms of the alternate stress invariants (p, q, and θ) given by where is the total stress tensor and is the deviatoric stress tensor.

Appendix A.2. Elastic Stress-Strain Relationship and Failure Surface
The incremental elastic constitutive law can be expressed in the following matrix form: σ (A2) where σ stress increment elastic stress strain matrix elastic strain increment The expression for in Equation (1)  where is a dilatancy parameter (0 1 . No plastic volume change for = 0 and associated flow for = 1.
Thus, in general: where During yielding, the consistency equation forces the stress to move along the failure surface Back-substituting Equation (A11) into Equation (A10), the stress increment is directly related to the total strain increment as follows: