Original Article

Tunnel and Underground Space. 31 December 2019. 542-557
https://doi.org/10.7474/TUS.2019.29.6.542

ABSTRACT


MAIN

  • 1. 서 론

  • 2. 3차원 X-ray CT 이미지 획득 및 보정

  •   2.1 X-ray CT 원리

  •   2.2 X-ray CT의 암석역학 분야 활용사례

  •   2.3 X-ray CT 이미지 보정 및 재구성

  • 3. 암석 종류에 따른 X-ray CT값의 분포특성

  • 4. X-ray CT 이미를 이용한 암석 내부구조 평가

  •   4.1 암석의 이방성 정량화

  •   4.2 대표단위길이(Representative Unit Length, LR)

  • 5. X-ray CT 이미지를 활용한 유체거동 해석

  •   5.1 공극 네트워크 모델 (Pore network model)

  •   5.2 격자볼츠만 방법 (lattice Boltzmann method)

  • 6. X-ray CT 이미지를 이용한 암석 불연속면 평가

  • 7. 결 론

1. 서 론

지구내부 지각(crust)의 대부분을 구성하고 있는 물질인 암석은 성인 및 오랜시간의 변성과정에서 결정되어지는 구성광물의 성분, 입자의 배열 및 크기, 공극 또는 균열의 크기 및 분포 양상 등에 의해 다양한 내부구조와 거동 특성을 나타낸다. 암반을 대상으로 하는 대표적인 전통적 산업분야인 석유와 가스 생산, 터널, 지하 유류 비축기지 건설 외에도 최근에 심부지열발전, 이산화탄소 지중저장, 방사성폐기물 처분장 등 환경과 에너지 측면에서 큰 관심을 받고 있는 다양한 산업분야에서 암반의 거동에 대한 평가는 매우 중요한 부분을 차지하고 있다.

암반의 거동 평가는 현지 암반을 대상으로 하는 평가(현지응력 평가, 현지 불연속면 평가 등)를 통해 진행되기도 하지만, 이러한 현지 암반 평가는 매우 제한적이기 때문에 주로 현장에서 채취한 암석 시편을 대상으로 하는 평가방법이 주로 이용된다. 이러한 실내 시험을 통한 암석의 물성 평가는 국제암반역학위원회(ISRM)에서 제안된 방법이 주로 활용된다. 기존에 제안되어 사용되고 있는 대부분의 실내실험을 통한 물성평가 방법은 암석의 내부구조를 별도로 고려하지 하지 않고 결과를 도출한다. 이 경우 불균질성과 이방성이 매우 강한 암석의 경우에는 실험값의 오차가 커지기 때문에 통계적으로 유의한 결과값을 확보하기 위해 많은 실험이 필요하다. 암석 내부의 마이크로 특성(구성광물의 특성 및 결합형태와 관련된 물성, 공극의 크기 분포 등)이 매크로 특성(강도특성, 수리특성 등)에 미치는 영향을 평가하기 위한 절차 및 방법은 아직까지 구체적으로 제안되지 않고 있다. 이는 암석 내부를 3차원으로 시각화하거나 광물단위의 마이크로 물성을 물리적으로 측정하는 것이 매우 어렵기 때문이다.

하지만, 앞서 언급하였듯이 암석의 물성 및 거동특성은 암석 내부구조에서 기인하기 때문에 이러한 내부구조를 평가하여 암석의 마이크로 물성을 산출 할 수 있다면, 매크로 물성 평가시 이러한 마이크로 물성을 적극적으로 활용하여 보다 신뢰도 있는 평가가 가능하다. 마이크로 스케일에서의 암석 특성이 매크로 스케일에서의 암석 거동에 미치는 영향이 지속적으로 보고되고 있으며(Brace, 1965, Lo et al., 1986, Lee et al., 2001), 이러한 특성을 다양한 미소파괴음 모니터링 등의 비파괴 방법을 이용해 평가하기 위한 시도가 계속 이루어지고 있다(Gueguen and Schubnel, 2003, Kemmerer and Oelze, 2012, Kim et al., 2014, Oelze et al., 2002, Sayers, 1988)). 매크로 스케일에서의 암석 물성이 동일하더라도 실제 변형 및 파괴 거동 양상은 마이크로 스케일에서의 암석 구조에 따라 달라질 수 있다는 것도 보고되고 있다(Takemura et al., 2003, Fujii et al., 2007).

특히, 석유 및 가스, 지열, 이산화탄소 지중 저장 등과 같이 암석내의 유체유동이 매우 중요한 산업분야에서는 마이크로스케일에서의 공극구조 특성(공극률, 공극크기분포, 공극 연결도 등)이 암석의 수리-역학적 거동 평가를 위한 핵심정보로 활용되어지고 있다. 특히, 최근 셰일가스 개발이 활발히 진행되면서 매장량 및 생산성 평가를 위한 암석의 평가시 셰일 내부 공극구조의 시각화와 이미지 기반의 유동해석 등이 보편적으로 진행되고 있다(Pini et al., 2012).

암석의 내부구조를 시각화하는 가장 일반적인 방법은 편광현미경, 주사전자현미경, 투과전자현미경 등을 활용한 2차원 단면 관찰이다. 이러한 2차원 단면 이미지를 통해 암석을 구성하고 있는 입자들의 크기, 형상, 미세균열 등의 관찰이 가능하다. 하지만, 최근에는 산업용 X선 단층촬영( computed tomography, 이하 X-ray CT) 기술을 이용하여 내부구조를 3차원으로 시각화하는 시도가 많이 이루어지고 있다. 또한, 암석의 내부구조를 3차원으로 시각화하여 정성적으로 평가하는 사례뿐만 아니라, 공극률 평가, 공극크기분포 평가, 암석내 3차원 균열네트워크 모델을 이용한 유동해석, 재료내부에 존재하는 이방성 및 불균질성 특성평가, 암석 내 불연속면 평가 등에 X-ray CT 이미지를 활용하는 연구가 증가하고 있는 추세이다. 본 연구에서는 X-ray CT 이미지를 활용하여 암석의 물리 역학적 특성 및 수리적 거동 특성을 평가하는 방법에 대하여 고찰하였다.

