Research Article

Tunnel and Underground Space. 31 December 2022. 568-585
https://doi.org/10.7474/TUS.2022.32.6.568

ABSTRACT


MAIN

  • 1. 서 론

  • 2. 연구 배경

  •   2.1 DECOVALEX-2023 TASK G

  •   2.2 한국건설기술연구원 Thermoshearing 실험

  • 3. 수치모델 및 해석개요

  • 4. 해석 결과

  •   4.1 온도 분포

  •   4.2 응력 분포

  •   4.3 변위 분포

  • 5. 요약 및 결론

1. 서 론

원자력 발전으로 인해 생성되는 고준위 방사성폐기물을 인간 생활권으로부터 안전하게 격리하기 위해, 현재 가장 유력하게 고려되는 방식은 심지층 처분방식이다(IAEA, 2011). 심지층 처분은 지하의 수백 미터 심도의 균열 밀도가 낮은 안정된 암반에 공학적방벽(engineered barrier system, EBS)과 자연방벽(natural barrier system, NBS)으로 구성된 다중방벽 시스템(multi-barrier system)을 구축하고 폐기물을 영구저장하는 방식으로서, 스웨덴의 Swedish Nuclear Fuel and Waste Management Company (SKB)가 제안한 심층처분시스템(Kärnbränslesäkerhet, KBS)이 대표적이다(SKB, 2010). 국내의 경우 한국원자력연구원이 위 처분개념을 토대로 한국형 기준 처분시스템(Korean Reference HLW Disposal System, KRS)을 제안한 바 있으며(Lee et al., 2007), 2021년에 과학기술정보통신부, 산업통상자원부, 원자력안전위원회가 공동으로 설립한 사용후핵연료관리핵심기술사업단(Institute for Korea Spent Nuclear Fuel, iKSNF)을 주축으로 사용후핵연료 저장 및 처분에 관한 활발한 연구개발이 진행되고 있다(iKSNF, 2022).

심부 지층은 고온, 고압, 고응력의 극한 조건 하에 있고, 제한된 접근성으로 인해 그 거동 특성을 평가하고 예측하는 것이 매우 어렵다. 또한, 암반의 열적, 수리적, 역학적 거동이 불연속면에 의해 지배적인 영향을 받으므로 그 불확실성은 더욱 증가한다. 따라서 10만년 이상의 장기적 안정성과 성능 확보를 목표로 하는 처분장 설계에 있어 신뢰성 높은 수치모델의 개발과 활용은 매우 중요한 기술요소 중 하나이다.

처분장에서 암반 응력 조건의 변화를 야기하는 주요 원인은 터널과 처분공의 굴착 시 발생하는 응력 재분포와 사용후핵연료의 붕괴열로 인한 열응력의 발달이다(Min et al., 2013). 처분장 근계영역(near-field area)에서는 굴착으로 인해 새로운 균열이 생성되거나 기존 균열이 성장하여 열, 수리, 역학적인 물성 변화가 발생할 수 있고, 그 영향은 심도가 깊고 현지응력이 클수록 크게 나타나는 것으로 알려져 있다(Malmgren et al., 2007). 암반 온도가 상승하면 수직방향으로는 열팽창이 가능하지만, 수평방향으로는 변형이 구속되므로 방향에 따라 상이한 열응력이 발달한다. 이에 따라 축차응력(differential stress)과 전단력 증가에 의한 균열 미끄러짐이 발생할 수 있다(Rutqvist, 2020). 암반 내 많은 단층이 임계응력상태(critical stress state)에 있음(Zoback and Gorelick, 2012)을 고려하면, 원계영역(far-field area)에서도 열전달에 의해 발생한 작은 크기의 응력 증분만으로 단층 재활성(fault reactivation)을 유발할 수 있다. 특히, 처분장 부지로서 결정질 암반을 고려하는 경우, 암반 불연속면 거동을 합리적으로 반영할 수 있는 열-수리-역학적 연계거동 해석기법의 개발이 필수적이다.

본 연구에서는 입자기반 개별요소모델(grain-based distinct element model, GBDEM)을 이용하여, 열에 의한 암석 균열의 미끄러짐 거동(fracture thermoshearing or thermally induced fracture slip)을 수치해석적으로 살펴보았다. 이는 DECOVALEX-2023 국제공동연구 Task G의 일환으로 수행된 연구로서, 해석대상은 한국건설기술연구원에서 수행된 암석 균열에 대한 열-역학적 하중 재하 실험(Sun et al., 2022)이다. 이하 본문에서는 위 실험을 ‘KICT-Thermoshearing 실험’으로 지칭하였다. 본 연구는 상기 해석모델을 통해 가열로 인한 암석 표면의 온도 분포, 열응력의 증가에 따른 주응력 변화, 균열의 전단변위와 수직변위 등을 모델링하고 이를 실험 결과와 비교함으로써, 해석모델의 암석 균열 thermoshearing 거동에 대한 적용성을 검토하고 모델을 개선하는 데에 주요 목적이 있다.

2. 연구 배경

2.1 DECOVALEX-2023 TASK G

DECOVALEX 프로젝트는 방사성폐기물의 지층처분 관련 수치해석 코드 개발 및 검증을 위한 국제공동연구로 1992년에 스웨덴의 Swedish Nuclear Power Inspectorate(SKI)와 미국의 Lawrence Berkeley National Laboratory(LBNL)을 중심으로 시작되어, 현재까지 지층처분을 고려하는 대부분의 국가가 참여하고 있다. DECOVALEX 프로젝트에서는 4년을 주기로 처분 관련 THMC 문제에 대한 새로운 연구 주제(Task)들을 선정하여 각 주제에 대한 공동연구가 진행되어 왔으며, 현재 2020년 4월에 착수된 DECOVALEX-2023(2020-2023)에서 7가지 연구 주제에 대해 공동연구를 수행하고 있다(Kim et al., 2021).

