Home > Vol. 16, No. 1


[ Research Articles ]
Journal of the Korea Institute of Ecological Architecture and Environment - Vol. 16, No. 1, pp. 81-87
ISSN: 2288-968X (Print) 2288-9698 (Online)
Print publication date Feb 2016
Received 04 Jan 2016 Revised 18 Feb 2016 Accepted 22 Feb 2016
DOI: https://doi.org/10.12813/kieae.2016.16.1.081

건물의 3차원 구조체에 대한 전열해석 프로그램 개발 중 서로 다른 열전도율을 갖는 복합재질 3차원 구조의 비정렬 격자에 대한 전산해석 방법
이주희* ; 장진우** ; 이현균** ; 이용준*** ; 이규성****

Numerical heat transfer analysis methodology for multiple materials with different heat transfer coefficient in unstructured grid for development of heat transfer analysis program for 3 dimensional structure of building
Lee, Juhee* ; Jang, Jinwoo** ; Lee, Hyeonkyun** ; Lee, Youngjun*** ; Lee, Kyusung****
*Corresponding author, Dept. of Computer Engineering, Hoseo Univ., South Korea (Juheelee@hoseo.edu)
**Graduate School of Mechanical Engineering, Hanyang Univ., Korea
***BEL Technology, Seoul, Korea
****Doall Tech, Seoul, Korea (kslee@doalltech.com)

© Copyright Korea Institute of Ecological Architecture and Environment
Funding Information ▼

Abstract
Purpose

Heat transfers phenomena are described by the second order partial differential equation and its boundary conditions. In a three-dimensional structure of a building, the heat transfer phenomena generally include more than one material, and thus, become complicate. The analytic solutions are useful to understand heat transfer phenomena, but they can hardly be applied in engineering or design problems. Engineers and designers have generally been forced to use numerical methods providing reliable results. Finite volume methods with the unstructured grid system is only the suitable means of the analysis for the complex and arbitrary domains.

Method

To obtain an numerical solution, a discretization method, which approximates the differential equations, and the interpolation methods for temperature and heat flux between two or more materials are required. The discretization methods are applied to small domains in space and time, and these numerical solutions form the descretized equations provide approximated solutions in both space and time. The accuracy of numerical solutions is dependent on the quality of discretizations and size of cells used. The higher accuracy, the higher numerical resources are required. The balance between the accuracy and difficulty of the numerical methods is critical for the success of the numerical analysis. A simple and easy interpolation methods among multiple materials are developed. The linear equations are solved with the BiCGSTAB being a effective matrix solver.

Result

This study provides an overview of discretization methods, boundary interface, and matrix solver for the 3-dimensional numerical heat transfer including two materials.


Keywords: Unstructured grid, Heat transfer analysis, Multiple materials, Boundary treatment
키워드: 비정렬격자, 전열해석, 복합소재, 경계면처리

1. 서론
1.1. 연구의 배경 및 목적

