1. 서 론

입상체 사이에서의 열전달은 다공질 재료나 지반공학 등 많은 공학 분야에 다양하게 사용될 수 있다. 입자들 사이의 에너지 흐름에 대한 관계를 이용하여 매질의 유효열전도도(effective thermal conductivity)를 파악하고자 했던 노력은 꾸준히 지속되어 왔다. 대표적으로 Maxwell 모델(Maxwell’s approximation based models; Maxwell, 1873)과 Self-consistent(Landauer, 1952)모델 그리고 Different effective medium(Bruggeman’s asymmetric model (BAM), Bruggeman, 1935; Landauer, 1952) 모델이 가장 많이 사용되었다. 기존 모델들은 주요한 매질에 대한 정보를 이미 알고 있는 상태에서 다른 고체 입자가 추가될 때 예상되는 합성물질의 유효 열전도도를 예측하는데 사용되었다. 하지만 불포화 상태가 존재하는 경우에는 고체입자 이외에도 물과 공기 같은 유체가 각각 매질내에 존재한다. 따라서 고체입자 사이의 관계가 중요한 분야에 대해서는 고체입자들 사이의 열전달 관계를 먼저 정의한 다음 고체입자 주변에 존재하는 유체가 열전달에 영향을 줄 수 있는 부분을 고려하여 이를 전체 유효열전도도 계산에 적용하는 모델이 필요하다.

Feng et al.(2008)은 흙 입자를 원형판(circular disk)으로 가정하고 주요 열전달 과정은 전도에 의해 발생한다는 가정 하에서 BEM(Boundary Element Method)법을 사용한 DTEM(Discrete Thermal Element Model)을 제안하였다. 하지만 Feng et al.(2008)의 모델에서도 불포화 상태에 대한 고려는 어렵고 입자와 입자 사이의 접촉각이나 접촉면적 등을 적절히 가정해야 하는 어려움이 있다. Kanuparthi et al.(2008), Yun and Evans(2010)는 임의의 입자 분포에서 나타난 입자들의 접촉에 대해 입자의 열전도도와 유체의 열전도도를 고려한 네트워크 모델을 제시하고 비교적 균일한 입자크기를 갖는 글라스 비즈(Glass Beads)를 이용하여 실내시험과 수치해석 결과를 비교하였다. 입자와 입자 사이에 접촉이 발생했을 때 접촉각을 고려한 입자 사이의 열전도력(thermal conductance)과 입자가 적당한 거리에 놓여 있을 때를 고려한 열전도력을 각각 정의하여 입자 사이에 존재하는 유체로 인해 유발되는 입자간 열전달을 고려하였다. 네트워크 모델은 입자 사이의 접촉정도에 따라 Batchelor and O’Brien (1977)가 제안한 입자와 입자 사이의 열 전도현상에 대한 유도식을 적용한다. 또한 네트워크 모델은 입자와 입자 사이의 접촉정도와 공극을 채우는 재료의 열전도도에 기반하기 때문에 공극률과 유효응력에 따라 증가되는 입자와 입자간의 접촉면적과 그에 따른 열전도도 변화를 평가한다.

본 연구에서는 3차원 개별요소 수치 모델을 이용한 프로그램인 PFC3D를 이용하여 시료의 입도분포를 기반으로 대상 시료의 입자를 생성하고 입자간 접촉을 유도하기 위해 낮은 수준의 구속압을 적용하였으며 이후 네트워크 모델을 사용하여 시료의 유효 열전도도를 산정하였다. 시료의 불포화 특성을 고려하기 위해 압력판 추출 시험을 수행하여 시료의 흙-수분 특성 곡선을 획득하고 이를 수치 모델에 적용하였다. 또한, 포화도를 고려한 유효 열전달 실린더 반경 예측식을 네트워크 모델 해석을 통해 제시하였다.

2. 실내 열전도도 측정 및 흙-수분 특성 곡선

2.1 시험 시료

본 연구에서는 표준 시료로서 3가지의 평균 입경이 다른 글라스 비즈와 주문진사를 이용하여 실내시험을 수행하였다. 실내시험을 통해 시료의 기본물성을 획득한 후 포화도에 따른 열전도도를 측정하였다. 포화도에 따른 열전도도 측정은 임의의 포화도에 따라 성형 한 시료에 대해 압력판 추출시험을 수행한 후 급속 열전도도 측정기(QTM-500, Kyoto Electronics)을 사용하여 측정하였다(Choi et al., 2008). 또한 네트워크 모델을 이용한 불포화 열전도 모델에 적용하기 위해 압력판 추출시험을 통한 흙-수분 특성 곡선(SWCC, Soil-Water Characteristic Curve)을 도출하고 평균 모세관 반경(meniscus radius)을 산정 하였다. 본 연구에 사용된 글라스 비즈와 주문진사의 기본물성과 입도분포는 Table 1과 Fig.1에 각각 나타냈다. 또한, 기존 연구결과와 비교하기 위해 Song et al.(2010)이 제시한 주문진사에 대한 입도분포를 Fig. 1에 첨부하였다. 본 논문에서 사용한 주문진사가 Song et al.(2010)의 주문진사 보다 세립분 함유량이 다소 많은 것으로 나타났다.