본 연구는 DECOVALEX-2023의 Task G에 참여하여 수행된 연구로서, Task G의 연구목표는 결정질 암석 내 균열의 열-수리-역학적 복합거동을 해석할 수 있는 수치해석기법을 개발, 비교·검증하는 것이다. 독일 Freiberg University of Mining and Technology(TUBAF), Helmholtz Centre for Environmental Research(UFZ), 스코틀랜드 Edinburgh 대학교, 독일 DyanFrax 사에서 공동 제안하였으며 총 9개의 연구팀이 참여하고 있다. Task G에 대한 보다 상세한 내용은 Park et al.(2020), Kim et al.(2021)에 보고된 바 있다.

2.2 한국건설기술연구원 Thermoshearing 실험

Task G에서는 Edinburgh 대학교에서 수행된 GREAT Cell 수리전단시험(McDermott et al., 2018)과 한국건설기술연구원에서 수행된 KICT-Thermoshearing 실험(Sun et al., 2021, Sun et al., 2022)을 대상으로 암석 균열의 열-수리-역학적 복합거동을 해석할 수 있는 수치모델을 개발하고, 벤치마크 모델링, 참여팀 간 결과 비교, 실험결과 비교 등을 통해 개발한 수치모델을 검증하게 된다.

KICT-Thermoshearing 실험은 암석 균열(saw-cut fracture, laser-marketing fracture, tensile fracture)에 하중 및 열을 가하고 열-역학적 요인에 의한 균열의 미끄러짐 거동을 살펴본 실험이다. 2022년 12월 현재, Task G에 참여하는 연구팀들은 벤치마크 모델링을 대부분 완료하고 수치모델의 적용성을 검토하는 단계에 있으며, 연구팀 간의 결과 비교 및 수치모델 검증을 위한 구체적 실험결과는 아직 배포되지 않은 상태이다.

본 연구에서는 그 예비연구로서 Sun et al.(2022)이 발표한 연구내용 중 saw-cut 균열에 대한 실험을 모사하고자 하였다. 실험에 사용된 암석은 국내에서 채집한 포천화강암으로 길이 100 mm의 정육면체 형상이며 시료 내부에 42° 경사의 saw-cut 균열을 포함하고 있다. Fig. 1은 시험편의 형상을 보여주며, Table 1은 시험편의 물성을 정리한 것이다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F1.jpg
Fig. 1

Test specimen: a single saw-cut fracture in Pocheon granite (modified from source: Sun et al., 2021, Sun et al., 2022)

Table 1.

Physical and mechanical properties of rock and fracture (Sun et al., 2022)

Properties Value Unit
Model size (length × height × width ) 0.1 × 0.1 × 0.1 m3
Rock density 2600 kg/m3
Rock Young’s modulus 55.09 GPa
Rock Poisson’s ratio 0.275 -
Rock internal friction angle 60.4 °
Rock cohesion 27.58 MPa
Rock thermal expansion coefficient 6.24 × 10-6 1/°C
Rock thermal conductivity 1.96 W/m/K
Rock specific heat 645 J/kg/K
Fracture tensile strength 0.1 MPa
Fracture dip angle, β 42 °
Fracture static friction coefficient, μs 0.58 -

위 실험은 크게 역학적 하중 재하와 열-역학적 하중 재하 단계로 구성된다. Fig. 2는 실험에 적용된 경계조건을 보여주는 것으로 역학적 하중 재하 단계에서는 수평방향과 수직방향으로 서서히 응력을 가하여 목표한 주응력 상태(σmin = 3 MPa, σmax = 9 MPa)를 재현한다. 시료에 가해진 주응력은 균열 전단강도와 경사를 고려하여 균열이 임계응력 상태에 놓이도록 설정된 값이다. Sun et al.(2022)은 saw-cut 균열 이외에도 거친 균열(laser marketing 균열, 인장 균열)을 사용하였는데, 임계응력 상태를 구현하기 위해 σmin은 3 MPa로 모두 동일하게 적용하되, 거칠기에 따라 laser marking 균열에 대한 σmax는 11.11 MPa로, 인장 균열에 대한 σmax는 24.61 MPa로 상이하게 적용하였다.

열-역학적 하중 재하 단계에서는 시료의 상부와 하부에 설치된 히터를 통해 시료를 가열하면서 균열의 전단거동을 모니터링한다. 가열 중에는 수직방향으로 일정 응력(σmin)을 유지하되 수평방향 변위를 구속하여 열응력(σT)으로 인한 수평 주응력 증가를 유도하게 된다. Fig. 2에서 시료의 전면(front)과 후면(back)은 응력이 작용하지 않는 자유면(free surface)이며, 모든 경계면에서 단열을 위한 별도의 장치는 고려되지 않았다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F2.jpg
Fig. 2

Boundary conditions during mechanical loading test and thermo-mechanical loading test (modified from source: Sun et al., 2022)

실험 중에는 균열의 전단 거동을 모니터링하기 위하여 시간에 따른 변위, 변형률, 온도, 미소파괴음(acoustic emission, AE) 등이 측정되었다. Fig. 3은 측정 장치 및 센서 위치를 보여주는 모식도이다. 상부와 하부의 가압판(loading plate)에는 가열을 위하여 4개의 히터가 설치되었다. 시험편 표면에 부착된 6개(하부면 1개, 전면 5개)의 센서를 통해 온도 변화를 측정하였고, 좌측면, 우측면, 전면에는 10개의 센서를 부착하여 AE를 모니터링하였다. 한편, 균열 변위를 측정하기 위해 2개의 변위 변환기(clip-on displacement transducer)와 디지털 이미지 상관관계(digital image correlation, DIC) 분석 기법이 활용되었다. 시험편 전면에 균열을 교차하는 2개의 변위 변환기를 부착하였으며, 후면에서는 실시간으로 이미지를 촬영하고 균열 주변 5개 지점에서의 변위를 DIC 기법을 통해 분석하였다. 실험장비 및 모니터링과 관련된 자세한 내용은 Sun et al.(2022)의 논문에 기술되어 있다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F3.jpg
Fig. 3

Schematic diagram showing the locations of heaters, sensors, and gauges on the six sides of the cubic granite specimen with a single fracture (Sun et al., 2022)