2. 3차원 X-ray CT 이미지 획득 및 보정

2.1 X-ray CT 원리

X-ray 재료를 투과하면서 에너지가 감쇄되는 X-ray 고유의 성질을 이용하여 재료 내부를 3차원 영상화 하고 특성을 평가하는데 활용되어지는 비파괴 이미징 기술이다. 일반적인 산업용 X-ray CT 장비는 Fig. 1과 같이 X-ray 발생장치, 영상 수집장치(detector), 대상체가 놓여지는 회전부(rotating table)로 구성된다. X-ray 발생장치는 대상체에 따라 투과되는 X-ray의 에너지를 조절하고 생성하는 역할을 한다. 발생한 X-ray는 대상체를 투과하면서 재료의 밀도에 따라 에너지가 감쇄되고, 남은 에너지를 영상수집장치에서 기록하게 된다. 영상수집장치의 각 픽셀에 도달하는 감쇄후 X-ray 양은 CT 이미지를 구성하는 직접적인 정보가 된다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F1.jpg
Fig. 1.

Schematic of typical X-ray CT imaging

CT 이미지는 대상체의 특정 단면 영상을 의미하며, 일반적인 디지털 이미지와의 차이점은 두께를 갖는 다는 것이다. 즉, 대상체를 일정한 두께를 가진 단면으로 잘라서 내부를 볼 수 있게 한 것이 CT 이미지라고 할 수 있다. 각각의 CT 이미지들은 다수의 복셀(voxel)로 구성되어지며, 이러한 복셀은 해당 위치에 있는 재료의 밀도 및 원자번호에 따라 결정되어지는 X-ray 선형감약계수(linear attenuation coefficient)라는 속성을 통해 변환과정을 통해 고유한 CT 값을 가지게 된다. 가장 널리 사용되는 CT값은 물과 공기의 선형감약계수와의 상대적 비율을 통해 변환된 값인 하운스필드 단위(Hounsfield unit)이다. 일반적으로 밀도와 원자번호가 높을수록 큰 CT 값을 가지게 된다. 이러한 CT 이미지 단면들을 연결하여 대상체 전체를 대표하는 3차원 객체를 만드는 작업을 재구성(reconstruction) 이라고 한다. 이러한 재구성 과정을 거쳐 CT 값들로 이루어진 3차원 객체는 재료 내부를 3차원으로 시각화할 수 있으며, 재료의 물성을 평가를 위한 핵심정보로 활용될 수 있다.

2.2 X-ray CT의 암석역학 분야 활용사례

암석의 종류에 따라 구성광물 입자의 공간적 분포 양상이 달라지며 이러한 공간적 분포는 재료의 물리-역학적 특성에 영향을 미치게 된다. 암석의 경우 밀도가 낮은 공극은 X-ray 이미지에서 어두운 색으로, 반면 밀도가 높은 광물은 밝은 색으로 표현되며 광물의 밀도 및 조성에 따라 서로 다른 값으로 표현된다. X-ray CT를 통해 암석을 구성하는 광물의 조성비, 입자의 크기 및 3차원 공간 분포를 정량화 할 수 있으며, CT 이미지가 검출 가능한 충분한 크기의 공극이 존재할 경우 공극의 특성화도 가능하다. X-ray 이미지를 활용하여 암석의 내부 구조 및 물리-역학적 특성화 및 공학적 응용분야에서 연구의 예는 Table 1과 같다.

Table 1. Representative applications in rock mechanics using X-ray CT images

Target subjects Applications
Characterization
of internal
structure
Heterogeneity/
anisotropy
- Quantification of rock heterogeneity (Semnani et al., 2017)
- Quantification of rock anisotropy (Yun et al., 2013)
Grain - Mineral size distribution/shape characterization (Lai et al., 2015, Cnudde et al., 2011, Kim et al., 2016b),
Fracture - Characterization of fracture size & pattern, discontinuity roughness(Diaz et al., 2017, Yang et al., 2016, Gouze et al., 2003)
Engineering
utilization
Physical and
mechanical behavior
- Evaluation of coupled hydro-mechanical fracturing behavior (Kim et al., 2016a, Li et al., 2012, Dann et al., 2011)
Single- phase/multi-
phase fluid flow
- Visualization of fluid flow through fractures (Hirono et al., 2003)
- Effect of irregular pore structure on residual saturation in multi-phase fluid flow (Pini et al., 2012)
Dynamic process
in rocks
- Characterization of faulting process(Martínez-Martínez et al., 2016, Kawakata et al., 1999)
CO2 geological
sequestration
- Evaluation of pore structure of porous rock and quantification of CO2 injection (Al-Menhali, A., 2015, Andrew et al., 2013)

2.3 X-ray CT 이미지 보정 및 재구성

X-ray 이미지는 에너지 빔이 재료를 투과하는 과정 및 영상수집장치에서 이미지를 재구성하는 과정에서 커핑결함(cupping artifact)과 링결함(ring artifact)가 나타날 수 있다. 커핑결함은 X-ray 발생장치로부터 조사된 고주파/저주파 X-ray가 시료를 통과하며 내부로 갈수록 커지는 에너지 감쇄현상으로 인해 발생된다. 일반적으로 시료 표면부보다 시료내부의 이미지 값이 낮게 나타나는 현상으로, Fig. 2(a) 이미지와 같이 중심부는 어둡고 주변부는 밝게 방사형으로 나타난다. 커핑결함이 존재하는 원본이미지에서는 중심부로 갈수록 어둡기 때문에 같은 재료라 하더라도 더 낮은 CT 값을 갖는 반면 주변부로 갈수록 더 높은 CT값을 갖는다. 따라서 원본 이미지의 2차원 곡면 이미지를 적합(fitting) 시킬 수 있는 함수를 지정한 후 특정 기준값과의 차이를 이용하여 보정하면 Fig. 2(b)와 같이 커핑형상이 제거된 이미지를 획득할 수 있다.