2.2 압력판 추출시험과 열전도도 측정시험

본 연구에서는 SWCC를 획득하기 위해 압력판 추출시험(PPE, Pressure Plate Extractor)을 수행하였다. 시험에 적용된 흡입력의 범위는 5～100kPa로 ASTM D2325- 68(1997)에 의거하여 압력판 추출시험을 수행하였다.

시료는 Fig. 2와 같이 황동 성형틀을 이용하여 모든 시료의 건조단위중량을 1.61t/m3으로 고정시킨 상태에서 정다짐을 통해 성형하였다. 이후, 열전도도 측정에는 다소 시간이 소요되기 때문에 압력판 추출시험을 수행한 후 황동 성형틀에 랩을 씌워 시료의 함수비를 유지한 상태에서 급속 열전도도 측정기를 이용하여 시료의 열전도도를 측정하였다.

2.3 시험 결과 분석

Fig. 3은 3가지 크기의 글라스 비즈와 주문진사의 압력판 추출시험 결과를 나타낸다. 글라스비즈의 크기가 증가할수록 AEV(Air Entry Value)가 감소하는 것을 볼 수 있다. 본 연구에 사용된 글라스 비즈의 경우, 크기에 따라 AEV는 약 0.8～5kPa 범위로 나타났다.

Table 2는 VG 모델(Van Genuchten, 1980)에 사용한 계수와 AEV을 나타낸다. 실험결과를 통해 VG모델의 계수를 산정하고 이후 수치해석에서는 VG모델을 이용하여 포화도에 따른 모세관 반경을 도출하였다. Fig. 4와 같이 Song et al.(2010)은 주문진사에 대한 SWCC에서 AEV를 건조과정과 습윤과정에서 각각 약 1.8kPa와 1kPa로 보고하였다(Fig. 4a). 또한, Moon and Kim(2011)은 낙동강 모래에 대한 SWCC 시험에서 1～10kPa의 범위에서 체적함수비가 급격히 감소한다고 보고한 바 있다(Fig. 4b). 본 연구에서 사용된 주문진사의 경우, AEV가 1～2kPa사이에서 나타났다. 이는 본 논문의 산정한 주문진사의 SWCC가 기존 문헌들과 유사한 경향을 보인다고 할 수 있다. 다만, Fig. 3에 제시한 압력판 추출시험으로부터 산정한 주문진사의 SWCC의 기울기와 잔류포화도가 Song et al.(2010)의 결과보다 다소 큰 것은 두 연구에서 사용한 주문진사의 세립분 함유량 차이(Fig. 1 참조)에서 온 것으로 추정된다. 실험에서 산정된 SWCC는 네트워크 모델에서 불포화상태의 입자간 접촉을 고려하여 열전도력을 산정하는데에 이용되었다.

Fig. 5는 3가지 크기의 글라스 비즈와 주문진사의 포화도에 따른 열전도도 측정결과와 대표적인 사질토의 유효열전도도를 표현하는 Kersten의 경험식(Kersten, 1949)을 비교한 그림이다. 실험에서는 3가지 크기의 글라스 비즈와 주문진사를 사용하여 건조 단위중량(1.61t/m3)과 간극률(0.36)을 고정 시킨 후 다음 2가지 방법을 이용하여 시료의 열전도도를 측정하였다. 첫 번째는 포화도(0, 5, 10, 20, 40, 60, 80, 100%)에 따라 시료를 성형한 후 측정하는 방법(그림에서 Preconditioned saturation test)이고, 두 번째는 압력판 추출시험 직후 열전도도를 측정하는 방법(그림에서 Pressure Plate Extractor test) 이다. Large size, Medium size, Small size 글라스비즈의 열전도도는 포화상태에서 순서대로 0.94, 0.91, 0.87W/m･K이고 건조상태에서는 0.20W/m･K으로 동일하게 나타났다. 모든 경우에서 Kersten의 경험식이 측정값에 비해 상대적으로 더 큰 값을 보인다.

주문진사의 경우 포화상태의 열전도도는 2.36W/m･K이고 건조상태의 열전도도는 0.28W/m･K이다. 건조 상태의 열전도도는 글라스 비즈와 큰 차이가 없었으나 포화 상태에서는 주문진사의 열전도도가 2배 이상 크게 측정 되었다. 이는 글라스 비즈 입자의 열전도도(1.13 W/m･K)가 주문진사 입자의 열전도도(3.0W/m･K) 보다 상대적으로 작고 입자의 모양이 다르기 때문에 나타난 결과로 추정된다. 주문진사의 경우에도 글라스 비즈와 같이 Kersten의 경험식이 상대적으로 실내 열전도도 측정값보다 큰 경향을 보인다.