Fig. 4는 열-역학적 하중 재하 단계에서 측정된 온도 변화 곡선을 보여주는 것으로 시험편에 표기된 숫자는 모니터링 지점을 나타낸 것이다. 가열은 상부와 하부의 가압판(로드셀)에 설치된 히터를 통해 약 6,000초 동안 이루어졌다. 초기 1,000초 동안 히터의 온도를 서서히 증가시켜 150°C에 도달하면 실험이 종료될 때까지 온도를 유지하였다. Fig. 4에서 1로 표현된 온도변화 곡선은 시험편 하부면과 가압판(loading plate) 사이에 부착된 1번 온도센서의 측정 결과로서 실험 종료 시 약 106°C에 수렴하는 것으로 나타났다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F4.jpg
Fig. 4

Temperature distributions and revolutions at six measuring points during thermoshearing test (modified from source: Sun et al., 2022)

Fig. 5는 시간에 따른 수평응력 증가량(stress increment)과 균열 변위를 보여주는 것으로, 균열 변위는 변위 변환기를 통해 측정된 값이다. Fig. 6은 DIC 기법을 통해 5개 지점에서 분석된 균열 전단변위를 변위 변환기를 통해 측정된 결과와 비교한 그림이다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F5.jpg
Fig. 5

Maximum principal stress increment and fracture shear and normal displacements (Sun et al., 2022)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F6.jpg
Fig. 6

Comparison of fracture shear displacements during the thermoshearing measured by the clip-on displacement transducer and the DIC technique (Sun et al., 2022)

3. 수치모델 및 해석개요

본 연구에서는 입자기반 개별요소모델(grain-based distinct element model, GBDEM)을 이용하여 제2.2절에서 요약한 saw-cut 균열 시료에 대한 실험 결과를 모사하였다. GBDEM은 암석의 구조를 다각형 또는 다면체 입자의 결합체로 가정하고, 입자와 입자 간 접촉에서 발생하는 역학적 거동을 개별요소법에 의해 계산하는 방법이다. GBDEM을 암석거동 해석에 활용한 다양한 연구들은 위 모델이 암석의 미시구조와 점진적 파괴 메커니즘을 직접적으로 모방함으로써 암석 재료의 불균질성과 취성거동, 이방성 등을 효과적으로 모사할 수 있음을 보고하였다(Lan et al., 2010, Ghazvinian et al., 2014, Park et al., 2017, Fu et al., 2020, Wang et al., 2021). 본 연구의 선행연구인 Park et al.(2020), Park et al.(2021a), Park et al.(2021b)의 연구에서는 GBDEM을 이용해 Task G의 역학적, 수리-역학적, 열-역학적 벤치마크 모델링 연구를 수행하고 위 수치모델의 적합성과 적용성을 제시한 바 있다.

본 연구에서는 수치모델링을 위해 3차원 Voronoi 다면체 요소와 개별요소법 상용코드인 3DEC(Itasca CGI, 2022)을 이용하였다. Fig. 7은 해석모델의 형상과 구성 요소를 보여주는 것이다. 생성된 Voronoi 다면체의 개수는 총 10,240개이며, 각 다면체는 3DEC의 block 요소에 일대일 대응한다(Fig. 7(a)). 시험편의 부피를 다면체 개수로 나누어 환산한, block 요소의 평균 등가직경은 약 5.7 mm이다. Voronoi 다면체를 생성한 이후에는 saw-cut 균열을 모사하기 위하여 42° 경사의 평면을 교차시켜 시험편을 크게 두 그룹으로 분할하였다(Fig. 7(b)). 3DEC에서 열-역학적 해석을 수행하기 위해서는 열팽창에 의한 시료의 변형을 계산하여야 하므로, 변형이 허용되는 block (deformable block) 요소를 사용하였고, 이 경우 각 block은 DEC에서 다시 사면체의 zone 요소로 분할된다(Fig. 7(c)). Zone 요소의 개수는 295,225개이며, 이는 한 개 다면체 block이 평균 28개 이상의 zone 요소로 이루어짐을 의미한다. Zone 요소의 평균 등가직경은 약 1.9 mm이다. Fig. 7(d)는 block 요소 간의 경계면(interface)을 모두 도시한 것이다. 각 경계면에서의 힘 전달과 변위는 실질적으로 contact 요소를 통해 계산되며(Fig. 7(e)), Fig. 7(f)는 saw-cut 균열 상에 존재하는 contact 요소들만을 별도로 표시한 그림이다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F7.jpg
Fig. 7

Numerical specimen for simulating thermo-mechanical behavior of a single saw-cut fracture embedded in rock

본 연구에서는 block 요소 자체는 파괴가 발생하지 않는 탄성체로 간주하고, contact 요소에는 Coulomb slip 모델을 적용하였다. Contact는 전단파괴나 인장파괴 이전에는 강성(수직강성, 전단강성)에 의해 탄성적으로 거동하게 된다. 응력 조건이 전단파괴 상태에 도달하면 미끄러짐(slip)이 발생하고, 전단응력은 전단강도로 재설정된다. 인장파괴가 발생한 이후에는 인장강도가 0이 되고 점착력에 의해서만 전단응력에 저항하게 된다.

