1. 서 론

가스 하이드레이트 퇴적층에서의 에너지 회수 효율 예측, 지하수 내 NAPL오염 정화 등과 같은 문제에 있어서 공극 내 모세관압(Capillary pressure), 상대투수율(Relative permeability)과 같은 불포화 지반물성치의 파악은 매우 중요하다(Valvante, 2004; Jang and Santamarina, 2014). 특히, 상대투수율은 경우에 따라서 실험적으로 획득하기 어렵기 때문에(Li and Horne, 2006), 함수특성곡선(SWCC) 실험 결과에 Mualem-van Genuchten식과 같은 수학적 모델을 적용하여 예측하기도 한다. 하지만 SWCC실험을 통한 함수특성곡선 획득 역시 오랜 시간을 필요로 하기 때문에 그 대안으로 공극스케일에서 다상 유체 거동을 모사할 수 있는 많은 수치해석 모델이 제안되었다(Raoof et al., 2013).

1956년 Fatt에 의하여 최초로 공극 네트워크 모델링을 통한 다상 유체 해석이 이루어진 후 지반공학분야 뿐만 아니라 석유공학, 환경공학분야에서도 공극 네트워크 모델링은 현재까지도 널리 이용되고 있다(Blunt, 2001; Dong et al., 2009). 공극 네트워크 모델링이란, 다공성 재료의 공극방(Pore chamber)과 공극목(Pore throat)의 연결도를 바탕으로 공극을 하나의 관망으로 대체하여 유체 해석을 수행하는 기법을 말한다. 1990년대 이전까지의 공극 네트워크 모델링은 주로 2차원 혹은 3차원 격자구조 관망에 공극의 특징을 반영하여 수행되어 왔으나, 1992년 Bryant 등의 연구 이후 X-ray CT 등을 통해 얻은 다공성 재료의 3차원 영상으로부터 직접 공극 네트워크 관망을 추출하여 모델링 하는 것이 일반화 되었다. 3차원 영상을 통한 공극 네트워크 추출은 다방면 검색 방법(Zhao, 1994), 보로노이 다이어그램 방법(Bryant and Blunt, 1992; Øren and Bakke, 2002), 중앙축(Medial axis) 방법(Baldwin et al., 1996; Liang et al., 1999; Prodanović et al., 2006), 최대내접구 방법(Silin and Patzek, 2006) 등 여러 방법이 존재하며 이 중에서 중앙축 방법은 공극의 위상학적 형태가 보존되는 장점이 있다(Dong, 2007).

이에 따라 본 연구에서는 세선화알고리즘(Thinning algorithm)을 이용한 중앙축 방법을 수정하여 주문진표준사의 공극 네트워크를 추출하고, 추출한 공극 네트워크 상에 SWCC실험을 모사하여 다상 유체 거동 해석을 실시하였다. 또한 수치해석 결과와 주문진표준사의 SWCC실험 결과를 비교하였다.

2. 해석 방법

2.1 X-ray CT 영상 처리

내경 11mm의 원통형 용기에 주문진표준사를 다지지 않고 담아 3차원 X-ray CT 촬영을 수행하였다. 촬영된 영상은 2048×2048 크기의 2차원 영상 2048장으로 구성되어 있으며, 복셀(Voxel)의 한 변의 길이는 0.0110mm이다. 2차원 영상을 수직방향으로 쌓아 구성한 2048× 2048×2048 크기의 3차원 도메인 가운데 용기 벽면의 영향 및 도메인 축소효과에 의한 영향이 없는 영상의 중심부 406×406×406를 임의로 정하여 수치해석에 사용하였다. 해석 도메인은 임계값(Threshold value) 3150을 기준으로 이진화(Binarization) 처리하여 Fig. 1(a), 1(b)와 같이 흙 입자 및 공극의 3차원 구조를 추출하였으며, 공극률 및 상대밀도는 0.4531과 15.37%로 계산되었다.

3차원 이진 흙 입자 영상에 플러드필알고리즘(Flood- fill algorithm)을 적용하여 개별 입자들을 분리한 뒤, 해당 입자들의 직경을 계산하여 입도 분포 곡선을 획득 하였으며 이를 주문진표준사의 체분석 시험 결과(Yang, 2014) 및 시험 결과에 수학적 모델인 Fredlund식(Fredlund et al., 1994)을 적용한 곡선과 비교하였다(Fig. 1(c)).

3차원 이진 공극 영상에 세선화알고리즘을 적용하여 얻은 결과물인 공극 스켈레톤(Pore skeleton)은 기하학적으로 공극목의 단면의 도심을 연결한 직선과 동일하다. 3차원 이진 공극영상을 유클리디언 거리 변환(Euclidean Distance Transform)한 뒤, 공극 스켈레톤에 해당하는 부분만을 추출하여 공극도 분포곡선을 획득하였다(Fig. 1(d)).