열전도도 측정결과, 임의의 포화도로 조절하여 성형한 시료로 측정한 열전도도(그림에서 Preconditioned saturation test)와 압력판 추출시험 직후 측정한 열전도도(그림에서 Pressure Plate Extractor test)가 다르게 나타났다. 이는 압력판 추출시험의 경우 입자 사이에 존재하는 물이 압력에 따라 입자간 거리가 큰 부분부터 빠져나가게 되므로 동일한 포화도라고 할지라도 임의로 성형한 시료의 입자내 물의 분포와 압력판 추출시험 직후 채취한 시료의 입자내 물의 분포가 서로 다를 수 있기 때문에 나타난 결과로 판단된다.

Fig. 6은 글라스 비즈 입자의 크기와 포화도에 따른 열전도도 측정 결과를 나타낸다. 본 연구의 경우 낮은 포화도 구간(0～20%)에서는 3가지 시료가 비슷한 열전도도를 보였으나 포화도가 증가할수록 Large size 글라스 비즈가 가장 큰 열전도도 값을 보이고 Small size 글라스 비즈가 가장 낮은 열전도도 값을 보였다.

각 크기별 글라스 비즈의 열전도도는 포화도가 증가할수록 열전도도가 증가하는 경향을 보인다. 이는 포화도가 증가할수록 간극사이에 존재하는 물의 부피비가 증가하여 입자 사이의 열 전달이 더 원할해지기 때문이다. 또한 글라스 비즈의 크기가 증가함에 따라 동일한 포화도에서 입자의 크기가 큰 시료의 열전도도가 상대적으로 큰 값을 보였다.

3. 네트워크 모델을 이용한 열전달 해석

3.1 개별 요소법을 이용한 입자 생성

최근 개별요소법(DEM, Discrete Element Method)을 이용하여 지반이나 조립재료 등의 열특성 예측을 위한 연구가 진행된 바 있다(Kanuparthi et al., 2008; Yun and Evans, 2010). 특히, 모래와 같은 조립재료를 대상으로 한 지반공학 분야에 개별요소법 적용이 증가하고 있다. 개별요소법을 통한 해석은 재료의 특성 중 입도분포와 간극률을 고려하여 실제 현상에 가까운 해석이 가능하다는 장점을 갖고 있다.

PFC(Particle Flow Code, ITASCA)는 2차원(disk)과 3차원(sphere) 해석이 가능하다. 3차원 프로그램인 PFC3D는 2차원 기반의 PFC2D와 기본개념은 유사하나 2차원 개별요소 형태가 원판(disk)형태인데 반해 구(sphere)의 형태를 띠고 있기 때문에 본 연구에서는 실제 재료의 간극률 조건에 맞게 모델링이 가능한 PFC3D를 사용하였다. Fig. 7은 개별요소법에서 생성한 입경에 따른 분포를 보여준다. 그림에서는 입경별로 다른 색으로 표현하였다. 생성된 입자의 위치는 PFC프로그램 내부의 랜덤함수를 이용하여 구성하였고 입자의 크기는 각 시료의 입도분포를 기준으로 생성하였다.

구 형태의 3차원 입자로 생성된 공간은 해석 모델의 외곽 경계면 형태에 따라 바뀌게 되는데 본 연구에서는 Fig. 8과 같이 원통형 경계면을 사용하여 모델링 하였다. 경계면 구성은 2개의 상하부 사각형판 경계면과 1개의 원통형 경계면을 이용하였다. 실험에 사용한 글라스 비즈의 크기가 3가지(large, medium, small size) 이고 각각의 크기가 다르므로 Table 3과 같이 원통 경계면의 크기를 다르게 하여 입자의 개수가 2500개 내외가 되도록 하였다. 개별요소법은 생성되는 입자들의 관계를 이용하여 매질의 특성을 모사하기 때문에, 응력-변형 모사의 경우 전체 모델의 크기와 입자 크기의 비가 유사한 경우 비슷한 해석 결과를 얻게 된다. 하지만, 본 연구에서는 개별요소법으로 생성한 입자의 개수가 유사하고 전체 모델의 크기와 입자 크기의 비가 비슷하더라도 각 모델에서 생성된 입자의 크기와 분포가 다르고 네트워크 모델에서는 입자간 열전도력 산정 시 입자간 거리와 접촉면적이 고려되기 때문에 각각의 모델이 다른 특성을 보인다고 할 수 있다.