GBDEM을 통해 해석을 수행하고자 할 때는 일반적으로 block과 contact의 물성을 조정하는 과정(micro-parameter calibration)이 필요하다. 이는 전체 시료가 매우 작은 block과 수많은 contact의 결합으로 이루어져 있어서 block(또는 zone) 요소의 물성뿐만 아니라 contact 요소의 물성에도 크게 영향을 받기 때문이다. 보통 해석모델이 기지의 실험실 시험(압축강도시험 등) 결과를 재현할 때까지 block과 contact 물성을 반복적으로 개선, 적용함으로써 시행착오를 통해 최적의 조합을 결정하게 된다. 여기에서는 block 요소와 contact 요소의 물성을 결정하기 위하여 시행착오법과 함께 Park et al.(2020)이 제안한 등가연속체 접근법(equivalent continuum approach)을 이용하였다. Block 요소의 탄성정수(탄성계수, 포아송비)와 contact 요소의 강성(수직강성, 전단강성)은 등가연속체 접근법을 통해 결정하였으며, contact의 요소의 전단응력(마찰각, 점착력)과 인장강도는 일축압축강도시험(균열 미포함 시료)을 수치적으로 모사하여 시행착오에 의해 결정하였다. 최종적으로 결정된 입력물성 적용 시, 모델의 탄성계수, 포아송비, 일축압축강도는 각각 54.2 GPa, 0.27, 215.2 MPa로 나타났다. 화강암 시료의 일축압축강도(Table 1의 내부마찰각과 점착력에 따른 계산치)는 208.78 MPa로서 해석모델의 강도와 약 3%의 상대오차를 보였다. Saw-cut 균열에 해당하는 contact의 경우(Fig. 7(f)), 실험실에서 측정한 정적 마찰계수(0.58)로부터 마찰각 30.11°을 계산하여 입력하였고, 점착력과 인장강도는 모두 0으로 설정하였다. 전단강성과 수직강성의 경우, 활용가능한 데이터의 부재로 임의의 값을 가정하였다.

열물성(비열, 열전도도, 선형 열팽창계수)의 경우에는 실험실 측정결과를 block 요소의 입력물성으로 그대로 적용하였다. 단, 열-역학적 커플링에 결정적인 영향을 미치는 선형 열팽창계수에 대해서는 수치적으로 열팽창시험을 모사하여 contact 요소의 강성이 전체 열팽창변형률에 미치는 영향이 있는지 검토하였다. 열팽창시험 모델링 시 초기온도를 25°C로 가정하였고, 50°C, 75°C, 100°C, 125°C에서 열팽창에 의한 방향별 변형률을 측정함으로써 전체 시료의 선형 열팽창계수를 계산하였다. 모든 온도 조건에서 계산된 선형열팽창계수는 실험실 측정치와 동일하게 나타났으며, contact 요소의 강성이 열변형률에 영향을 미치지 않는 것을 확인하였다. Table 2는 해석모델의 입력물성과 이를 적용한 모델에 대한 수치적 일축압축강도시험 결과와 열팽창시험 결과를 요약한 표이다.

Table 2.

Calibrated input parameters, and results of numerical uniaxial compressive strength test and linear thermal expansion test

Parameters or results Value Unit
Model size (length × height × width ) 0.1 × 0.1 × 0.1 m3
Block (zone) density 2600 kg/m3
Block (zone) Young’s modulus 57.98 GPa
Block (zone) Poisson’s ratio 0.275 -
Contact friction coefficient rock 1.07 -
fracture 0.58
Contact cohesion rock 65 MPa
fracture 0
Contact tensile strength rock 15 MPa
fracture 0
Contact normal stiffness rock 1.93 × 1014 Pa/m
fracture 5.00 × 1011 (assumed)
Contact shear stiffness rock 7.57 × 1013 Pa/m
fracture 5.00 × 1011 (assumed)
Block (zone) thermal expansion coefficient 6.24 × 10-6 1/°C
Block (zone) thermal conductivity 1.96 W/m/K
Block (zone) specific heat 645 J/kg/K
Uniaxial compressive strength of calibrated specimen 215.2 MPa
Young’s modulus of calibrated specimen 54.2 GPa
Poisson’s ratio of calibrated specimen 0.27 -
Linear thermal expansion coefficient of calibrated specimen 6.24 × 10-6 1/°C

본 연구의 해석 과정은 실내실험과 동일하게 두 단계로 구분된다. 먼저 초기응력이 0, 초기온도가 25°C인 시료에 경계응력(σxx=0 MPa, σyy=3 MPa, σzz=9 MPa)을 가한 뒤 준정적 평형해석을 실시하여 시료와 균열면에 목표한 응력 상태를 구현하였다. Fig. 2의 실험 조건에서 σmin과 σmax는 수치모델의 σyy와 σzz에 대응한다. 엄밀한 의미에서 σxx가 σmin에 해당되지만, 균열의 주향이 X축과 평행하므로 이론적으로 σxx가 균열의 역학적 거동에 영향을 미치지 않게 되며, 2차원 평면응력 조건의 문제로 가정하여도 무방하다. 시료가 평형상태에 도달한 이후에는 Y축 방향의 변위를 구속한 상태에서 상하부의 가열에 의한 열전달(transient heat conduction) 해석을 실시하였다. 이때 열-역학적 커플링에 의해 암석 블록(zone 요소)의 열팽창 및 열응력의 증가가 계산되며, 열해석과 역학적 해석이 반복된다.

해석모델에서는 상부와 하부에 설치된 가열판과 히터를 별도의 해석요소로 고려하지 않았다. 따라서 히터의 온도를 경계조건으로 이용하는 대신, 시료의 하부면에서 실제로 측정된 온도변화 곡선(Fig. 4의 1번 센서 측정결과)을 수치모델의 상부, 하부 경계조건으로 적용하였다. 한편, 실험실에서 별도의 단열재는 설치되지 않았으므로, 시료 표면을 통한 열손실을 모델링에 반영할 필요가 있었다. 시험편의 전면과 후면에서는 공기의 대류에 의한 열손실이 발생할 수 있으며, 좌측면과 우측면의 경우 가압판으로의 전도와 가압판에서 공기로의 대류에 의한 열손실이 발생하였을 것이다. 본 연구에서는 열손실을 고려하기 위하여 전면과 후면, 좌측면과 우측면에 각각 30 W/m2K, 15 W/m2K의 열전달계수를 적용하였다. 참고로 기체의 자연대류(natural convection)에 의한 열손실계수는 2–25 W/m2K, 강제대류(forced convection)에 의한 열손실계수는 25–250 W/m2K 범위를 갖는 것으로 알려져 있다(Bergman et al., 2011). Fig. 8은 열해석 시 적용한 경계조건을 보여주는 모식도이다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F8.jpg
Fig. 8