2.2 공극 네트워크 추출

공극 네트워크는 공극방을 나타내는 노드와 그들을 연결하는 공극목을 나타내는 튜브로 이상화되며 네트워크 추출 기법에 따라서 노드 및 튜브의 개수가 달라질 수 있다. 본 연구에서는 해당 복셀과 연결된 인접복셀수를 평가하여 그 개수가 3개 이상일 경우 채널의 합류점으로 판단하여 노드로 결정하였고, 노드와 노드를 연결하는 복셀의 단위집합을 채널로 정의하였다. 이 때, 노드는 공극 내 압력 분포를 대변하며 부피는 존재하지 않는다고 가정하였으며, 튜브는 원통형으로 공극의 부피를 대변하며 압력은 존재하지 않는다고 가정하였다.

하나의 채널을 이루는 복셀들은 튜브 후보가 되며, 아래 Fig. 2(a)의 가장 왼쪽 그림과 같은 형상을 갖는다. 이 때 튜브의 길이는 양단의 노드와 노드 사이의 거리, 튜브의 반경은 튜브 후보를 중심으로 하는 내접구의 반경으로 정의하였다(Fig. 2(a)). 이는 다상 유체 거동 해석 시 공극목의 모세관압을 실제에 가깝게 반영하기 위함이다.

모세관압은 공극목에 작용하는 힘에 의해 결정되므로 튜브의 반경은 채널 내 존재하는 최소내접구의 반지름으로 정의하였다. 이를 통하여 해석에 사용된 도메인은 9094개의 노드와 10946개의 튜브를 갖는 관망으로 결정되었다(Fig. 2(b)).

2.3 다상 유체 거동 해석

SWCC실험을 모사하기 위하여 모든 튜브가 어떠한 유체로 포화 된 육면체 공극 네트워크에서 배수 방향을 결정하고, 그에 수직인 하나의 평면(침투 평면)에 다른 침투 유체의 압력을 점진적으로 증가시켜 다른 하나의 평면(배수 평면)으로 기존 유체를 배수시키는 침투 시뮬레이션을 실시하였다. 이 때 흐름방향에 수직인 침투 평면과 배수 평면은 배수경계조건으로, 평행한 나머지 평면은 불투수경계조건으로 설정하였다(Fig. 3(a)). 또한 모든 유체는 비압축성이며 공극 내 유속은 충분히 느려서 항상 층류(Laminar flow)의 흐름을 갖는 것으로 가정하였다.

공극 네트워크 추출 단계에서 자동으로 결정되는 튜브의 모세관압은 아래의 Young-Laplace식(식 (1))으로부터 계산된다.

 (1)

P_c,ij와 R_ij는 각각 i와 j번째 노드를 연결하는 튜브의 모세관압[Pa]과 반경[m]이며, T_s는 표면장력[N/m] 그리고 φ는 서로 다른 유체 사이의 접촉각[deg]을 나타낸다. 튜브 양단의 노드에 서로 다른 유체가 점유하고 있고, 양단의 압력 차이가 모세관압보다 클 경우 흐름이 발생한다. 이 때 침투를 허용한 튜브와 노드는 모두 침투 유체가 점유하게 되고 튜브 양단 노드의 압력은 침투유압으로 같아진다고 가정하였다. 양단의 압력차이가 모세관압보다 작고, 다른 경로로 해당 양단 노드에 침투 유체가 점유하고 있는 경우, 튜브는 고립되어 침투되지 않는다고 가정하였다(Fig. 3(b)).

이를 통해 공극 네트워크의 침투평면과 배수평면에 특정한 압력차를 가했을 때의 정상상태 침투 유체 점유 양상을 파악하였으며, 압력 차 ΔP에 따른 두 유체의 포화도를 계산하였다.

특정한 압력차에 의하여 침투 경로와 포화도가 결정되면, 침투 유체의 점유 부분과 기존 유체의 점유 부분을 서로 다른 두개의 관망으로 나누어 유동 해석을 실시하였다. 아래 Hagen-Poiseuille식(식 (2))과 노드에서의 연속방정식(식 (3))을 통하여 침투평면과 배수평면의 압력차에 의한 유량 및 내부압력을 계산하였다.

 (2)

(3)

q_ij는 i와 j번째 노드를 연결하는 튜브의 유량[m³/s], μ_ij는 튜브를 점유하고 있는 유체의 점성[Pa･s], R_ij와 L_ij는 각각 튜브의 반경[m]과 길이[m]를 나타내며, P_i와 P_j는 i와 j번째 노드의 유체 압력을 나타낸다. 위 두 식을 연립하여 행렬식 형태로 풀어내면 아래 식 (4)와 같다.