접촉상태의 개별요소 모델을 생성하기 위해 사용된 모델은 Hertz-Mindlin 접촉모델(Hertz, 1895; Mindlin and Deresiewicz, 1953)이다. Hertz-Mindlin 모델은 두 물체가 접촉하여 힘을 받으면 접촉부분이 변형하여 접촉면이 생기는 동시에 접촉압력을 발생시킨다. 이 접촉부분은 물체에 대해 지극히 작은 경우를 말한다(Hertz, 1895). 실내시험에서 사용한 간극률과 동일하게 개별요소법에서 입자를 생성한 후 중력과 5kPa 내외의 낮은 수준의 구속압을 적용하여 입자간 접촉을 유지하도록 하였다. PFC에서는 경계면을 통한 직접적인 구속압을 적용시키는 것이 불가능하기 때문에 경계면을 특정 속도로 원통의 중심방향으로 이동시키고 그 과정에서 경계면에 작용하는 반력을 측정하는 방식으로 구속압을 적용시킨다(PFC3D, 2006). 따라서 해석과정 중 일정한 구속압 조건을 유지시키기 위해 측부 경계면은 반력과 경계면 속도의 자동조절을 반복하게 하는 Servo- mechanism을 적용하였다. Fig. 8은 구속압이 작용된 상태에서 입자간 접촉상태를 보여준다.

개별요소법에 의해 구성된 입자간의 접촉관계를 이용하여 네트워크 모델에 적용한 해석 결과는 Fig. 9와 같은 온도분포를 나타낸다. 본 연구에서는 경계면 온도를 상부에서 5K, 하부에서 1K를 적용하였다. 그림에서는 온도구간을 7단계로 나눠 표시하였다. 그림에 나타난 바와 같이 열전달 네트워크 모델은 온도 경계조건을 잘 반영할 수 있는 것으로 나타난다.

3.2 입자간 열전도력 산정

네트워크 모델을 적용하기 위해서는 특정 재료를 대표하는 3차원 입자구조 시스템 내에서의 열흐름을 수치적으로 표현해야 한다. 입자구조 내 열흐름의 가장 중요한 수단은 입자 또는 입자 사이의 접촉에서의 전도(conductance) 이다. 이때, 입자 사이의 공극을 채우는 유체(fluid)의 종류(공기, 물)에 따라 열전달이 좌우된다. 자연상태에서 건조조건은 네트워크 모델에서 공극을 채우는 유체의 종류가 공기인 경우이고, 자연상태에서 포화조건은 유체의 종류가 물인 경우이다. 본 연구에서 사용된 네트워크 모델에서는 포화상태에서 유체가 균등하게 분포하는 것으로 가정하였고 불포화 상태에서는 함수비와 SWCC에 나타난 흡입력에 따라 입자와 입자 사이에 균등하게 분포하는 것으로 가정하였다.

입자의 유효 열전도력은 두 입자 자체의 열전도력과 입자 사이의 열전도력을 직렬 형태로 합하여 계산할 수 있다(Kanuparthi et al., 2008; Yun and Evans, 2010). 입자의 접촉은 입자 사이에 겹침(overlap)이 있는 경우와 분리(gap)되어 있는 경우로 나눌 수 있다. 본 연구에서는 Yun and Evans(2010)에 의해 제시된 열전도도 산출을 위한 열 네트워크 모델이 사용되었다. 즉, 입자간 겹침이 있는 경우(입자간 중심거리가 반지름의 합보다 작은 경우) 이를 겹침 조건으로 간주하여 열전도력을 계산하고, 입자가 분리되어 있는 경우 문헌(Yun and Evans, 2010)에서 제시한 글라스 비즈의 유효거리 계수()인 0.5를 사용하여 계산하였다.

3.3 네트워크 모델을 이용한 입자간 열전달 해석 알고리즘

Fig. 10은 본 논문에서 적용한 해석 알고리즘을 보여준다. 수립된 알고리즘은 Fortran90 프로그램을 사용하여 작성되었다. 다음은 해석순서를 요약한 것이다.

(1)개별요소법에서 생성된 각 입자의 3차원 위치와 반지름을 획득한다.

(2)각 입자에 대해 인접하는 입자가 접촉하는 조건인지 또는 유효 거리 내에 존재하는 입자인지 판별하고 해당 조건에서 열전도력을 계산한다. 이는 입자의 위치를 고려하여 입자간 거리로 판별한다.

(3)계산된 열전도력을 이용하여 입자간 접촉조건과 열전도력을 정리하고 경계조건을 고려하여 각 입자의 온도를 계산한다.

(4)계산된 각입자의 온도 변화의 합이 허용 오차범위 내에 만족할 때까지 계산을 반복하고 각 입자의 온도를 얻는다(Gauss-Seidel 반복법 이용).

(5)정리된 각 입자의 온도 중 경계면과 접촉 조건 상태에 있는 입자의 온도와 열전도력을 고려하여 최종적으로 매질의 열전도도를 계산한다.