Thermal boundary condition applied to numerical model for simulating thermoshearing test (modified from source: Zhuang et al., 2022)

3DEC의 열-역학 연계거동 해석에서는 기본적으로 일방향(one-way) 커플링을 가정한다. 다시 말해, 가열에 의한 열응력의 증가는 계산되지만, 역학적 변화가 열 유동 특성에 영향을 미치지는 않는다. 온도의 증가로 인해 zone 요소에 발달하는 열응력은 열팽창으로 인한 변형률에 기인하며, 이를 관계식으로 표현하면 식 (1), (2)와 같다.

(1)
ij=αtTδij

여기서 ϵij는 열에 의한 변형률, αt는 단위가 [1/°C]인 선형열팽창계수, δij는 Kronecker delta이다.

(2)
σij=-3Kij

여기서 σij는 열로 인한 응력이며, K는 체적탄성계수로서 응력의 단위를 갖는다.

4. 해석 결과

4.1 온도 분포

Fig. 9는 열-역학 해석이 완료된 시점(가열 6,000초 경과)에서 해석모델의 온도 분포를 보여준다. Fig. 10은 시간의 흐름에 따른 모니터링 지점에서의 온도 변화를 나타낸 것으로 수치해석 결과를 실험 결과(Fig. 4)와 비교하여 도시한 것이다. 그림에서 T1, T2, T3, T6은 온도센서 1, 2, 3, 6번이 부착된 지점과 동일한 위치에서의 결과이며, 온도센서 4, 5번의 경우 각각 6, 2번과 동일한 조건에 있어 별도로 고려하지 않았다. P0, P1, P2는 해석모델에서만 정의한 모니터링 지점으로 모두 균열면 상에 위치하며, P0은 균열 중심에, P1과 P2는 각각 균열의 경사방향과 주향방향을 따라 중심에서 0.25 m 거리에 놓인다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F9.jpg
Fig. 9

Contour plot of rock temperature after thermoshearing test

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F10.jpg
Fig. 10

Temperature distribution and evolution during the thermoshearing: (a) experimental and numerical results of surface temperatures at T1, T2, T3 and T4, and (b) numerical results of fracture temperatures at P0, P1 and P2

Fig. 9Fig. 10의 결과를 종합적으로 살펴보면, 수치모델의 온도 변화 곡선은 실험 결과와 거의 일치하는 것으로 계산되었다. 가열 초반, T2와 T3 지점에서 수치해석 결과가 실험 결과에 비해 다소 높게 나타나는 것을 제외하면, 수치모델이 실험의 열전달 조건과 온도 분포를 적절히 모사하고 있음을 확인할 수 있었다. 가열 후 약 2,000초까지는 열전도에 의해 모든 지점에서 온도가 서서히 증가하였고, 해석 종료 시(6,000 초)에 일정한 값으로 수렴하면서 열적 평형상태에 도달하는 경향을 보였다. P0, P1, T6의 결과를 함께 살펴보면, 동일한 높이에서 시료 표면에서 내부로 갈수록 높은 온도를 보이는 것으로 나타났는데, 이를 통해 시료 표면에서 발생하는 열손실이 해석 결과에 적절히 반영된 것을 확인할 수 있다. 또한, Fig. 10에서 시료의 좌우 측면부에 비해 전면부의 온도가 더 낮은 것을 확인할 수 있는데, 이는 더 큰 값의 열손실계수를 적용하였기 때문이다.

4.2 응력 분포

Fig. 11은 시간의 흐름에 따른 최대주응력(σyy)의 변화를 보여주며, 해석 결과와 실험 결과를 비교한 그림이다. σzz의 경우, 수치해석과 실험에서 모두 3.0 MPa로 일정하게 유지되었으므로 그림에 포함하지 않았다. 수치해석에서 계산된 주응력 변화는 실험 결과와 다소 상이한 경향을 나타내었다. 실험 결과에서는 가열이 종료될 때까지 σyy가 지속적으로 증가하여 6,000초 경과 시에는 약 12.3 MPa로 측정되었으나, 수치해석 모델은 초기 가파른 열응력 증가를 보이며 최대 11.2 MPa로 증가하였다가 약 800초를 전후하여 다소 감소하거나 수렴하는 경향을 보였다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F11.jpg
Fig. 11

Comparison of experimentally measured and numerically calculated maximum principal stress change

본 연구에서는 수치해석에서 열응력이 감소하는 원인이 균열의 전단파괴와 수직방향으로의 열팽창에 있는 것으로 판단하였다. 수치해석상 균열 전체에서 전단파괴가 발생하는 시점은 약 800초로 이 시점 이후로 시료의 수직방향 변위가 증가하게 된다. 시료의 가열 시 변위가 구속되어 있는 경우, 열응력은 지속적으로 증가하게 되며, 열팽창에 상응하는 변위가 허용되는 경우 응력 증가는 발생하지 않았다. 균열의 마찰계수를 0.58 대신 0.71로 상향조정한 경우, 전단파괴 개시 시점이 지연되고 더 작은 변위가 발생하였으며 실험 결과와 같이 σyy가 지속적으로 증가하는 경향을 보임을 확인하였다. 균열의 변위에 대한 해석 결과는 제4.3절에서 자세히 논의하였다.

한편, 수치해석 모델에서는 가열 초기 열응력의 증가율이 상대적으로 크게 나타났는데, 이는 경계조건의 영향 때문인 것으로 판단된다. 수치해석에서는 수평방향의 변위를 완전히 구속한 롤러(roller) 경계조건이 사용되었으나, 실험실에서 시편의 경계는 가압판과 맞닿아 있으며 가압판의 강성과 열팽창에 따라 실제로는 변위가 일부 허용되는 조건이다. Park et al.(2021b)의 연구에서는 변위 구속 조건에 따라 열응력의 발달이 상이하게 나타나게 됨을 수치해석적으로 제시한 바 있다. 향후 후속 연구를 통해 변위구속 조건이 아닌 일정수직강성(constant normal stiffness) 조건 등을 적용함으로써 보다 현실적인 실험실 조건을 수치해석 모델에 반영하여야 할 것으로 판단된다.

