1. 서 론
2. DECOVALEX-2023 TASK G
2.1 DECOVALEX 프로젝트 및 Task G 개요
2.2 Task G 2차원 Benchmark Exercise
3. GBM-3DEC 모델을 통한 Benchmark 해석 모델링
3.1 해석모델 개요
3.2 무결함 시험편의 모델링
3.3 평면 균열의 모델링
4. 해석결과 및 토의
4.1 해석 결과
4.2 해석해 비교 검증
4.3 도메인 크기의 영향
5. 요약 및 결론
1. 서 론
암반은 여러 광물들의 조합인 암석과 다양한 스케일의 불연속면들로 이루어진 불균질 재료로서, 암반의 특성을 이해하고 예측하는 데 있어 불연속면에 대한 이해는 매우 중요하다. 미시적 관점에서, 광물 스케일의 여러 요인들로 인해 발생하는 미세균열의 개시와 전파는 실험을 통해 그 메커니즘을 규명하기 매우 어렵고, 암석 물성에 불확실성을 야기한다. 암반에 분포하는 단층, 절리와 같은 거시적 규모의 불연속면 역시 암반의 열적, 수리적, 역학적 특성에 지배적인 영향을 미치며, 암반 및 암반구조물의 불안정성을 초래할 수 있다. 최근 심부 암반 또는 심부 지하공간을 활용한 에너지 관련기술, CO2 지중저장, 방사성폐기물 지층처분 등과 관련된 기술 수요가 증대하면서 다양한 압력과 온도 조건에서 불연속면의 거동에 대한 연구가 증가하고 있으며, 심부 암반의 극한 조건과 제한된 접근성을 고려하여 수치해석적 방법을 통해 그 거동을 평가하려는 시도도 계속되고 있다.
암반공학 분야의 여러 수치해석기법 중 불연속면의 역학적 거동을 가장 직접적인 방식으로 모사하는 기법은 Cundall과 Strack(1979)이 제안한 개별요소법(distinct element method)이라고 할 수 있다. 이 해석기법에서는 불연속암반을 여러 개별 블록(block)의 집합체로 표현하고, 블록 및 블록간 경계(joint)에 구성모델과 입력물성을 부여하여 암반의 거동을 해석하게 된다. 이 모델은 주로 암반 구조물이나 사면의 거동 해석에 적용되어 공학적 스케일에서 불연속면의 영향과 대변형을 해석하기 위한 도구로 널리 사용되어 왔다(Cundall, 1988, Gu and Huang, 2016, Li et al., 2019, Hu et al., 2020). 보다 최근에는 미시적 관점에서 광물입자간 경계(또는 접촉)에서 발생하는 상호작용과 미세균열을 위 모델을 통해 재현하고자 하는 연구들이 많이 보고되고 있다(Lan et al., 2010, Ghazvinian et al., 2014, Park et al., 2017, Fu et al., 2020). 입자기반 개별요소모델(Grain-based distinct element model)로 지칭되는 위 모델링 기법에서는 보통 실험실 스케일에서 암석의 광물학적, 구조적 특징을 작은 개별 블록의 집합체를 통해 재구성하고 그 역학적 거동을 UDEC(Itasca, 2014)이나 3DEC(Itasca, 2017)과 같은 개별요소코드를 통해 계산한다. 위 연구들에서는 모두, 다각형이나 다면체 입자를 사용하는 개별요소모델이 입자 맞물림(interlocking), 입자의 회전과 변위에 대한 저항, 미세균열의 개시와 성장 등을 직접적으로 재현하여 암석의 파괴 메커니즘과 취성거동을 효과적으로 재현할 수 있음을 보고하였다.
본 연구에서는 개별요소코드인 3DEC과 입자기반모델(Grain-Based Model combined with 3DEC, GBM-3DEC)을 이용하여 암석 내 존재하는 단일 균열의 응력과 변위를 분석할 수 있는 모델링 기법을 제시하였다. 이는 DECOVALEX - 2023(DEvelopment of COupled models and their VALidation against EXperiments 2020-2023) 프로젝트 Task G의 일환으로 한국지질자원연구원이 수행한 연구 내용으로, 여기에서는 Task G의 간략한 소개와 함께 현재까지의 해석 결과를 소개하였다.
2. DECOVALEX-2023 TASK G
2.1 DECOVALEX 프로젝트 및 Task G 개요
DECOVALEX 프로젝트는 1992년부터 수행되고 있는 국제공동연구로서 방사성폐기물의 지층처분과 관련하여 공학적 방벽과 자연 방벽에서의 열-수리-역학-화학적(thermal-hydrological-mechanical-chemical, THMC) 연계거동을 평가할 수 있는 수치해석기법의 개발을 목표로 현재까지 다양한 연구가 이루어져 왔다(Birkholzer et al., 2019). DECOVALEX 프로젝트에서는 4년을 주기로 다양한 Task에 대한 공동연구가 수행되어 왔으며, 현재 2020년 4월에 DECOVALEX-2023이 착수되어 2023년까지 7가지 Task에 대해 공동연구를 수행할 예정이다.
DECOVALEX-2023의 정식회원기관(약칭, 국가)은 Agence nationale pour la gestion des déchets radioactifs(ANDRA, 프랑스), Federal Office for the Safety of Nuclear Waste Management (BASE, 독일), Federal Company for Radioactive Waste Disposal (BGE, 독일), Federal Institute for Geosciences and Natural Resources(BGR, 독일), Chinese Academy of Science (CAS, 중국), Canadian Nuclear Safety Commision (CNSC, 캐나다), Center Organisation for Radioactive Waste (COVRA, 네덜란드), Department of Energy(DOE, 미국), Empresa Nacional de Residuos Radioactivos(ENRESA, 스페인), Swiss Federal Nuclear Safety Inspectorate(ENSI, 스위스), Japan Atomic Energy Agency(JAEA, 일본), 한국원자력연구원(KAERI, 한국), Nuclear Waste Management Organization(NWMO, 캐나다), Radioactive Waste Management(RWM, 영국), Radioactive Waste Repositories Administration(SÚRAO, 체코), Swedish Radiation Safety Authority(SSM, 스웨덴), Taiwan Power Company(Taipower, 대만)의 17기관(14개국)이며, 관련 연구기관은 총 50기관이다(Birkholzer and Bond, 2020). 국내에서는 한국원자력연구원이 DECOVALEX-2011(2008-2011)부터 정식회원기관과 연구기관으로 참여하고 있으며, 한국지질자원연구원은 DECOVALEX-2019(2016-2019)부터 한국원자력연구원의 협력연구기관으로 참여하고 있다(Lee et al., 2020).
Task G는 독일 Freiberg University of Mining and Technology(TUBAF)/Helmholtz Centre for Environmental Research (UFZ), 스코틀랜드 Edinburgh 대학교, 독일 DynaFrax에서 공동 제안한 Task로 결정질 암반 내 균열의 열-수리-역학적 복합거동을 해석할 수 있는 수치해석기법을 개발, 비교, 검증하는 데에 목표가 있다. Task G는 기본적으로 실내실험 결과를 모사하고 각 연구팀들의 수치해석 결과를 상호비교 및 토의하는 형식으로 진행되며, 필요에 따라 해석해(analytical solution)와 비교하는 benchmark 모델링과 blind prediction 모델링을 병행한다. Task G의 주요 연구 테마는 다음과 같이 요약할 수 있다.
•Step 1: Mechanical (M) process in rock fracture – TUBAF에서 수행된 일정수직하중(constant normal load, CNL) 및 일정수직강성(constant normal stiffness, CNS) 조건의 직접전단시험 결과를 바탕으로 암석 균열의 역학적 거동 모델링 기법 수립
•Step 2: Coupled hydro-mechanical (HM) process in rock fracture – Edinburgh 대학교에서 수행된 GREAT cell 수리전단시험 결과(McDermott et al., 2018)를 바탕으로 다양한 응력 조건에서 발생하는 수리역학적 연계거동 모델링 기법 수립
•Step 3: Coupled thermo-mechanical (TM) process in rock fracture – 한국건설기술연구원에서 수행중인 고온조건 진삼축압축시험 결과를 바탕으로 열응력으로 인한 균열 미끄러짐 모델링 기법 수립
•Step 4: Combining and upscaling near-field approaches for THM analysis
Task G의 연구 테마는 크게 M-HM-THM 연구단계와 M-TM-THM 연구단계로 구성되며(Fig. 1), 각 연구는 benchmark exercise(BE), experimental analysis(EA), blind prediction(BP), up-scaling(UP)의 순으로 진행될 예정이다. 각 참여팀들은 관심도나 필요에 따라 연구 테마를 선택하고, 개별적 또는 통합적으로 연구를 진행하게 된다. 현재 BE와 EA에 대한 상세 계획이 수립된 상태이며, BP와 UP의 연구내용은 추후 논의를 거쳐 구체화될 예정이다. Table 1은 Task G의 참여기관과 해석기법, 해석코드를 보여준다.
Table 1.
DECOVALEX-2023 Task G participating research teams, numerical methods and codes
BE 연구단계에서는 참여연구팀들이 단순화된 문제에 대해 모델링을 실시하고, 연구팀 간 해석결과 비교 또는 해석해와의 비교를 통해 각자의 해석기법을 개선하고 검증하게 된다. 이를 위해 2차원(평면변형률 조건) 모델과 3차원 모델에 대한 benchmark 해석이 고려되며, 3차원 모델은 GREAT cell 시험(HM)과 고온 조건 진삼축압축시험(TM)에 사용된 응력조건과 시편의 형태를 고려해 결정되었다(Fig. 2).
Fig. 2
Benchmark exercises of Task G - Step 1 (M), Step 2 (HM), Step 3 (TM) and Step 4 (THM): (a) cubic model for triaxial test and (b) cylinderical model for polyaxial test (Kolditz et al., 2020)
2.2 Task G 2차원 Benchmark Exercise
Task G 연구는 2020년 9월에 본격적으로 착수되어 11월에 개최된 제2차 워크숍에서 2차원 benchmark 해석에 대한 결과 보고와 논의가 진행되었다. 2차원 benchmark 모델링의 주요내용은 암석 내부에 단일한 균열이 포함되어 있을 때 압축력에 의해 암석과 균열면에서 발생하는 역학적 거동을 해석하는 것으로, 모델링 결과를 Pollard와 Segall(1987)의 해석해와 비교하거나 연구팀 간 상호 비교함으로써 각 연구팀의 수치모델링 기법을 검증 개선하게 된다. 본 논문의 주요내용은 제2차 워크숍에서 발표된 한국지질자원연구원의 해석 결과를 정리한 것이다.
Fig. 3은 해석모델의 형상과 경계조건(σx = 5 MPa, σy = 10 MPa), 모니터링 지점을 보여준다. 3차원 benchmark 모델에서는 암석 시험편을 분할하는 균열(infinite fracture)이 고려되는 반면(Fig. 2), 2차원 모델에서는 암석 내부에 포함된 균열(finite fracture)을 가정한다. 해석 도메인은 한 변의 길이가 0.2 m인 정사각형이며, 내포된 균열은 길이가 0.17 m이다. 균열의 거칠기에 따라 매끄러운 평면(planar fracture)과 거친 면(rough fracture) 두 가지가 고려되며, 거친 면의 모델링을 위해서는 TUBAF의 직접전단실험에 사용된 암석 절리면의 3차원 거칠기 측정 자료(좌표 정보)가 이용된다. 균열의 경사는 수평면으로부터 0°, 30°, 45°, 90°로 기울어져 있다고 가정하므로 거칠기와 경사에 따라 총 8개의 해석 케이스가 고려된다. Fig. 3의 숫자는 모니터링 지점을 보여주는 것으로 지점 1 ~ 7(P1 ~ P7)은 균열면 상에, 8과 9(P8 ~ P9)는 시료 내부에 위치한다. P1은 균열의 중심에 있으며, P2과 P5, P3과 P6, P4와 P7은 각각 중심으로부터 0.02, 0.05, 0.085 m의 이격거리를 갖는다. P8과 P9는 좌표축을 기준으로 1사분면과 3사분면의 중심에 위치한다. P1 ~ P7에서는 균열의 전단변위를, P8, P9에서는 시험편의 응력(σxx, σyy, σxy) 상태를 측정한다.
Pollard와 Segall(1987)의 해석해는 균열을 포함한 탄성재료의 응력과 변위에 대한 이론해로서 Fig. 4는 위 연구에서 사용된 좌표계와 균열의 개념도를 보여준다. 좌표계의 원점은 균열 중심에 위치하며, 1, 2, 3 좌표축은 각각 균열의 인장(normal), 평면내 전단(in-plane shear), 평면외 전단(out-of plane shear) 방향을 가리킨다. 균열은 1-2 평면상에서 길이가 2a이며, 3축 방향으로 무한대의 길이를 갖는 매우 얇은 두께의 절단면 또는 공극(cut or void)으로 가정된다. 위 이론에 따르면, 균열에 Δσ1, Δσ2, Δσ3의 응력이 가해질 때, 각 방향으로 발생하는 균열의 변위 Δu1, Δu2, Δu3는 식 (1)과 같이 표현된다. 이때 변위는 정확히 균열을 구성하는 두 면의 상대변위를 의미한다.
여기서 υ는 포아송비, G는 전단변형계수, a는 균열길이의 절반, x2는 2축 방향 좌표를 의미한다.
2차원 평면변형률 조건에서 균열면에 작용하는 수직응력과 전단응력을 각각 Δσn과 Δτ이라고 하면, 이로 인해 발생하는 수직변위(Δun)와 전단변위(Δus)는 식 (1)을 변형하여 식 (2)와 (3)으로 표현된다.
여기서 E는 탄성계수이며, r은 균열 중심으로부터의 거리를 의미한다.
3. GBM-3DEC 모델을 통한 Benchmark 해석 모델링
3.1 해석모델 개요
본 연구에서는 3차원 개별요소코드인 3DEC으로 구현한 입자기반모델(GBM-3DEC)을 이용하여 benchmark 해석을 수행하였다. 3DEC은 Itasca의 상용해석프로그램으로 지반, 암반공학 분야에서 불연속 암반의 거동 해석을 위해 널리 활용되고 있다. 3DEC에서는 해석 대상의 거동을 블록(block)과 불연속면(joint)의 조합을 통해 모델링한다. 블록은 강체 블록(rigid block) 또는 변형가능한 블록(deformable block)으로 설정할 수 있으며, 변형가능한 블록이 사용되는 경우, 이는 더 작은 사면체 요소로 분할되어 응력과 변형이 계산된다. 두 블록 사이에 평면 형태의 joint 요소가 정의되면, 평면상에 수많은 접촉점(contact)들이 생성된다. 접촉점에서는 힘-변위 관계, 구성모델, 입력물성에 의해 접촉력(contact force)과 joint의 개폐(separation)와 미끄러짐(slip) 여부가 결정된다. 3DEC의 연산은 외재적 시간 전진 알고리즘(explicit time marching algorithm)을 기반으로 하며, 접촉점에서의 힘-변위 법칙과 개별요소의 운동방정식으로 구성된다. 매 timestep마다 개별블록에 작용하는 합력으로부터 블록의 운동이 결정되고, 접촉점의 변위로부터 새로운 접촉력이 계산된다. 새로운 접촉력은 다음 timestep에서의 연산에 반영되고 동일한 연산과정이 반복된다.
입자기반 개별요소모델의 입력물성 산정방식은 재료의 물성을 직접 입력자료로 활용하는 다른 모델링 기법과는 다소 상이하다. 일반적인 암석역학 문제에서 입자 스케일에서의 물성은 미지의 값이거나 측정이 불가하기 때문에, 시뮬레이션을 수행하기 전 해당 모델이 기지의 거시적 거동(일반적으로 실험실 스케일)을 재현할 수 있도록 미시물성(입자나 접촉의 물성)을 조정하는 과정(micro-parameter calibration)이 필요하다. 암석의 경우 보통 압축강도나 인장강도 시험이 활용되며, 내부에 균열이 포함되어 있다면 전단시험을 고려할 수도 있다. 이 같은 이유로 본 연구의 benchmark 해석을 위해서는 크게 두 가지 시뮬레이션이 수행되었다. 먼저 무결함 시료에 대한 압축시험을 수행하여 미시물성을 결정하였으며, 이를 바탕으로 응력 재하에 따른 균열의 거동을 모델링하였다. 균열 미시물성의 경우, 조정과정을 별도로 수행하지 않고 Pollard와 Segall (1987)의 이론을 재해석함으로써 결정하였다. 이와 관련된 자세한 내용은 3.3절에 기술하였다.
3.2 무결함 시험편의 모델링
본 연구에서는 사면체 강체 블록의 조합을 통하여 균열을 포함한 암석 재료를 모델링하였다. 여기에서는 해석도메인의 크기가 0.2 m×0.2 m×0.005 m인 얇은 두께의 3차원 모델을 생성하였으며, 수치해석 시 평면변형률 조건을 설정하여 0.2 m×0.2 m인 평면 내에서 발생하는 2차원 거동을 살펴보았다.
일반적으로 3DEC에서 균열을 모델링할 때에는 해석 모델의 경계까지 연장된 무한한 균열를 가정하지만, 작은 개별 블록들을 사용하면 일련의 블록 경계에만 균열 물성을 부여함으로써 도메인 내부 포함된 균열을 모델링하는 것이 가능하다. 본 연구에서는 개별 블록의 생성을 위해 오픈소스 소프트웨어인 Gmsh(Geuzaine and Remacle, 2009)를 이용하였다. Fig. 5는 Gmsh에서 구성된 45° 경사의 평면 균열 모델 생성과정을 보여주는 것으로 모니터링 지점의 위치와 균열의 위치를 반영해 사면체 요소들이 생성될 수 있도록 하였다. 그림에서 P8과 P9 지점을 살펴보면, 서로 교차하는 두 평면이 그리드 경계로 설정된 것을 확인할 수 있는데, 이는 두 지점에서 작용하는 응력(σxx, σyy, σxy)을 산정하기 위해 정의된 것이다. 강체 블록을 사용하는 경우 블록 내부의 응력을 계산하지 않는다. 따라서 그 대안으로서 X, Y 방향에 수직한 두 평면을 삽입하여, 이 평면상에 작용하는 contact 수직응력과 contact 전단응력을 통해 위 응력값들을 계산하고자 하였다. 한편, Gmsh에서 생성된 사면체 그리드의 좌표정보를 3DEC 블록(edge, vertex)의 정보로 변환하여 해석모델을 생성하였으며, Fig. 6은 45° 경사의 평면 균열을 포함한 3DEC 해석모델을 보여준다. 생성된 모델의 블록의 개수는 약 12,500개로 평균 부피를 통해 계산된 블록의 등가 직경(equivalent spherical diameter)은 약 3.0 mm이다.
등가 연속체 접근법(Equivalent continuum approach)은 암반 내 절리의 밀도가 높고, 절리 간격이 구조물의 크기에 비해 상당히 적을 때 불연속 암반의 변형거동을 등가의 연속체로 간주하여 해석하는 방법이다. 등가 연속체 개념에서 암반의 탄성변형은 다음 식과 같이 기술된다.
여기서 Em과 Gm은 각각 암반의 등가탄성계수와 등가전단탄성계수, Er과 Gr은 각각 무결암의 탄성계수와 전단탄성계수, kn과 ks는 절리의 수직강성과 전단강성, s는 절리의 간격이다.
입자기반모델에 등가 연속체 개념을 도입하면, 모델 전체의 변형은 입자의 변형과 접촉점의 변형에 의해 결정되며, 두 인자의 기여도는 입자의 간격(크기)에 따라 좌우된다. 본 연구에서와 같은 강체 블록 모델의 경우, 전체 시료의 변형은 오직 입자 크기와 접촉점의 강성에 의해 제어된다.
본 연구에서는 입자기반 모델에 대한 압축시험을 수치적으로 모델링하여, 암석의 탄성거동을 재현할 수 있도록 접촉점(joint 또는 contact)의 수직강성(kn)과 전단강성(ks)을 반복 조정하였다. 압축시험에서는 시료 상하단에 재하판을 설치하고, Y축 방향으로 일정 속도를 부여하여 σy = 30 MPa이 될 때까지 가압을 진행하였고, 응력-변형률 곡선의 기울기로부터 탄성계수와 포아송비를 측정하였다. Benchmark 해석에 고려되는 암석의 일축압축강도는 약 120 MPa로 주어진 응력조건 하에서 탄성거동을 하는 것으로 가정하였다. 따라서 시료 내부의 입자간 접촉점의 인장강도, 점착력, 마찰각을 각각 50 MPa, 50 MPa, 70°로 높게 설정하여, 압축시험 중 시료에 파괴가 유도되지 않도록 하였다. 암석 시험편의 탄성계수와 포아송비는 각각 49.75 GPa, 0.26이다.
Fig. 7은 블록 간 접촉의 수직강성과 전단강성이 4.0×104 GPa/m와 1.27×104 GPa/m일 때의 압축시험 결과를 보여주는 것으로 (a)는 블록의 변위벡터를, (b)는 응력-변형률 곡선을 보여준다. 응력-변형률 곡선으로부터 구한 GBM-3DEC 모델의 탄성계수와 포아송비는 49.27 GPa과 0.26로 나타났으며, 이는 암석 시험편의 탄성정수와 거의 일치하는 값이다.
3.3 평면 균열의 모델링
본 해석모델에서 암석 시료 내에 포함된 단일 균열은 3.2절에서 기술한 암석시험편과는 별개의 물성을 갖는 불연속면으로 가정된다. 균열의 물성을 결정하기 위해서는 전단시험을 통해 얻어진 강성이나 강도를 이용하여 미시물성 조정과정을 수행하는 것이 일반적이나, 본 연구에서는 Pollard와 Segall (1987)의 이론을 통해 균열의 물성을 결정하고자 하였다.
식 (2), (3)을 변형하면, 응력과 변위의 관계를 식 (6), (7)과 같이 표현할 수 있으며, 이들은 수직강성과 전단강성의 정의와 정확히 일치하는 개념이다.
식 (6)과 (7)을 살펴보면, 전단강성과 수직강성은 동일한 형태로서 균열중심과의 거리의 함수로 표현된다. 암석 절리면에서 측정되는 수직강성은 일반적으로 전단강성에 비해 현저히 클 뿐만 아니라, 동일한 균열의 물성을 위치의 함수로 표현하는 것은 다소 부적절한 측면이 있다. 위 이론에서는 균열을 탄성재료 내 얇은 결함이나 공극으로 가정한다. 따라서 해석해에서는 균열 자체의 강성이나 강도가 고려되지 않으며, 균열이 미끄러진 이후에도 탄성정수에 의해 변위가 제어되는 것으로 간주한다.
사실상 이러한 방식은 일반적인 암석역학에서 균열의 거동을 기술하는 방식과 크게 다르다. 그러나 본 연구에서는 해석해의 재현에 목표를 두고, 탄성재료 일부로서 균열의 거동을 모델링하고자 하였다. 균열의 인장강도와 전단강도는 모두 0으로 설정하였으며, 식 (6)과 (7)을 통해 계산된 수치를 수직강성과 전단강성으로 적용하였다(kn = ks = 1.57×102 GPa/m). 3DEC에서 전단강도와 인장강도가 0인 균열은 아주 작은 응력을 적용하더라도 재하와 동시에 파괴상태로 인식된다. 즉 탄성영역을 벗어나게 되므로, 상기한 강성을 통해서는 균열의 변위 제어가 불가능하게 된다. 따라서 본 연구에서는 잔류전단강도와 잔류인장강도를 임의로 부여하여 파괴상태에 이른 이후에도 수직강성과 전단강성을 통해 선형적으로 변위가 제어될 수 있도록 하였다. 해석모델의 초기조건은 시료 내부에 응력이 작용하지 않는 상태이며, XY 평면의 경계에 일정 응력조건(σx = 5 MPa, σy = 10 MPa)을 설정하였다. Z축 방향으로는 변위를 구속하여 평면변형률 조건에서 정해석을 통해 균열의 변위를 유도하였다.
4. 해석결과 및 토의
4.1 해석 결과
Figs. 8~11은 각 균열면의 경사에 따른 해석결과를 보여주는 것으로 각 그림에서 (a)는 블록 절점(vertex)에서의 변위 contour를, (b)와 (c)는 contact에서의 수직방향과 전단방향 변위벡터를 보여준다.
전반적인 결과를 살펴보면, 해석모델의 변위는 대체로 균열을 따라 발생하며, 암석의 변형 자체는 미미한 수준이다. 앞서 기술한 바와 같이 개별 블록에서는 변형이 허용되지 않으므로 모델 전체의 변위는 입자간 접촉에서의 변형 특성에 지배된다. 따라서 무결함 시료 내부에 비해 상대적으로 낮은 강성을 부여한 균열면에서 변위가 지배적으로 발생하는 것은 사실상 당연한 해석결과이다. 모든 해석 케이스에서 균열은 압축하중 하에 놓이게 되며, 균열면 경사가 작을수록 더 큰 압축 변형을 보였다. 또한 균열의 중심에서 최대 수직변위를 보였다. 전단변위의 경우, 균열 경사가 0°와 90°인 경우에는 거의 발생하지 않았고 균열의 양 끝단에 위치한 접촉점들에서 응력집중에 의해 상대적으로 큰 전단변위를 보였다. 단, 위 접촉점들은 균열면이 아닌 시료 내부에 위치하며, 전단방향이 균열 전단방향과 상이하므로 균열 양끝단에서 최대 전단변위가 유도되었음을 의미하지는 않는다. 균열의 경사가 45°인 경우, 30°인 경우에 비해 더 큰 전단응력과 전단변위가 유도되었다.
Table 2는 모니터링 지점 P8과 P9에서의 응력상태를 나타낸다. 두 지점에서는 해석모델의 대칭성으로 인해 대체로 유사한 결과가 계산되었다. 모델 경계에 적용된 응력 수준(σx = 5 MPa, σy = 10 MPa)을 고려해 보면, 균열의 변위로 인해 시험편 내부에서 응력의 재배치가 일어난 것을 확인할 수 있다. 예를 들면, 경사가 0°일 때와 90°일 때의 응력을 살펴보면, 균열과 평행한 방향의 응력은 경계 응력과 거의 유사한 반면, 균열에 수직한 방향으로는 경계 응력보다 크게 낮은 수준을 보이는 것으로 나타났다. 모든 해석 케이스에서 균열의 수직변위가 전단변위에 비해 크게 발생하였음을 감안할 때, 균열의 수직변위로 인해 시료 내부의 응력이 크게 영향을 받는 것으로 해석할 수 있다.
Table 2.
Stresses calculated at P8 and P9
4.2 해석해 비교 검증
Task G의 benchmark 해석에서는 균열의 전단변위를 모사하는 데에 우선순위를 두었으나, 균열의 거동 해석을 위해서는 전단변위는 물론 수직변위의 적절한 모델링이 매우 중요하다. 균열면의 역학적 수직변위는 수리간극의 크기와 직접적인 연관이 있어 수리유동을 고려하는 경우, 매우 중요한 고려사항이다. 따라서 본 연구에서는 수직변위와 전단변위의 수치모델링 결과를 상기한 해석해와 비교함으로써 제시한 모델링 기법의 타당성을 검증하고자 하였다.
식 (2)와 (3)에서 확인할 수 있듯이, 균열 변위에 대한 해석해를 구하기 위해서는 우선 균열에 작용하는 응력 조건을 알아야 한다. 2차원 문제에서 임의 방향으로 기울어진 평면에 작용하는 수직응력과 전단응력은 힘과 모멘트 평형해석을 통해 식 (8)과 (9)로 표현된다.
여기서 σn와 τ는 기울어진 평면에 작용하는 수직응력과 전단응력, σx, σy, τxy는 요소 경계에 작용하는 수직응력과 전단응력, θ는 평면의 법선벡터가 X축과 반시계방향으로 이루는 각도이다.
위 식들의 부호 규약은 Fig. 12(a)를 따른다. Fig. 12(b)는 본 연구의 benchmark와 같이 연속체 모델 내에 유한한 길이의 균열이 포함되어 있는 경우를 보여준다. Table 3은 균열의 기울기(α = 90 - θ)에 따라 식 (7)과 (8)에 의해 계산된 수직응력과 전단응력을 수치해석에서 얻어진 값과 비교한 표이다. 여기서, 수치해석 결과는 균열 중심의 접촉점에서 획득한 값이다. 표에서 확인할 수 있듯이, 균열면의 방향이 주응력 방향과 일치하는 경우(α = 0°, α = 90°), 이론적으로 균열에 전단응력이 작용하지 않아야 한다. 그러나 수치해석모델에서는 하중 재하에 따른 블록들의 변위로 인해 균열을 따라 전단응력이 유도되었을 뿐만 아니라, 이론적으로 계산된 응력에 비해 낮은 수준의 응력이 유도되었음을 확인할 수 있다. 이는 수치모델에서 균열과 무결함 시험편의 상이한 변형 특성으로 인해 응력이 재배치되었음 의미한다. 이러한 결과는 이미 Table 2를 통해서도 확인한 바 있다.
Table 3.
Stress acting on the embedded fracture
Table 3의 응력 값과 식 (2), (3)을 이용하면, 균열면을 따르는 수직변위와 전단변위에 대한 해석해를 계산할 수 있으며, Fig. 13과 Fig. 14은 수치해석 결과와 함께 해석해를 비교한 그림이다. 그림에서 수치해석 결과(scatter plot)는 균열 상의 접촉점들에서 계산된 결과이며, 실선은 이론적인 응력을 사용한 경우의 해석해, 점선은 균열 중심에 실제 작용하는 응력을 사용한 경우의 해석해를 보여준다. 각 그림에서 r은 균열 중심으로부터의 거리를 나타내며, 균열 좌표가 1, 4분면에 있는 경우 양의 부호를, 2, 3분면에 있는 경우 음의 부호를 부여하였다(Fig. 3 참고).
Fig. 14(a)와 (d)를 살펴보면, 균열 경사가 0°와 90°인 경우 균열을 따라 전단변위가 무작위적인 분포를 보이는 것을 확인할 수 있는데, 이는 균열을 따르는 블록간 접촉점에서 미세한 변위가 발생하기 때문이다. 이러한 결과는 여러 입자의 접촉을 통해 단일한 균열을 표현하는 본 해석모델의 기하학적 특성에 기인한 것으로, 대체로 0.5 μm 이하의 작은 수치를 고려할 때 균열을 따라 전단변위가 거의 발생하지 않는 것으로 해석할 수 있다 (Fig. 8(c), Fig. 11(c) 참고). 또한, 균열을 따라 산발적으로 발생하는 변위와 응력 분포를 고려할 때 해석해와의 비교도 무의미하다.
위 해석결과를 제외한 나머지 결과들을 살펴보면, 수치모델에서 균열에 실제 유도되는 응력이 이론적 응력에 비해 작기 때문에(Table 3), 두 응력을 이용한 해석해가 서로 큰 차이를 보였다. 그러나 실제 균열 중심에 유도된 응력을 사용하는 경우, 모든 해석케이스에서 모델링 결과와 해석해가 거의 일치하였다. Pollard와 Segall(1987)의 이론은 기본적으로 균열면에 작용하는 응력을 전제한다. 따라서 본 연구의 해석모델이 위 해석해를 효과적으로 재현하고 있는 것으로 판단할 수 있다. 다만, Table 3에 보인 응력의 불일치에 대해서는 추가적인 논의가 필요하다.
4.3 도메인 크기의 영향
4.2절에서 살펴본 바와 같이, 본 연구의 수치해석모델에서는 암석 내 균열에서 예상보다 낮은 응력수준이 유도되었다. 따라서 도메인 크기에 따른 수치모델링을 추가적으로 수행하여, 도메인과 균열의 상대적 크기가 응력에 미치는 영향을 검토하였다.
Fig. 15는 도메인 크기에 따른 해석모델을 보여주는 것으로, 해석 도메인 한 변의 길이는 각각 0.12 m(Case 1), 0.2 m(Case 2, benchmark와 동일), 0.5 m(Case 3)이다. 균열 기울기가 45º인 경우를 대상으로 하였으며, 도메인 크기를 제외한 나머지 조건은 benchmark 해석조건과 동일하다. Case 1의 경우, 해석모델을 관통하는 조건이며, Case 3은 benchmark 해석조건에 비해 더 큰 도메인 사이즈를 설정하였다.
Fig. 16은 각 해석케이스의 수치해석 결과로서 균열을 따라 유도된 수직응력과 전단응력을 보여준다. 균열이 도메인을 완전히 가로지르는 Case 1의 경우, 균열에 작용하는 수직응력과 전단응력은 응력해석으로부터 계산된 이론치와 거의 동일하게 나타났다. 또한, Case 2와 Case 3의 해석결과에서는 도메인의 크기가 커질수록 균열에 유도되는 응력 수준이 낮아지는 것을 확인할 수 있다. 이러한 결과는 본 연구의 해석모델에서 균열에 유도되는 응력은 결국 균열 변위의 구속 여부와 관련이 있음을 시사한다. Case 1의 경우에는 균열의 변위가 모든 방향으로 자유롭지만, Case 2와 Case 3에서는 암석 내부의 높은 강성으로 인해 균열의 변위가 구속된다. 또한, 도메인의 크기가 클수록 균열 변위는 더욱 제한적이게 되며, 변위를 기반으로 힘을 계산하는 개별요소모델의 연산과정을 고려할 때 균열에 낮은 응력이 유도되는 것은 자연스러운 해석 현상이다. 이러한 결과가 실제 암석 균열에서 발생하는 메커니즘을 반영하는 것인지, 또는 본 수치해석모델의 개선점일지에 대해서는 향후 추가적인 연구가 수반되어야 할 것으로 판단된다.
5. 요약 및 결론
본 논문에서는 국제공동연구인 DECOVALEX-2023 프로젝트 Task G의 현황과 현재까지의 연구내용을 소개하였다. Task G는 결정질 암반 내 균열의 생성과 성장 메커니즘 및 균열에서 발생하는 열-수리-역학적인 복합거동을 해석하기 위한 수치해석기법을 개발하고 이를 실내실험에 적용하여 타당성을 검증하는 데에 목적이 있다.
본 연구에서는 Task G의 1단계 연구로서 암석 내 포함된 단일 균열의 역학적 거동을 3차원 입자기반 개별요소모델을 통해 모사하고, 그 결과를 해석해와 비교함으로써 모델링 기법의 타당성을 살펴보았다. 해석 결과, 압축하중 하에서 균열과 암석의 불균등한 변위로 인해 응력 재배치가 발생하였고, 균열에 유도되는 응력 수준은 응력해석을 통해 이론적으로 계산된 수치보다 낮게 나타났다. 하지만 수치모델을 통해 계산된 균열의 수직변위와 전단변위는 유도 응력을 바탕으로 계산된 해석해와 거의 동일하게 나타났으며, 본 연구의 해석모델이 탄성재료 내 포함된 균열의 응력-변위 거동을 합리적으로 재현할 수 있음을 확인하였다. 본 연구의 수치모델링 기법은 Task G에 참여하는 국외 연구팀들과의 의견 교류와 워크숍을 통해 지속적으로 개선하는 한편, 향후 다양한 조건의 실내시험에 적용하여 타당성을 검증할 예정이다.