링결함은 재료를 투과한 X-ray를 받아들이는 디텍터 패널에서의 결함 때문에 발생하며, 360도 회전하며 이미지를 촬영하면 ‘원형’의 띠가 방사형으로 규칙적으로 나타나게 된다. 동심원으로 발생하는 링결함은 원본 이미지를 극좌표로 변환하여 방사형 띠를 1차원 ‘선’ 형태로 전환시킨 후 2차원 푸리에 변환(Fourier Transform)을 수행하면 규칙적으로 반복되는 결함을 대표하는 특정 주파수 영역이 추출된다. 해당 주파수 영역을 제거한 후 역푸리에 변환(Inverse Fourier Transform)을 하면 Fig 2(c)와 같이 링결함의 제거가 가능하다. 앞서 소개한 2가지 결함요소의 제거를 위한 상세한 알고리즘 및 적용예는 Jung et al.(2011)에 서술된 바와 같다. 노이즈가 보정된 2차원 이미지는 Fig. 3과 같이 3차원에서 재구성되어 암석의 내부구조평가 및 추가 수치해석을 위한 기본 이미지로 활용이 가능하다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F2.jpg
Fig. 2.

Image calibration of raw 3D X-ray Computed Tomography images

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F3.jpg
Fig. 3.

3D reconstruction of calibrated images

3. 암석 종류에 따른 X-ray CT값의 분포특성

암석을 구성하는 다양한 광물의 특징에 따라 X-ray 에너지의 감쇄정도가 달라지므로 X-ray 이미지는 그 자체로 암석 내부 구성 물질에 대한 정보를 담고 있다. 암석의 종류에 따른 2차원 X-ray 이미지 및 해당 이미지의 CT값 히스토그램은 Fig. 4와 같다. Fig. 4(a)의 다공성 암석인 현무암의 경우 공극(어두운 색, P1)과 암질(밝은색, P1)로 이루어진 구조가 이미지에서 명확히 보이며, 히스토그램에서도 약 1500 (P1)과 4500 (P2)에서 최대값을 갖는 가우시안 함수의 합으로 표현된다 (Fig.4의 실선). 반면, 각섬암의 경우 가시적으로 명확한 공극이 존재하지 않고 주로 높은 밀도를 갖는 각섬석과 다소 낮은 밀도를 갖는 사장석으로 구성되어 있어 CT값 영역이 중첩되어 히스토그램 상에 표현되고 (Fig. 4의 점선) P4에 해당하는 CT값에서의 최대값이 P3에 비해 더 높은 확률로 존재한다. 공극이 많이 존재하지 않는 암석들의 경우 구성 물질간의 상대적인 밀도 차이 및 양에 따라 CT값의 히스토그램은 달라진다. Fig. 4(b)의 대리석과 편암의 경우 히스토그램 상의 최대값은 약 4250으로 같으나, 구성 물질의 분포에 따라 표준편차가 달라진다. 즉, 석영, 사장석, 녹염석, 각섬석과 같이 다양한 광물이 혼재된 편암은 큰 표준편차를 갖는 반명, 대리석을 구성하는 광물의 밀도 차는 상대적으로 크지 않으므로 작은 표준편차를 보인다. 화강암은 석영과 장석계열로 주로 구성되어 해당 광물간의 밀도차가 작으며 이로 인해 대리암과 유사한 표준편차 값을 보이는 반면, CT값의 최대값은 대리석보다 작은 약 3700에 존재한다. Fig. 4의 CT값 히스토그램에서 보이듯 다양한 광물이 존재하는 암석의 경우 각 구성 물질을 반영하는 CT값의 차이가 크지 않아 개별 광물을 완전히 분리하는 과정이 쉽지 않다. 따라서 X-ray를 이용한 암석의 내부 구조 연구는 획득된 2차원 이미지에 대한 정성적인 기술, 혹은 명백한 형태로 공극이 존재하는 다공성 암석에서의 이진화를 통해 공극 구조를 추출한 후 내부 구조를 분석하고 추가적인 수치해석기법을 적용하는 연구가 주로 이루어져 왔다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F4.jpg
Fig. 4.

CTnumber histograms of XCT images of various rock types

4. X-ray CT 이미를 이용한 암석 내부구조 평가

X-ray 이미지는 암석 내부 구조를 가시화하여 보여줄 수 있다는 매우 큰 장점이 있으나, 불균질하며 이방성이 존재하는 암석 구조의 특성화를 정량적인 방법으로 수행하기 위해서는 이미지에 기반한 별도의 알고리즘 개발이 필요하다. 본 장에서는 암석에서 흔히 관찰되는 이방성 평가법과 내부 구성물질의 대표 크기를 추출하는 알고리즘을 소개한다.

4.1 암석의 이방성 정량화

Slicing Plane Method(SPM)은 3차원 이미지 정보에 특정 단위벡터(unit vector)를 법선벡터(normal vector)로 갖는 평면들을 이용해 통계계수를 획득하는 방법으로 이방성이 존재하는 면의 방향에 대한 정보를 추출할 수 있으며 그 과정은 다음과 같다 (Yun).

동일한 normal vector를 갖는 평면들을 3차원 이미지에 가상으로 끼워 넣은 후, 서로 다른 면적을 갖는 해당 면이 교차하는 이미지의 CT값들의 평균값을 계산한다. 이후 해당 평면의 중첩면적을 구한 후 식 (1)과 같이 가중평균(Weighted arithmetic mean)과 가중표준편차(Weighted standard deviation), 그리고 변동계수(Coefficient of variation)를 산출한다.