(4)

C는 각 노드에서의 유량의 합계 벡터, P는 각 노드의 압력 벡터, K는 유량과 압력과의 관계를 나타내는 전도 행렬이다. 이를 통하여 배수경로에 수직인 임의의 평면을 통과하는 유량을 계산하였고, 포화 상태의 유량과의 비율인 상대투수율을 계산하였다.

SWCC실험 모사를 위하여 사용한 두 유체의 물성치는 아래 표와 같으며, 각각 상온에서의 물과 공기에 가까운 값을 사용하였다. 두 유체 및 튜브의 재료에 따라서 달라지는 접촉각은, 완전 습윤 조건(Fully wetting condition)을 가정하여 0으로 설정하였다(Jang and Santamarina, 2011).

2.4 함수특성곡선 및 상대투수율의 수학적 모델

공극 네트워크 모델을 통한 SWCC를 평가할 경우 공기함입에 의해 결정되는 특정 포화도에서의 공극구조 산출이 가능하므로 투수계수의 직접적인 계산이 가능하다. 따라서 본 연구에서는 SWCC실험으로부터 모델에 기반한 상대투수율곡선을 공극 네트워크 모델로부터 직접 산출이 가능한 상대투수율곡선과 비교하였다.

SWCC실험에 의하여 모관흡수력에 따른 체적함수비의 변화를 나타내는 함수특성곡선을 예측하는 많은 수학적인 모델이 존재한다. 그 중 본 연구에서는 van Genuchten모델(van Genuchten, 1980)을 적용하였으며, 이는 다음과 같다.

(5)

ψ는 모관흡수력(Soil suction)[kPa]으로, 시뮬레이션에서의 침투평면과 배수평면의 압력차 ΔP와 동일하다. θ는 체적함수비이며, θ_s와 θ_r는 각각 포화 및 잔류체적함수비를 나타낸다. a_p, n_p는 공기함입치(Air entry value), 배수와 관련된 계수이다.

상대투수율은 모관흡수력에 따른 포화도에 관한 식으로 표현되며, 본 연구에 사용한 Corey모델(식 (6))(Corey, 1954)과 Brooks-Corey모델(식 (7))(Brooks and Corey, 1964)는 각각 아래와 같다.

 (6)

 (7)

k_rw는 모관흡수력에 따른 유체(물)의 상대투수율을 나타내며, S_w와 S_wr은 각각 유체(물)의 모관흡수력에 포화도와 잔류포화도를 나타낸다. λ는 공극도 분포와 관련된 계수로서 λ=2/3 일 경우, 두 모델은 동일하다. 또한, 상대투수율의 실험적 획득의 어려움으로 위의 식 (5)로부터 계산된 계수로부터 SWCC실험에 따른 상대투수율을 아래 식 (8)을 통하여 예측하였으며 이를 Mualem- van Genuchten모델이라 한다(van Genuchten, 1980).

 (8)

3. 결과 및 분석

건조과정(Drying path) SWCC실험을 모사한 공극 네트워크 모델 시뮬레이션을 서로 다른 여섯개의 침투 평면과 경계조건에 따라 실시하였다(Fig. 4). 이 때 시뮬레이션은 배수평면의 압력을 고정하고, 침투평면의 압력을 10Pa씩 증가시켜가며 수행되었다.

그 결과 주문진표준사는 침투평면의 면적과 배수거리가 일정할 때, 공기함입치 및 잔류체적함수비는 침투방향에 무관한 등방성을 나타내었다. 이 중 +x 방향의 침투 시뮬레이션 결과와 건조과정 SWCC실험 결과를 비교하였으며 이 때 물의 포화도 및 체적함수비에 따른 +x 방향 침투 양상은 아래와 같다(Fig. 5).

건조과정 SWCC실험(Yang, 2014)은 3차원 X-ray CT 영상과 동일한 시료로 수행되었으며 공극률과 상대밀도가 각각 0.4574와 10.00%로 계산되었다. 이를 공극 네트워크 모델의 시뮬레이션 결과와 비교하였으며 각각의 결과에 식 (5)를 적용하여 함수특성곡선을 획득하였다(Fig. 6 및 Table 2).