Fig. 12는 가열시간이 0초(가열 시작 단계), 500, 6000초일 때의 해석 결과로서 균열상 contact 요소들의 응력 상태(연두색 원, 노란색 원, 분홍색 원)와 함께 파괴포락선과 주응력 Mohr 원을 도시한 것이다. 그림에서 0초의 Mohr 원 위에 있는 빨간색 원은 균열의 이론적 응력상태를 표시한 것이다. 응력 해석(Jaeger et al., 2007)에 따르면, 주어진 응력 조건에서 42°의 경사를 갖는 균열에 작용하는 수직응력(σn)과 전단응력(τ)은 각각 5.69 MPa과 2.98 MPa이다. 마찰계수를 고려하면 상기한 수직응력하에서 균열에 작용하는 전단강도는 3.30 MPa로서 가열 전 균열면의 응력상태는 포락선에 근접한 임계응력(critical stress) 상태에 있음을 확인할 수 있다. 가열 시작 단계(0초)에서 균열의 contact 요소들의 응력 상태는 이론치와 거의 근접하게 해석되었으며, 몇몇 contact 요소의 응력상태는 이미 전단파괴 포락선 위에 있음을 확인할 수 있었다. 시간이 경과함에 따라 수직응력과 전단응력이 모두 증가하였고(500초), 가열 종료 시점에는 균열 상의 모든 contact 요소가 전단파괴 상태에 놓였다(6000초).

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F12.jpg
Fig. 12

Changes in normal and shear stresses of fracture contacts with increasing temperature and Mohr’s circle presentation of boundary stresses (numerical results of 0, 500 and 6000 seconds of heating)

Fig. 13은 가열 시간(300, 400, 500, 600, 700, 800초)에 따라 전단파괴가 발생한 contact 요소(빨간색 +)의 위치를 보여준다. 균열을 따라서 약 300초 이후 서서히 국부적 전단파괴가 발생하였으며, 시간의 경과에 따라 중심에서부터 파괴영역이 점점 확대되어 약 800초 이후에는 균열면상의 거의 모든 contact 요소가 전단파괴 상태에 놓임을 확인할 수 있었다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F13.jpg
Fig. 13

Failure states of contacts existing on fracture plane; a cross denotes shear failure state

4.3 변위 분포

Fig. 14는 가열 500, 1000, 1500초일 때의 시료 변위(displacement magnitude)를 보여주는 것으로 변위 벡터(scaled displacement vector)와 함께 표시하였다. 수치해석에서 좌측면과 우측면(Y축 방향)에 변위 구속조건을 부여하였으므로, 전면과 후면(X축 방향) 방향과 상부면, 하부면(Z축 방향) 방향으로만, 변위가 허용된다. 가열에 의해 수직방향 열응력이 증가하게 되면 서보컨트롤에 의해 상부경계면은 위로, 하부 경계면은 아래로 이동시킴으로써 일정응력을 유지하게 된다. 따라서 균열의 전단방향은 그림의 화살표 방향으로 발생하는 것으로 나타났다. 시료의 온도가 비교적 낮은 수준인 500초일 때에는, 가열이 이루어지는 상부면과 하부면에서 열팽창으로 인해 X축 방향(자유면 방향)으로도 우세한 변위가 관찰된다. 그러나 시간이 경과함에 따라 전체 시료의 변위는 균열 전단변위의 방향에 지배적인 영향을 받는 것으로 나타났다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F14.jpg
Fig. 14

Contour plots of displacement magnitude at 500, 1000, 6000 seconds of heating (unit: m)

Fig. 15는 균열의 전단변위를 보여주는 것으로 해석 결과와 실험 결과를 비교한 그림이다. 그림에서 붉은 실선으로 표시된 데이터는 변위 변환기를 통해 실험실에서 측정된 전단변위이며, 원으로 표시된 데이터는 DIC 기법을 통해 균열면을 따라 5개 지점에서 분석된 전단변위이다. 실험 결과(붉은 실선)를 살펴보면, 약 2,000초까지 전단변위가 음의 값을 보이다가 다시 증가하여 가열 종료 시 약 80 µm의 전단변위를 보였다. 이러한 결과는 균열이 Fig. 14의 화살표의 역방향으로 미끄러지다가 2,000초 이후 다시 화살표 방향으로 전단변위를 발생시킴을 의미한다. 사실상 실험에서 측정된 데이터만으로 그 원인과 메커니즘을 추론하기는 쉽지 않으며, 이에 대한 수치모델링 역시 본 논문의 연구범위를 벗어난다. 단, DIC 분석을 통해 얻어진 균열의 전단변위가 2,000초까지 경미한 수준임을 고려할 때, 하중 재하 초기 단계에서 균열면의 접촉상태나 매칭정도에 따라 균열면의 역방향 미끄러짐이 측정결과에 반영된 것으로 추정해 볼 수 있다. 따라서 균열의 전단변위는 2,000초 전후부터 증가하여 해석 종료 시까지 약 100 µm 발생한 것으로 해석할 수 있다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F15.jpg
Fig. 15

Comparison of experimentally measured and numerically calculated fracture shear displacement

수치해석 결과에서는 실험 결과보다 더 이른 시점(낮은 온도)에서 전단변위가 개시되는 것으로 나타났으며, 전반적으로 더 낮은 수치의 전단변위를 보였다. 약 800초 경과 시 전단변위가 크게 증가하기 시작하여 3,000초 이후 일정한 값으로 수렴하는 경향을 보였으며, 해석종료 시점(6,000초 경과 시)에 총 42 µm의 전단변위를 보였다. Fig. 12에서 살펴본 바와 같이 본 연구의 해석모델에서는 가열이 시작되기 전 이미 전단파괴 상태에 있는 contact 요소들이 관찰되었는데, 균열의 마찰계수를 0.58 대신 0.71로 상향조정하는 경우, 초기상태에서 contact 요소들이 모두 역학적 안정 상태에 있고, 가열에 따른 전단파괴 시점이 지연되는 결과를 보였다. Sun et al.(2022)의 연구에서는 연마하지 않은 saw-cut 균열이 거칠기(JRC > 2)를 가지고 있음을 보고한 바 있으며, 향후 균열의 마찰계수의 적용과 관련된 추가 연구가 필요할 것으로 판단된다.