불포화 조건일 경우, 포화조건에서 계산된 유효실린더 반경과 입자간 유효거리를 AEV에서의 meniscus 반지름으로 보정하는 과정을 거친다. 포화상태의 네트워크 모델에서와 마찬가지로 접촉 조건인지 유효 거리 내에 존재하는 조건인지 판별하는 과정을 거치지만 열전도력을 계산하는 과정에서는 불포화 조건을 고려한 유효 실린더 반경과 유효거리가 적용되어 계산된다. 이후 각 입자의 온도를 계산하고 반복하는 과정은 포화상태와 같다.

3.4 불포화 상태를 고려한 네트워크 모델

3.4.1 포화상태의 로 정규화된 불포화 열전달 해석

본 연구에서 사용된 글라스 비즈의 포화도에 따른 열전도도를 기존 네트워크 모델(Yun and Evans, 2010; Kanuparthi et al., 2008)을 사용하여 각각의 포화도에 따른 값의 변화를 Fig. 11에 나타냈다. 즉, 입자의 열전도력(Cg)과 입자간 열전도력(Ccontact)의 접촉면적은 값에 따라 동시에 감소하거나 증가하는 경우이다. Fig. 11은 기존 네트워크 모델에서 유체(fluid)는 물을 적용하고 실내실험에서 측정한 포화도에 따른 열전도도에서 획득한 값이다. 그림에서 볼 수 있듯이 글라스 비즈의 크기에 따라 값이 변화 하는데, 건조상태에서의 값은 동일하고 포화도가 증가할수록, 글라스 비즈 크기가 증가할수록 값도 커짐을 알 수 있다.

Fig. 12에서 보이는 바와 같이 기존 네트워크 모델을 이용한 해석 시 건조상태에서는 동일한 값을 적용하였을 때 글라스 비즈와 주문진사의 열전도도가 같게 산정되었으나 포화상태일 경우에는 동일한 값을 적용하였을 때 실내실험으로 측정한 열전도도보다 네트워크 모델을 통해 예측한 열전도도 값이 더 큰 값을 보인다.

따라서 본 연구에서는 네트워크 모델을 사용하여 산정한 포화상태의 (sat)를 사용하여 각각의 포화도에 따라 네트워크 모델에서 산정한 를 정규화하고 이를 Fig. 13에 나타냈다. 정규화된 값()은 서로 유사한 경향을 보였으며 그림에 나타난 바와 같이 타원의 방정식을 변형하여 다음 식 (1)과 같이 표현할 수 있다. 여기서 는 포화도(%)를 나타낸다.

                                           (1)

3.4.2 SWCC를 이용한 불포화 열전달 해석

자연 상태에서 입자간 간격이 존재할 때 간극에 존재하는 물의 형태를 다음 그림과 같이 meniscus 반지름을 사용하여 표현할 수 있다. 본 연구에서는 다음 4가지 경우에 대해 접촉에 대한 유효 열전달 실린더의 반경을 고려하였다(Fig. 14). 즉, 입자의 크기 같고 적당히 떨어져 있는 경우(a), 입자의 크기가 다르고 적당히 떨어져 있는 경우(b), 입자의 크기가 같고 접촉되어 있는 경우(c) 그리고 입자의 크기가 다르고 접촉되어 있는 경우(d)이다.

meniscus 반지름(rp)을 고려한 유효 실린더 반지름()은 입자의 크기가 같고 두 입자가 일정 간격을 유지하는 경우에 대해 다음 식 (2)로 표현된다.

        (2)

여기서, 와 는 각각 i와 j번째 입자의 반지름(m), 는 i번째 입자와 j번째 입자의 떨어진 거리 또는 간격(m)으로 Dij - (Ri + Rj)로 표현된다. Dij는 입자의 중심간 거리(m)이다. 입자의 크기가 다르고 적당히 떨어져 있는 경우는 다음 식 (3)으로 표현된다.

            (3)

여기서, . 입자의 크기가 같고 접촉되어 있는 경우는 다음 식 (4)로 표현된다.

      (4)

여기서, 는 i번째 입자와 j번째 입자의 입자간 겹친거리(m)이다. 입자의 크기가 다르고 접촉되어 있는 경우는 다음 식 (5)로 표현된다.

         (5)

여기서 는 meniscus 반지름(m)으로 모든 경우에서 다음 식 (6)을 사용하여 산정한다.

                (6)

여기서, 는 계면장력(N/m), 는 모세관압(kPa), 는 계면각 또는 접촉각이다. 본 연구에서는 를 0으로 가정하고(cos=1) 모세관 압은 압력판 추출시험을 통해 산정한 SWCC곡선 상의 포화도에 따른 모세관압을 사용한다. 위 4가지 방법을 통해 계산된 유효 실린더 반경은 접촉면에 해당하는 것으로 각각의 입자에 적용되는 유효 실린더 반경은 포화상태에서 보정한 값을 그대로 사용하였다. 즉, 포화상태에서는 기존 네트워크 모델(Yun and Evans, 2010)에서 사용한 보정계수 를 입자의 열전도력을 고려할 때에는 변하지 않는 것으로 간주하였다. 따라서 포화상태에서 고려한 입자의 열전달 반경은 모든 포화도 구간에서 일정한 것으로 간주하였다. 그러므로 기존 Yun and Evans(2010)이나 Kanuparthi et al.(2008)에서 고려한 입자간 접촉에서의 열전도력은 포화도를 고려하여 다음 식 (7)과 식 (8)로 바꾸어 표현할 수 있다.

                                           (7)

 (8)

