1. 서 론
2. 연구 배경
2.1 DECOVALEX-2023 TASK G
2.2 한국건설기술연구원 Thermoshearing 실험 개요
3. 수치모델 및 해석개요
4. 해석결과
4.1 온도 분포
4.2 열응력의 증가
4.3 균열 변위
5. 요약 및 결론
1. 서 론
단층재활성(fault reactivation)은 지각활동(지진, 지체응력 등)이나 인간활동(암반, 지하공간 활용 등)에 의해 기존에 비활성화되어 있던 단층이 불안정한 응력 상태에 놓이면서 다시 활성화되는 현상을 의미한다. 이와 관련된 연구는 주로 이산화탄소 지중저장, 인공저류층 지열발전, 원유 회수증진 등과 같이 고압 유체의 주입을 고려하는 기술 분야에서 많이 이루어져 있다(Gudmundsson, 2004, Vidal-Gilbert et al., 2009, Morris et al., 2011, Cappa and Rutqvist, 2011, Rinaldi et al., 2014). 유체의 주입은 인접한 단층에 작용하는 유효 수직응력을 감소시켜 이에 따른 전단저항력의 저하와 전단파괴(shear failure, slip)를 야기할 수 있다. 또한, 이로 인해 인공지진을 유발할 수 있으므로 주변 지역 및 시설물의 안정성과 프로젝트의 성패를 가르는 중요한 요소로 인식되어 왔다. 최근에는 방사성폐기물 지층처분과 관련하여서도, 처분장 주변 단층(또는 균열)의 열에 의한 재활성 문제에 관한 연구가 증가하고 있는 추세이다(Min et al., 2013). Urpi et al.(2019)은 저투수성 지층에서의 폐기물 처분을 모사한 수치해석 연구를 통해, 열응력(thermal stress)과 열압력(thermal pressurization or thermally driven fluid pressure)의 증가가 최대 500미터 거리에 있는 단층에 재활성을 유도할 수 있음을 보고하였다. 처분장에서의 단층재활성 문제는 조사 및 시공단계에서 저투수성 또는 불투수성으로 분류되었던 단층이 운영단계에서 영구적인 핵종 이동 경로(flow path)로 작용할 수 있다는 점에서 매우 중요한 고려사항이다(Rutqvist, 2020, Rutqvist et al., 2020).
단층재활성 문제를 평가하기 위한 기법은 해석적 방법(analytic method)과 수치적 방법(numerical method)으로 나뉜다(Bohloli et al., 2015). 해석적 방법은 단층의 방향과 물성, 압력, 현지응력 등을 바탕으로 Mohr-Coulomb 파괴기준에 의한 전단파괴 가능성을 개략적으로 평가하는 방법(Wiprut and Zoback, 2002, Vidal-Gilbert et al., 2010)이다. 그러나 현장 응력조건에 대한 여러 가정과 단순화를 전제하며, 수리-역학적 일방향 커플링 문제에만 적용이 국한된다. 반면 수치적 접근법을 통하면, 균열과 암반의 상호작용과 열-수리-역학 커플링에 의한 복합거동을 고려할 수 있고, 초기응력과 유도응력, 시간 경과에 따른 물성 변화 및 점진적 파괴과정 등을 예측할 수 있다는 장점이 있다. 따라서 상기한 여러 암반공학 분야에서 불연속면의 특성을 합리적으로 반영할 수 있는 열-수리-역학 복합거동 수치해석 기법의 개발은 단층 재활성 위험을 평가하기 위해 필수적인 연구과제라고 할 수 있다. 방사성폐기물 지층처분과 관련된 대표적인 국제공동연구인 DECOVALEX (DEvelopment of COupled models and their VALidation against EXperiments) 프로젝트에서도 DECOVALEX-2019 Task B와 DECOVALEX-2023 Task G를 통해 단층재활성 문제에 대한 연구가 진행 중에 있다(Birkholzer et al., 2019, Kim et al., 2021).
Cundall and Strack(1979)이 제안한 개별요소법(distinct element method, DEM)은 단층 등 암반 불연속면의 거동을 가장 직접적인 방식으로 모사하는 대표적인 해석기법이다. 개별요소법은 해석대상을 개별적, 독립적으로 거동하는 요소들의 집합으로 표현하고 요소들 간 상호작용과 경계면에서의 거동을 계산함으로써 불연속 암반체의 거동을 평가한다. 암반공학 분야에서 대부분의 개별요소법 적용사례에서는 단층과 절리와 같은 공학적 스케일의 불연속면이 경계면으로 고려되지만, 최근의 연구들에서는 미시적 스케일에서 광물이나 입자의 경계면에서 발생하는 거동을 재현하고자 하는 연구들도 많이 수행되고 있다(Lan et al., 2010, Ghazvinian et al., 2014, Park et al., 2017, Fu et al., 2020). 상기한 연구들에서는 다면체나 다각형 입자를 사용하는 개별요소모델이 입자간 맞물림(interlocking), 회전 모멘트에 대한 저항, 암석 재료의 불균질성 및 이방성, 미세균열의 개시와 성장에 따른 취성파괴 메커니즘 등을 효과적으로 재현할 수 있음을 보고하였다(Park et al., 2022).
본 연구는 DECOVALEX-2023 국제공동연구 Task G의 일환으로 수행된 연구로서, 입자기반 개별요소모델(grain-based distinct element model, GBDEM)을 이용하여, 암석 균열의 열에 의한 미끄러짐 거동(fracture thermoshearing or thermally induced fracture slip)을 모사하였다. Task G는 결정질 암석 내 균열의 열-수리-역학적 복합거동을 해석할 수 있는 수치해석 코드를 개발하는 데에 최종 목표가 있으며, 여기에서는 실내실험 모델링 단계에서 수행한 연구 결과를 제시하였다. 해석대상은 한국건설기술연구원의 ‘Thermoshearing’ 실험(Sun et al., 2022)으로 거친 암석 균열에 대한 열-역학 재하실험이다.
2. 연구 배경
2.1 DECOVALEX-2023 TASK G
DECOVALEX 프로젝트는 1992년부터 현재까지 진행되고 있는 국제공동연구로서 방사성폐기물의 지층처분과 관련하여 처분장의 열-수리-역학-화학적(Thermal-Hydro-Mechanical-Chemical, THMC) 복합거동을 평가할 수 있는 수치해석 코드를 개발, 검증하는 데에 목표가 있다(Birkholzer et al., 2019). DECOVALEX 프로젝트에서는 4년을 주기로 새로운 주제(Task)들에 대한 연구가 이루어지며, 현재 DECOVLAEX-2023(2020-2023)에 14개국에서 50여개의 관련기관이 참여하여 7개 Task에 대한 공동연구를 수행하고 있다(Kim et al., 2021). 국내에서는 한국원자력연구원이 정식회원기관과 연구기관으로, 한국지질자원연구원이 한국원자력연구원의 협력 연구기관으로 참여하고 있다.
한국지질자원연구원이 참여하고 있는 Task G는 독일의 Freiberg University of Mining and Technology (TUBAF)와 Helmholtz Centre for Environmental Research (UFZ), 스코틀랜드의 Edinburgh 대학교, 독일의 DyanFrax 사에서 공동 제안한 연구주제로, 총 9개의 연구팀이 참여하고 있다. Task G의 주요 연구목표는 결정질 암반 내 균열의 열-수리-역학적 복합거동을 해석할 수 있는 수치해석 코드를 개발, 검증하는 데에 있으며, 정해(analytical solution)에 기반한 벤치마크 모델링, 참여팀간의 연구결과 상호비교, 실험 결과 모델링을 통한 검증(validation) 등을 주요 연구내용으로 한다. 실험 결과 모델링은 크게 균열에서의 수리-역학적 커플링과 열-역학적 커플링 문제로 구분되며, 각각 Edinburgh 대학교에서 수행된 GREAT Cell 실험(McDermott et al., 2018)과 한국건설기술연구원에서 수행된 Thermoshearing 실험(Sun et al., 2021, Sun et al., 2022)을 대상으로 한다. 본 연구의 선행연구인 Park et al.(2020), Park et al.(2021a), Park et al.(2021b), Park et al.(2022)에 Task G의 역학적, 수리-역학적, 열-역학적 벤치마크 모델링, saw-cut 균열에 대한 실험 모델링 결과가 소개된 바 있다. Table 1은 Task G 참여팀과 해석 코드, 참여중인 연구테마를 보여준다.
Table 1.
Task G modeling teams and approaches
SSM: Swedish Radiation Safety Authority, UFZ: Helmholtz Centre for Environmental Research, LBNL: Lawrence Berkeley National Laboratory, DOE: Department of Energy, CNSC: Canadian Nuclear Safety Commision, SNL: Sandia National Laboratories, CAS: Chinese Academy of Science, KIGAM: Korea Institute of Geoscience and Mineral Resources, KAERI: Korea Atomic Energy Research Institute, UoE: University of Edinburgh, NWA: Nuclear Wasted Service, DEM: distinct element method, FEM: finite element method, CEM: cohesive element method and FVM: finite volume method
2.2 한국건설기술연구원 Thermoshearing 실험 개요
Thermoshearing 실험은 균열(saw-cut fracture, laser-marketing fracture, tensile fracture)을 포함한 암석 시료에 하중을 가하여 균열면이 임계응력 상태에 놓이도록 유도하고 가열하여 열에 의한 균열의 미끄러짐 거동을 살펴본 실험이다(Sun et al., 2021, Sun et al., 2022). Park et al.(2022)은 입자기반 개별요소모델을 통해 saw-cut 균열에 대한 실험 결과를 수치적으로 모사한 바 있으며, 본 논문에서는 그 후속 연구로서 거친 암석 균열(tensile-splitting fracture)에 대한 모델링 결과를 제시하였다. Saw-cut 균열과 거친 균열에 대한 실험은 실험장치와 방법이 거의 동일하여, 본 절의 내용은 Park et al.(2022)의 실험 개요를 상당 부분 인용하여 작성하였다.
실험에 사용된 시료는 국내 포천화강암으로 길이 100 mm의 정육면체 형상이며, 시료 내부에 거친 인장균열(tensile-splitting fracture)을 포함하고 있다. 균열의 경사는 수평면으로부터 약 42°로서, 상부와 하부 시료가 균열의 전단방향으로 약 3 mm의 변위량(offset)을 갖는 상태를 초기 조건으로 가정한다. 즉 실험 시 상부와 하부의 균열면이 맞물리지 않은 상태(unmated fracture)에서 가압과 가열이 수행되었다. 3차원 레이저 스캐닝 기법을 통해 상부와 하부 균열면의 프로파일이 측정되었으며, 이를 바탕으로 접촉면적을 분석한 결과, 실험 시 총면적의 약 10%만이 접촉상태에 놓일 것으로 추정되었다(Zhuang et al., 2022). Fig. 1은 시험편의 형상을 보여주며, Table 2는 시험편의 물성을 정리한 것이다.
Table 2.
Physical and mechanical properties of rock and fracture (Sun et al., 2022)
Thermoshearing 실험은 크게 역학적 하중과 열하중을 재하하는 단계로 구분할 수 있다. Fig. 2는 역학적 하중 재하 단계의 경계응력 조건을 보여주는 것으로 X축 방향의 전면과 후면을 자유면으로 하는 2축 응력조건이다. 이 단계에서는 수평(Y축) 방향과 수직(Z축)방향의 응력이 각각 24.6 MPa()과 3.0 MPa()에 이를 때까지 점진적으로 하중을 재하한다. 응력해석(Jaeger et al., 2007)에 따르면, 위 경계응력 조건에서 경사각이 42°인 균열면에 작용하는 수직응력()과 전단응력()은 각각 12.67 MPa과 10.74 MPa이다. 수직응력과 마찰계수(0.86)로부터 균열의 전단강도를 구하면 10.90 MPa로서, 위 조건에서 균열은 이론적으로 전단파괴에 대한 임계응력(critical stress) 상태에 놓인다.
Fig. 2.
Boundary conditions during mechanical loading test and thermo-mechanical loading test (modified from source: Sun et al., 2022)
열하중 재하단계에서는 시료의 상부와 하부 로드셀 외곽에서 히터를 통해 시료를 가열한다. 가열 중에는 를 일정하게 유지하되, 수평방향으로 변위(열팽창)를 구속하여 열응력에 의한 의 증가를 유도한다. Fig. 3은 실험에서 사용된 히터와 모니터링 장치의 위치를 보여준다. 그림에서 시료의 상하부 경계면에는 가압판(loading plate) 외곽에 각각 4개의 히터가 설치되었으며, 시험편 표면에 부착된 6개(하부면 1개, 전면 5개)의 센서(thermometer)를 통해 시간의 경과에 따른 표면 온도 변화가 측정되었다. 균열 변위를 측정하기 위해 2개의 변위 변환기(clip-on displacement transducer)와 디지털 이미지 상관관계(digital image correlation, DIC) 분석 기법이 활용되었다. 시험편 전면에 균열을 교차하는 2개의 변위 변환기를 부착하였으며, 후면에서는 실시간으로 촬영된 이미지를 바탕으로 균열 주변 5개 지점에서의 변위를 DIC 기법을 통해 분석하였다. 한편, 좌측면, 우측면, 전면에 10개의 센서를 부착하여 시험 중 발생하는 미소파괴음(acoustic emisson, AE)을 모니터링하였다. 실험장비 및 모니터링과 관련된 자세한 내용은 Sun et al.(2022)의 논문에 기술되어 있다.
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)
본 수치모델링에서는 시험편의 온도 분포와 균열의 변위를 모사하는 데에 목표를 두었으며, AE 측정 자료는 비교 분석에 고려하지 않았다. 인장균열 시료의 경우, 시험 중 표면에서 온도가 측정되지 않았으므로, 수치모델링에서는 saw-cut 균열 시료의 온도 측정 자료(Sun et al., 2022)를 비교 분석에 이용하였다. Fig. 4는 가열 중 모니터링 지점에서 측정된 온도 변화 곡선을 보여주는 것으로 시험편에 표기된 숫자는 각 지점의 식별번호이다. 히터는 상부와 하부의 가압판 외곽에 설치되었으며, 초기 1,000초 동안 히터 온도를 서서히 증가시켜 150°C에 도달하면 실험이 종료될 때까지 온도를 유지하였다. Fig. 4에서 1로 표시된 곡선은 시험편의 하부면과 가압판 사이에 부착된 1번 온도센서에서 측정된 온도로, 약 200초부터 서서히 증가하여 실험 종료 시 약 106°C에 수렴하였다.
Fig. 4.
Temperature distributions and revolutions at six measuring points during thermoshearing test (modified from source: Sun et al., 2022)
Fig. 5는 시간에 따른 최대주응력()의 증가량과 변위 변환기를 통해 측정된 균열의 전단변위, 수직변위 곡선이다. Fig. 6은 균열 전단변위 곡선으로 DIC 기법을 통해 5개 지점에서 분석된 결과와 변위 변환기를 통해 측정된 결과를 함께 도시한 것이다.
Fig. 5.
Maximum principal stress increment and fracture shear and normal displacements (Zhuang et al., 2022)
Fig. 6.
Fracture shear displacements during the thermoshearing measured by the clip-on displacement transducer and the DIC technique (Zhuang et al., 2022)
3. 수치모델 및 해석개요
본 연구에서는 입자기반 개별요소모델(GBDEM)을 이용하여 거친 균열에 대한 thermoshearing 실험 결과를 모사하였다. GBDEM은 해석대상을 다각형(2차원)이나 다면체(3차원) 입자의 결합체로 모델링하고, 입자와 입자간 접촉에서 발생하는 거동을 개별요소법에 의해 계산하는 방법이다. 여기에서는 사면체 형상의 입자모델과 Itasca 사의 개별요소 코드인 3DEC(Itasca CGI, 2022)을 이용하였다.
균열의 거칠기를 정의하는 스케일은 공학적 스케일에서 비해 매우 작기 때문에, 수치모델링 시 직선(2차원)이나 평면(3차원) 요소를 이용해 균열의 거동을 단순화하는 것이 일반적이다. 이 경우, 거칠기에 의해 발현되는 역학적 거동 특성은 전단강성, 정점마찰각(peak friction angle), 팽창각(dilation angle)과 같이 수치화된 물성이나 구성모델을 통해 해석에 반영된다. 그러나 본 연구에서는 균열의 전단거동에 거칠기가 미치는 효과를 직접적으로 모사하는 데에 목표를 두고 거친 균열면의 형상을 수치모델에 반영하고자 하였다.
Fig. 7은 실험실 시험편 균열의 거칠기 프로파일 측정 자료를 도시한 것이다. 거칠기 프로파일은 레이저 스캐닝 기법에 의해 상하부 균열면에서 각각 측정되었으며, 측정방향에 의해 상부 표면은 하부 표면과 서로 반전된 이미지를 갖게 된다. 즉, 상부 균열면에서 함몰부에 해당하는 지점은 하부 균열면의 돌출부에 상응하게 되므로, 두 그림에서 A와 B는 동일한 지점을 가리킨다. 인장균열의 생성 과정에서 발생하는 표면의 부분적 탈락과 레이저 스캐닝 과정의 측정오차 등으로 인해 두 표면의 좌표가 정확하게 일치하지 않았으나, 두 면이 유사한 특징을 보임을 확인할 수 있다.
본 연구에서는 해석모델의 생성을 위하여 하부 균열의 프로파일 자료만을 이용하였다. 주어진 자료는 불규칙한 간격의 점 클라우드(point cloud) 자료 형식을 가졌으므로 수치모델의 메쉬(mesh) 작성을 위한 데이터 변환이 요구되었다. 이를 위해, 해석모델의 원점에 균열 중심이 위치하도록 데이터를 평행이동하고 YZ 평면상에서 42° 회전이동하였으며, XY 그리드상에서 2 mm 간격의 좌표를 생성하기 위한 2차원 내삽(interpolation)을 실시하였다.
3DEC에서 제공하는 균열 생성 기능을 사용하면, 해석도메인을 가로지르는 무한한 평면 균열(infinite planer fracture)만을 모델링할 수 있으나, 개별 블록의 집합체를 사용하면 인접한 블록간의 경계면(contact 요소)에 균열 물성을 부여함으로써 거친 균열면의 생성이 가능하다. 여기에서는 그리드 작성을 위하여 오픈소스 코드인 GMSH (Geuzaine and Remacle, 2009)를 사용하였다. 균열의 3차원 프로파일을 이용해 상부와 하부 시료의 도메인 경계를 정의하고, 각 시료 내부에 사면체 요소를 발생시켰다. 3DEC에서 GMSH 그리드의 좌표 정보를 변환하고, 개별블록을 생성하기 위해서 3DEC의 스크립트 언어인 FISH를 이용하였다. GMSH의 각 사면체 요소는 3DEC의 개별 block에 해당되며, 프로파일 좌표에 해당하는 경계면(contact)에는 균열 물성을 부여하였다.
Fig. 8은 위 과정을 통해 3DEC에서 작성한 해석모델과 실제 시험편의 형상을 보여준다. 실험실 시험에서는 두 표면이 맞물리지 않은 시료(unmated fracture)가 사용되었으나(Fig. 8(a)), 수치모델링에서는 균열면의 맞물림(matedness) 정도에 따라 두 가지 모델을 구성하여(Fig. 8(b), Fig. 8(c)), 그 거동 특성을 비교, 분석하고자 하였다. 아래 본문에서는 실험실 시험편과 동일하게 상하부의 균열면이 맞물려 있지 않은 모델(offset 3 mm)을 UF (unmated fracture) 모델로, 균열면이 완전히 일치하여 맞물려 있는 모델을 MF (mated fracture) 모델로 표기하였다.
Fig. 8.
Grain-based model for simulating thermo-mechanical test on rough fracture: (a) experimental specimen (modified from source: Zhuang et al., 2022), (b) unmated fracture (UF) model and (c) mated fracture (MF) model
MF 모델 생성 시에는 GMSH에서 생성된 사면체 요소 정보를 그대로 사용하여 각 사면체에 일대일 대응하는 block 요소를 정의하였고, UF 모델의 경우 상부 블록의 좌료를 수평이동함으로써 두 균열면이 전단방향으로 3 mm의 offset을 갖도록 하였다. Block 요소의 총 개수는 30,230개로서 균열면에 가까울수록 조밀한 그리드 밀도를 갖는다. 열팽창을 모사하기 위해 각 block을 ‘deformable block’ 요소로 설정하였으며, 해석의 정확도를 위하여 각 block을 더 작은 크기의 사면체 zone 요소로 분할하였다. 3DEC에서 block 요소 간의 경계(interface)에서 발생하는 힘과 변위, 그리고 전단파괴(sliding)와 인장파괴(seperation) 여부는 접촉의 형태(vertex-to-face, edge-to-ege 등)에 따라 contact 요소와 subcontact 요소를 통해 계산된다.
Fig. 8(b)와 Fig. 8(c)의 가장 오른쪽 그림은 균열의 접촉면 위치를 붉은색으로 표시한 것으로, 경계응력( = 0 MPa, = 24.6 MPa, = 3 MPa)을 적용하고 평형해석을 종료하였을 때(가열전), 접촉상태에 놓이는 것으로 판별된 위치를 보여준다. UF 모델의 경우 일부 영역만이 접촉상태에 놓이며, MF 모델에서는 상부와 하부 균열이 대부분 접촉상태에 있는 것을 확인할 수 있다. Fig. 8(a)의 실험실 시험편과 UF 모델을 비교하여 살펴보면, 균열의 형상(trace)이나 접촉면적 분포가 상당한 유사성을 보였다.
본 해석모델에서 사면체 입자(block 요소)는 탄성모델을, 입자간 경계(contact 요소)는 Coulomb slip 모델을 따르는 것으로 가정하였다. Coulomb slip 모델에서는 contact에 작용하는 힘이 강도를 초과하기 전까지 탄성적으로 거동하며, 전단응력이 강도를 초과해 전단파괴에 도달하면 전단응력이 파괴 시의 전단강도로 재설정된다. 인장파괴가 발생하는 경우, 인장강도가 0이 되며 점착력에 의해서만 전단응력에 저항한다. Contact 요소들은 크게 시료 내부를 구성하는 그룹과 균열을 구성하는 그룹으로 나눌 수 있으며, 전자의 경우 암석에 해당하는 그룹으로 상대적으로 높은 전단강도와 인장강도를 가진다. 후자의 경우 주어진 물성에 상응하는 전단강도를 가지지만, 인장강도는 0인 불연속면 경계에 해당한다.
입자 스케일의 미시물성(micro-parameter)을 사용하는 GBDEM에서는 block과 contact의 물성을 조정하는 과정(calibration)을 필요로 한다. 일반적으로 해석모델이 기지의 실험실 시험(압축강도시험 등) 결과를 재현할 때까지 미시물성을 반복적으로 개선함(시행착오법)으로써 최적의 조합을 결정하게 된다. 본 연구에서는 시행착오법과 함께 Park et al.(2020)이 제안한 등가연속체 접근법(equivalent continuum approach)을 이용하여 block과 contact의 물성을 결정하였으며, 이 과정에서 일축압축강도시험(균열 미포함 시료)과 열팽창계수 측정시험 결과를 사용하였다. 암석의 물성을 재현하기 위한 calibration 과정은 Park et al.(2020)과 Park et al.(2022)에 상세히 기재되어 있다.
실험실에서 측정한 인장균열의 마찰계수는 0.86(마찰각 40.70°)으로, 이는 거칠기의 맞물림에 의해 정점응력(peak stress stage) 단계에서 발현된 마찰계수에 해당한다. 균열을 평면으로 모사하는 일반적인 모델링에서는 정점 마찰계수를 통해 거칠기 효과를 반영하게 된다. 그러나 거칠기 형상을 직접적으로 재현한 본 연구의 경우, 균열 물성으로서 정점 마찰계수(peak friction coefficient)가 아닌 기본 마찰계수(basic friction coefficient)나 잔류 마찰계수(residual friction coefficient)를 부여하는 것이 합리적이다. 따라서 균열 contact 요소들에는 saw-cut 균열의 마찰계수인 0.58(마찰각 30.11°)을 적용하였으며, 점착력과 인장강도는 모두 0으로 설정하였다. 팽창각(dilation angle) 역시 거칠기에 의해 발현되는 수직변위를 고려하기 위한 물성이므로 해석에서 별도로 고려하지 않았다. 전단강성과 수직강성의 경우에는 활용가능한 데이터의 부재로 임의의 값을 가정하였다. Table 3은 해석모델의 입력물성과 이를 적용하여 수치적으로 수행한 일축압축강도시험과 열팽창시험의 결과를 요약한 표이다.
Table 3.
Calibrated input parameters, and results of numerical uniaxial compressive strength test and linear thermal expansion test (Park et al., 2022)
수치모델링의 과정은 실내실험과 동일하게 두 단계로 구분된다. 먼저 초기응력이 0, 초기온도가 25°C인 시료에 경계응력을 가한 뒤 준정적 평형해석을 실시하여 시료와 균열면에 목표한 응력 상태를 구현하였다. 시료가 평형상태에 도달한 뒤 Z축 방향으로는 서보컨트롤에 의해 일정응력을 유지하고 Y축 방향의 변위를 구속한 상태에서 상하부의 가열에 의한 과도 열전달(transient heat conduction) 해석을 실시하였다. 이 과정에서 열-역학적 커플링에 의해 암석 블록(zone 요소)의 열팽창 및 열응력의 증가가 계산되며, 열해석과 역학적 해석이 순차적으로 반복된다.
실내시험에서는 히터와 암석 시료 사이에 가압판(로드셀)이 존재하므로, 정확한 해석을 위해 히터, 가압판, 암석 간의 열전달 특성을 정확하게 모사할 필요가 있었다. 그러나 히터와 가압판의 물성 및 사양 등에 대한 정보가 부재하였으므로, 본 연구에서는 이들을 별도의 해석요소로 고려하지 않았다. 히터 온도를 경계조건으로 이용하는 대신, 실험실에서 측정된 시료 하부면의 온도 변화 곡선(Fig. 4의 1번 센서 측정결과)을 수치모델의 상하부 경계면의 온도 조건으로 설정하였다.
한편, 실험실에서 단열재가 설치되지 않았으므로, 시험편의 표면을 통한 열손실을 해석에 고려하였다. 시험편의 전면과 후면에서는 공기의 대류에 의한 열손실이, 좌측면과 우측면에서는 가압판으로의 전도와 가압판에서 공기로의 대류에 의한 열손실이 발생하였을 것으로 추정하였다. 따라서 본 연구에서는 전면과 후면, 좌측면과 우측면에 각각 30 W/m2K, 15 W/m2K의 열손실계수를 적용하였으며, 각 수치는 모니터링 지점(Fig. 4의 2, 3, 6번 센서)에서 측정된 온도 분포를 바탕으로 시행착오법을 통해 결정하였다. Fig. 9는 열해석 시 적용한 경계조건을 보여주는 모식도이다.
Fig. 9.
Thermal boundary condition applied to numerical model for simulating thermoshearing test (modified from source: Zhuang et al., 2022)
4. 해석결과
4.1 온도 분포
해석결과의 분석에 앞서, 실험실 온도 측정자료인 Fig. 4를 살펴보면, 시료의 대칭조건에 의해 4번과 6번, 2번과 5번 센서에서 유사한 온도분포를 보일 것으로 예상할 수 있다. 그러나 4번과 2번의 측정결과에 비해 6번과 5번에서 구한 온도는 낮은 온도 분포를 보였다. 특히, 5번 센서의 경우 히터와 매우 인접한 위치에도 불구하고 더 많은 열손실이 예상되는 3번 센서의 값에 비해 낮은 온도분포를 보였다. 앞서 제3절에서 기술한 바와 같이, 표면에서의 열손실을 반영하기 위해 수치모델의 온도분포가 시료의 온도분포를 재현할 때까지 열손실계수를 시행착오법에 의해 결정하였는데, 5번과 6번 센서의 측정데이터를 함께 고려하는 경우 실험 데이터를 재현하는데에 많은 제약이 있었다. 따라서 2, 3, 6번 센서의 측정데이터만을 이용하여 수치모델 시료 표면의 열손실계수를 결정하였다. 수치해석에서는 4번과 6번, 2번과 5번 센서에 동일한 온도가 계산되었으며, 아래 논의는 2, 3, 6번 센서 위치에서의 온도분포만 비교 분석하였다. 단, 상기한 열손실계수는 Park et al.(2022)의 연구를 통해 saw-cut 균열 시료의 실험실 측정 결과와 수치해석 결과를 바탕으로 결정된 것이다.
Fig. 10은 시간의 흐름에 따른 모니터링 지점에서의 온도 변화를 나타낸 것으로 해석결과를 실험 결과와 비교하여 도시한 것이다. T1, T2, T3, T6은 1, 2, 3, 6번 온도센서(Fig. 4) 위치에서 해석된 결과이다. 두 해석모델의 결과를 실험 결과와 비교하여 살펴보면, 정량적인 수치에는 다소 차이가 있으나, 대체로 유사한 경향을 보였다. UF 모델의 경우, T2 지점에서 실험 결과와 거의 일치하는 결과를 보였으며, 가열 종료 시 T2과 T6 지점에서 실험에서보다 더 낮은 온도에 수렴하였다. MF 모델의 경우, UF 모델과 전반적인 경향은 유사하나, 시료 전체에서 더 높은 온도 분포를 보였다.
MF 모델이 UF 모델보다 더 높은 온도분포를 보이는 이유는, 본 해석모델에서 상하부 시료간의 열전달이 오직 균열의 접촉면을 통한 전도에 의해서만 이루어지는 것으로 가정하기 때문이다. Fig. 8(b)와 Fig. 8(c)에서 살펴본 바와 같이, UF 모델에 비해 접촉면적이 넓은 MF 모델에서 균열을 통한 열전달이 더욱 용이하게 발생하게 된다. Fig. 11은 두 해석모델에서 가열 종료 시(가열 6000초) 해석된 온도 분포와 균열면 접촉상태를 보여준다. 초기 접촉면적이 작은 UF 모델의 경우, 시간 경과에 따라 접촉면의 변화가 거의 없고, 상부 시료와 하부 시료 내에서의 열전달이 거의 독립적으로 진행되는 것을 확인할 수 있다. 그림에서 A 지점과 B 지점의 온도를 살펴보면, 두 지점이 모두 열원과 인접하지만 A 지점이 B 지점에 비해 크게 낮은 온도를 나타낸다. 이는 A 지점에서는 상부경계로부터의 열전달이 거의 이루어지지 않고 하부경계에 의해서만 온도가 상승했음을 의미한다. MF 모델의 경우에는, 가열 초기에는 균열 전체에서 접촉면을 통한 열전달이 발생하다가 전단파괴를 따르는 수직변위 발생과 함께 접촉면적이 감소하면서 상하부 시료간의 열전달이 단절되는 현상을 보였다.
Fig. 12는 MF 모델에서 가열 1570초, 1590초, 3000초 경과 시 하부 시료의 온도 분포와 균열 접촉면적을 보여준다. MF 모델의 초기 접촉면적(Fig. 8(c))과 6000초 경과 시의 해석결과(Fig. 11)를 종합적으로 살펴보면, 1570초 경과 시까지 접촉면적의 변화를 거의 보이지 않다가, 1570 ~ 1590초 사이에 접촉면적이 급격히 감소하여 해석종료 시까지 큰 변화를 보이지 않았다. 1570초에는 A과 B 지점에서 유사한 온도 수준을 보이는데, 이는 상부 경계의 열원에 의한 열전달이 접촉면을 통해 하부 시료의 A 지점까지 이루어지고 있음을 의미한다. 그러나 1590초 이후 접촉면적이 감소하면서 상부 시료에 의한 가열보다는 하부 열원에 의한 열전달이 우세하게 나타났으며, A 지점의 온도는 표면에서의 열손실에 의해 오히려 시간에 따라 감소하는 경향을 보였다. 1590초와 3000초일 때 A와 B 지점에서의 온도를 비교하면, 3000초일 때 B 지점에서의 온도가 증가함에도 불구하고, A 지점에서 온도가 감소함을 확인할 수 있다.
균열을 통한 열전달을 오직 접촉면에서의 열전도로 모사하는 상기 모델의 접근방식은 매우 직관적이며 합리적이다. 현실에서 균열이 포화되어 있지 않다면 비접촉면에서 복사나 공기 대류에 의한 열전달은 미미할 것이며, 접촉상태에 따라 불균등한 열전도가 이루어질 것이다. 앞서 언급한 바와 같이, Fig. 10의 실험 결과는 인장균열 시료가 아닌 saw-cut 균열 시료에서 측정한 결과이다. 만약 인장균열 시료에서 온도가 측정되었다면 saw-cut 균열 시료보다 더 낮은 온도 분포를 보였을 것으로 판단되며, 이에 따라 UF 모델이 실험실 측정데이터에 더 부합하는 결과를 보였을 것으로 추정된다.
4.2 열응력의 증가
Fig. 13은 시간의 흐름에 따른 최대주응력()의 변화를 보여주는 것으로 각 해석모델의 결과와 실험 결과를 비교한 그림이다. ()의 경우, 실험과 수치모델에서 모두 3.0 MPa로 일정하게 유지되었으므로 별도로 도시하지 않았다.
UF 모델의 결과를 살펴보면, 가열이 진행됨에 따라 응력이 24.6 MPa에서 27.1 MPa로 약 2.5 MPa 증가하는 것으로 나타났으며, 3300초를 전후하여 다시 감소하는 경향을 보였다. 실험 결과에서 측정된 최대 응력은 26.0 MPa로 수치해석 모델에 비해 다소 낮은 수치를 보였다. 실험 초반 24.6 MPa의 목표 응력에 비해 더 높은 응력이 측정되었고, 1500초까지 응력이 감소하다가 다시 증가하는 경향을 보인다. 이러한 결과는 실험실과 시험편 조건에 따른 것으로 사료되며, 측정된 데이터만으로 그 원인을 추정하기는 쉽지 않으며 본 수치모델링의 연구범위를 벗어난다. 다만, 1500초 이후 측정된 응력 변화를 해석 결과와 비교해 볼 때, 상대적으로 응력 증가량이 적고 증가율이 작은데 이는 수치해석에 적용한 변위 구속 경계조건의 영향 때문일 것으로 판단된다. 수치해석에서는 수평방향(Y축 방향)의 변위를 완전히 구속한 롤러(roller) 경계조건이 사용되었으나, 실험실에서 시편의 경계는 가압판과 맞닿아 있으며 가압판의 강성과 열팽창에 따라 실제로는 변위가 일부 허용되는 조건이다. 향후, 더 합리적인 모델링을 위해서 가압판과 암석의 상이한 열팽창 특성 및 수평방향 변위 구속조건을 현실적으로 반영할 수 있는 경계조건(일정 수직강성 조건 등)이 적용되어야 할 것으로 판단된다.
MF 모델에서는 UF 모델의 결과와는 달리, 응력이 최대 34 MPa까지 증가하였으며, 정점응력(peak stress)과 잔류응력(residual stress) 단계가 뚜렷이 나타났다. 가열 초반에는 수평응력이 선형적으로 증가하다가 소규모 asperity의 파괴(contact failure)와 전단탈락(shear-off)에 의한 급격한 응력연화(unstable stress sofening)와 수직팽창(normal dilation)이 관찰되었다. 이러한 결과는 MF 모델의 거동이 실험 결과에 부합하지는 않지만, 실제 암석 절리면(rock joint)의 전단과정(Jaeger et al., 2007)에서 발현되는 거칠기 효과를 매우 직접적인 방식으로 모방하고 있음을 보여준다.
4.3 균열 변위
Fig. 14와 Fig. 15는 각각 두 해석모델의 균열 수직변위와 전단변위를 실험 결과와 비교한 그림이다. Fig. 14에서 음의 수직변위는 균열의 닫힘(fracture closure)을, 양의 수직변위는 열림(opening)을 의미하며, 변위량은 두 균열면의 상대변위이다. Fig. 15에서 양의 전단변위는 붉은 화살표 방향을 의미한다.
수직변위의 실험 결과에서는 DIC 분석 결과와 변위 변환기(Transducer) 측정결과가 다소 다른 경향을 나타내었다. Thermoshearing 실험의 특성상 온도와 위치에 따라 균열을 따르는 응력과 변위 분포가 불균등하고, 특히 거친 균열면의 경우 균열의 aspertity 특성와 접촉여부도 수직변위에 영향을 미칠 수 있으므로 한 지점에서의 측정 결과가 대표성을 갖기 어렵다. 따라서 본 연구에서는 실험실에서 구한 DIC 분석 결과만을 수치해석 결과와 비교하였다.
전반적인 결과를 살펴보면, UF 모델이 MF 모델에 비해 더 실험 결과를 잘 재현하고 있음을 확인할 수 있다. 가열 종료 시 UF 모델의 수직변위는 약 10 µm 이내로 실험 결과보다 다소 작은 변위량을 보였으며, 대체로 일치하는 경향을 나타내었다. 전단변위의 경우, 최대 25 µm로 실험 결과보다 현저히 낮은 수준을 보였으나, 변위의 개시 시점과 정성적인 경향성은 대체로 일치하였다. UF 모델에서는 접촉면에서 균열면이 중첩이 발생하고 이에 따라 수직응력의 집중되는 현상이 관찰되었는데, 이에 따라 접촉면에서의 전단파괴가 지연되고 낮은 수준의 변위를 동반하는 것으로 판단할 수 있었다. 이는 근본적으로 자연적 연속성을 갖는 암석 균열의 형상을 이산화된 모델을 통해 단순화하려는 접근법에서 기인한다. 균열면의 형상을 2 mm 간격의 그리드로 모델링하는 경우에는 2 mm보다 작은 규모의 asperity가 갖는 기하학적, 역학적 특성을 배제한다는 것을 전제한다. 즉, 암석 균열의 offset에 따른 접촉면적이나 중첩량을 완벽히 재현하는 것은 사실상 불가능하며, 실험 결과를 바탕으로 수치모델의 반복적 보정이 요구된다. 만약 UF 모델의 offset 모델링 시 수직방향으로 더 큰 이격거리를 가정하였다면, 더 작은 수평주응력 증분과 더 큰 균열 변위를 유도할 수 있을 것이다.
한편, MF 모델의 변위 결과는 수평응력 결과와 유사하게, 실험 결과와는 전혀 다른 거동을 나타내었다. 가열 초반 거칠기의 맞물림에 의한 균열 닫힘 현상이 관찰되었고, 정점응력 발현 시(가열 후 1587초) 급격한 수직팽창과 전단변위를 보였다. 열에 의해 전단응력이 증가하면서 고경사의 균열면에서 contact 요소의 전단파괴가 발생하고, 한 표면이 다른 표면을 타고 이동하면서 전단변위와 수직변위가 급격하게 발생하는 현상이 나타났다.
5. 요약 및 결론
본 논문에서는 3차원 입자기반 개별요소모델(grain-based distinct element model, GBDEM)을 통해 임계응력 조건에 놓인 거친 암석 균열의 열-역학적 커플링과 열에 의한 미끄러짐 거동을 모사하였다. 해석대상은 한국건설기술연구원에서 수행된 ‘Thermoshearing 실험’으로 임계응력 상태의 균열에 열을 가하고 열응력의 발생에 따른 점진적 전단파괴 거동을 살펴본 실험이다. 균열면의 레이저 스캐닝 자료를 바탕으로 균열의 형상을 직접적으로 모사하여 균열면의 맞물림(matedness) 정도와 접촉면적에 따른 전단거동을 비교, 분석하였으며, 불균등한 열팽창 조건에서 가열에 따른 균열의 미끄러짐 거동을 살펴보았다.
시험편의 형상을 동일하게 재현한 UF 모델의 경우, 실험에서 측정된 열전달 특성, 가열에 의한 응력의 증가와 균열의 점진적 전단파괴 프로세스 등을 합리적인 수준에서 모사하고 있는 것으로 나타났다. 균열면이 완전히 맞물린 MF 모델은 비록 실험 결과를 적절히 모사하지는 못하였지만, 일반적인 암석 절리면에서 관찰되는 전단 프로세스(정점응력 단계의 급격한 수직팽창과 전단 미끄러짐, 응력연화, 잔류응력 상태로의 전이 등)를 현실적으로 재현하였다. 또한, 개별요소법의 장점을 활용하여 균열 접촉면적의 변화를 직접적으로 모사함으로써 시료 내 불균등한 열전달 특성과 거친 균열면의 전단거동 메커니즘 등을 정량적으로 분석하였다. 본 연구는 국제공동연구인 DECOVALEX-2023 프로젝트 Task G의 일환으로 수행된 연구로서 향후 Task G에 참여하는 국외 연구팀들과의 의견 교류 및 협력을 통해 지속적으로 개선, 검증할 예정이다.