실험결과로부터 예측한 함수특성곡선과 공극 네트워크 모델로부터 예측한 함수특성곡선은 모두 건조과정의 함수특성곡선이므로 이력현상(Hysterisis)의 영향은 존재하지 않으며, 유사한 포화체적함수비와 잔류체적함수비를 나타냈다. 두 결과에 적용한 van Genuchten 모델의 계수 a_p는 각각 0.4453과 0.4657로 계산되었으며, 이를 통해 공극 네트워크 모델 시뮬레이션의 공기함입치가 실험 결과에 비하여 과대평가 되었음을 알 수 있다. 이는 최소내접구로 정의된 튜브의 반경이 실제보다 과소평가 된 것으로 판단된다. 계수 n_p는 각각 5.6807과 11.6826으로 계산되었으며, 이를 통해 시뮬레이션으로부터 얻은 함수특성곡선의 기울기가 실험 결과에 비하여 급하다는 것을 알 수 있다. 그 이유는 최소내접구로 계산된 전체 튜브 반경의 분포가 실제 분포보다 균일하게 계산되었기 때문으로 판단된다.

공극 네트워크 모델로부터 물의 상대투수율을 직접 계산한 결과로부터 Corey모델과 Brooks-Corey모델을 이용하여 상대투수율곡선을 획득하였다. 또한 이를 SWCC실험 결과로부터 도출된 van Genuchten모델의 계수를 이용하여 예측한 상대투수율곡선과 함께 비교하였다(Fig. 7).

시뮬레이션 결과를 이용한 Corey모델과 실험결과를 이용한 Mualem-van Genuchten모델은 매우 유사한 경향을 보였으며, 시뮬레이션을 결과를 이용한 Brooks-Corey모델은 그보다 낮은 상대투수율 변화를 보였다. 실제 공극 네트워크 모델 시뮬레이션 결과값은 포화도 0.8 이상일 경우, Brooks-Corey모델과 비슷하거나 낮게 나타났으며, 그보다 작은 포화도에서는 Brooks-Corey모델보다는 크게, Corey모델보다는 비슷하거나 작게 나타났다. 또한 시뮬레이션 결과값은 0.6 이상의 포화도에서 실험으로부터 도출된 상대투수율 곡선보다 낮게 나타났으며, 그보다 낮은 포화도에서는 실험으로부터 도출된 상대투수율 곡선과 유사하게 나타났다. 이는 공기함입에 따른 물의 유동 경로를 완벽히 차단하는 공극 네트워크 모델의 특성 때문으로 판단된다. 그러나 실험적으로 획득한 함수특성곡선에 Mualem-van Genuchten모델을 적용하여 획득한 상대투수율 역시 간접적인 방법이기 때문에 실험으로부터 도출된 상대투수율과 공극 네트워크 모델을 이용하여 계산한 상대투수율의 차이를 물리적으로 설명하는 데에는 한계가 존재한다고 사료된다.

4. 결 론

본 연구에서는 주문진표준사를 용기에 담아 3차원 X-ray CT 촬영하여 공극 영상을 얻은 후, 공극 네트워크를 추출하여 건조과정 SWCC실험을 모사한 다상 유체 거동 시뮬레이션을 실시하였으며 시뮬레이션 결과를 실제 건조과정 SWCC실험 결과와 비교하였다. 이로부터 도출된 결론은 다음과 같다.

(1)다지지 않은 주문진표준사 시편으로부터 동일한 조건 하에서 배수 방향에 관계없이 동일한 함수특성곡선을 획득할 수 있다.

(2)공극 네트워크 추출 시, 최소내접구의 반지름을 이용한 방법은 공극목의 단면을 과소평가하며, 그로 인해 시뮬레이션 결과는 실험 결과에 비하여 큰 공기함입치를 나타냈다. 또한 같은 이유로 높은 포화도에서는 상대투수율이 실험으로부터 도출 된 결과에 비하여 낮게 평가되었다.

(3)또한 최소내접구의 반지름을 이용한 방법은, 공극목의 단면적 분포를 실제보다 균일하게 평가하며, 그로 인해 시뮬레이션 결과는 실험 결과에 비하여 급한 함수특성곡선 기울기를 나타냈다.

(4)3차원 X-ray CT 영상과 공극 네트워크 모델을 통하여 조립토의 함수특성곡선과 상대투수율곡선을 획득 할 수 있으며, 이들은 실험결과로부터 도출된 결과와 유사한 양상을 보이고 시료 내부의 유체 점유 변화 양상을 확인 가능하다는 장점을 가지므로, 시간 및 환경적 제약이 따르는 실험의 대안으로써 활용될 수 있을 것으로 사료된다.

다상유체의 거동 해석은 매우 복잡하나 공극 네트워크 모델을 통해 이를 상대적으로 수월하게 수행할 수 있으며, X-ray CT의 발달이 그 기반을 제공하였다. 또한 위에서 언급한 내용을 토대로 공극목을 이상화한 튜브의 반지름을 결정하는 방법에 대한 추가적인 연구의 필요성을 확인하였다.