1. 서 론
2. 시험 방법
2.1 절리 시험편 준비
2.2 절리 표면 좌표 측정 및 후처리 작업
2.3 수리시험 수행
3. 시험 결과 및 논의
3.1 절리 표면의 기하학적 특성
3.2 수리시험 결과
4. 수치해석
4.1 모델 구성
4.2 해석 결과
5. 결 론
1. 서 론
암반은 크게 무결암과 불연속면으로 구성된다. 무결암은 확인 가능한 균열이 없거나 매우 적은 여러 광물의 집합체를 의미하며, 불연속면은 암반의 물리적 연속성을 끊는 모든 면 구조를 의미한다. 불연속면은 규모와 성인(genesis)에 따라 단층, 층리, 엽리 등으로 구분할 수 있으며, 이중 관찰 빈도 측면에서 가장 대표적인 암반 불연속면은 절리이다. 절리는 암반의 전반적인 열-수리-역학적(Thermal-Hydraulic-Mechanical, THM) 거동에 지대한 영향을 미친다. 역학적 관점에서 절리는 연약면(weak plane)에 해당하며, 수리적으로는 유체의 주요 유동 통로 역할을 한다. 특히, 매질 자체의 투수성이 매우 낮은 결정질 암반의 경우, 대부분의 지하수 유동은 절리를 통해 발생한다. 고준위방사성폐기물처분장에서 핵종이 유출된 경우, 방사성 물질의 위해성이 인간 생활권까지 도달하는 주요 경로는 지하수를 통한 이동이며, 결정질 암반에서 지하수는 상술한 것처럼 대개 절리를 통해 유동한다. 따라서 절리를 통한 수리 유동 특성을 정확하게 파악하는 것은 처분장 부지 특성화, 공학적방벽 설계 및 장기 안전성 확보 측면에서 매우 중요한 사항이다.
암석 절리를 통한 수리 유동은 암반공학 및 토목, 지질학 분야의 오랜 연구 주제 중 하나이기 때문에 이에 대한 이론적, 실험적, 수치해석적 연구는 광범위하고 지속적으로 수행되었다(Witherspoon et al., 1980, Brown, 1987, Zimmerman and Bodvarsson, 1996, Barton and de Quadros, 1997, Yeo et al., 1998, Lee, 2002, Singh et al., 2015, Yu et al., 2022). 절리를 통한 비압축성 유체의 유동은 관성력(유체 가속도), 점성력(내부 마찰), 압력 구배, 외력 등의 항을 고려한 Navier-Stokes 방정식으로 표현될 수 있다. 여기서 비선형 항의 영향을 무시하면(즉, 관성력의 영향을 무시) Stokes 방정식이 유도되며, 더 나아가 평판 조건(smooth parallel plate)을 가정하면 잘 알려진 삼승 법칙(cubic law)이 유도된다. 삼승 법칙은 적용의 편이성으로 인해 널리 활용되었으나, 실제 절리를 통한 수리 유동과는 많은 차이를 보일 수 있다(Witherspoon et al., 1980). 이러한 오차의 원인으로는 절리 거칠기, 접촉 면적, 간극 분포 등을 들 수 있으며, 이를 반영하기 위해 기존 식에 보정 계수를 추가하는 형태의 연구가 다수 수행된 바 있다.
절리를 통한 유체 유동에서 가해진 수압 구배(pressure gradient)와 배출 유량 사이에는 선형적 혹은 비선형적 관계가 존재한다. 유동이 층류 영역(laminar flow regime)에 속하는 경우, 압력과 유량 사이에는 선형적 관계가 존재하며 이는 Darcy law로 표현할 수 있다. 한편, 비선형 관계를 보이는 경우, Forchheimer 방정식 혹은 Izbash 방정식 등이 적용될 수 있으며, 이러한 비선형성의 원인으로는 관성 효과, 유체 주입 조건 그리고 절리 표면의 기하학적 특성 등을 들 수 있다(Zeng and Grigg, 2006, Javadi et al., 2014, Zoorabadi et al., 2015, Rong et al., 2020). 절리 표면 형상에 의해 결정되는 기하학적 물성은 거칠기(roughness), 간극(aperture), 접촉 면적(contact), 맞물림(matedness), 공간 상관성(spatial correlation) 등 다양하며(Hakami and Larsson, 1996), 이러한 기하학적 특성은 절리의 전반적인 수리역학적 거동에 영향을 미친다(Witherspoon et al., 1980, Brown, 1987, Zou et al., 2015). 예를 들어, 결정질 암석 절리의 간극은 상술한 것처럼 대부분의 유체가 실제로 유동하는 공간에 해당하며, 또한 절리에 하중이 가해질 경우, 초기 변형의 상당 부분을 발생시키는 역할을 한다. 표면 거칠기 수준에 따라 하중에 의한 절리의 수직, 전단 변형이 결정되고, 이는 다시 잔류 간극 수준과 그에 따른 수리 유동에 영향을 준다. 또한 절리 맞물림 혹은 접촉 면적 수준에 따라 수리 유동 경로와 형태가 결정되며, 이는 유동의 채널링(channeling) 및 국부적인 와류(eddy flow) 형성에 영향을 미친다. 이처럼 절리 표면의 기하학적 특성과 수리역학적 거동은 밀접하게 연계되기 때문에 절리를 통한 수리 유동을 해석하는 경우, 시험편의 기하학적 특성을 종합적으로 고려할 필요가 있다.
본 논문은 결정질 암석 절리를 통한 유체 유동의 선형, 비선형 거동을 확인하고, 이에 대한 절리 표면 기하학적 물성의 영향을 파악하기 위해 작성되었다. 이를 위해 먼저 KURT (KAERI Underground Research Tunnel) 인근에서 취득한 결정질 암석을 활용하여 일련의 실내시험을 수행하였다. 무결암 시험편을 인장 파괴하여 인공 절리 시험편을 제작했으며, 그 표면 좌표를 정밀하게 스캔하고 처리하여 거칠기, 간극, 공간 상관성 등의 기하학적 물성을 산정하였다. 이후 일정 압력 주입 형태의 수리시험을 수행하여 압력-유량 관계를 관찰하고 선형, 비선형 유동 여부를 판별함과 동시에 절리 시편의 수리적 물성을 측정하였다. 다수의 절리 시편에 대한 시험을 수행하였으나, 연구 목적에 적합한, 즉 선형, 비선형 유동 차이가 확연한 두 시험편을 선정하여 그 결과를 비교 평가하는 방식으로 논문을 구성하였다. 마지막으로 측정된 간극 공간에 대한 정보를 활용하여 수치해석을 수행하였으며, 이를 통해 기하학적 특성에 의한 두 시험편의 수리 유동 차이를 채널링, 와류 발생 등 유동 형태 관찰을 중심으로 비교하였다. 본 논문에서는 상기 과정을 통해 단일 암석 절리의 수리 유동을 해석할 수 있는 기초 자료를 생산하였다. 개별 절리에 대한 정밀한 유동 평가는 암반 단위의 수리역학적 거동 해석을 위한 입력자료에 해당하기 때문에 본 논문의 결과는 향후 처분장 주변 근계암반 해석 및 안전성 평가 시 참조할 수 있을 것으로 판단된다.
2. 시험 방법
2.1 절리 시험편 준비
본 연구에서 사용한 암석 시험편은 대전시 유성구 지역에서 취득하였으며, 정확히는 KURT 인근에 위치한 1 km 심도 시추공 DB-2에서 확보하였다. 연구 지역의 지질학적 특성을 간략히 기재하면 다음과 같다. 연구 지역은 경기변성암 복합체 내에 위치하며, 주로 선캄브리아기 편마암류와 중생대의 심성암과 관입맥암류로 구성된다. 선캄브리아기의 변성암류는 흑운모편마암 및 편암으로 나누어지며, 이들은 KURT 북서부에 주로 분포한다. 심성암류는 크게 시대 미상의 편상화강암과 연구 지역 전 범위에 걸쳐 광범위하게 분포하는 복운모화강암으로 나뉜다. 이 중, 복운모화강암이 편상화강암을 관입하고 있는 것으로 알려져 있다. 대전도폭에서는 편상화강암을 쥐라기의 편마상 화강암으로 기재하였고, 복운모화강암과는 동일한 마그마에서 유래된 것으로 기재하였다(Park et al., 1977, Lee et al., 1980).
DB-2 시추 코어 박스를 면밀히 확인했으나 본 논문의 실험 조건에 적합한 자연 절리를 충분히 확보할 수 없었다. 특히, 후술할 수리시험 장비의 구조를 고려하면(Fig. 3) 절리가 시추 코어의 중심을 지나고, 경사가 수직에 가까우며, 절삭 및 연마 등 성형이 가능한 수준의 건전성을 보이는 시험편이 요구된다. 이러한 조건을 지닌 시험편을 충분히 찾을 수 없었기 때문에 무결암 시험편을 인장 파괴한 인공 절리 시험편을 대안으로 사용하였다. 무결암 시험편은 DB-2 시추공의 심도 395~430 m 구간에서 취득하였다. 해당 구간의 암종은 대부분 중/조립질 화강암으로 분류되어 있으며, 국부적으로 세립질 화강암 혹은 화강섬록암이 분포하고 있다. 후술할 절리면 수리시험 외에도 해당 암석 시험편의 기본 물성을 측정하기 위한 실내시험을 별도 수행하였으며, 그 결과를 간략히 정리하면 Table 1과 같다.
Table 1.
Fundamental physical and mechanical properties of DB-2 intact specimens
| S.G. (-) | 2.61 | E (GPa) | 50.27 | |
| Porosity (%) | 1.01 | v (-) | 0.26 | |
| UCS (MPa) | 191.26 | BTS (MPa) | 8.96 |
DB-2에서 확보된 화강암 시험편은 단축압축강도 기준 경암(high strength, 100~250 MPa)으로 분류할 수 있으며(ISRM, 1978b), 비중이 높고 공극률이 작아 조직이 치밀한 결정질 무결암으로 확인되었다. 따라서 절리 시험편에 대한 수리시험 수행 시, 대부분의 유체는 매질이 아닌 절리를 통해 유동할 것으로 판단하였다.
실내시험을 위한 모든 시험편의 크기 및 형상은 단축압축시험용 시험편을 기준으로 제작하였다. 이후, 인공 절리 시험편을 생성하기 위해 간접인장시험과 같은 형태의 파괴 시험을 수행하였다. 가압 프레임은 길이 약 100 mm인 NX 크기 원주형 시험편의 크기 및 형상에 적합하도록 제작했으며, 가압부 만곡은 시험편 직경의 약 1.5배 수준을 갖도록 제작하여 일반적인 간접인장시험의 접촉 조건과 동일하게 가공하였다(ISRM, 1978a). 시험편 파괴를 위한 하중은 MTS 816 rock testing system을 사용하여 제어하였으며, 접촉부의 과도한 파쇄를 막기 위해 일정 변위 조건으로 가압하였다. 인공 절리 생성을 위한 전체적인 시험 전경은 Fig. 1과 같다. 총 20개 시험편에 대해 시험을 수행하였으나, 본 논문의 연구 목적을 위해 선정된 두 시험편(HMN2와 HMN9)의 결과에 대해서만 기술하도록 하겠다.
2.2 절리 표면 좌표 측정 및 후처리 작업
절리면의 기하학적 특성을 정량적으로 산정하기 위해서는 표면의 3차원 좌표가 요구된다. 이를 측정하기 위한 다양한 방법들이 제시되었으며, 이는 크게 접촉식과 비접촉식으로 분류할 수 있다. 측정 기술의 발전으로 인해 최근에는 비접촉식 측정이 보편화되었으며, 본 연구에서도 비접촉 방식의 일종인 레이저 센서(laser sensor)를 기반으로 한 3D 스캐너를 사용하였다(Fig. 2(a)). Fig. 2의 스캐너는 Gocator 2530 레이저 라인 센서를 기반으로 제작되었다. 이 센서는 라인 형태의 레이저를 시험편 표면에 조사하고 반사된 레이저를 측정하여 Z축 데이터를 취득한다. 총 3개의 센서를 동시에 적용하면 최대 300 mm 너비의 표면 좌표를 한 번에 측정할 수 있다. 절리 시험편은 좌우로 이동 가능한 측정 테이블 위에 거치되어 이동 정밀도 0.005 mm 수준으로 이동한다. 센서의 해상도는 28~54 μm, Z축 반복성(repeatability)은 0.5 μm 수준이다. 인장 절리면의 상/하부 좌표를 모두 측정하였으며, 이때 X, Y 방향의 측정 점 간격은 0.2 mm로 설정하였다.
한편, Fig. 2(a)를 통해 측정한 절리 표면 자료를 활용하기 위해서는 적절한 후처리 작업이 요구된다. 본 연구에서 적용한 후처리 작업 과정을 간략히 정리하면 다음과 같다. 먼저, 절리면을 V-block 위에 최대한 수평하게 거치한 후 표면을 스캔한다. 이때, 실제 스캔 영역(FOV)은 절리면 넓이보다 약간 크게 설정하여 표면 정보 손실을 최소화하였다(Fig. 2(a)). 따라서 경계부의 불필요한 부분과 측정 시 일부 발생한 노이즈를 제거할 필요가 있으며, 이는 DBSCAN (density based spatial clustering of applications with noise)를 통해 수행하였다. DBSCAN은 데이터 클러스터링 및 이상치(outlier) 제거를 위해 활용되는 알고리즘으로, 미리 설정한 데이터 간 거리 및 최소 데이터 수를 바탕으로 클러스터링을 수행한다(Ester et al., 1996). 이를 Fig. 2(a)와 같이 측정한 자료에 적용할 경우, 소수의 노이즈와 경계부의 불필요 부분은 특정 클러스터에 포함되지 않거나 혹은 상대적으로 작은 클러스터를 생성한다. 측정 영역 설정에 따라 실제 필요한 관심 절리면 부분이 최대 크기의 클러스터를 생성하기 때문에, 이를 제외한 나머지 부분을 모두 제거하였다(Fig. 2(b)). 이후, 간접인장시험 시 접촉부에 생겼을 가능성이 있는 2차 균열의 영향, 즉 경계 효과를 제거하기 위해 네 변의 외측 길이 2.5% 수준을 강제로 제거하였다. 또한 잠재적인 추세의 영향을 제거하기 위해 회귀 평면(regression plane)을 생성하고 이를 제거하는 detrending 작업을 수행하였다(Fig. 2(c)).
역학적 간극(mechanical aperture)은 절리 상/하부면 상 대응되는 두 지점 사이의 물리적 거리를 의미한다. 이를 계산하기 위해서는 상/하부 한 쌍의 절리면 좌표를 동일한 좌표계 내에서 수직 방향으로 정렬할 필요가 있다. 본 연구에서는 점 군(point cloud) 간 좌표를 정렬하기 위해 ICP (iterative closest point) 알고리즘을 적용하였다(Fig. 2(d)). ICP는 최근접점 대응 탐색과 변환 행렬 계산을 반복하여 점 군 간 좌표를 정렬하는 알고리즘이다(Zhang, 1994). 정렬이 완료되면 상/하부면을 공통 좌표계의 XY 평면에 투영하고, 중첩되지 않는 부분이 존재하면 이를 삭제하였다. 마지막으로 일정한 점 간격을 갖도록 내삽을 수행하여 (x, y, z) 좌표를 정리하였으며, 이때 적용한 점 간격은 측정 시의 간격과 동일한 0.2 mm로 설정하였다.
2.3 수리시험 수행
표면 좌표 측정이 완료된 절리 시험편은 수침시킨 후 진공 데시케이터에 24시간 이상 보관하여 포화시켰다. 포화된 시험편의 상/하부에 수리 조건을 부여할 수 있는 엔드캡(endcap)을 거치하고, 그 외부를 열수축 튜브와 방수 테이프를 사용하여 밀봉하였다. 이후 이 조합체를 삼축압축챔버 내 로드셀 위에 거치하고 외부의 MTS 공극수압 가압장치(pore pressure intensifier)와 연결되도록 SUS 튜브 및 피팅으로 체결하였다. 상술한 챔버 내부 세팅(in-vessel setup)은 Fig. 3의 좌측 상부와 같다.
이후 삼축압축챔버는 MTS 816 main frame 위에 거치되고 MTS 구속압 제어 유닛(confining pressure unit)에 연결되어 시험편에 구속압을 가하였다. 구속압 제어 유닛은 작동유의 압력을 서보 제어하여 시험 시간 동안 일정한 압력을 유지시키며, 이때 가한 구속압은 절리면에 대한 수직응력으로 작용한다. 시험편 축 방향 압력은 MTS main frame을 통해 가해지며 마찬가지로 서보 제어하였다. 공극수압 가압장치는 시험편을 통과하는 유체의 압력을 서보 제어하며, 본 논문에서 사용한 유체는 상온의 증류수이다. 시험편 하부 엔드캡을 통해 주입된 증류수는 절리 간극과 상부 엔드캡을 거쳐 챔버 외부로 배출되며, 배출된 증류수는 전자 저울에 집수되어 시간에 따른 유출량을 측정하였다.
절리 시험편에 가해진 수직응력(구속압)과 수압은 실제 처분장의 암반 환경을 반영하여 설정하였다. 먼저, 처분장의 심도는 500 m로 가정하였고, 암반의 비중은 국내 화강암의 일반적인 비중인 2.5를 가정하였다. 이 경우, 암반 상재 하중으로 인한 수직응력은 12.5 MPa로 산정된다(). 또한, 국내 원위치 응력 측정 결과를 참조하여(Kim et al., 2021) 500 m 심도에서 측압비()를 1.5 정도로 가정하면 약 20 MPa 수준의 원위치 응력이 예상된다. 이를 반영하여 수리시험 시 가한 수직응력의 최대값은 20 MPa 수준에서 결정되었다. 한편, 지표에서 심도 500 m의 처분장까지 모든 위치 수두가 작용할 경우 예상되는 수압은 약 5 MPa이다. 따라서 수리시험 시 적용한 수압 역시 이를 반영하여 최대 5 MPa까지 적용하였고, 모든 경우에서 일정 압력 조건에 따라 유체를 주입하였다. 상술한 수직응력 - 수압 범위에 따라 설정한 시험 조합은 총 22가지이며(수압이 너무 낮거나 혹은 수압이 수직응력 보다 높아서 유체 유출 가능성이 있는 조합은 제외), 각 단계의 구체적 시험 조건은 Table 2와 같다.
Table 2.
Test cases defined by the combination of normal stress and water pressure (22 cases)
| Normal stress (MPa) | |||||
| 3.0 | 5.0 | 10.0 | 15.0 | 20.0 | |
|
Water pressure (MPa) | 0.5 | 0.5 | - | - | - |
| 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | |
| 2.0 | 2.0 | 2.0 | 2.0 | 2.0 | |
| - | 3.0 | 3.0 | 3.0 | 3.0 | |
| - | - | 4.0 | 4.0 | 4.0 | |
| - | - | 5.0 | 5.0 | 5.0 | |
전체적인 시험 과정은 다음과 같다. 먼저, 구속압(수직응력)과 축 방향 응력을 동시에 증가시켜 정수압 상태를 형성하고, 이를 시험이 진행되는 동안 유지시켰다. 수직응력이 일정하게 유지되는 것을 확인한 후, 유체 주입을 개시하였다. 주입은 각 단계에서 설정한 수압이 일정하게 유지되며, 동시에 전자 저울에서 측정되는 배출 유량이 일정해질 때까지 수행하였다. 이때, 측정된 배출수의 무게는 시간에 따라 측정하였다. 한 구속압 수준에서 모든 수압 조건을 수행한 후, 다음 구속압 수준으로 이행하여 동일 절차를 반복하였다. 이 과정을 Table 2의 전체 조합에 대해 순차적으로 적용하여 시험을 수행하였다.
3. 시험 결과 및 논의
3.1 절리 표면의 기하학적 특성
3.1.1 절리 거칠기
본 연구에서는 절리면의 여러 가지 기하학적 특성 중, 거칠기, 초기 역학적 간극, 간극의 공간 상관성을 측정하였으며, 그 결과는 3.1 절을 통해 순차적으로 기술하였다. 먼저, 거칠기는 절리면의 중요 기하학적 물성 중 하나이며, 이를 평가하는 가장 대표적인 지수는 JRC (joint roughness coefficient)를 들 수 있다(Barton and Choubey, 1977). JRC는 그 활용성과 대표성에도 불구하고 몇 가지 단점이 지적되고 있는데 예를 들면, JRC는 정량성이 다소 부족하며 육안으로 거칠기 평가 시, 사용자의 주관과 경험이 반영될 수 있다. 이러한 단점을 보완하기 위해 다양한 정량적 거칠기 지수들이 제안되었다. 각각의 지수는 고유의 정의를 지니고 그에 따른 장/단점이 존재하지만, 대부분은 표면 형상의 한 가지 특성만을 반영하는 경우가 많다. 절리면의 복잡한 기하학적 성질을 한 가지 특성만으로 평가하는 것은 한계가 있기 때문에, 서로 다른 정의에서 유도된 복수의 거칠기 지수를 사용하는 것이 거칠기의 종합적인 평가를 위해 바람직하다. 따라서 본 연구에서는 세 가지 서로 다른 정의를 지닌 거칠기 지수를 적용하여 두 절리 시험편의 거칠기를 평가하였다.
본 논문에서 적용한 첫 번째 거칠기 지수는 프랙탈 차원(fractal dimension)이다. Mandelbrot(1967) 이후, 암석 절리 형상을 프랙탈 개념에서 해석하는 연구는 꾸준히 진행되어 왔다. 절리 형상이 self-affine 특성을 지닌다는 가정 하에, 절리면의 프랙탈 물성을 측정하거나(Brown and Scholz, 1985b, Kulatilake and Um, 1999, Jacob et al., 2017), 혹은 특정 물성을 따르는 프로파일/표면을 수치적으로 합성하는 연구는 매우 다양하게 수행되었다(Patir, 1978, Candela et al., 2009, Stigsson and Ivars, 2019). 이 중 본 연구에서는 절리 프로파일에 대한 FFT (fast Fourier transform) 기반의 PSD (power spectral density)를 이용하여 프랙탈 지수를 산정하였다. 즉, 신호(절리 프로파일)를 주파수 영역으로 변환한 후 직선 회귀하였으며, 회귀 직선의 기울기가 계산되면 식 (1)과 같이 Hurst exponent()와 프랙탈 차원 를 산정할 수 있다().
여기서 는 유클리디안 차원, 𝛽는 PSD에서 회귀한 직선의 기울기를 의미한다.
두 번째, 세 번째 거칠기 지수는 (RMS of the first derivative of the profile)와 (center line average)를 선정하였다(Tse and Cruden, 1979). 는 절리면의 미세 거칠기(micro-asperity)를 바탕으로 산정되고, 는 절리면의 전반적인 진폭(amplitude)을 강조하여 계산되기 때문에 서로 다른 관점에서 절리면의 기하학적 특성을 평가할 수 있을 것으로 판단하였다. 두 지수의 정의는 각각 식 (2), (3)과 같다.
한편, 상술한 세 가지 거칠기 지수는 모두 2차원 지수이다. 절리면은 기본적으로 3차원 객체이기 때문에 2차원 지수를 적용하면 거칠기 방향성 등의 특성화가 제한된다. 이를 보완하기 위해 3차원 거칠기 지수들이 제안되었으나(Belem et al., 2000, Tatone and Grasselli, 2009, Park and Song, 2013), 아직까지 2차원 지수들의 활용성을 넘어서지 못하는 실정이다. 본 연구에서는 2차원 거칠기 지수를 사용하되, Fig. 4와 같은 방법을 적용하여 거칠기의 이방성을 특성화하였다.
먼저 Fig. 4(a)와 같이, 절리면 장축이자 수리 유동 방향을 y축으로 설정하고, x축(𝜃=0°)을 기준으로 반시계 방향으로 15° 만큼 회전(0°~165°까지 15° 간격)하며 여러 단면을 생성하였다(Fig. 4(b)). 이 단면들과 교차하는 복수의 2차원 프로파일 그룹을 생성하고(Fig. 4(c)), 이 프로파일 그룹에 상술한 세 가지 거칠기 지수를 각각 적용하였다. 이를 통해 Fig. 4(d)와 같이 지수별 분포를 얻은 후, 그 평균을 해당 방향의 거칠기 대표값으로 산정하였다. 이 방식을 통해 계산된 HMN2, HMN9 시험편의 거칠기는 Table 3과 같다. 이때, 거칠기 최대/최소값에 대응하는 방향과 이방성 지수(최대/최소 비율)를 함께 기재하였다.
Table 3.
Evaluation of roughness parameters and its directionality
두 시험편의 거칠기 평가 결과를 비교하면 다음과 같다. 먼저 의 경우, 두 절리면 모두 유사한 수준의 거칠기가 측정되었다. 최대, 최소값이 각각 약 1.5, 1.3 수준으로 평가되었고, 이에 따른 이방성 수준 역시 1.1 정도의 낮은 값을 보였다. 의 최대값 방향은 대체적으로 절리면 장축 방향이 우세한 것으로 확인되었고, 최소 방향은 특정한 경향을 보이지 않았다. 절리면의 미세 거칠기로 계산되는 역시 두 시험편에서 유사하게 평가되었다. 최대값은 약 0.26, 최소값은 0.19로 측정되었고, 이방성 비는 1.4로 에 비해 높은 거칠기 이방성이 평가되었다. 의 최대 방향은 시험편 장축 방향으로 산정되었는데, 이는 Fig. 1에서 설정한 것과 같이 인장 절리 생성 시 가한 가압축에 수직한 방향이며, 반면에 최소 방향은 가압 방향으로 측정되었다. 이는 간접인장시험 시 미세 거칠기의 형성은 가압 방향의 영향을 받을 수 있음을 의미한다. 한편, 표면의 전반적인 진폭에 따라 결정되는 는 세 거칠기 지수 중 가장 대조적인 결과를 보였다. HMN9가 HMN2에 비해 최대, 최소 모두에서 더 큰 진폭 수준을 보였으며, 이방성 비 역시 1.2로 HMN2 보다 크게 측정되었다. 특히, HMN9의 최대값은 135° 방향에서 측정되어 전체적으로 대각선 방향의 1차 거칠기(undulation)가 발달하였다. 종합하면, HMN2과 HMN9 두 시험편의 거칠기는 프랙탈 차원과 미세 거칠기는 유사한 수준이며 전반적인 진폭은 HMN9 시험편이 다소 높게 측정되었다. 즉, 국부적인 기울기와 스펙트럼 성분은 유사하지만, 저주파 성분은 HMN9가 더 뚜렷한 것으로 확인되었다.
3.1.2 초기 역학적 간극
역학적 간극은 절리면 상/하부 대응되는 두 지점 사이의 물리적 거리를 의미한다. 본 연구에서는 ICP 기법을 통해 동일한 좌표계 내에서 두 절리면을 정렬했기 때문에 상부 절리 좌표에서 하부를 빼주는 방식으로 간극을 계산하였다. 한편, 초기(initial) 역학적 간극()은 하중이 없거나 작은 초기 상태에서의 간극을 의미한다. 이를 정의하는 방식은 다양할 수 있으나, 본 연구에서는 절리 시험편 자중에 의해 접촉된 상태를 초기 상태로 가정하였다. Esaki et al.(1998)은 자중에 의한 절리면 접촉 시, 전체 면적의 약 1% 수준이 접촉한다고 보고하였다. 본 연구에서는 이를 반영하여 상부 절리면을 하부로 평행이동 시키며 중첩 상태를 확인하고, 중첩된 표면 노드(node)가 전체 노드의 1%에 도달하면 초기 상태인 것으로 간주하였다. 이러한 방식에 따라 두 절리 시험편의 초기 역학적 간극을 산정했으며, 그 결과에 대한 여러 기술통계량을 기재하면 Table 4와 같다. 또한 초기 역학적 간극 분포의 히스토그램을 도시하면 Fig. 5과 같다.
Table 4.
Descriptive statistics of initial mechanical aperture
| No | Ave. | STD | Median | IQR | Skewness | Kurtosis |
| HMN2 | 0.4185 | 0.1890 | 0.4082 | 0.2327 | 0.6979 | 2.1447 |
| HMN9 | 0.4700 | 0.1913 | 0.4684 | 0.2398 | 0.2239 | 0.5609 |
Table 4와 Fig. 5에서 확인할 수 있듯, 초기 역학적 간극의 평균은 HMN9가 더 크게 측정되었다. 그러나 상대적인 분산은 HMN2가 크기 때문에 더 넓은 분포를 보이고 있다. 두 시험편 모두 간극의 왜도가 양수로 측정되어 오른쪽 꼬리가 생성되었으며, HMN2의 왜도가 더 크기 때문에 중심이 왼쪽으로 치우친 분포를 보였다. 또한 간극 분포는 두 시험편 모두 단봉(unimodal) 형태를 보였으나 KS (Kolmogorov – Smirnov) 검정을 수행한 결과, 유의수준 0.05를 기준으로 특정 분포(정규, 대수정규, Weibull, Gamma)를 따르지 않는 것으로 판정되었다.
3.1.3 간극에 대한 공간 상관성
공간 상관성은 분리 거리(lag distance)에 따른 물성의 유사성을 나타내며 베리오그램(variogram), 매도그램(madogram), 자기 상관함수(autocorrelation) 등의 지구통계학적 기법을 통해 정량화할 수 있다. 본 연구에서는 이중 가장 대표적인 베리오그램을 적용했으며, 수리 유동의 통로에 해당하는 간극의 공간 상관성을 분석하여 수리시험 결과와의 연관성을 확인하고자 하였다. 활용 가능한 간극 정보는 초기 상태를 가정하고 산정하였기 때문에, 본 논문에서의 공간 상관성은 모두 초기 역학적 간극에 대한 분석 결과이다. 베리오그램은 분리 거리에 대한 자료들의 산포도를 의미하며 식 (4)와 같이 정의된다.
여기서 는 특정 위치, 는 분리 거리, 는 해당 위치에서의 관심 물성을 의미하며, 본 연구에서는 초기 역학적 간극에 해당한다. 일반적인 베리오그램 분석 절차는 다음과 같다. 먼저, 주어진 자료에 대해 식 (4)를 적용하여 실험적 베리오그램을 구성한다. 다음 이를 일반화하기 위해 이론적 베리오그램 수식을 적용하여 모델링한다. 이론적 모델은 실험적 베리오그램 형태에 따라 결정되며 크게 문턱값(sill)이 있는 모델, 없는 모델, 주기 모델 등으로 분류할 수 있다. 이후 보다 정교한 모델링을 위해 상관거리(range), 문턱값(sill), 너깃(nugget) 등의 베리오그램 인자들을 세부적으로 보정한다. 본 연구에서는 두 시험편 간극의 실험적 베리오그램 형태를 확인한 후, 구형 모델(spherical model)의 적합성이 높은 것으로 판단하여 이를 기본 모델로 선정하였다(식 (5)).
여기서 는 너깃, 는 부분 문턱값(partial sill), 는 상관거리를 의미한다. 너깃은 자료 자체에 내재된 불확실성을 의미하며 이론적으로 베리오그램의 y 절편에 해당한다. ()는 자료의 분산 문턱값이며 특정 거리 이상에서 일정하게 유지되거나 혹은 공간 상관성이 불분명한 경우 발산할 수 있다. 상관거리는 공간 상관성이 유지되는 거리를 의미하며, 이 이상의 거리에서는 자료들 사이의 공간적 독립을 간주할 수 있다.
경우에 따라 단일 모델로는 이론적 베리오그램 모델링에 한계가 있을 수 있으며, 이때 여러 개의 기본 모델을 선형 결합하는 모델(nested model)을 사용할 수 있다. 본 논문에서는 두 개의 구형 모델을 결합하여 모델링을 수행하였다. 두 기본 모델 중, 상관거리가 짧은 성분의 베리오그램 인자는 작은 범위(scale)에서의 공간적 특성을 나타내며, 상관거리가 긴 성분은 보다 넓은 범위에서의 공간적 특성을 나타낸다. 본 연구에서는 범위가 다른 두 모델 각각의 상관거리와 문턱값을 대표 인자로 선정하여 계산하였다. 한편, 3차원 자료에 대한 베리오그램은 특정 방향에 따라 계산하거나, 혹은 전방향(omni-direction) 형태로 계산하여 방향 구분 없이 자료 전체의 상관성을 확인할 수 있다. 본 연구의 수리시험은 절리의 장축 방향을 따라 유체가 유동하도록 설계되었기 때문에 전방향과 𝜃=90°, 두 경우에 대해 모델링을 수행하였다. 또한 이론적 베리오그램 피팅(fitting) 시, 첫 번째 정점(peak) 중심으로 모델링을 수행하였다. 방향과 범위에 따라 계산한 베리오그램 대표 인자는 Table 5와 같고 이를 두 시험편 초기 역학적 간극의 공간적 분포와 함께 도시하면 Fig. 6과 같다.
Table 5.
Variogram parameters according to the direction and scale
여기서 는 상관거리, 아래 첨자 , 90, , 는 각각 omni-direction, 𝜃=90°방향, 작은 범위(short scale)과 넓은 범위(large scale)을 의미한다.
Fig. 6(a)와 같이 두 시험편의 공간적 간극 분포는 차이를 보였다. 여기서 파란색이 간극이 가장 큰 지역을 의미하며 초록색이 가장 작은 지역을 의미한다. HMN2의 경우, 유동 방향(장축 방향)으로 간극의 연결성이 발달하였으며, 간극이 작은 잠재적 접촉 부분은 절리 상부에 집중적으로 분포하는 것을 확인할 수 있다. 반면 HMN9의 경우, 간극은 유동 방향이 아닌 좌측 상부 방향으로 우세하게 연결되었으며, 잠재적 접촉부 역시 상대적으로 넓게 분포하였다. 이러한 분포는 Table 3의 거칠기 방향성과 정성적으로 일치한다. HMN2의 거칠기는 지수에 관계없이 𝜃=90°방향으로 최대값을 보였는데, 이는 간극의 주요 연결 방향과 일치한다. 또한 HMN9는 90°~135° 방향으로 거칠기 최대값을 보였는데, 이 역시 간극의 발달 방향과 대체적으로 유사하다. 간극의 크기가 거칠기만으로 결정되지 않지만 Fig. 6(a)와 같이 어느 정도의 상관성이 존재함이 관찰되었다. 두 시험편의 간극 분포 특성은 Fig. 6(b)의 베리오그램 인자를 통해 정량적으로 확인할 수 있다. HMN2 시험편의 간극 상관거리는 방향과 규모에 관계없이 HMN9에 비해 크게 측정되었다. 특히, 유동 방향을 따른 에서 큰 차이를 보였는데, HMN2의 상관거리는 HMN9의 약 1.6배로 계산되어 Fig. 6(a)에서 육안으로 관찰된 분포의 연속성이 정량적으로 확인되었다.
절리면의 기하학적 성질에 대한 분석 결과를 종합하면 다음과 같다. 먼저, 프랙탈과 미세 거칠기 관점에서 두 시험편의 거칠기는 유사한 수준이었으나 전체적인 진폭에 의한 거칠기는 HMN9가 더 큰 값을 보였다. 이는 적용하는 지수의 특성에 따라 서로 다른 거칠기 산정 결과를 보일 수 있음을 의미한다. 초기 역학적 간극의 경우, HMN9 시험편의 평균 간극이 HMN2에 비해 크게 측정되었다. 그러나 단순 평균이 아닌 공간 상관성을 고려하면 HMN2의 상관거리가 상대적으로 크게 측정되었으며, 동시에 HMN2 시험편은 후술한 수리시험 시 유체의 유동 방향으로 더 우세한 간극 연결성을 보였다.
3.2 수리시험 결과
Fig. 3의 시험 시스템과 Table 2의 시험 조건을 적용하여 일정 압력 조건에서 절리 시험편에 대한 수리시험을 수행하였다. 수직응력, 수압 수준에 따른 배출량(discharged flow rate)의 변화는 Fig. 7과 같다. 여기서 닫힌 도형이 HMN2, 열린 도형이 HMN9의 결과에 해당한다. Fig. 7(a)와 같이, 주입 압력이 일정할 때, 수직응력 증가에 따라 배출량은 지수적으로 감소하는 것을 확인할 수 있다. 이는 응력 증가에 의해 간극이 감소하기 때문으로 암석 절리를 통한 수리시험에서 일반적으로 확인되는 경향이다. 또한 Fig. 7(b)와 같이 일정 수직응력 조건에서 수압이 증가할수록 유량이 증가하는 경향 역시 일반적으로 확인되는 현상이다(Zhao and Brown, 1992, Singh et al., 2015, Kulatilake et al., 2020, Ye et al., 2022).
서론에 기술한 것처럼, 수직응력, 수압 구배, 기하학적 조건 등에 따라 압력 구배()와 배출 유량 사이에 특정한 관계가 존재할 수 있다. 만약 선형적인 관계가 확인되면 유동은 층류 영역에서 발생한 것으로 볼 수 있으며, 이 경우 압력과 유량 사이 관계는 식 (6)과 같이 Darcy law 혹은 삼승 법칙(cubic law)으로 표현할 수 있다.
여기서 는 배출 유량, 𝜇는 유체의 동점성 계수(dynamic viscosity), 는 유동 단면적, 는 절리 폭, 는 수리적 간극(hydraulic aperture), 는 수압 구배, 는 절리의 투수계수(permeability)로 와 같다. 삼승 법칙은 층류 조건과 평판 모델을 가정했을 때 성립한다. 그러나 실제 암석 절리는 이상적인 평판 모델과 큰 차이가 있으며(Witherspoon et al., 1980, Zimmerman and Bodvarsson, 1996), 층류 조건에서 삼승 법칙이 성립함을 가정했을 때, 이에 대한 등가 간극에 해당하는 것이 수리적 간극이다.
한편, 절리면의 기하학적 형상이나 압력 조건에 따라 주입 압력과 배출 유량 사이에 비선형적 관계가 존재할 수 있다. 이러한 비선형, non-Darcy 유동은 Forchheimer 방정식 혹은 Izbash 방정식을 통해 표현할 수 있다. 본 연구에서는 Forchheimer 방정식을 적용했으며 이는 식 (7)과 같다.
여기서, 는 Forchheimer 상수로 각각 , 와 같으며, 𝜌는 유체의 밀도, 𝛽는 표면 형상과 관련된 Forchheimer coefficient에 해당한다. Fig. 7(b)를 수압 구배와 유량의 관계로 다시 도시하면 Fig. 8(a)과 같다. 여기서 닫힌 도형이 HMN2, 열린 도형이 HMN9 결과에 해당하며 실선 및 점선으로 표현된 곡선은 Forchheimer 식에 대한 회귀 곡선을 의미한다.
Fig. 8(a)에서 확인할 수 있듯, Forchheimer 방정식은 본 연구에서 수행된 수리시험 결과와 그 비선형성을 양호하게 모사하는 것으로 확인되었다. 특히, HMN9의 비선형성이 더욱 확연한 것으로 보이며, 동시에 같은 압력 구배에서 HMN9 시험편의 배출 유량이 더 크게 측정되었다. Table 4에서 확인할 수 있듯, 두 시험편 중 HMN9의 초기 역학적 간극이 상대적으로 크게 측정되었다. 본 연구에서는 특정 수직응력 하에서 간극의 변화를 측정하지 않았지만, 초기 역학적 간극의 차이를 바탕으로 판단할 때 유사한 수준의 수직응력 하에서는 HMN9의 (잔류)간극이 더 클 가능성이 높다. 따라서 동일 수압 조건에서 HMN9의 배출 유량이 상대적으로 큰 것으로 판단되었다. 다만, 비선형성 증가에는 간극의 크기뿐 아니라 공간 상관성 차이에 따른 현상도 기여할 수 있으므로, 해석 시 이를 함께 고려할 필요가 있다.
한편, 전체적인 선형, 비선형 유동 차이를 확인하기 위해 Fig. 8(b)를 추가 도시하였다. 명확한 대비를 위해 HMN9의 결과를 사용했으며, 여기서 실선은 식 (7)에 따른 Forchheimer 곡선이고 점선은 전체 유동 중 선형 성분, 즉 식 (7)에서 만 도시한 결과이다. Fig. 8(b)에서 확인할 수 있듯 수압 구배가 증가함에 따라 직선에서의 편차가 증가하여 비선형성이 강화되는 것을 볼 수 있다.
Fig. 8의 시험 결과에 식 (7)을 적용하여 계산한 각종 수리 상수(식 (7)의 , , , 𝛽 등)는 Table 6과 같다. 시험 시 유체는 증류수를 사용하였기 때문에 상온 기준으로 밀도는 1000 kg/m3, 점성계수 1.002×10-3 Pa·s를 적용하였다.
Table 6.
Hydraulic properties calculated from Forchheimer equation
Table 6의 수리 상수들은 수직응력에 따라 변하는 것을 확인할 수 있으며, 이를 수직응력에 따라 도시하면 Fig. 9와 같다. 먼저, 수직응력 증가에 따라 식 (7)의 두 Forchheimer 계수 , 는 모두 증가하는 경향을 보였다. 3~20 MPa의 수직응력 범위에서 HMN2, HMN9 시험편의 계수 는 각각 16.1, 18.8배 증가하였으며, 동일 범위에서 는 각각 99.9, 39.4배 증가하였다. 이는 암석 절리를 통한 유체의 비선형 유동에서 일반적으로 확인되는 경향이다(Zoorabadi, 2015, Rong et al., 2020, Yu et al., 2022). 식 (7)의 Forchheimer 방정식은 두 항의 결합으로 구성되는데, 이중 첫 번째 항(즉, 계수가 인 항)은 선형 거동에서 유동에 대한 저항을 나타내며, 두 번째 항(계수가 인 항)은 마찬가지로 비선형 거동에서의 추가 저항을 의미한다. 수직응력이 증가함에 따라 절리의 간극은 수축하고, 접촉 면적, 채널링 등의 효과는 증가하는 것이 일반적이다. 이 경우, 유동에 대한 저항이 증가하는 것으로 해석할 수 있기 때문에 두 항의 계수는 모두 수직응력에 따라 증가한다. 한편, 수직응력 증가에 따라 배출 유량이 감소하기 때문에 선형 유동 항에서 계산되는 수리적 간극()은 감소했으며, 또 다른 계수 𝛽는 수직응력과 일정한 상관관계를 보이지 않았다.
상술한 것처럼 식 (7)의 두 항 , 는 각각 선형, 비선형 유동과 관련된다. 두 항을 식 (8)과 같이 전체 압력 구배에 대한 비율로 정의하면, 이는 nonlinear effect factor에 해당한다(Zeng and Grigg, 2006). 공학적으로 식 (8)의 E가 10%에 해당하면 선형-비선형 유동 간 전이가 발생한 것으로 판단할 수 있다(Zimmerman et al., 2004, Javadi et al., 2014, Xiong et al., 2018).
한편, Reynolds number()는 유체 유동에서 관성과 점성의 비율로 정의되는 무차원 수이며, 층류와 난류를 구분하는 중요한 기준이 된다(식 (9)). 즉, 가 작을수록 유체 분자 간 마찰에 의한 점성이 높아 층류가 발생하여 반대의 경우 관성력이 우세하여 난류로 전이된다. 이는 유동 해석을 단순화하여 예측할 수 있어 다양한 공학적 문제에서 유동 특성을 예측하는데 활용된다.
식 (8)을 에 대해 정리하여 식 (9)에 대입하면 최종적으로 식 (10)과 같이 정리할 수 있다. 또한 상술한 것처럼 E = 10%에서 선형-비선형 전이가 발생하는 것으로 가정하면 식 (10)과 같은 critical Reynolds number()가 정의된다(Javadi et al., 2014, Xiong et al., 2018).
두 시험편에서 측정된 결과를 바탕으로 수직응력과 수압 조건에 따른 , 를 계산하면 Table 7과 같다. 여기서 굵은 글자는 임계값을 초과하여 비선형 유동이 발생한 경우를 의미한다. Table 7에서 확인할 수 있듯, 동일 수직응력 조건에서 수압이 증가할수록 가 증가하였는데, 이는 수압 구배 증가에 따라 평균 유속이 증가하기 때문이다. 따라서 동일 수직응력에서 주입 수압이 증가할수록 비선형 유동 발생 가능성이 커짐을 알 수 있다. 반대로 동일 수압 조건에서 수직응력이 증가할수록 가 감소하였다. 이는 수직응력 증가에 따른 간극 감소와 이에 따른 전도도 감소로 평균 유속이 저하되기 때문이다. 따라서 가 임계값을 초과하기 어려워져 층류가 유지될 가능성이 상대적으로 커짐을 의미한다.
Table 7.
Comparison of Re and Rec according to the normal and water pressure
한편, 와 이에 따른 비선형 유동 발생 빈도를 확인하면 HMN2와 HMN9 두 시험편의 유동은 큰 차이가 있음을 알 수 있다. HMN2의 경우, 총 22개의 수직응력-수압 조건 중에서 비선형 유동을 보인 것은 한 가지 경우(수직응력 5 MPa, 주입 수압 3 MPa)에 불과하지만, HMN9는 14 경우에서 비선형 유동이 확인되었다. 이러한 유동 차이의 원인은 두 시험편의 기하학적 물성, 그 중에서도 간극과 그 공간 상관성에서 기인한 것으로 판단하였다. Table 5와 Fig. 6을 통하여 두 시험편의 간극 연결성에 차이가 있는 점은 정성적, 정량적으로 확인되었다. HMN2의 간극 분포는 시험편의 장축 방향, 즉 수리 유동 방향을 따라 지배적으로 분포하고 있으며, 이에 비해 HMN9는 단축 방향 혹은 좌상부 방향으로의 연결성이 상대적으로 높게 형성되었다. HMN2의 우세한 간극 연결성(특히, 유동 방향)은 그만큼 수리 유동이 원활하게 발생하여 선형적인 층류 거동이 유지되기 용이함을 의미한다. 반면에 압력 구배에 대각인 방향으로 간극이 발달했으며 상관거리가 짧은 HMN9 경우, 복잡하고 짧은 연결성으로 인해 유체 유동 시 채널링의 영향을 크게 받았을 가능성이 있고, 이는 다시 HMN9 시험편의 비선형적 수리 거동에 영향을 미쳤을 가능성이 있다. 이처럼 암석 절리를 통한 수리 유동은 역학적, 수리적 조건뿐 아니라 절리면의 다양한 기하학적 특성을 함께 고려하여 종합적으로 판단할 필요가 있다.
4. 수치해석
4.1 모델 구성
상술한 실내시험 결과 분석을 통해 절리를 통한 수리 유동은 표면의 기하학적 물성의 영향을 받는다는 점을 확인하였다. 그러나 이는 압력과 유량 측정값에 따른 거시적 해석 결과이기 때문에, 간극 공간을 따른 미세한 유동의 차이를 보다 명확하게 확인하고자 수치해석을 수행하였다. 다만 이는 정량적인 분석이 아닌 정성적인 비교에 목적을 두었으며, 그 이유는 다음과 같다. 시험 결과와 해석 결과를 정량적으로 비교하기 위해서는 하중 조건에 따른 간극 공간의 변화를 반영하여야 한다. 즉, 초기 역학적 간극() 뿐 아니라 하중에 따른 수직 변위(절리 닫힘, )를 측정하여 특정 응력 조건에서의 간극 상태를 계산하여야 한다(). 본 논문에서 사용한 절리 시험편은 단축압축시험편을 직경 방향으로 가압하여 파쇄한 인장 절리이다. 이러한 형태의 절리면을 수직 방향(직경 방향)으로 가압하면 가압 도중 시편 파괴의 위험이 높기 때문에 하중으로 인한 변위와 그에 따른 잔류 간극을 실험적으로 측정할 수 없었다. 이에 대한 대안으로 절리면의 형상을 바탕으로 Hertzian 접촉 이론 등에 따른 절리 압축을 모사하거나(Brown and Scholz, 1985a, Matsuki et al., 2008, Ma et al., 2024) 혹은 수치해석적 접근을 적용할 수 있다. 그러나 이는 본 논문의 연구 범위를 벗어나며, 이에 대한 연구는 추후 수행하여 보다 정확한 모델 구성에 반영할 계획이다. 따라서 본 절의 수치해석은 간극 공간 내의 유동 차이를 가시적으로 대비하고 비선형 유동의 원인을 파악하기 위한 보조적 수단으로 수행되었다. 해석 모델은 주어진 정보(상/하부 절리 좌표)를 바탕으로 구성하되, 그 외의 조건은 HMN2, HMN9 두 시험편의 수리 유동에 대한 정성적 비교가 용이하도록 적절히 가정하여 적용하였다.
수리 해석은 COMSOL multiphysics을 통해 수행하였으며, 다음과 같은 과정에 따라 모델을 구성하였다. 먼저, 2.2절에서 측정한 두 절리면의 상/하부 절리 좌표를 활용하였으며, 각각의 절리가 접촉 면적 0%, 1%, 4%인 상태의 간극 공간을 가정하였다. 원하는 접촉 면적에 도달하도록 상/하부 절리 좌표를 수직 방향으로 평행 이동시킨 후, 이를 만족하는 상태의 절리 좌표를 COMSOL에 입력하여 상/하부 면을 구성하였다. 따라서 수직 하중에 의한 절리 표면 변형 등은 고려하지 않았다. 이후, 간극 공간의 네 측면을 생성하고 접촉부(contact area)를 모두 삭제하여 간극을 모사하는 닫힌 공간을 구성하였다(Fig. 10). 앞선 실내시험과 같이 유동은 y축 방향을 따라 발생하는 것으로 설정하였으며, 이에 따른 inlet과 outlet 경계를 지정하였다. Inlet에서는 일정 속도의 주입 조건을, outlet은 대기압 상태의 압력 조건을, 그 외의 모든 면은 ‘no slip’ 경계 조건을 설정하였다.
서론에 기재한 것과 같이 절리를 통한 비압축성, 뉴턴 유체의 유동을 해석하기 위해 일반적으로 사용되는 지배 방정식은 Navier-Stokes 방정식과 질량 보존 방정식이며, 이는 각각 식 (11), (12)와 같다(Zimmerman and Bodvarsson, 1996).
여기서 𝜌는 유체 밀도, 는 속도 벡터, 는 압력, 𝜇는 점성계수, 는 외력을 의미한다. 식 (11)은 운동량 보존을 의미하며 식 (12)는 비압축성 조건을 나타내는 연속 방정식(continuity equation)에 해당한다. Navier-Stokes 방정식에서 관성력이 점성력에 비해 무시할 수 있는 수준으로 작은 경우, 식 (11)의 비선형 항을 무시하여 Stokes 방정식으로 단순화할 수 있으며, 다시 평판 조건을 가정하면 식 (6)과 같은 삼승 법칙이 유도된다. 본 논문에서는 Navier-Stokes 방정식을 적용하여 절리 간극 공간을 통한 유동 해석을 수행하였고, 해석은 정상상태, 비압축성, 뉴턴 유체 가정하에 수행하였다. 이때 유체의 물성은 𝜌 = 1,000 kg/m3, 𝜇 = 1.002 × 10-3 Pa·s를 적용하였다.
4.2 해석 결과
두 절리면 정보에 0%, 1%, 4%인 접촉 면적을 구현하여 간극 공간을 각각 구성하고 해석을 수행하였다. 이는 실내시험 결과 중, 수직응력 증가에 따른 유동 변화와 비교하기 위해 계획되었다. 접촉 면적의 수준에 따라 HMN2 모델의 평균 역학적 간극은 약 0.92, 0.42, 0.29 mm 수준이었고, HMN9는 약 0.97, 0.47, 0.32 mm이 적용되었다. Inlet을 통한 주입 속도는 = 100에 해당하는 유속을 입력하였으며, 이는 Table 7의 critical Reynolds number를 참고하여 모든 경우에서 비선형적 유동이 충분히 발생할 수 있도록 설정하였다. 접촉 면적(수직응력)에 따른 두 절리 시험편의 해석 결과는 Fig. 11과 같으며, 이는 유속에 대한 contour와 유선(streamline)으로 도시하였다.
Fig. 11(a), (d)에서 확인할 수 있듯, 접촉 면적을 고려하지 않는 상태에서도 두 절리면을 통한 유체의 유동은 약간의 차이를 보였다. 두 경우 모두에서 유선은 압력 구배 방향으로 비교적 직선적으로 나타났으나 HMN9 시험편의 경우, 절리 중앙부에 위치한 영역에서 유속이 비교적 작게 측정되었다. Fig. 11(d) 결과를 Fig. 6(a)의 HMN9 초기 역학적 간극 분포와 비교하면 해당 영역(Fig. 11(d)에서 흰색 영역)은 간극이 작은 지역과 비교적 일치하며, 이를 통해 본 해석 모델의 간극 공간이 실제 절리 정보를 비교적 잘 모사하였음을 확인할 수 있다. 한편, 수직응력 증가와 그에 따른 접촉 면적 증가로 인해 유동 경로가 휘어지고 유로가 제한되는 채널링 현상이 확인되었다. 이는 접촉 면적이 증가함에 따라 더욱 확연히 구분되었는데, Fig. 11(c), (f)와 같이 접촉 면적이 4%에 달한 경우, 두 절리면 모두에서 뚜렷한 채널링 현상이 확인되었다. 이 과정에서 간극이 작은 지역의 국부 가속 구간과 저속 구간이 모두 확인되며, 접촉부 주변의 유선 굴곡과 분기가 관찰되었다. 이러한 현상은 관성항의 영향을 증가시켜 전체 유동의 비선형성 강화로 연결될 수 있다.
두 시험편의 해석 결과를 상대적으로 비교하면 다음과 같다. HMN2의 경우, 주요 유선이 y축을 따라 비교적 잘 유지된 반면, HMN9은 접촉 면적이 증가함에 따라 대각선 방향의 유선이 지배적인 형상을 보였다. 이는 베리오그램 모델링과 실내 수리시험 결과에서 확인된 방향성, 연결성 차이와 일관된 결과이다. 즉, 간극의 공간적 배열과 이방성이 강한 채널링을 유발하고, 그 결과 비선형 유동 경향을 강화시킬 수 있음을 정성적으로 재확인하였다.
유동 경로의 채널링 외에도 간극 공간에 국부적으로 발생할 수 있는 와류(eddy flow) 역시 절리를 통한 비선형 유동의 원인이 될 수 있다. 이를 확인하기 위해 HMN9 모델을 접촉 면적 1%, = 100 조건에서 해석했으며, 그 결과를 수압 구배 방향에서 등간격으로 취한 10개의 단면에서의 유동장(velocity field)으로 도시하였다(Fig. 12).
Fig. 12에서 국부적으로 높은 유속(붉은색 영역)이 관찰되었으며, 이는 주로 간극이 급격히 좁아지는 지역에서 발생하였다. 이를 보다 자세히 확인하기 위해 HMN9 모델의 절리면에서 단면 프로파일을 발췌하여 2D 상태에서의 유동 해석을 수행하였다. YZ 평면에 평행한 단면을 생성하였고, 이 단면이 접촉부를 포함하지 않지만 그에 매우 근접하도록 평행 이동하여 2D 해석 도메인을 생성하였다. 정확히는 HMN9 모델에서 x = 11.5 mm인 위치의 YZ 단면을 취득했으며, 이 2D 모델에 대해 수리 해석을 수행하였다. 기본적인 모델 구성 조건은 3D 모델과 동일하며 2D 해석이기 때문에 접촉 면적은 고려하지 않았다. Inlet을 따라 일정 속도로 경계 조건을 부여하였으며, 이때 속도는 = 1, 50, 100에 해당하는 유속을 가하였다. 전체 2D 모델 결과에서 y = 16~21 mm인 부분만 발췌하여 도시하면 Fig. 13과 같다.
Fig. 13(a)에서 확인할 수 있듯, y = 16 mm 부근에서 간극이 국부적으로 작아지며, 이후 다시 넓어지는 형태의 공간이 형성되었다. 따라서 해당 구역의 유속이 상대적으로 높고 이 구역을 통과하면 유속의 재분배가 발생한다. 이처럼 간극의 변화와 그에 따른 속도 변화가 큰 구역에서는 유동 조건에 따라 와류가 발생할 수 있다. Fig. 13(a)는 = 1 조건이므로 전체적인 유동은 층류인 것을 확인할 수 있다. 그러나 Fig. 13(b)처럼 = 50으로 증가하면 절리면의 굴곡이 큰 지역에 와류가 생성된 것을 볼 수 있으며, Fig. 13(c)처럼 가 더욱 증가하면( = 100) 와류가 더 넓은 지역에서 발생하는 것을 확인할 수 있다. 이러한 와류의 생성은 유체가 실질적으로 이동하는 유효 간극(effective aperture)의 크기를 감소시키기 때문에 전체적인 유동의 비선형성이 강화된다(Zou et al., 2015, Wang et al., 2016, Dou et al., 2018). 한편, Table 5에서 HMN9의 상관거리는 규모에 관계 없이 HMN2 보다 작다는 점을 확인하였다. 베리오그램 상관거리가 작다는 것은 짧은 거리에서 공간적 변화가 크다는 것을 의미한다. 이는 다시 Fig. 13과 같은 와류 발생 구역이 존재할 가능성이 상대적으로 높다는 것을 의미하여 결과적으로 HMN9 시험편의 비선형성이 보다 높게 측정된 것으로 해석할 수 있다.
5. 결 론
국내 지질학적 조건을 고려했을 때, 향후 예정된 고준위방사성폐기물처분장은 화강암과 같은 결정질 암반에 건설될 가능성이 높다. 처분장에서 핵종이 유출된 경우, 그 위해성이 인간 생활권에 도달하는 주요 경로는 지하수를 통한 이동이며, 매질 자체의 투수성이 매우 낮은 결정질 암반에서 지하수 유동은 대부분 절리를 통해 발생한다. 따라서 절리를 통한 유체 유동 특성을 정확히 파악하는 것은 처분장 전체의 장기 안전성을 확보하는데 매우 중요한 부분을 차지한다. 본 논문에서는 결정질 암석 절리를 통한 수리 유동의 형태, 즉 선형적, 비선형적 거동 특성을 파악하고 절리 표면의 기하학적 물성이 유동 특성에 미치는 영향을 파악하고자 실내시험과 수치해석을 수행하였다. 수리적으로 가장 대조적인 거동을 보인 두 절리 시험편(HMN2, HMN9)의 결과를 비교하는 형태로 논문을 작성했으며, 이를 통해 도출한 주요 결론은 아래와 같다.
∙간접인장시험과 같은 형태로 생성한 두 인공 절리 시험편 HMN2, HMN9의 표면 기하학적 물성을 측정한 결과, 와 는 대체로 유사하게 측정되었고 는 HMN9가 더 크게 측정되었다. 초기 역학적 간극의 평균은 HMN9 시험편이 더 크게 측정되었으나, 유체 유동의 통로인 간극의 공간 상관성은 HMN2이 높게 측정되었으며, 유동 방향으로 양호한 간극 연결성을 보였다.
∙유체 주입을 통한 수리시험 결과, 두 시험편 모두 수직응력에 의한 유량 감소, 수압 증가에 따른 유량 증가 등 일반적인 수리 거동이 확인되었다. 또한 Forchheimer 방정식을 통해 시험 결과의 비선형을 적절히 모사할 수 있었으며, 수직응력 증가에 따른 Forchheimer coefficients와 수리적 간극의 감소 현상이 확인되었다.
∙두 절리 시험편에 대해 critical Reynolds number를 적용한 결과, 선형, 비선형 유동 관점에서 대조적인 결과가 확인되었다. HMN9 시험편은 수직응력과 수압 조건에 따라 훨씬 빈번한 비선형 유동을 보였으며, 본 논문에서 적용한 수리시험 조건에서 HMN9의 critical Reynolds number는 1.2~4.7 수준으로 평가되었다. 이는 평균 간극이 큰 HMN9에서 critical Reynolds number가 상대적으로 작아지고, 유동에 대한 관성의 영향이 이른 단계에서 유의미하게 발생하기 때문으로 판단된다.
∙두 시험편의 표면 기하학적 물성을 분석한 결과, 전반적인 수리 유동에 영향을 준 요인은 간극의 공간 상관성인 것으로 판단되었다. 유동 방향으로 상관성이 높은 HMN2 시험편은 많은 경우에서 층류 형태를 보였으며, 이에 비해 연결성이 짧은 HMN9 시험편은 채널링 영향을 크게 받아 비선형성이 증가한 것으로 판단되었다.
∙수리 유동 차이를 보다 명확하게 확인하고 비교하고자 간극 공간에 대한 수치 모델을 구성하고 해석을 수행하였다. 수치해석 결과, 수직응력(접촉 면적) 증가에 따라 HMN9 시험편의 채널링 현상이 뚜렷하게 확인되어 간극의 공간적 배열이 전반적인 유동에 영향을 미침을 확인하였다.
∙간극의 공간적 변화가 심한 구역은 국부적으로 와류가 발생할 수 있으며, 이는 유동의 유효 간극을 감소시켜 전반적인 비선형성을 증가시키는 요인이 되었다. 상관거리가 짧은 HMN9 시험편은 작은 규모에서 간극의 공간적 변화가 심하며 이는 와류 발생 확률이 상대적으로 높음을 의미한다. 결과적으로 상기 채널링과 와류 발생 가능성은 HMN9 시험편의 비선형 유동의 빈도와 강도를 강화하는데 기여한 것으로 판단되었다.
∙절리 간극 공간의 구성은 수리 해석 결과에 직접적 영향을 미치나, 본 연구에서는 이를 직접 측정할 수 없었기 때문에 가정하였다. 이는 향후 수치해석, 접촉 이론 도입 등을 통해 개선할 계획이다. 또한 본 연구에서 사용한 시험편은 자연 절리가 아닌 인공적으로 생성한 인장 절리이다. 인공 절리는 표면 거칠기, 충진물 여부 등에서 자연 절리와 수리적 거동의 차이를 보일 수 있다. 따라서 향후 이에 대한 보완이 바람직할 것으로 판단된다.
