1. 서 론
2. 열-수리 파이프 유한요소 및 역해석 개발
2.1 순환 파이프에서의 열전달
2.2 열-수리 파이프 요소의 유한요소 수식화
2.3 유체 대류에 의한 인접한 연속체의 열흐름 방정식의 변경
2.4 비선형 역해석
3. 수치 예제 해석
3.1 파이프 요소의 일차원 응답해석
3.2 열응답 실험 수치해석
3.3 열응답 실험에 대한 역해석
4. 결 론
1. 서 론
지중 열교환 시스템(ground heat exchange system)은 친환경적 재생 에너지원으로 지속적으로 에너지 효율이 개선되고 있다. 열교환 시스템에서 중요한 열적 상호작용(thermal interaction)은 열교환 파이프(heat exchange pipe)와 주위지반과의 접촉면에서 이루어지게 된다. 이러한 접촉면의 열적 상호작용은 양방향으로 가능한데, 우선 파이프 내의 유체흐름에 의한 강제 대류(forced convection)는 주위지반의 역학적 변형이나 수리학적 특성의 변화를 초래한다. 역으로, 다공질 재료인 주위지반에서 지하수 흐름은 열교환 시스템의 주위 지반과의 열교환율(heat exchange rate)을 변화시킬 수 있다. 열 교환 파이프와 주위 지반과의 이러한 열적 상호작용은 지중 열교환기의 장-단기적 거동특징에 지대한 영향을 미치게 된다.
지중 열교환 시스템에 대한 수치해석은 시간경과에 따른 열성능의 변화, 주위지반의 열적 변화, 지반 유효 열물성치의 추정, 또는 열성능 평가 등 다양한 열-유체 흐름과 변형에 관하여 수행되고 있다. 국내의 경우 대부분 상용프로그램을 이용하여 열적 해석 혹은 열-역학적 해석을 수행하였으며, 파이프 내의 유체흐름은 CFD 기법에 의하여 열-수리해석을 수행하였다(Choi, 2011; Lee et al., 2013; Park et al., 2012; Woo et al., 2007).
지중 열교환에 대한 수치해석에서 가장 중요한 부분은 파이프내의 유체 흐름과 이에 접한 파이프 그리고 주위재료에 대한 수치해석 모델링이다. 특히, 유체가 순환하는 파이프를 적절히 모사하는 것은 매우 중요하다. 하지만, 난류 흐름의 열-수리 현상에 대한 이론적 불확실성(Wilcox, 2006)과 안정된 해석결과를 얻기 위해 요구되는 요소수의 급증으로 인한 과다한 컴퓨터 연산 시간(Marcotte and Pasquier, 2008)의 어려움이 있다. 그리고 유체가 순환하는 파이프가 주위 해석영역에 비하여 기하학적으로 큰 세장비를 갖는 것은 유체흐름에 대한 직접적인 수치해석을 더욱 어렵게 한다(Al-Khoury and Bonnier, 2006). 따라서 지금까지의 많은 수치해석들이 깊이 방향의 열적 유동을 무시한 2차원 수치해석을 수행하거나(Yavuzturk et al., 1999), 순환수 흐름에 대한 모델링의 어려움으로 인하여 상용프로그램에서 순환수의 유입-유출 온도를 사전에 지정하여 해석모델에 적용하기도 한다(Schiavi, 2009). 이러한 문제점을 개선할 수 있는 보다 효율적이고 개선된 수치해석 접근법이 필요하다(Diersch et al., 2011).
따라서 본 연구에서는 원형 단면에 대한 파이프내 난류의 열-유체 흐름에 대하여 직접적인 수치해석을 수행하지 않고, 파이프를 일차원 요소로 가정하고 실험결과에서 얻어진 대류 관계식으로부터 주위재료과 유체의 상호 열교환을 모사하였다. 그리고 개발된 일차원 열-수리 파이프 요소를 다공질 재료의 THM(Thermo-Hydro-Mechanical) 현상에 관한 유한요소 프로그램에 결합하여, 지반의 열전도도를 평가하기 위해 주로 사용하는 TRT(Thermal Response Test)를 수행하고 효과적인 유효 열물성치를 평가하기 위한 역해석의 적용성에 대하여 논의하였다.
2. 열-수리 파이프 유한요소 및 역해석 개발
2.1 순환 파이프에서의 열전달
지중 열교환 시스템에서 열전달은 전도(conduction)와 대류(convection)에 의하여 이루어 진다. 전도는 연속체 재료들(예, 다공질 지반, 파이프, 그라우트재 등) 사이에서 가장 중요한 열전달 수단이다. 반면 대류에 의한 열전달은 파이프 내부의 순환 유체에 의하여 발생한다. 지중 열교환기에서 전도와 대류의 상대적인 중요성을 평가하기 위하여 다음의 무차원 계수를 사용한다(Table 1).
파이프내의 유동은 유체가 벽면에 의해 둘러 싸여진 내부유동이다. 이러한 내부 유체의 열흐름은 유체의 이송(advection)과 전도에 의하여 발생하며, 고정된 자유유동온도가 없고 단면내의 온도 분포가 일정하지 않기 때문에 유체의 평균온도를 사용한다. 그리고 위치에 따른 평균온도의 변화는 에너지 평형을 적용하여 산정할 수 있다. 주어진 속도로 이동하는 관내의 유체와 관을 둘러싼 주위재료 사이에서 발생하는 대류에 의한 열전달률은 일반적으로 Newton의 냉각법칙으로 간단히 표현된다.
(1)
여기서, 
는 관내 유체의 온도, 
는 관을 둘러싼 주위재료의 온도이다. 관내 유체와 벽체 사이의 열교환량을 표현하는 상수를 대류 열전달 계수(convective heat transfer coefficient, 
)라고 하며 원형단면에서 
의 관계가 성립한다.
본 연구에서는 순환 유체와 관벽사이의 열대류를 평가하기 위한 Nusselt number(
)는 Dittus-Boelter식을 사용하였다.
(heating,
)
(cooling,
) (2)
2.2 열-수리 파이프 요소의 유한요소 수식화
열교환 파이프에서 순환하는 유체의 열-수리적 거동을 모사하기 위하여 일차원 유한요소를 개발하였다. 외부 열유입이 일정한 파이프 내부의 온도분포는 수리학적 완전히 발달된 영역에서 온도분포가 일정하게 유지된다(Cengel and Ghajar, 2012). 그리고 전체해석 영역에 비하여 파이프와 내부유체의 세장비가 극히 크므로, 파이프의 단면방향의 온도분포는 동일하다고 가정하였다. 따라서 개발된 파이프 요소는 길이방향으로 온도의 변화를 모사할 수 있다.
파이프내의 유체흐름과 관련되는 열에너지의 흐름은 유체의 전도(
), 변형과 유체의 이송(advection)에 의한 energy flux(
), 파이프내의 유체 흐름에 의한 대류(
), 그리고 내부 source(
, 예, TRT 시험)이 있다. 이들에 대한 주어진 시간 
동안 에너지 보존의 법칙을 적용하면 다음과 같다.
(3)
여기서, 파이프의 단면적은 
이고 주면장 
이다. 내부 유체의 물성치는 밀도(density, 
), 비열(heat capacity, 
), 내부에너지(internal energy, 
), 그리고 유속(fluid velocity, 
)이 있다. 그리고 열에너지의 흐름과 관련된 항목들은 
이다.
제안된 에너지 방정식을 정리하면 다음과 같다.
(4)
열-수리 파이프 요소의 에너지 보존의 법칙(식 (4))을 Galerkin 수식화하고 시간 interval에 대하여 시간적분을 수행하였으며, 이를 유한요소 코딩을 수행하였다. 개발된 파이프요소는 파이프 온도(
)와 이에 접한 주위 지반의 재료의 온도(
)가 결합되어 있으며, 적용물성치는 파이프의 단면적, 유체의 속도, 비중, 비열, 그리고 단열여부 등이다.
파이프 요소는 내부 유체 온도에 대한 2절점 요소형태이며, 유체의 흐름방향은 요소의 node numbering순서를 따른다. 주위 재료와의 convective boundary 조건은 단열(insulated), 대류(convecting), 그리고 열량 주입(heat sourve/sink)의 3가지 종류가 있어 주위 재료와의 다양한 경계조건에 대한 해석을 수행할 수 있다. 따라서, 개발된 일차원 파이프 요소를 순차적으로 연결함으로써 다양한 지중 열파이프 배치(U-type, W-type, Coil-type)를 모사할 수 있으며, 파이프 요소를 closed-loop 형태로 모델링하여 TRT와 같이 지속적인 열유입에 의한 파이프와 주위 지반의 온도변화를 모사할수도 있다.
2.3 유체 대류에 의한 인접한 연속체의 열흐름 방정식의 변경
다양한 종류의 물질이동에 대한 거시적인 평형방정식은 다음과 같이 일반화된 편미분 방정식으로 정리할 수 있다(Shin, 2011).
(5)
~ ~
여기서, 
는 다공질 재료의 종류(예로, 흙입자, 물, 공기)을 나타내며, 
는 그 종류의 단위체적당 질량을 나타낸다. 그리고 
는 해당 종류의 질량 흐름을 나타내고, 
는 source/sink term이다.
에너지 보존의 법칙은 내부 에너지(internal energy)를 이용하여 표현하였으며, 고체, 액체, 그리고 기체의 다공질 재료에 대한 에너지 보존 방정식을 일반화된 평형 방정식을 이용하여 유도하였다. 연속체의 열에너지 보존법칙에서는 전도(conduction)에 의한 열전달만이 고려되므로, 연속체내에 삽입된 유체와의 열적 상호작용인 대류(convection) 현상을 고려하기 위하여 열에너지에 관한 지배방정식을 수정하였다.
~ ~
(6)
여기서 
와 
는 흙입자의 단위중량과 내부에너지이
며, 
는 공극의 포화도를 나타낸다. 
과 
는 액체 상태
~ ~
와 기체 상태의 단위면적당 유량(flow rate)이며, 
, 
~ ~
는 물의 액체와 기체 상태에서의 분자확산(molecular diffusion)과 분산(mechanical dispersion)을 고려할 수 있는 비이류 유량(non-advective flow rate) 항목이다. 그리
고, 다공질 재료에서 열전도(heat conduction)에 의한 에
너지 전달은 Fourier의 법칙(
)을 이용하였으
~ ~
며, 유체 흐름에 의한 대류(
)는 에너지가 유출되므로 음의 부호를 갖게 된다.
2.4 비선형 역해석
지중 열교환 시스템에서 주위 지반의 유효 열물성치를 평가하거나(TRT시험), 추후 계측 data를 이용한 열적 분포특성 평가를 위하여 역해석 기법을 적용하였다. 비선형 역해석은 계측된 값과 예측값의 차이로 정의되는 목적함수(objective function)를 최소화하여 최적의 물성치를 추정하는 방법이다. 이때 목적함수는 값들의 차이의 제곱(L2 norm)으로 정의되고 다음과 같다.
=― = = = =
(7)
=― = = = =
여기서, 
는 미지의 물성치이고, 
는 물성치의 예측 초
기값, 
는 계측된 값(변위, 온도, 간극수압 등), 
는
=
물성치-예측값의 변환 행렬, 그리고 
는 가중치이다.
식 (7)을 
에 관하여 편미분한 후 목적함수를 최소화 하는 
를 산정하면, 새로운 미지수 
는 다음과 같이 산정된다(Santamarina and Dante, 2005).
= = = =
(8)
= = =― = =
비선형 역해석에서 Input-Output의 상관성을 결정하는 변환행렬 
는 해의 수렴성과 밀접한 관련이 있다.
예로, 선형탄성재료의 일차원 변형에 대한 응력과 변형
량에 관한 식은 
과 같다. 이때 계측값이 변위
(
)이고 미지의 변수가 탄성계수(
)일 때, 
이 된다. 일반적으로 사용하는 
와 같은 선형적 변형행렬(Fig. 2)은 해의 수렴성을 크게 저하시킨다. 하지만, 변위(
)와 탄성계수(
)의 역의 관계를 적용하는 
의 변환행렬 설정은 가장 높은 수렴성을 나타낸다(Fig. 2). 그러나, 탄소성재료에서는 역의 변환행렬이 아닌 다른 형태에서 최적의 수렴성을 보일 수도 있다. 따라서 주어진 문제에 따라 다양한 형태의 변환행렬에 대한 수렴성 확인이 필요하다.
역해석에서 개선된 알고리즘을 사용할수록 수렴성이 향상되나, 필요한 순차해석(forward simulation) 횟수가 증가하여 총 연산 소요시간이 증가하게 된다. 본 연구에서는 최소제곱법에 의한 역해석 모듈을 추가하고 계측값과 미지의 변수의 상관관계에 따른 수렴성을 확인하였다. 특히, 지중열교환 시스템에 대한 역해석을 수행하여 계측값과 구하고자 하는 변수의 상관관계를 확인하고 수렴성을 개선하였다.
. (a) Young’s modulus (Etrue=1000), (b) Poisson’s ratio (νtrue=0.3)
3. 수치 예제 해석
개발된 열-수리 파이프 요소에 대한 유한요소 프로그램을 이용하여 수치해석을 수행하였다. 단열된 파이프의 일차원 흐름 해석, 열응답 해석과 이에 대한 역해석에 의한 열물성치 추정에 대한 수치해석을 수행하여 프로그램의 안정성과 적용성을 살펴보았다.
3.1 파이프 요소의 일차원 응답해석
개발된 열-수리 파이프 요소의 수치적인 안정성을 평가하기 위하여 파이프 외부로의 단열상태에서 일차원 흐름에 대한 위치별 온도변화에 대한 수치해석을 수행하였다. 파이프의 직경은 0.1m이고, 내부 유체의 초기온도는 10℃이며, a지점에서 c지점으로 유체의 흐름속도는 0.2m/s이다. Fig. 3은 a지점에 유입되는 유체의 온도를 5℃로 변경하였을 때, 위치별 시간경과에 따른 온도변화를 나타내었다. 산술적으로 5℃의 유체가 b지점에 도달하는 시간은 2.5sec이고, c지점에 도달하는 시간은 5.0sec이다. 개발된 일차원 열-유동 요소가 파이프내의 열흐름을 효과적으로 모사할 있음을 알 수 있다. 다만, 본 예제와 같이 유체 흐름에 의한 열의 이동이 전도에 의한 양보다 큰 경우는 온도구배(thermal gradient)가 매우 큰 경계면 부근에서 의사 진동(spurious oscillations)이 발생할 수 있다(Idelsohn et al., 1996).
= 1.0 w/m･K3.2 열응답 실험 수치해석
지중열교환기의 성능은 지중 열교환기 파이프 내를 순환하는 유체와 파이프 주위 지반의 구성, 그라우팅 재료 등과 밀접한 관련이 있다. 하지만 보어홀의 열저항, 그라우팅 재료의 열물성, 파이프와 보어홀 벽면 사이의 간격 등 열교환기 성능에 미치는 변수들을 모두 개별적으로 고려하기 어려우므로 현장 열응답 시험(TRT, Thermal Response Test)을 통해 지중의 구성 인자들의 특성을 종합하여 유효 열전도도라는 개념을 설계에 적용한다.
현재 현장 열응답실험 결과로부터 지반의 대표 열전도도를 산정하는 방법은 대부분 이론식인 선형 열원 모델(Line-source model)을 적용하고 있다(Carslaw and Jasger, 1959). 지중 열교환기를 무한길이의 직선 열원으로 가정한 선형 열원모델은 열원의 반경 방향만으로 전도에 의한 열전달이 이루어진다고 가정하였다.
(9)
여기서, 
는 보어홀의 길이당 주입열량, 
는 유효 열전도도, 
는 열확산계수, 
는 초기 온도, 그리고 
는 보어홀의 열저항이다
지중에 설치된 열교환 파이프에 지표면의 히트펌프에서 일정한 열량(Q)을 주입하였을 때, 파이프와 지중의 온도 변화에 대한 수치해석을 수행하였다. 해석에 사용된 기하학적 조건, 적용 물성치 및 초기조건은 Table 2에 정리하였다.
Fig. 4는 해석에 사용한 메쉬와 시간경과에 따른 주위지반의 온도변화를 보여주고 있다. 파이프의 유입구 쪽에 일정된 열량이 지속적으로 공급되면서, 파이프내 유체의 온도는 초기에 선형적으로 증가하다가, 시간경과후 지반내부로의 열확산으로 인하여 유체 온도의 증가경향이 급속히 저하된다(Fig. 5a). 순환 초기를 제외하고 유입구와 유출구는 온도차이는 거의 일정하였으며(Fig. 5b), 유출구와 유입구의 평균온도로부터 선형열원모델 이론식(Eq. 9)에 의한 유효 열전도도를 추정하였다(Fig. 5c). Fig. 5d에서 선형열원모델 이론식에 의한 지반의 열전도도 추정방법은 주어진 원지반 물성치에 비하여 10%이상 과다 평가함을 알 수 있다. 또한, 원지반의 다양한 열전도도와 설치심도에 대한 추가적인 수치해석들에도 이론식에 의한 지반 열전도도 평가는 10%내외의 오차가 유발하는 것으로 확인할 수 있었다. 이러한 지반 열전도도의 과대평가는 이론식에서 고려되지 않은 인접 파이프간의 열적상호작용과 파이프 단부효과에 의한 추가적인 열에너지의 외부 방출에 의한 것으로 판단된다. 즉, 유한의 길이를 갖는 파이프의 단부에서 지중으로 열에너지의 지속적 이동과 상-하 방향 열교환 파이프간의 온도차에 의한 열이동이 이론식에서 고려되지 않아, 이러한 에너지 손실을 이론식에서 주위 지반의 열전도도가 큰 것으로 잘못 평가하게 된다.
현재 원지반의 유효 열전도도를 평가하기 위하여 대부분 이론식인 선형 열원 모델(Line-source model)을 적용하고 있지만, 이는 원지반의 열전도도를 과다 산정하여 결과적으로 장기적인 지중 열교환기의 효율을 저하시킬 수 있다.
3.3 열응답 실험에 대한 역해석
TRT 실험결과로부터 지반의 유효 열전도도를 평가하기 위하여 이론식을 적용하는 것은 이론식의 기본가정의 문제로 인하여 지반물성치를 과대평가하는 경향이 있으므로, 열응답 수치해석 결과에 대한 역해석을 수행하여 유효 열전도도를 추정하였다. 이를 통하여 실제 원지반에 대한 열응답시험 자료를 이용하여 안정적이고 빠른 유효 열전도도 도출방법을 제안하였다.
비선형 역해석에서 Input unknown variables(
)과 Output measurement(
) 사이의 상관성을 확인하기 위하여, 다양한 변환행렬 
에 관하여 지반의 유효 열전도도를 평가하기 위한 역해석을 수행하였다. 지반의 열전도도가 높아질수록 주위지반으로 전달되는 열의 증가로 인하여 파이프내의 유체는 느린 속도로 상승하게 된다. 따라서 간략화된 Eq. 9와 같이 유체온도와 지반 열전도도를 역의 관계로 변형행렬 
을 수립할 때 가장 높은 수렴성을 나타내었다(Fig. 6). 변환행렬 
도 높은 수렴성을 나타냈으며, 변환행렬 
을 사용시 해가 수렴하지 못하고 발산함을 알 수 있다. 따라서, TRT 계측자료에 대한 역해석을 수행하여 지반의 유효 열전도도를 평가할 때, 변환행렬 
나 
를 사용하는 것이 바람직한 것으로 판단된다.
in TRT (
true = 1.0 W/m･K)4. 결 론
지중 열교환 시스템에서 유체가 순환하는 파이프와 이에 인접한 지반의 열적 상호작용은 매우 중요한 열에너지 통로 역할을 한다. 하지만, 파이프의 수치모델링에서 열-수리가 연관된 난류해석에 대한 어려움과 파이프의 긴 세장비에 의한 메쉬 사이즈의 부적합성은 열교환 시스템의 적절한 수치해석을 복잡하게 하고 있다. 따라서, 본 논문에서는 원형 단면 파이프내의 난류흐름에 대하여 열-유체가 연동된 일차원 파이프 요소로 개발하고, 이를 기 개발된 다공질 재료를 위한 THM 유한요소 프로그램과 결합하였다. 개발된 프로그램을 이용하여 지반의 열전도도 평가를 위한 TRT해석과 유효 열물성치 평가하기 위한 역해석을 수행하였다.
(1)파이프요소의 에너지 흐름은 유체의 전도, 변형과 유체의 이송에 의한 energy flux, 파이프내의 유체 흐름에 의한 대류, 그리고 내부 source로 구성되며, 이에 대한 에너지 보존의 법칙을 적용하였다. 제안된 에너지 방정식에 Galerkin수식화와 시간적분을 수행하고 이를 바탕으로 유한요소 프로그램 코딩을 수행하였다. 그리고 연속체의 열에너지 지배방정식에 파이프 요소와 주변 연속체 사이의 대류에 의한 열적 상호작용을 추가 수정하였다.
(2)개발된 요소를 이용한 수치해석 결과에서 TRT 실험결과로부터 주위지반의 유효 열전도도를 추정하기 위하여 사용하는 선형 열원 모델이 지반의 열전도도를 과다 평가하는 것으로 보여주었다. 이는 이론식에서 고려되지 않은 인접 파이프간의 열적 상호작용과 파이프의 단부효과에 의하여 추가적인 열에너지 방출에 기인한 것으로 사료된다.
(3)TRT 실험결과로부터 유효 열전도도를 평가하기 위하여 역해석을 수행하였으며, 역해석시 
나 
의 변환행렬이 높은 수렴성을 나타내었다.
본 연구에서 개발된 열-수리 연동 파이프요소를 이용하여 지중열교환 시스템의 다양한 열-유체 흐름과 변형에 대한 수치해석을 수행할 수 있다. 그리고 파이프의 배치(U, W, Coil type 등)에 따른 열성능 평가, 지하수 흐름에 대한 등가 열전도도의 평가와 장기 열거동 예측 등에 활용할 수 있다. 또한 향후 실내모형이나 현장계측 결과를 이용하여 수치해석 결과에 대한 비교 검증을 수행하고자 한다.