여기서, 는 입자간 간격이 존재할 때의 열전도력(W/m), 는 입자간 겹침이 존재할 때의 열전도력(W/m), kf는 유체의 열전도도(W/mK), 는 kg/kf, kg는 입자의 열전도도(W/mK), 는 포화도를 고려한 입자간 접촉 거리(m), Kc와는 네트워크 모델에서 산정된 계수(Yun and Evans, 2010), l는 입자의 간격에 대한 계수(Yun and Evans, 2010)이다. 따라서 포화도가 감소 할수록 유효 실린더 반경은 다음 Fig. 15와 같이 감소하게 된다.

본 연구에서는 실내실험을 통해 각 시료의 SWCC를 획득하고 입자간 접촉면에서의 유효 접촉면적을 각 포화도에 해당하는 모세관 흡입력을 사용하여 산정하였다. 이를 위해 불포화 상태를 고려한 네트워크 모델에서는 AEV에 해당하는 상태를 포화상태로 가정하고 기존 네트워크 모델에서 산정한 포화상태의 열전도도 값을 갖도록 보정하는 과정이 필요하다. 즉, AEV에서 계산한 유효 실린더 반경과 유효 거리가 포화상태에서 고려한 유효 실린더 반경과 유효 거리와 같도록 접촉면의 관계를 보정한 후 이를 이용하여 포화도에 따른 접촉면의 유효면적을 다시 계산하는 방법을 사용하였다. 이는 개별요소법으로 구성한 입자의 구성 모델이 실제 상태가 아니므로 이를 보정하기 위해 포화상태에 대한 네트워크 모델의 유효 실런더 반경과 실험결과 나타난 AEV의 모세관 흡입력을 meniscus 반경으로 환산하여 식 (2)～(5)에서 계산한 유효 실린더 반경과 서로 일치시키도록 한 것이다. 그러므로 AEV에서의 열전도도는 기존 네트워크 모델에서 계산한 매질의 포화상태의 열전도도와 동일한 값을 갖게 되고 모세관압이 증가할수록(포화도가 감소할수록) meniscus 반경이 감소하기 때문에 입자간 접촉면적과 계산에 고려되는 유효거리는 감소하게 된다. 또한, 기존 네트워크 모델에서, 입자간 접촉이 없으나 열전도 영향에 대한 유효거리 계산에 사용되는 계수 은 meniscus 직경에 따라 감소하도록 하되 유효 실린더 반경에서와 마찬가지로 AEV 상태에서 계산된 값과 기존 네트워크 모델의 포화상태로 사용된 값으로 보정하는 과정을 거친다. 따라서 포화도가 감소할수록 입자간 영향거리는 감소하게 되므로 계산에 고려되는(입자간 열전달이 발생하는) 입자의 개수가 SWCC의 형태에 따라 감소하게 되고 포화도가 감소함에 따라 전체 계산횟수도 감소하게 된다. 먼저 유효 실린더 반경에 관련된 는 AEV에서 계산된 와 기존 네트워크 모델에서 계산한 를 사용하여 다음 식 (9)에서 보정계수()를 산정하였다.

              (9)

계산된 보정계수는 AEV에 해당하는 모세관압 보다 큰 구간(S<100%)에서 모두 적용되었다. 그러므로 SWCC를 적용하여 계산한 유효 실린더 반경이 포화상태를 모사한 네트워크 모델에서의 유효 실린더 반경보다 클 경우 모든 구간에서 계산한 유효 실린더 반경은 보정계수를 통해 감소되어 적용되고 반대의 경우 계산한 유효 실린더 반경이 증가되어 계산에 고려되었다. 다음으로 입자간 간격에 고려되는 계수인 의 경우는 포화상태를 모사한 네트워크 모델에서 AEV에 해당하는 meniscus 반지름을 고려하여 다음 식 (10)에서 보정계수()를 산정하였다.

             (10)

하지만 입자간 거리가 meniscus 직경()과 같은 경우는 계산에 고려하기 어렵기 때문에 입자간 거리가 meniscus 직경의 80%이하일 때 계산에 고려하였다. 즉, 입자간 거리가 meniscus 직경과 같을 경우 식 (2)와 (3)에서 계산된 유효 실린더 반경이 0이 되고 meniscus 직경보다 입자간 간격이 작더라도 간격이 상대적으로 작은 경우 입자간 유효 실린더 반경이 매우 작아지기 때문에 80% 이하의 간격에서만 계산에 고려하고 그 이상일 때에는 입자간 영향이 없는 것으로 간주하였다. 따라서, meniscus 반지름은 모세관압이 증가함(포화도 감소)에 따라 감소하게 되고 입자간 영향을 고려하는 거리(･)도 SWCC의 형태에 따라 함께 감소한다.