한편, 앞서 Fig. 13에서 800초 이전에도 균열 내부에 비교적 큰 규모의 국부적 전단파괴가 발생하는 것을 관찰하였는데, Fig. 15에서는 800초 이전의 전단변위는 미미한 수준이다. 이를 통해 해석모델에서 균열의 모든 contact 요소가 전단파괴 상태에 이른 이후에 거시적 규모의 균열 미끄러짐이 발생함을 확인할 수 있었다.

Fig. 16은 균열의 수직변위를 보여주는 것으로 해석 결과와 실험 결과를 비교한 그림이다. 여기서 음의 수직변위는 균열의 닫힘(fracture closure)을 의미한다. 수치해석 결과가 실험 결과에 비하여 전반적으로 작은 수직변위량(magnitude)을 보였지만, 두 데이터는 매우 유사한 경향을 나타내고 있음을 확인할 수 있다. 가열 초기에는 균열의 압축 방향으로 수직변위량이 증가하다가 급격한 전단변위가 발생하기 시작하면서(수치해석: 약 800 초, 실험: 약 2000 초) 수직변위량이 감소하였다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F16.jpg
Fig. 16

Comparison of experimentally measured and numerically calculated fracture normal displacement

Fig. 17Fig. 18는 각각 500, 1000, 6000초 경과 시 균열의 전단변위와 수직변위를 나타낸 그림이다. 그림에서 A와 B는 실험에서 변위 변환기가 설치된 지점을 표시한 것으로, 상기한 Fig. 15Fig. 16의 수치해석 결과는 A 지점과 B 지점에서 구한 변위의 평균값이다. Theromoshearing 실험의 경우, 모니터링 지점과 열전달 정도에 따라 균열의 응력과 변위 분포가 상이하므로 A와 B 지점에서 측정한 변위가 균열면의 변위를 대표한다고 보기 어렵다. 전단 변위의 경우 측정 지점보다 시료 내부에서 더 큰 변위를 보이는 것으로 나타났고, 온도가 높은 상부 경계면과 하부 경계면에서 가장 큰 전단변위를 보였다. 수직 변위의 경우에도 측정지점에서 균열이 압축 상태(fracture closure)에 있는 것과는 달리, 균열면 중심에서는 균열의 열림(fracture opening) 현상에 의한 인장파괴와 전단파괴가 동시에 관찰되었다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F17.jpg
Fig. 17

Contour plots of fracture shear displacements at 500, 1000 and 6000 seconds of heating (unit: m)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-06/N0120320618/images/ksrm_32_06_18_F18.jpg
Fig. 18

Contour plots of fracture normal displacements at 500, 1000, and 6000 seconds of heating (unit: m)

5. 요약 및 결론

본 논문에서는 국제공동연구인 DECOVALEX-2023 프로젝트 Task G의 일환으로 수행된 암석 균열의 열-역학적 연계거동 수치모델링 연구 결과를 소개하였다. 해석대상은 한국건설기술연구원에서 수행된 ‘KICT-Thermoshearing 실험’ 결과로서, 3차원 Voronoi 입자기반 개별요소모델을 통해 열적, 역학적 하중을 받는 saw-cut 암석 균열의 전단거동을 살펴보았다. 주요 해석내용은 가열로 인한 암석 표면의 온도 분포, 열응력의 증가에 따른 주응력 변화, 균열의 전단변위와 수직변위로서 실험결과와의 비교를 통해 해석모델의 적용성을 검토하였다. 수치해석의 장점을 활용하여, 실험실에서 직접적인 측정이 어려운 시료 내부의 온도 분포, 균열 응력과 변위의 시공간적 발달 양상, 점진적인 전단파괴 거동 등을 추가적으로 논의하였다. 온도 분포의 경우, 수치모델이 실험 결과를 상당히 유사한 수준으로 재현하였으나 주응력 변화와 균열 변위는 다소 상이한 결과를 나타내었다. 가열이 진행됨에 따라서 수치모델이 암석 균열에 비해 수평주응력의 변화율이 더 크게 나타났으며, 이에 따라 더 이른 시점에서 전단파괴에 도달하였다. 그러나 수치해석모델은 실험에서 측정된 열응력 증가, 전단파괴로 인한 변위 발생 경향 등을 합리적으로 재현하였으며, 전반적으로 해석결과와 실험 결과는 유사한 범위(order of magnitude) 내에 분포하였다.

본 논문에서 기술한 내용은 상기 수치해석모델의 KICT-thermoshearing 실험에 대한 적용성을 검토하기 위한 예비해석 결과로서, 향후 실험 결과를 보다 현실적으로 재현하기 위해 균열의 마찰계수와 응력상태, 암석 시편과 가압판 사이의 강성, 가압판의 열팽창 등이 해석 결과에 미치는 영향을 면밀히 검토하고, 해석모델을 개선할 예정이다.

Acknowledgements

The authors appreciate and thank the DECOVALEX-2023 Funding Organizations, Andra, BASE, BGE, BGR, CAS, CNSC, COVRA, US DOE, ENRESA, ENSI, JAEA, KAERI, NWMO, RWM, SÚRAO, SSM and Taipower for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the Funding Organizations. This research was supported by the Basic Research Project of the Korea Institute of Geoscience and Mineral Resources (GP2020-010) funded by the Ministry of Science and ICT, Korea.

References