건물의 실내외 온도차이로 인하여 벽체나 창호를 통해 열전달이 발생하게 된다. 열전달현상은 대표적인 포아송방정식(Poisson's equation)으로 중첩원리를 이용하여 비교적 복잡한 형상으로 이루어진 영역에 대해서도 해석해(analytic solution)를 구할 수 있다. 그러나 대상이 되는 계의 형상이 매우 복잡해짐에 따라 해석적인 방법보다는 수치적인 방법에 의존해 근사적인 해를 구하게 된다. 수치적 방법으로 해를 구하기 위해서는 전체 계산 영역을 여러 개의 작은 제어체적(control volume) 혹은 격자(cell)로 나누고 작은 제어체적에 대하여 지배방정식(governing equation)인 포아송방정식을 적용한다. 사용하는 격자의 모양과 이를 다루는 자료구조(data structure)에 따라 크게 정렬격자(structured grid)와 비정렬격자(unstructured grid)로 나누어진다. 정렬격자는 3차원의 경우 6면체, 2차원의 경우 4면체격자로 구성이 되며 격자선(grid line)을 따라 인덱스(일반적으로 i, j, k)를 이용하여 형상, 온도 등 계산에 필요한 정보에 직접 접근할 수 있다. 가장 대표적인 카트시안 격자(Cartesian grid)의 경우 각 격자는 좌표축을 따라 형성되어 있어 일차미분(대류항)과 이차미분(확산항)의 차분화가 매우 직관적이다. 또한 Taylor series에 의한 분석으로 해의 정밀도 분석이 가능하다. 복잡한 형상에 대하여 물리적공간(physical domain)과 계산공간(computational domain)을 서로 맵핑(mapping)하는 방법을 사용하거나 다중블록(multi-block)방법을 이용하여 축에 나란하지 않은 임의의 형상에 제한적으로 적용이 가능하다[1]. 그러나 이러한 정렬격자는 해석하고자 하는 형상이 복잡해지고 영역이 단순연결(simple connected domain)이 아니라 다중연결(multiple connected domain)의 경우 격자생성이 어렵거나 불가능해지는 단점을 가지고 있다. 이에 반해 비정렬격자는 복잡한 형상에 적용이 가능하며 필요한 부분에 쉽게 격자를 밀집시킬 수 있다는 장점을 가지고 있다[2]. 그러므로 대부분의 상용해석프로그램은 비정렬격자를 사용하고 있다. 비정렬격자는 이산화된 대수방정식은 성긴 행렬식(sparse matrix)이기되기 때문에 메모리의 낭비와 성긴 행렬을 다루기 위한 기술이 필요하다.

비정렬격자를 사용하는 유한체적법(finite volume method, FVM) 중에서 해를 저장하는 방식에 따라 셀중심법(cell-centerd method), 격자중심법(vertex-centered method), 셀-격자점법(cell-vertex method)로 나눌 수 있으며 이 중 셀 중심방법은 기하학적으로 매우 단순하며 검사체적에 대하여 특별한 보간방법을 사용하지 않고 보존법칙을 만족하는 장점을 가지고 있어 가장 널리 사용되고 있다[2]. 제어체적이 임의의 형상을 가지고 있으므로 이산화한 방정식이 2차이상의 정밀도를 가지도록 하는데 어려움이 있다. 확산플럭스의 수치근사에서 격자면을 공유하는 인접한 두 셀의 중심을 잇는 선과 셀면의 수직한 두 성분에 대하여 이산화하는 방법을 사용하여 정확성과 안정성은 높이는 시도가 있었다[3].

건축물의 벽체와 같은 구조체에서는 단열을 위해 여러 개의 재질로 이루어진 다층구조를 가지게 되며 이들 사이의 열전달은 동일한 재질로 이루어진 구조체의 열전달과는 상이한 열전달 특성을 가진다. 동일 재료로 이루어진 두 인접한 격자의 공유면에서는 물성치가 분명하게 정의되나 복합재질로 이루어진 구조체의 접합면의 경우 서로 다른 재질로 이루어진 두 격자가 공유하고 있는 면에서 물성치는 물리적으로 명확하게 정의 되지 않는다. 이러한 면에서의 물성치는 실제 물성치가 아니며 인접한 두 재료간의 열전달을 정확하게 모사할 수 있는 방법으로 정의되어야 한다. 격자 표면에서의 물성치를 적용이 쉽고 단순한 산술평균(arithmetic mean)을 사용하는 경우, 이 면을 통한 열유속이 적절히 모사되지 않으며 불합리한 결과를 얻을 수 없다[4]. 또한 다른 재료 간에 접촉면에서 격자점이 서로 일치하지 않는 경우 서로 만나는 면적에 관한 추가적인연산이 필요하게 된다. 본 연구에서는 이전 연구인 단일재료로 이루어진 전열해석[5]에 이어 복잡한 형상에 적용이 가능하도록 3차원 비정렬격자를 이용하여 복합재료로 이루어진 영역의 전열해석을 하기 위하여 이산화, 행렬해석기(matrix solver), 복합재료간의 경계조건을 다루기위한 수치적방법에 관한 연구를 수행하였다.


2. 전산응용 전열해석
2.1. 적분형 수송방정식

유체 및 열유동의 경우, 지배방정식을 유도하기 위하여 개개의 입자를 고려한 접근보다는 공간내의 임의의 영역에서의 검사체적을 이용하여 에너지(온도), 운동량, 질량의 연속된 흐름을 관찰하게 된다. 이러한 관찰은 최종적으로 다음과 같은 일반화된 적분형의 Reynolds수송방정식(Reynolds transport equation)을 얻게 된다[3].

tVρdV+Sρυ-ΓdS=VQdV(1) 

식(1)은 왼쪽부터 첫 번째 항은 시간에 따라 변화하는 양인 비선형항(unsteady), 두 번째 항의 첫 번째 항은 제어체적의 면을 통과하는 유체에 의해 전달이 일어나는 대류(convection), 두 번째항은 확산항(diffusion)이며, 오른쪽항은 생성항(source term)을 나타낸다. 식(1)에서 t , ρ, V , υ, S 는 각각 시간(sec), 밀도(kg/m3), 제어체적(m3), 속도(m/s), 제어체적의(m2) 표면을 나타낸다. 또 Γ는 종속변수(Φ)에 따라 달라지며 전열해석의 경우 Γ는 열전도 계수(k)와 관련된 값이 된다. 본 연구의 전열해석의 경우 대류항이 없는 수송방정식이 된다. 전열해석의 경우 확산항과 생성항이 차분화의 대상이 되며 비정렬격자를 사용하기 때문에 확산항과 격자의 표면에서의 종속변수 값에 대한 적절한 모델링하거나 근사화가 필요하다.


Fig. 1. 
Control volume of unstructured cell and its surfaces.


Fig. 2. 
Gradient calculation on cell face

본 연구의 대상은 열전도에 의한 전열해석을 위하여 대류항을 제외하였으며, Φ = h = cT , Γ = k/c을 취하고, 확산항에 대하여 Gauss의 발산정리를 이용하면 식(1)은 미분형태의 일반적인 비정상 열전도 방정식이 된다.

ρcTt=k2T+Q(2) 

여기서c는 비열(kJ/kg), k열전도계수(W/m), Q 는 생성항을 나타낸다.

2.2. 셀중심-비정렬격자 이산화

비정렬격자에서 사용하는 격자는 사면체, 육면체, 피라미드, 프리즘, 다면체등 다양한 형상의 격자가 사용이 가능하다. 이를 위해서 이산화는 좌표와 무관하게 이루어져야 한다. 가장 일반적인 형상의 격자를 fig 1에 나타내었다. 본 연구에서는 명현국(2013)이 사용한 방법과 같은 방법을 이용하여 셀중심-비정렬 격자에 적용할 수 있도록 이산화를 수행하였다[2][3]. 격자의 중심에 종속변수를 저장하며, Sj는 이웃한 미소제어체적과 공유하는 면(surface)을 나타내며 P0 , Pj는 계산하고자하는 격자와 이웃한 격자의 중심을 각각 나타낸다. dsj는 이웃한 격자까지의 거리를 나타내는 벡터이다. 식(2)의 확산한은 가우스 이론(Gauss theorem)을 이용하면 제어체적의 표면적분을 이용하여 다음과 같이 이산화가 가능하다.

SΓϕS~jΓϕj*Sj(3) 

식(3)의 적분은 제어체적의 표면에서 이루어지는 적분이므로 이산화를 완성하기 위해서는 제어체적의 표면에서 확산항을 적절히 모델링하여야 한다. 제어체적의 표면에서 값인 (∇Φ)j*을 구하기 위하여 Fig 2에서와 같이 cell의 중심을 연결하는 방향과, 표면벡터와 직각을 이루는 2개의 방향으로 나누고 s에 대해서는 두 점사이의 값을 직접 이산화하고 나머지 한 방향에 대해서는 모델을 사용하게 된다. n-αs가 제어체적의 표면과 나란하다면 n-αsn=0이 되어야함으로 α=1/sn이 된다. 그러므로 Fig 2에서 보는 것과 같이 면에서의 ϕjn는 두 벡터의 합으로 나타낼 수 있다.

ϕjn=ϕjαs+ϕjn-αs~~ϕjαs+¯ϕjn-αs(4) 

여기서 n는 표면에 수직인 단위벡터를 나타내며, srPj-rP0방향의 단위벡터를 나타낸다.

첫째 항은 두 셀의 중심방향을 따른 미분성분을 나타냄으로 차분방법을 사용하여 근사화 할 수 있다. 또 두 번째 항은 셀의 중심에서 값을 이용해서 내삽하여 사용할 수 있다.

ϕj~αs=αϕPj-ϕPodsj(5) 
ϕj¯=ωjϕP0+1-ωjϕPj(6) 

α을 계산하여 대입하면 다음과 같이 된다.

α=1sn=nnsn=dsjSjdsjSjn(7) 
ϕj*=ϕPj-ϕP0SjdsjSj+ϕj¯-ϕj¯dsjSjdsjSj(8) 

만약 ns가 서로 나란하다면 즉 정렬격자와 같은 격자가 된다면, 우변의 두 번째 항과 세 번째 항은 서로 상쇄되어 없어진다. 그러므로 다음과 같이 FDM(finite Difference Method)와 같이 면을 중심으로 양쪽 격자에서의 값이 된다.

ϕj*=ϕPi-ϕPjdsj(9) 

식(8)을 이용하여 확산항을 계산하면 다음과 같이 된다.

Dj=ΓϕjSj=ϕPj-ϕPoΓSjdsjn+Γϕj¯Sj-ϕj¯dsjSjdsjn(10) 

식(10)에서 오른쪽 두 항은 내재적으로 다루며 나머지 항은 외재적으로 다루며 생성항으로 처리하게 된다[6].

2.3. 복합재질 경계면처리

서로 다른 두 물질이 만나는 경계에서 격자 생성방법에 따라 격자면이 서로 일치하지 않을 수 있다. 이러한 경계를 지나는 플럭스를 계산하기 위해서는 추가적인 형상에 대한 정보가 필요하다. 한 면이 만나는 반대쪽 경계면의 개수, 만나서 나누어지는 작은 면들의 크기, 중심위치 등의 정보가 필요하다. 본 연구에서는 이러한 추가적인 형상정보를 얻기 위하여 그림 4와 같이 두 면사이체 추가적인 점을 넣고 이를 이용해서 구하는 방법을 이용하였다.


Fig. 3. 
Neighboring cells and virtual points for interface

먼저 경계면에 일정한 간격으로 가상의 점을 분포시킨다. 이 점을 포함하고 있는 면을 그림에서 보는 것과 같이 검색한다. 그 다음 양쪽 격자의 개수와 같은 크기의 2차원 배열을 정의하고 가상의 점들을 각 격자에 해당하는 배열 안에 가상의 점들을 저장한다. 간단한 2차원의 예로서 Fig 3에서 가상의 점을 위에서부터 차례대로 1~9번의 점이라고 하면, 왼쪽은 3개 오른쪽2개의 격자로 구성됨으로 필요한 배열은 3x2크기의 2차원 배열이 된다. 1번점의 왼쪽은 FL1 오른쪽은 FR2이기 때문에 (FL1,FR2)에 1번이 저장이 된다. 마찬가지로 2, 3번 점은 1번과 같은 (FL1,FR2)에 저장이 되고 나머지 점들은 Fig. 4와 같이 저장된다. 이 저장된 값을 바탕으로 양쪽격자가 서로 만나 생기는 접면의 크기와 그 면의 중심 좌표를 구할 수 있다. 배열에서 비어 있는 공간에 해당하는 한쪽 격자와 다른 한쪽 격자는 실제로 만나지 않게 되는 격자에 해당된다. 이 방법의 특징은 유한체적법에서 중요한 경계면에서 열유속값이 서로 일치하며 두 격자사이에 만나는 면이 경미하거나 세장비가 긴 경우에 의한 문제가 발생하지 않는다. 또한 정밀도를 높이기 위해서는 단순하게 점의 개수를 늘림으로 이루어진다. 단점으로는 점의 개수를 늘림에 따라 계산시간이 기하급수적으로 늘어나게 된다. 이러한 부분에 대한 추가적인 연구를 수행할 예정이다.

차분화한 선형방정식인 식(7)의 해를 풀기 위해서는 인접 면에서의 온도를 구해야 한다. Fig. 4과 같이 열전도 계수가 k1k2로 다른 두 셀 사이에서의 온도는 다음과 같은 방법으로 구할 수 있다.


Fig. 4. 
Heat flux of interface at the material interface

q1=k1A1T1-TwL1(11) 
q2=k1A2T2-TwL2(12) 
qi=k2A1+A2+Tw-TiLi(13) 
qi=q1+q2+(14) 

식(11)은 1번 셀 중심에서 인접 면까지의 플럭스를 구한 식이다. k1은 열전도 계수, A1은 면적, T1은 1번 셀 중심온도, Tω 은 구하고자 하는 인접 면에서의 온도이다. 식(12), (13)식(11)과 마찬가지로 2번, i번 셀에서의 플럭스를 구한식이다. 식(11), (12), (13)에서 구한 플럭스의 총합이 같으므로 Tω를 미지수로 놓고 인접면의 왼쪽에 위치한 셀의 개수를 무한하다고 가정하고 풀게 되면 인접 면에서의 온도 Tω은 다음 식 (15)과 같다.

Tω=k1A1T1L1+A2T2L2++k2A1+A2+TiLik2A1+A2+Li+k1 A1L1+A2L2+(15) 
2.4. 반복계산 행렬해석

식(1)을 이산화하면 격자개수와 같은 크기의 행렬이 구성된다. 전체 계산시간 중 선형방정식의 해를 구하는데 많은 시간이 소모된다. 그렇기 때문에 행렬의 수렴성을 높이는 것은 수치해석 의 성능을 높이는데 중요한 요인이 된다.

Ax=b(16) 

An × n의 정방행렬이며, xb는 차원이 n인 벡터가 된다. 식(16)에서 해를 구하는 방법으로 직접적인 방법과 반복적인 방법이 있다. 직접적인 방법은 소거에 의한 Gauss-Jordan법, Gauss법이 있으며, 반복적인 방법에는 Jacobi법, Gauss-Seidel법이 있다. 이러한 방법은 연산이 직관적이고 간단하다는 장점이 있지만 수렴성이 떨어지는 치명적 단점이 있다. 반면에 반복법 중 공액구배법(CGM, Conjugate Gradient Method)은 A 행렬이 대칭성을 가지고 양의 부호(positive definite)라면 효율적으로 해를 구할 수 있다[7].

식(16)은 다음과 같이 목적함수를 가지며 제한(constraint)이 없는 최적해 문제로 변환이 가능하다.

min Xx=12xTAx+xTb(17) 

식(17)를 미분하면 식(16)이 얻어지며 식(17)는 볼록함수(convex function)이기 때문에 최소값은 식(16)의 해가 된다. 식(17)의 최적해를 얻기 위하여 최급강하방향(steepest decent direction)을 따라 이동하며 최소값을 검색하고 이 최솟값을 갖는 위치에서 새롭게 검색방향을 찾아 지역적인 최소해(local optimum)를 찾게 된다. 최급강하 방향은 다음과 같이 구할 수 있다.

pn=Xn=b-Axn(18) 

윗첨자 nn번째 계산에서의 값을 나타낸다. pn 방향으로 최솟값은 다음과 같이 구할 수 있다.

xn+1=xn+ηpn(19) 

여기서, ηp방향으로 이동거리를 나타내는 값으로 p방향으로 최솟값이 되는 곳의 위치를 나타내며 다음과 같이 얻어진다.

η=pnTpnpnTApn(20) 

이러한 방법을 최급강하법(steepest decent method)이라 하며 매 반복마다 최급강하방향을 따라가게 됨으로 해 근방에서 수렴성이 급격히 떨어진다. 이렇게 반복하여 얻어진 최종해는 xn + 1xn 이 되며 식(17)의 최솟값이 되며 동시에 식(16)의 해가 된다. 공액구배법은 최급강하방향이 항시 서로 직교한다는 원리를 이용하여 수렴속도를 향상시키는 방법이다. 다음과 같은 식을 이용하여 해를 개선해 나간다.

η=rnTrnpnTApn(21) 
rn+1=rn-ηnApn(22) 
pn+1=rn+1+βnpn(23) 
βn=rn+1Trn+1rnTrn(24) 

식(22)은 나머지 벡터(residual vector)를 의미한다. 공액구배법은 대칭행렬에서만 쓰일 수 있다는 단점 때문에 BiCG(Bi-Conjugate Gradient Method) 또는 BiCGSTAB(BiCG STABlized)방법이 많이 사용된다. BiCG법은 비대칭행렬에 적용이 가능하나 전치행렬을 구해야 하는 단점이 있다. BiCGSTAB는 BiCG법 보다는 기울기를 구할 때 좀 더 안정적이고 빠르게 수렴하기 위한 방법이 사용된다[8]. 최적의 기울기를 찾기 위한 식은 다음과 같다.

s=r-αAp(25) 
ω=As, sAs, As(26) 

식(25)식(26)과 같이 최적의 기울기를 찾기 위한 새로운 행렬과 값을 정의한다. 여기서, α는 나머지 행렬이 직교한다는 가정으로부터 나온다. 최종적으로 BiCGSTAB에서 해 벡터를 다음과 같이 구할 수 있다.

xN+1=xN+αpN+ωsN(27) 

또한, 전치행렬을 구할 필요가 없어서 계산에 필요한 메모리를 줄일 수 있다는 장점이 있다. 본 연구에서는 BiCGSTAB방법을 사용하여 식(16)의 해를 구하였다. 알고리즘은 다음과 같다[9].



step01 : r0 = b -Ax0 , r0^=r0
step02 : ρ0 = α = ω0 = 1
step03 : υ0 = p0 = 0
DO i = 1,2,3 ⋯
step04 : ρi =<r0^,ri - 1>
step05 : pi = ri - 1 + β (pi - 1 - ωi - 1υi - 1)
step06 : υi = Api
step07 : α = pi/ <r0^,υi >
step08 : s = ri - 1 - αυi
step09 : t = As
step10 : ωi =< t,s > / < t,t >
step11 : xi = xi - 1 + αpi + ωis
step12 : ri = s - ωit
END

여기서 <,>은 내적연산을 나타낸다.


3. 결과 및 토의
3.1. 행렬 계산기 검증

미분방정식을 차분화하면 선형방정식이 얻어지며 컴퓨터의 빠른 계산 능력을 이용한 반복 계산이 가능하다 . 이러한 행렬식을 풀기 위하여 본 연구에서는 반복법 중 하나인 BiCGSTAB를 사용하였으며 코드의 정확성을 확인해보기 위하여 시험모델에 적용하였다. 시험모델의 행렬은 다음과 같다.

Ax=b(28) 
A-12I, b=-sinx(29) 
A=-0.50.0051160.005121-0.50.0053160.0051340.0050710.005278  -0.5  (30) 
x2sinx(31) 

식 (28)는 일반적인 선형시스템의 방정식을 나타내며, 이 때 행렬 A식 (29)(30)과 같은 값을 가지도록 하였다. 행렬A의 대각성분은 모두 -0.5을 가지며 나머지 성분들은 대각성분의 절대값 보다 충분히 작은 임의의 값을 가진다. 벡터b은 -sinx의 값을 가지도록 하였다. 이렇게 함으로 식(28)의 해는 근사적으로 2sinx의 값을 가지게 된다. 행렬A 가 정확하게 단위행렬을 나타내고 있지 않으며 대각성분을 제외한 값들은 완전한 “0”이 아닌 “0”에 가까운 작은 값들을 나타내기 때문에 벡터b 역시 정확한 2sinx값을 나타내지는 않는다. 또한 행렬들의 크기가 커지면서 상태수(condition number)가 증가하고, 크기가 100에 가까울 때는 상태수가 크게 증가하여 수치적으로 해를 구하기 어려운 불량조건(ill condition)상태가 된다. 이러한 방정식을 BiCGSTAB방법을 이용하여 해를 구하는데, 식(30)을 보게 되면 행렬A은 대칭행렬이 아니기 때문에 CG법으로는 계산할 수 없다. Fig 5에서 2sinx의 값과 비교하였다. 그림에서는 잘 나타나지 않으나 미세한 차이가 있으며 전체적인 모양은 2sinx와 같은 모양을 나타낸다. 이러한 미소한 차이는 행렬A가 정확한 단위행렬이 아니기 때문에 발생하는 노이즈(noise)이며 이 노이즈는 같은 행렬크기에서 상태수를 제어하는데 사용된다.

계산시간을 확인하기 위해서 일반적으로 해를 구하는데 많이 사용되는 LU분해법(lower upper decomposition)과 비교하였다. 행렬A의 크기를 500 × 500, 1000 × 1000, 2000 × 2000, 5000 × 5000으로 바꿔가며 계산시간을 비교하였으며 결과는 Fig 6에 나타내었다. 격자의 개수가 많아짐에 따라 행렬의 크기가 증가하게 됨으로 계산 시간은 모두 지수함수적으로 증가하고 있다. 그러나 LU방법과 BiCGSTAB를 비교하면 BiCGSTAB가 우수한 성능을 보이고 있으며 격자수가 많아짐에 따라 BiCGSTAB가 더욱 유리해짐을 알 수 있다.

행렬해석기의 정확성을 확인하기 위하여 실제적인 유한체적법과 유사하면서 단순화하여 쉽게 해를 구할 수 있는 정렬격자의 유한차분법(FDM, Finite Difference Method)을 이용한 2차원 정상상태의 열전달 해석을 수행 하였다. Fig 7에 보는 것과 같이 계산영역의 아래 면에 일정한 온도(200K)를 주었으며 나머지는 면은 0K으로 고정을 시켰다. 정상상태 2차원 문제로 가정하였으며 사용한 사각형 격자를 사용했으며 격자의 개수는 800개(20 × 40)이다. 결과는 기존 연구[5]와 비교하였다. 유한차분화법에 대한 자세한 내용은 본 연구의 범위를 벗어남으로 생략했으며 최종적으로 얻어진 방정식의 행렬A는 폭이 3인 삼대각행렬(tri-diagonal matrix)이다. 정상상태의 온도는 이전 연구결과와 함께 Table 1에 나타내었으며 소숫점 2자리까지 서로 일치하고 있으며 오차율은 1 × 10-2%이하를 보이고 있다.

Table 1. 
Temperature along x at y=10
x y Temperature[K]
TDMA[5] CGM[5] BiCGSTAB
0 10 0.00000 0.00000 0.00000
5 10 38.02782 38.02782 38.02812
15 10 52.20471 52.20471 52.20377
15 10 38.02782 38.02782 38.02812
20 10 0.00000 0.00000 0.00000


Fig. 5. 
Solution comparison between analytic function and matrix solver


Fig. 6. 
Computational time between LU decomposition and BiCGSTAB


Fig. 7. 
Boundary conditions and temperature distribution at y=10

3.2. 프로그램 검증

FVM방법은 복잡한 형상에 쉽게 적용할 수 있다는 장점을 가지고 있는 반면 자료구조가 복잡하고 해를 안정시키기 위한 다양한 기술을 필요로 한다. 격자의 형태 중 육면체 격자가 가장 열의 전달방향과 물리적으로 비슷할 뿐 아니라 차분화에 의한 수치적 오차가 가장 작다. 본 연구에서는 복잡한 형태의 열전도 해석을 위하여 사면체 격자를 사용하여 구현하였으며 이를 검증하기 위하여 기존연구[5]의 육면체 격자를 사용한 해석 결과와 본 연구에서 4면체 격자를 사용한 계산 결과와 비교하였다. 또한 해석해(analytical solution)과 비교하여 수치적방법의 타당성을 확인하였다.

Fig 8은 단일재료로 길이는 5이며 가로 세로는 각각 1, 2인 직육면체이다. 직육면체의 양쪽 끝단의 온도는 각각 250K, 300K으로 고정되어 있으며 그 외의 면은 반복경계를 사용하였다. 정상상태로 가정하였으므로 단일재료에서의 온도는 선형으로 나타난다. 온도는 계산 영역내의 모든 점의 온도를 출력하고 길이 방향(z)으로 정렬하여 나타낸 것이다. 육면체의 경우 격자는 640개, 사면체의 경우 800개를 사용하였다. 사면체 격자의 경우 위치가 미미하게 차이를 보이고 있는데 이는 여러 개의 점을 한꺼번에 출력했기 때문이다. 육면체, 사면체 격자 모두 해석해와 잘 일치하고 있으며 모두 상대오차는 최대가 0.1%로 해석해와 잘 일치하고 있다.

복합재질에 의한 경계면을 가지는 수치해석방법을 검증하기 위하여 단순하면서 해석해를 구할 수 있는 열전도문제에 적용하였다. 계산 영역은 Fig 9에서 보는 것과 같이 z방향으로 긴 막대형상이며 직육면체의 양 끝단은 일정한 온도 (250K, 300K)의 유지된다고 가정하여 온도 경계 조건을 사용하고 나머지 면은 반복경계(periodic boundary condition)를 사용하였다. 2절에서 제시한 수치적방법을 이용해서 3차원 열전도 해석을 수행하였다. 3차원계산을 수행하나 양 끝단을 제외하고 모두 무한 경계가 됨으로 1차원형태의 해를 얻게 된다. 복합재질을 분석하기 위하여 z방향으로 두 영역으로 나누고 각기 다른 열전도계수를 갖는다고 가정하였다. Fig 9에서 보는 바와 같이 왼쪽 부분의 열전도 계수는 1로 하였으며 오른쪽 부분의 열전도 계수는 4로 설정하였다. 수렴은 상대오차가 1 × 10-4이하가 되면 수렴한 것으로 판별하였다. 정상상태 온도분포임으로 z축의 위치에 따른 온도는 선형분포를 가지게 된다. 온도분포는 두 재질이 서로 만나는 면을 기준으로 각기 다른 온도 분포를 가진다. 해석해와 온도 분포를 비교하였으며 결과는 해석해와 잘 일치하는 것을 알 수 있다. 이로서 본 연구의 4면체 비정렬격자 방법을 이용한 3차원 전열해석 코드는 복합재료에서의 온도전달 현상을 잘 모사하고 있음을 알 수 있다.


Fig. 8. 
Temperature distribution comparison between tetrahedral, hexahedral meshes and exact solution


Fig. 9 
Computational domain


Fig. 10. 
Comparison temperature distributions of multiple materials


4. 결론

본 연구는 3차원이면서 다양한 재질로 구성된 건축물 내의 전열해석을 수행하기 위하여 복잡한 형상으로 이루어졌을 뿐 아니라 상이한 열전달계수를 갖는 복합재질로 이루어진 계산 영역의 전열해석에 적용할 수 있는 비정렬 격자계를 사용한 정상 3차원 전열해석 프로그램을 구현하였다. 열전달 방정식을 차분화하여 얻은 선형방정식을 해석하기 위한 행렬 계산기, 차분화 방법 및 복합재질을 위한 열전달에 대한 검증을 수행하였다. 먼저 행렬해석기의 정확성을 확인하기 위하여 다양한 크기의 행렬에 대하여 계산을 수행하였다. 행렬해석에 있어서 BiCGSTAB는 반복법을 사용함으로 해가 안정되며 성긴행렬의 연산에 적합한 방법임을 알 수 있었다. 또 4면체 비정렬격자에 적용이 가능한 FVM 이산화방법의 정확성을 확인하기 위하여 간단한 3차원형상에 대해 전열해석을 수행하였다. 수치해석 결과는 정상해와 잘 일치하는 선형온도분포를 얻었으며 이를 통하여 이산화방법이 정확하게 구현되었음을 알 수 있었다. 또 2개의 서로 다른 전열계수를 가지는 재질의 경계면을 다룰 수 있는 간단하면서 구현이 쉬운 방법을 제안하였으며 그 결과를 해석해와 비교를 통하여 복합재료에 대한 수치기법 및 타당성을 검증하였다.

향후 3차원 정상/비정상 전열해석에 있어 사용자 숙련정도와 무관하게 해석이 가능하도록 수치안정성, 2개 이상의 복합재질, 완화법, 격자생성법, 행렬해석기 전처리의 연구를 수행할 예정이다.


Acknowledgments

본 연구는 국토교통부 주거환경연구사업의 연구비지원(15RERP-B082204-02)에 의해 수행 되었습니다.


References
1. THOMPSON, Joe F., WARSI, Zahir UA, MASTIN, C. Wayne, Numerical grid generation foundations and applications, Amsterdam North-holland, (1985).
2. 명현국, “비정렬 셀 중심 방법에서 대류플럭스의 수치근사방법 평가”, 한국전산유체공학회지, 11(1), (2006).
Myong, H.K., Evaluation of Numerical Approximations of Convection Flux in Unstructured Cell-centered Method, 11(1), (2006).
3. 명현국, “비정렬 셀 중심 방법에서 확산플럭스의 새로운 수치근사방법”, 한국전산유체공학회지, 11(1), (2006).
Myong, H.K., A New Numerical Approximatiojn of Diffusion Flux in Unstructured Cell-centerd Method, 11(1), (2006).
4. 명현국, 전산열유체공학, 문운당, (2009).
Myong, H.K., Numerical heat and fluid flow engineering, Munundang, (2009).
5. 이주희, 장진우, 이용준, 최준혁, 이상환, “3차원 비정렬격자를 이용한 전열해석 방법론”, 대한건축친환경설비학회 논문집, 8(6), (2014).
Lee, Juhee, Jang, Jinwoo, Lee, Yongjun, Choi, Junhyuck, Lee, Sanghwan, Methodology for numerical heat transfer considering 3-dimensional unstructured grid, Journal of KIAEBS, 8(6), (2014).
6. Hadzic, Hidajet, Development and Application of a Finite Volume Method for the Computation of Flows Around Moving Bodies on Unstructured, Overlapping Grids, Ph.D thesis, Hamburg University of Technology, (2005).
7. Golub, Gene H., Loan, Charles F. Van., Matrix computations, (4th ed), Baltimore, The John Hopkins University Press, (2013).
8. Sickel, S., Yeung, M. C., Held, M. J., A comparison of some iterative methods in scientific computing, Summer Research Apprentice Program, (2005).
9. Saad, Yousef, Iterative Methods for Sparse Linear Systems, 2nd Ed., Society for Industrial and Applied Mathematics, (2003).