본 연구에서 사용하는 불포화를 고려한 해석 알고리즘은 입자간 유효거리가 포화도가 감소함에 따라 함께 감소하기 때문에, 건조과정에서 발생하는 입자 사이의 잉크병 현상을 고려하게 된다. 즉, 포화도가 감소할수록 meniscus 직경이 감소하게 되고 주어진 알고리즘에 의해 입자간 영향을 주는 거리가 계산되므로 자연적으로 입자간 거리가 먼 경우에서부터 계산에서 제외되게 된다. 또한, 시료의 SWCC 곡선상에 나타난 AEV에 기반하여 해석 초기에 포화상태의 열전도도와 비교하여 유효 실린더 반경을 수정하게 되므로 대상 시료의 SWCC 기울기가 급하고 모세관 흡입력의 차가 심할수록 해석 결과로 나타나는 함수비에 따른 열전도도의 차가 커지게 된다.

4. 네트워크 모델을 이용한 열전달 해석 결과

4.1 포화상태의 로 정규화된 불포화 열전달 해석

식 (1)을 이용하여 글라스 비즈의 포화도에 따른 시료의 열전도도를 예측하고 실내실험에서 측정한 열전도도와 비교하여 Fig. 16에 나타냈다. 비교 결과에서 나타난 바와 같이 정규화를 이용한 열전도도 예측식은 열전도도 측정결과와 매우 유사한 경향을 보인다. 건조상태의 경우, 네트워크 모델에서는 매질에 채워진 유체가 공기로 가정되기 때문에 포화상태의 경우와 조건이 다르게 구성된다. 따라서 정규화를 이용한 비교에서는 건조상태의 결과는 제외하고 표시하였다. 또한, 예측식을 사용하여 주문진사에 적용한 결과를 Fig. 17에 나타냈다. 주문진사의 경우에서도 정규화를 이용한 예측식이 실내 측정결과를 잘 반영할 수 있는 것으로 나타났다. 포화상태에 대해 네트워크 모델을 사용하여 sat를 산정하고 본 연구에서 제시하는 예측식을 적용하여 포화도에 따른 를 계산한 후 조립토의 포화도에 따른 열전도도를 산정할 수 있을 것으로 판단된다.

Fig. 18은 포화도에 따른 실내 열전도도 측정결과, 예측식을 이용한 주문진사의 열전도도 해석 결과, 네트워크 모델을 사용한 예측결과와 Kersten 경험식을 통한 열전도도를 비교하여 나타냈다. 본 연구에서 제시한 예측식을 주문진사에 적용한 결과에서 볼 수 있듯이 불포화 구간의 열전도도 경향을 비교적 유사하게 예측할 수 있는 것으로 나타났다. Kersten 모델의 경우 글라스 비즈의 경우와 마찬가지로 실내 측정 열전도도 측정결과보다 다소 과대평가 되는 경향이 나타났다.

4.2 SWCC를 이용한 불포화 열전달 모델과 실내시험 결과 비교

Fig. 19는 3가지 크기의 글라스 비즈와 주문진사의 실내 열전도도 측정결과와 네트워크 모델의 해석결과를 비교한 그림이다. Large size 글라스 비즈에 대한 해석 결과는 포화도가 낮은구간(5～20%)에서 실내 열전도도 측정결과에 비해 상대적으로 크게 나타났으나 포화도가 40% 이상인 구간에서는 실내 열전도도 측정결과와 유사한 경향을 보였다. Medium size 글라스 비즈에 대한 해석 결과는 포화도가 5～10%에서 실내 열전도도 측정결과에 비해 상대적으로 크게 나타났으나 포화도가 20% 이상인 구간에서는 실내 열전도도 측정결과와 유사한 경향을 보였다. Small size 글라스 비즈에 대한 해석 결과는 포화도가 5～20%에서 실내 열전도도 측정결과에 비해 상대적으로 크게 나타났으나 포화도가 40% 이상인 구간에서는 실내 열전도도 측정결과와 유사한 경향을 보였다. 주문진사에 대한 수치해석 결과는 포화도가 5～20%에서 실내 열전도도 측정결과에 비해 상대적으로 크게 나타났으나 포화도가 40% 이상인 구간에서는 글라스 비즈의 경우와 같이 실내 열전도도 측정결과와 유사한 경향을 보였다.