$$s_w^{SPM,l}=\sqrt{\frac{\displaystyle\sum_{i=1}^NW_{i`}{(x_i-\overline{x_W}\;)}^2}{\displaystyle\frac{(N'-1)\sum_{i=1}^NW_i}{N'}}},\;\;s_w^{SPM,l}\;:\;weighted\;s\tan dard\;deviation$$ (1)

여기서, Wi : 가중치, N′: 0이 아닌 가중치, xW : 가중평균

33차원에서 존재할 수 있는 모든 가능한 법선벡터 방향으로 위 과정을 반복하면, 각 법선벡터마다 특정한 변동계수 획득이 가능하다. Fig. 5의 가상 도메인 내부에는 평행한 층리구조가 생성되어 있으며, 이방성이 우세한 면과 수직한 방향으로 법선벡터가 존재할 경우, 가상 절단면(virtual slicing plane)에서 계산된 CT값의 평균값이 크게 변하여 변동계수 값이 커지는 반면, 이방성 면과 평행한 방향으로 존재하는 법선벡터의 경우 CT값 평균값의 변화가 작다. 따라서 각 법선벡터에서 산출한 변동계수값을 평사투영법으로 표현하면 Fig. 5와 같이 특정방향으로 발달된 이방성을 정량적으로 평가할 수 있다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F5.jpg
Fig. 5.

Schematic illustration of Slicing Plane Method (Yun et al., 2013)

Fig. 6(a)는 사암의 3D X-ray CT 이미지를 분석한 결과로, 이미지상에 수직방향으로 지배적으로 존재하는 층리구조의 특성을 반영하여, Fig. 5의 예에서와 같이 정확한 방향을 잡아내고 있다. 또한 Fig. 6(b)의 화강암과 같이 여러 방향으로 중첩되어있는 복잡한 구조에서도 각 방향으로 발달한 이방성을 평가 가능하다. 본 분석기법은 재료와 스케일에 관계없이 확보한 이미지에 대한 내부 이방성의 방향을 평가할 수 있어 폭넓은 적용성을 갖고 있다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F6.jpg
Fig. 6.

Stereonet results from SPM (Yun et al., 2013)

4.2 대표단위길이(Representative Unit Length, LR)

Two–point correlation 기법은 일반적으로 이진화된 재료(즉, 암석의 경우에는 광물과 공극이 명확히 분리되는 경우)에 널리 이용되어 왔으나(Torquato, 1988), 공극률이 낮고 광물간의 분리가 어려운 경우에는 적용이 어렵다. 이를 위해 재료 분리(segmentation) 과정 없이, 3차원 X-ray 원본 이미지를 그대로 이용하여 내부 구조의 대표크기를 산출하는 분석기법인 Modified two-point correlation을 소개한다. Fig. 7과 같이 길이 L을 갖는 2개의 쌍으로 된 n개의 점들(P1n-P2n)을 무작위로 분포시킨 후 P1과 P2에 해당하는 CT값의 차이를 계산하고 n개의 차이값들로부터 산술평균(μ)을 구한다. P1과 P2간의 거리 L값을 0부터 도메인의 최대 크기까지 각 L에 대한 산술평균을 동일한 과정을 통해 계산하여 거리 L에 대한 산술평균의 변화를 평가하면 Fig.7의 결과와 같다. L이 0에 가까울수록 P1과 P2가 동일한 구성 물질에 위치할 확률이 증가하므로 차이값의 산출평균은 0에 가까운 값을 갖는다. L이 증가할수록 서로 다른 물질에 P1, P2가 놓일 확률이 증가하게 되며, 특정 L값에서 최대 변곡점을 갖게 된다. 최대 변곡점에서의 L값을 대표단위길이(LR)로 정의한다. 제안된 방법이 평가하는 크기(LR)는 통계적인 평균값을 이용하기 때문에 해당 암석에서의 광물 크기뿐만 아니라 공극 및 미세구조의 크기 모두를 반영한 대푯값이다. 대표적인 예로 Fig. 7과 같이 사암과 화강암에 적용할 경우, 획득한 대푯값이 화강암(LR=0.85 mm)에 비해 상대적으로 미세한 입자들로 구성되어 있는 사암(LR=0.15 mm)에서 더 작아 암석의 지배적인 구조의 크기를 잘 잡아냄을 확인할 수 있다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F7.jpg
Fig. 7.

Schematic image of the distribution of randomly spread two point sets throughout the image and results from the granite and sandstone images

5. X-ray CT 이미지를 활용한 유체거동 해석

5.1 공극 네트워크 모델 (Pore network model)

암석내부의 공극은 유체가 흐르는 주된 도메인으로써, 3차원에서 서로 복잡한 유동 경로를 갖는 공극구조 및 이를 통한 유체거동 평가는 많은 연구가 이루어져 왔다. 특히 공극 네트워크 모델은 공극 구조를 공극방(Pore chamber)과 공극목(Pore throat)으로 이루어진 관망으로써 이상화하여 유체 거동을 수행하는 해석 기법으로써 지반공학, 석유공학, 환경공학분야 등에서 널리 이용되고 있다. 1990년대 이전까지 공극 네트워크 모델링은 주로 2차원 혹은 3차원 격자구조 관망 도메인에 공극의 특징을 반영하여 수행되어 왔으나(Lenormand et al., 1983), 이후 3차원 이미징 기법의 발달로 재료의 3차원 영상으로부터 공극의 위상학적 특징이 보존된 관망을 추출하여 모델링하는 것이 일반화되었다(Suh et al., 2016). 공극 네트워크 모델은 공극스케일에서의 재료 내 다상 유체 해석 기법 중 가장 빠르고 간편한 기법으로 알려져 있으며(Blunt, 2001, Joekar-Niasar and Hassanizadeh, 2012), 이를 통한 유동 해석으로부터 획득할 수 있는 물성으로는 공극도분포(Pore size distribution), 함수특성곡선(Water retention curve), 상대투수율곡선(Relative permeability curve) 등이 있다. 유동 해석 시 공극목의 크기와 형상은 공극 내 모세관압을 결정하는 매우 중요한 인자이며, 공극 네트워크 모델은 이를 원 혹은 다각형으로 단순화하여 해석을 수행한다는 내재적 한계점을 지닌다. 그러나 최근 3차원 X-ray 영상으로부터 개별 공극목의 실제 단면 형상정보를 추출하여 보다 정확한 유체 해석을 수행할 수 있는 방법이 제안된 바 있다(Suh et al., 2017, Suh and Yun, 2017). 3차원 X-ray 영상으로부터 얻은 공극 네트워크를 통한 시뮬레이션의 개략적 순서는 Fig. 8과 같다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F8.jpg
Fig. 8.

Image-based pore network modeling sequence

재료의 3차원 X-ray CT 이미지로부터 해석을 수행하고자 하는 관심영역(Region of interest, ROI)을 추출하고, 이를 이진화(Binarization) 처리하여 공극 구조를 추출한다. 이진 공극 영상에 세선화알고리즘(Thinning algorithm)을 적용하여 공극 스켈레톤(Pore skeleton)을 획득할 수 있으며, 이는 위상학적으로 연속된 공극 단면의 도심을 연결한 곡선과 동일하다. 일반적으로 공극 스켈레톤을 구성하는 하나의 복셀과 연결된 인접복셀수를 평가하여, 그 개수가 3개 이상일 경우 이를 공극 관망의 합류점으로 판단해 공극방으로써 정의한다. 이 때, 공극방은 공극 내 압력분포 및 공극의 부피를 대변하게 된다. 공극방과 공극방 사이의 복셀의 단위집합은 공극 네트워크의 하나의 관을 형성하며, 이 중 관벽까지의 유클리디언 거리(Euclidean distance)가 가장 작은 영역 (공극 내 최소내접구의 중심)을 모세관압이 작용하는 공극목의 중심으로써 정의한다. 본 과정을 통해 공극 크기분포를 정량적으로 평가할 수 있으며, 전체 도메인의 공극도분포 뿐만 아니라, 공극방과 공극목의 크기 분포 및 공극의 연결도 분석이 가능하다(Suh, 2017).

공극 네트워크 모델을 통한 다상 유체 시뮬레이션은 비압축성 유체를 대상으로 하며, 정상 상태 하에서의 침투해석이 주를 이룬다. 이는 특정 유체로 포화된 육면체 형태의 공극 네트워크 내에서 침투 방향을 결정하고, 그에 수직인 하나의 평면(침투 평면)에 침투 유체를 주입하여 기존 유체를 그에 평행한 나머지 평면(배수 평면)을 통해 배수시키는 해석 방법을 의미한다. 이 때, 유체의 침투방향에 수직인 침투평면과 배수평면은 배수경계조건으로, 침투방향과 평행한 나머지 4개의 평면은 불투수 경계조건으로 선정하는 것이 일반적인 방법이다. 공극 네트워크 추출 단계에서 자동적으로 결정되는 공극목의 모세관압은 Young-Laplace식으로부터 계산되며, 공극목 양단의 공극방을 서로 다른 유체가 점유하고 있는 경우, 양단의 압력 차이가 공극목의 모세관압을 초과하게 되면 침투가 발생하게 된다. 이러한 방법을 도메인 전반에 적용하면 공극 네트워크의 침투평면과 배수평면에 특정한 압력차를 가했을 때의 정상상태 침투 유체의 점유 양상 분석이 가능하다.

5.2 격자볼츠만 방법 (lattice Boltzmann method)

공극네트워크 모델에서 이진화된 공극 이미지를 구성하여 유체 거동을 해석하는 것과 달리 복잡한 공극 구조를 직접 사용하여 유동 해석을 수행할 수 있다. 이를 위해 일반적인 CFD 방법을 사용할 수도 있으나 최근에는 격자 볼츠만 방법 (lattice Boltzmann method) 을 많이 사용하고 있다. 격자 볼츠만 방법은 일반적인 단상 유동 뿐 아니라 오염 물질 전달, 열 유동, 다상 유동 등과 같은 복잡한 유동 해석과 포와송 방정식과 같은 편미분 방정식을 푸는데도 적용되고 있다. 격자 볼츠만 방법은 기존 연속체 기반의 CFD 방법에 비해 계산량이 상대적으로 많으나 내재적인 병렬성으로 인해 연산 효율이 높다. 격자 볼츠만 방법에서 다상 유동 해석을 하기 위해 제안된 모델로는 Rothman-Keller(RK) 모델, Shan-Chen(SC) 모델, free energy 모델, mean-field 모델 등이 있으며, 다공성 암석과 같은 복잡한 공극 구조에서는 RK 모델과 SC 모델이 주로 사용되고 있다. 다공성 암석에서 격자 볼츠만 방법을 이용한 유동 해석으로부터 획득할 수 있는 물성으로는 절대 투수율 (Absolute Permeability), 상대 투수율(Relative Permeability), 함수 특성 곡선(Water Retention Curve), 잔여 포화도 (Residual Saturation) 등이 있다. Fig. 9는 3차원 X-ray 이미지를 이용하여 유체해석을 할 수 있는 격자 볼츠만 해석의 순서를 개략적으로 보여준다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F9.jpg
Fig. 9.

Lattice Boltzmann simulation based on X-ray CT image

다공성 암석 시편의 CT 이미지에서 해석을 수행하고자 하는 관심영역(Volume of interest)를 추출하고 공극복셀(Pore voxel)과 고체복셀(Solid voxel)로 이진화 한다. 단상 유동 해석의 경우 압력 차이를 이용한 유동 해석과 중력을 이용한 유동 해석을 주로 수행하며 압력 차이를 이용한 해석은 유입부와 유출부에 각각 압력 경계를 적용하고 유동 방향과 평행한 경계 면들은 고체 벽으로 감싸 불투수영역으로 모사한다. 중력을 이용한 해석의 경우 유입부와 유출부에 압력 경계를 적용하지 않고 주기적 경계를 적용하게 된다. 이때 유입부와 유출부의 형상이 동일해야 하므로 미러링(mirroring)을 이용한 대칭적인 도메인 구조를 생성하고 유동 방향과 평행한 경계면은 동일하게 고체 벽으로 감싼다. 압력 차이를 이용한 해석은 적용한 압력 차이에 따라 밀도 차이가 비례해서 발생하게 되는 문제가 있으며 압력 경계의 수치 안정성에 따라 문제가 생길 수도 있다.

다상 유동 해석의 경우 함수 특성 곡선을 얻기 위한 해석과 비정상, 정상 상태의 코어 주입 해석이 있다. 함수 특성 곡선 해석의 경우 실험과 같이 강한 습윤성을 가진 다공성 플레이트를 유출 부에 적용하고 양단에 일정한 압력을 가하면 적용된 압력 차에 의해 비습윤성 유체가 침투하며 평형 상태에 도달할 때까지 해석을 수행하고 이후 순차적으로 압력을 바꿔 가면서 동일한 과정을 반복한다. 비정상 상태의 코어 주입 해석은 일정한 유량으로 비습윤성 유체 또는 습윤성 유체를 유입 부를 통해 주입하며 유출부의 경우 압력 경계를 적용한다. 비정상 상태의 코어 주입 해석은 두 유체의 경계가 유출부 경계 조건에 도달할 시 수치적 불안정성으로 인해 해석이 중단될 가능성이 크다. 정상 상태의 코어 주입 해석은 중력을 이용한 단상 유동 해석과 동일하게 이루어지며 유일한 차이점은 초기에 두개의 유체가 정해진 포화도에 따라 배치된다는 점이다. 따라서 상태 투수율을 얻기 위해서는 포화도를 다르게 맞추어 여러 번의 해석을 수행해야 한다.

도메인의 구성, 경계 조건, 초기 조건 적용이 완료되면 격자 볼츠만 해석을 수행하게 되며 정해진 반복 횟수에서 모든 공극복셀에서의 밀도, 속도, 압력 등을 저장하게 된다. 이후 후처리 과정을 통해 결과들을 시각화하고 분석한다.

6. X-ray CT 이미지를 이용한 암석 불연속면 평가

암석 내부의 절리면과 같은 불연속면의 평가는 암반의 역학적 거동 평가뿐만 아니라 수리거동 평가시에도 매우 중요하다. 일반적으로 암반의 불연속해석을 위해 불연속면에 대한 거칠기, 간극 등과 같은 형태적 인자들에 대한 정보가 필요하다. 일반적으로 불연속면의 거칠기의 평가시 노출되어 있는 절리면에 대해서는 물리적, 광학적 프로파일러를 이용하여 측정이 가능하다. X-ray CT 이미지를 이용하게 되면, 노출되어 있는 절리면 뿐만 아니라 노출되어 있지 않은 암석 내부의 절리면에 대한 거칠기와 거칠기 이방성의 평가가 가능하다(Diaz et al., 2017). 노출되지 않은 불연속면을 CT 기술을 이용해 가시화 할 수 있는 것은 여러모로 장점이 있다. 실제 심부의 암석 코어에 대한 분석시, 노출되어 있는 절리면에 대한 분석은 매우 제한적이다. 하지만 코어 내부에 존재하는 절리는 상대적으로 수가 많으며, 방향성 및 거칠기 정보를 그대로 유지한 상태로 존재하기 때문에 많은 정보를 제공할 수 있다. 거칠기 정보뿐만 아니라 절리간극(aperture), 충진물 상태 등 암반의 불연속 거동에 영향을 미치는 중요한 요소들도 동시에 평가가 가능하게 된다.

CT를 이용하여 암석 코어시료 내부의 절리면을 평가하는 과정을 Fig. 10에 나타내었다. 암석에 무수히 발달되어 있는 미세균열과는 달리 절리와 같은 불연속면은 규모가 상대적으로 크고 CT 이미지에서 충분히 검출이 가능하기 때문에 분석이 용이하다. 암석의 CT 이미지를 통해 불연속면만을 효율적으로 추출해 내는 과정을 거치게 되면, 해당 면에 대한 3차원 프로파일을 얻을 수 있다. 이 프로파일을 3차원 매트릭스로 구성하게 되면 각 방향에 따른 거칠기 지수 값을 평가할 수 있다. Diaz et al.(2017)은 심부에서 생성된 절리면에 거칠기 이방성이 존재하며, 이러한 이방성이 해당 구간의 여러 절리면에서 일관성 있는 방향을 나타낸다는 것을 확인하였다. 이러한 절리면 이방성은 실제로 암반의 수리-역학 전단 거동에 영향을 미칠 뿐만 아니라 절리가 생성될 당시의 현지응력(in-situ stress)과도 상관성이 높을 것으로 보이기 때문에 매우 중요한 정보가 된다.

http://static.apub.kr/journalsite/sites/ksrm/2019-029-06/N0120290613/images/ksrm_29_06_13_F10.jpg
Fig. 10.

Joint surface roughness characterization using X-ray CT

7. 결 론

암석의 물리적 특성 및 다양한 거동 특성은 내부의 구성물질의 분포 특성에 지배적인 영향을 받는다. 암석의 성인과 변성에 따른 내부구조의 마이크로 특성이 실제 암석의 거시적 거동 특성에 미치는 영향이 매우 크다는 것은 자명하지만, 지금까지의 시험법은 이러한 마이크로 특성을 고려하여 암석을 평가하기 어려웠다. 이러한 점 때문에 다양한 암반공학적 설계 문제에서 주로 활용하는 암석의 여러 가지 기초물성들은 암반의 잠재적이고 복잡한 거동 양상을 반영할 수 없는 한계가 있다. 예를 들면, 동일한 강도를 갖는 사암과 화강암을 강도만 가지고 평가한다면 두 암종이 나타내는 매우 상이한 수리-역학적 거동특성은 고려될 수 없게 된다.

본 연구에서는 전통적인 암석의 평가 방법의 한계를 극복하기 위한 방안으로 암석 내부의 특성화에 기반한 새로운 암석의 평가방안을 고찰하였다. 암석의 불균질 특성 및 이방성 특성의 정량화, 암석의 구성광물 입자의 크기분포 및 형상특성, 공극이미지를 이용한 유동해석, 암석내부의 노출되지 않은 절리면 거칠기 평가 등 전통적인 암석의 시험법으로 측정하기 어려웠던 중요한 암석의 특성들이 X-ray CT를 이용하여 평가 될 수 있음을 확인하였다. 현재 의료분야 이외의 다른 산업분야에서 X-ray 활용이 급증하고 있으며, CT 장비의 제조사들은 여러 산업분야의 요구도에 맞추어 장비를 개발하고 있다. 암석, 콘크리트 등 건설산업분야에서도 이러한 활용도가 급증하고 있으며, X-ray 기반의 재료 평가방안을 표준화하려는 노력이 계속되고 있다. 암석역학 분야에서도 향후 추가적인 관련 연구들을 통해 X-ray 기반 분석의 활용성과 신뢰도가 확보된다면 새로운 표준시험법으로 제안될 수 있을 것이다.

Acknowledgements

이 논문은 2019년도 한국해양대학교 신진교수정착연구사업 및 2015년도 한국연구재단의 지원(No.2011-0030040, 2016 R1A2B4011292)을 받아 수행된 연구임.

References

1
Al-Menhali, A., Niu, B., Krevor, S., 2015, Capillarity and wetting of carbon dioxide and brine during drainage in Berea sandstone at reservoir conditions, Water Resources Research, 51, pp. 7895-7914.
10.1002/2015WR016947
2
Andrew, M., Bijeljic,B., Blunt, M.J., 2013, Pore-scale imaging of geological carbon dioxide storage under in situ conditions, Geophysical Research Letters, 40, pp. 3915-3918.
10.1002/grl.50771
3
Blunt, M.J., 2001, Flow in porous media-pore-network models and multiphase flow, Current opinion in colloid & interface science, 6, pp. 197-207.
10.1016/S1359-0294(01)00084-X
4
Brace, W.F., 1965, Some new measurements of linear compressibility of rocks, J. Geophys. Res., 70, pp. 391-398
10.1029/JZ070i002p00391
5
Cnudde, V., Boone, M.N., Dewanckele, J., Dierick, M., Hoorebeke, V., Jacobs, P., 2011, 3D characterization of sandstone by means of X-ray CT., Geosphere, 7, pp. 1-8.
10.1130/GES00563.1
6
Dann, R., Turner, M., Close, M., Knackstedt, M., 2011, Multi-scale characterisation of coastal sand aquifer media for contaminant transport using X-ray computed tomography, Environmental Earth Science, 63, pp. 1125-1137.
10.1007/s12665-010-0788-8
7
Diaz, K., Kim, K.Y., Yeom, S., Zhuang, L., Park, S., Min, K-B., 2017, Surface roughness characterization of open and closed rock joints in deep cores using X-ray computed tomography, International Journal of Rock Mechanics & Mining Sciences, 98, pp. 10-19.
10.1016/j.ijrmms.2017.07.001
8
Fujii, Y., Takemura, T., Takahashi, M., Lin, W., 2007, Surface features of uniaxial tensile fractures and their relation to rock anisotropy in Inada granite, Int. J. Rock Mech. Min. Sci., 44, pp. 98-107
10.1016/j.ijrmms.2006.05.001
9
Gouze, P., Noiriel, C., Bruderer, C., Loggia, D., Leprovost, R., 2003, X-ray tomography characterization of fracture surfaces during dissolution, Geophysical Research Letters, 30(5), 1267.
10.1029/2002GL016755
10
Gueguen, Y., Schubnel, A., 2003, Elastic wave velocities and permeability of cracked rocks, Tectonophysics, 370, pp. 163-176
10.1016/S0040-1951(03)00184-7
11
Hirono, T., Takahashi, M., Nakashima, S., 2003, In situ visualization of fluid flow image within deformed rock by X-ray CT, Engineering Geology, 70, pp. 37-46.
10.1016/S0013-7952(03)00074-7
12
Joekar-Niasar, V., Hassanizadeh, S.M., 2012, Analysis of Fundamentals of Two-Phase Flow in Porous Media using Dynamic Pore-Network Models: A Review, Critical Reviews in Environmental Science and Technology, 42, pp. 1895-1976.
10.1080/10643389.2011.574101
13
Jung, Y.J., Yun, T.S., Kim, K.Y., Choo. J., 2011, Image Calibration Techniques for Removing Cupping and Ring Artifacts in X-ray Micro-CT Images, Journal of the Korean geotechnical society, 27(11), pp. 93-101.
10.7843/kgs.2011.27.11.093
14
Kawakata, H., Cho, A., Kiyama, T., Yanagidani, T., Kusunose, K., Shimada, M., 1999, Three-dimensional observations of faulting process in Westerly granite under uniaxial and triaxial conditions by X-ray CT scan, Tectonophysics, 313, pp. 293-305.
10.1016/S0040-1951(99)00205-X
15
Kemmerer, J.P., Oelze, M.L., 2012, Ultrasonic assessment of thermal therapy in rat liver, Ultrasound Med. Biol., 38, pp. pp. 2130-2137
10.1016/j.ultrasmedbio.2012.07.02423062365PMC3511621
16
Kim K.Y., Zhuang, L., Yang, H., Kim, H., Min, K-B., 2016, Strength Anisotropy of Berea Sandstone: Results of X-Ray Computed Tomography, Compression Tests, and Discrete Modeling, Rock Mech Rock Eng, 49. pp. 1201-1210.
10.1007/s00603-015-0820-0
17
Kim, G., In, C.-W., Kim, J.-Y., Kurtis, K.E., Jacobs, L.J., 2014, Air-coupled detection of nonlinear Rayleigh surface waves in concrete-Application to microcracking detection, Ndt E Int., 67, pp. 64-70
10.1016/j.ndteint.2014.07.004
18
Kim, K.Y, Suh, H.S., Yun, T.S., Moon, S.-W., Seo, Y.S., 2016, Effect of particle shape on the shear strength of fault gouge, Geosciences Journal, 20(3), pp. 351-359
10.1007/s12303-015-0051-0
19
Lai, P., Moulton, K., Krevor, S., 2015, Pore-scale heterogeneity in the mineral distribution and reactive surface area of porous rocks., Chemical Geology, 411, pp. 260-273.
10.1016/j.chemgeo.2015.07.010
20
Lee, S.-E., Cho, S.-H., Seo, Y.-S., Yang, H.-S., Park, H.-M., 2001, The effect of microcracks on the mechanical anisotropy of granite, J. Soc. Mater. Sci., 50, pp. 7-13
10.2472/jsms.50.3Appendix_7
21
Lenormand, R., Zarcone, C., Sarr, A., 1983, Mechanisms of the displacement of one fluid by another in network of capillary ducts, Journal of Fluid Mechanics, 135, pp. 337-353.
10.1017/S0022112083003110
22
Li, J., Liu, J., Trefry, M.G., Liu, K., Park, J., Haq, B., Johnston, C.D., Clennell, M.B., Volk, H., 2012, Impact of Rock Heterogeneity on Interactions of Microbial-Enhanced Oil Recovery Processes., Trasport in Porous Media, 92, pp. 373-396.
10.1007/s11242-011-9908-5
23
Lo, T.-W., Coyner, K.B., Toksöz, M.N., 1986, Experimental determination of elastic anisotropy of Berea sandstone, Chicopee shale, and Chelmsford granite, Geophysics, 51, pp. 164-171
10.1190/1.1442029
24
Martínez-Martínez, J., Fusi, N., Galiana-Merino, J.J., Benavente, D., Crosta, G.B., 2016, Ultrasonic and X-ray computed tomography characterization of progressive fracture damage in low-porous carbonate rocks, Engineering Geology, 200, pp. 47-57.
10.1016/j.enggeo.2015.11.009
25
Oelze, M.L., O'Brien, W.D., Darmody, R.G., 2002, Measurement of attenuation and speed of sound in soils, Soil Sci. Soc. Am. J., 66, pp. 788-796
10.2136/sssaj2002.7880
26
Pini, R., Krevor, S.C.M., Benson, S.M., 2012, Capillary pressure and heterogeneity for the CO2/water system in sandstone rocks at reservoir conditions, Advances in Water Resources, 38, pp. 48-59.
10.1016/j.advwatres.2011.12.007
27
Sayers, C., 1988, Stress-induced ultrasonic wave velocity anisotropy in fractured rock, Ultrasonics, 26, pp. 311-317
10.1016/0041-624X(88)90028-5
28
Semnani, S.J., Borja, R.I., 2017, Quantifying the heterogeneity of shale through statistical combination of imaging across scales, Acta Geotechnica, DOI 10.1007/s11440-017-0576-7.
10.1007/s11440-017-0576-7
29
Suh, H.S., 2017, Estimation of water retention characteristics of geomaterials by pore network simulation, M.S. thesis, Yonsei University.
30
Suh, H.S., Kang, D.H., Jang, J., Kim, K.Y., Yun, T.S., 2017, Capillary pressure at irregularly shaped pore throats: Implications for water retention characteristics, Advances in Water Resources, doi: 10.1016/j.advwatres.2017.09.025.
10.1016/j.advwatres.2017.09.025
31
Suh, H.S., Yun, T.S., 2017, Modification of capillary pressure by considering pore throat geometry with the effects of particle shape and packing features on water retention curves for uniformly graded sands, Computers and Geotechnics, doi: 10.1016/j.compgeo.2017.10.007.
10.1016/j.compgeo.2017.10.007
32
Suh, H.S., Yun, T.S., Kim, K.Y., 2016, Prediction of Soil-Water Characteristic Curve and Relative Permeability of Jumunjin Sand Using Pore Network Model, Journal of the Korean geotechnical society, 32(1), pp. 55-62.
10.7843/kgs.2016.32.1.55
33
Takemura, T., Golshani, A., Oda, M., Suzuki, K., 2003, Preferred orientations of open microcracks in granite and their relation with anisotropic elasticity, Int. J. Rock Mech. Min. Sci., 40, pp. 443-454
10.1016/S1365-1609(03)00014-5
34
Torquato, S,, Beasley, J.D., Chlew, Y.C., 1988, Two-point cluster function for continuum percolation, J. Chem. Phys., 88(10), pp. 6540-6547.
10.1063/1.454440
35
Yang, S-Q., Ju, Y., Gao, F., Gui, Y-L., 2016, Strength, Deformability and X-ray Micro-CT Observations of Deeply Buried Marble Under Different Confining Pressures, Rock Mech Rock Eng, 49, pp. 4227-4244.
10.1007/s00603-016-1040-y
36
Yun, T.S., Jeong, Y.J., Kim, K.Y., Min, K-B., 2013, Evaluation of rock anisotropy using 3D X-ray Computed Tomography, Engineering Geology, 163, pp. 11-19.
10.1016/j.enggeo.2013.05.017
페이지 상단으로 이동하기