1
Bergman, T.L., Lavine, A.S., Incropera, F.P., and DeWitt, D.P., 2011, Fundamentals of Heat and Mass Transfer, 7th edition. Wiley, New York.
2
Fu, T. F., Xu, T., Heap, M. J., Meredith, P. G., and Mitchell, T. M., 2020, Mesoscopic time-dependent behavior of rocks based on three-dimensional discrete element grain-based model. Computers and Geotechnics, 121, 103472. 10.1016/j.compgeo.2020.103472
3
Ghazvinian, E., Diederichs, M. S., and Quey, R., 2014, 3D random Voronoi grain-based models for simulation of brittle rock damage and fabric-guided micro-fracturing. Journal of Rock Mechanics and Geotechnical Engineering, 6(6), 506-521. 10.1016/j.jrmge.2014.09.001
4
IAEA, 2011, Geological disposal facilities for radioactive waste. Specific safety guide No. SSG-14, IAEA, Vienna, Austria.
5
Institute for Korea Spent Nuclear Fuel, 2022, https://iksnf.or.kr/ accessed on 1 December 2022.
6
Itasca Consulting Group Inc., 2022. 3DEC (3 Dimensional Distinct Element Code) online manual. https://www.itascacg.com/software/3DEC accessed on 1 December 2022.
7
Jaeger, J.C., Cook, N.G.W., and Zimmerman, R.W., 2007, Fundamentals in Rock Mechanics. fourth ed. Oxford: Blackwell publishing.
8
Kim, T., Lee, C., Kim, J. W., Kang, S., Kwon, S., Kim, K. I., Park, J.W., Park, C.H., and Kim, J. S., 2021, Introduction to Tasks in the International Cooperation Project, DECOVALEX-2023 for the Simulation of Coupled Thermohydro-mechanical- chemical Behavior in a Deep Geological Disposal of High-level Radioactive Waste. Tunnel and Underground Space, 31(3), 167-183.
9
Lan, H.X., Martin, C.D., and Hu, B., 2010, Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. Journal of Geophysical Research: Solid Earth 115, B1. 10.1029/2009JB006496
10
Lee, J., Cho, D., Choi, H., and Choi, J., 2007, Concept of a Korean reference disposal system for spent fuels. Journal of Nuclear Science and Technology, 44(12), 1565-1573. 10.1080/18811248.2007.9711407
11
Malmgren, L., Saiang, D., Töyrä, J., and Bodare, A., 2007, The excavation disturbed zone (EDZ) at Kiirunavaara mine, Sweden - by seismic measurements. Journal of Applied Geophysics, 61, 1-15. 10.1016/j.jappgeo.2006.04.004
12
McDermott, C. I., Fraser-Harris, A., Sauter, M., Couples, G. D., Edlmann, K., Kolditz, O., Lightbody, A., Somerville, J., and Wang, W., 2018, New experimental equipment recreating geo-reservoir conditions in large, fractured, porous samples to investigate coupled thermal, hydraulic and polyaxial stress processes. Scientific Reports, 8(1), 1-12. 10.1038/s41598-018-32753-z30266937PMC6162306
13
Min, K., B., Lee, J., and Stephansson, O., 2013, Implications of thermally-induced fracture slip and permeability change on the long-term performance of a deep geological repository, International Journal of Rock Mechanics and Mining Sciences, 61, 275-288. 10.1016/j.ijrmms.2013.03.009
14
Park, J.W., Park, C., Song, J. W., Park, E. S., and Song, J. J., 2017, Polygonal grain-based distinct element modeling for mechanical behavior of brittle rock. International Journal for Numerical and Analytical Methods in Geomechanics, 41(6), 880-898. 10.1002/nag.2634
15
Park, J.W., Park, C.H., and Lee, C., 2021a, Hydro-Mechanical Modeling of Fracture Opening and Slip using Grain-Based Distinct Element Model: DECOVALEX-2023 Task G (Benchmark Simulation). Tunnel and Underground Space, 31(4), 270-288.
16
Park, J.W., Park, C.H., and Lee, C., 2021b, Voronoi Grain-Based Distinct Element Modeling of Thermally Induced Fracture Slip: DECOVALEX-2023 Task G (Benchmark Simulation). Tunnel and Underground Space, 31(6), 593-609.
17
Park, J.W., Park, C.H., Yoon, J.S., and Lee, C., 2020, Grain-Based Distinct Element Modelling of the Mechanical Behavior of a Single Fracture Embedded in Rock: DECOVALEX-2023 Task G (Benchmark Simulation). Tunnel and Underground Space, 30(6), 573-590.
18
Rutqvist, J., 2020. Thermal management associated with geologic disposal of large spent nuclear fuel canisters in tunnels with thermally engineered backfill. Tunnelling and Underground Space Technology, 102, 103454. 10.1016/j.tust.2020.103454
19
Sun, C., Zhuang, L., Jung, S., Lee, J., and Yoon, J. S., 2021, Thermally induced slip of a single sawcut granite fracture under biaxial loading. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 7(4), 1-13. 10.1007/s40948-021-00293-y
20
Sun, C., Zhuang, L., Yoon, J.S., and Min, K.B., 2022, Thermally induced shear reactivation of critically-stressed smooth and rough granite fractures, Proceedings of EUROCK 2022, https://www.ril.fi/media/2022-eurock/eurock-2022-programme.pdf accessed on 1 December 2022.
21
Swedish Nuclear Fuel and Waste Management Company, 2010, Choice of Method - Evaluation of Strategies and Systems for Disposal of Spent Nuclear Fuel, SKB report, SKB P-10-47.
22
Wang, Z., Wang, T., Wu, S., and Hao, Y., 2021, Investigation of microcracking behaviors in brittle rock using polygonal grain‐based distinct method. International Journal for Numerical and Analytical Methods in Geomechanics, 45(13), 1871-1899. 10.1002/nag.3246
23
Zhuang, L., Sun, C., and Yoon. J., 2022, Laboratory investigation of thermoshearing in critically stressed sawcut and unmated rough granite fractures, DECOVALEX-2023 5th Workshop, Virtual conference.
24
Zoback, M.D. and Gorelick, S.M., 2012, Earthquake triggering and large-scale geologic storage of carbon dioxide, Proceedings of the National Academy of Sciences of the United States of America, 109(26), 10164-10168. 10.1073/pnas.120247310922711814PMC3387039
페이지 상단으로 이동하기