본 연구에서 제시하는 수정된 네트워크 모델로 산정한 유효 열전도도와 실내실험을 통해 측정한 시료의 열전도도를 비교한 결과, 수정된 네트워크 모델의 경우, 포화도가 낮은구간(0～20%) 에서 유효열전도도가 실험결과에 비해 다소 과대평가되는 경향을 보였으나 이후에 구간에서는 실내실험 측정결과와 유사한 경향을 보였다. Kersten 모델의 경우 모든 시료에 대해서 실내실험 결과에 비해 다소 과대평과 되는 경향을 보인다.

Yun and Evans(2010)은 글라스 비즈와 주문진사의 열전도도를 네트워크 모델을 통해 모사한 결과 포화상태에서 값은 0.8, 건조상태에서 0.5 라고 보고하였다. 또한 이런 경향에 대해 시료가 포화된 경우, 입자와 간극수(물)의 열전도도 차이가 건조상태(공기)에 비해 줄어들기 때문에 입자의 열전도력(Cg)이 포화상태에서 더 큰 역할을 한다고 설명하였다.

본 연구에서는 포화도를 고려하기 위해 AEV 상태에서 기존 네트워크 모델의 포화상태 유효 열전도도와 보정하는 과정을 거쳐 입자 자체의 유효 반경은 고정하고 해석을 수행하였다. 따라서 포화도가 낮아지더라도 입자의 열전도력은 포화상태와 같아지고 다만 입자간 접촉(contact)에서의 열전도력만 감소하기 때문에 건조상태와 가까워 질수록 유효 열전도도가 다소 과대평가되는 경향을 보이는 것으로 판단된다.

5. 결 론

본 연구에서는 불포화토의 열전도도를 수치해석을 통해 모사하기 위해 네트워크 모델과 흙의 SWCC를 이용하는 방법에 대해 연구하였다. 본 연구를 통해 도출한 결론은 다음과 같다.

(1)3가지 크기의 글라스 비즈와 주문진사의 포화도에 따른 열전도도 측정결과, 포화도가 증가할수록 시료의 열전도도가 증가하는 경향을 보였다. 포화상태에서 열전도도는 글라스 비즈 입자의 크기 순서대로 0.94, 0.91, 0.87W/m･K이고, 건조상태 에서는 0.2W/m･K로 동일한 값을 보였다. 주문진사의 실내 열전도도 측정 결과, 포화 상태에서 2.36W/m･K, 건조상태에서 0.28W/m･K 으로 나타났다.

(2)개별요소법을 사용하여 입도분포에 따른 시료를 수치해석으로 모사할 수 있었다. 개별요소법을 사용하여 입자를 구성한 후 열전도도를 모사하는 방법을 사용할 때에는 Hertz-Mindlin 모델을 사용하여 입자간 접촉을 유지시키는 방법이 유효하였다. 또한 구성된 입자모델에 구속압을 가하여 입자 사이의 겹침(overlap)을 유도할 수 있었다. 본 연구에서는 낮은 수준의 구속압(5kPa 내외)을 가하여 입자사이의 겹침을 생성하였다 (normalized contact radius ≈0.006).

(3)네트워크 모델을 이용한 불포화 열전도도 예측식을 사용하여 포화도에 따른 시료의 열전도도를 산정한 결과 실내시험 결과와 유사한 경향을 보였다. 포화상태의 sat값은 기존 네트워크 모델을 사용하여 산정하고 본 연구에서 제시하는 예측식을 적용하여 포화도에 따른 를 산정할 수 있다. 제시한 예측식을 글라스 비즈와 주문진사에 적용한 결과, 실내시험 결과와 유사한 경향을 얻을 수 있었다. 따라서, 포화도에 따라 정규화된 값을 사용하여 네트워크 모델에 적용하여 임의의 포화도에서의 열전도도를 예측할 수 있을 것으로 판단된다.

(4)본 연구에서 제시하는 수정된 네트워크 모델로 산정한 유효 열전도도와 실내실험을 통해 측정한 시료의 열전도도를 비교한 결과, 수정된 네트워크 모델의 경우, 포화도가 낮은구간(0～20%) 에서 유효열전도도가 실험결과에 비해 다소 과대평가되는 경향을 보였으나 이후에 구간에서는 실내실험 측정결과와 유사한 경향을 보였다.

(5)본 연구에서는 포화도를 고려하기 위해 네트워크 모델 상에서 AEV 상태의 열전도도를 포화상태의 유효 열전도도와 보정하는 과정을 거쳐 입자(grain)의 유효 열전달 반경을 고정하고 입자의 접촉면에서 유효 열전달 반경을 수정하여 해석을 수행하였다. 따라서 포화도가 낮아지더라도 입자 자체의 열전도력은 포화상태와 같고 다만 입자사이의 접촉면에서 열전도력(Ccontact)만 감소하기 때문에 건조상태에 가까워 질수록 유효열전도도가 다소 과대평가되는 경향을 보이는 것으로 판단된다.