KIEAE Journal

한국 현대사회에 있어서 공동주택은 삶의 질 향상과 공간 활용의 극대화라는 장점으로 인하여 수요가 계속적으로 증가해 왔다. 2006년 1월부터 공동주택 발코니 확장이 법적으로 가능해지게 됨에 따라 활용할 수 있는 공간이 늘어났다는 장점이 있는 반면 열적완충공간이 사라져 단열 성능의 저하라는 문제를 야기하였다.[1] 그 결과 단열성능 향상 혹은 에너지의 손실을 막기 위해 전열에 관한 실험연구들[2, 3]이 수행되었다. 그러나 실험에 의한 연구는 실험대상구조물을 직접 만들고 건물의 제어된 외부와 실내조건을 구현해야 한다는 어려운 점을 가지고 있다. 즉, 관심대상이 되는 건축물이 실제 혹은 법규에 정한 것과 같은 조건이 되도록 조건을 조절해야 한다. 또한, 온도의 변화로 인한 불확실성과 비용, 시간 등 많은 현실적인 어려움이 있다. 그러므로 많은 연구자들은 외부의 조건을 정형화 할 수 있는 컴퓨터 시뮬레이션을 통한 연구를 진행해왔다. 단열에 취약한 창호를 보완하기 위한 발열유리의 열해석[4], 건물외피의 투과형 단열 벽체의 열성능 해석[5]의 연구에서 시뮬레이션은 실험을 대신하여 유용한 정보를 제공해줬다.

최근에는 복합재질로 구성된 구조체의 열해석[6, 7]과 같이 임의의 면에서 다른 물질이 만나는 복잡한 문제로 확장되었다. 여러 개의 층으로 이루어진 복잡한 형상에 대한 전열해석이 일반화 되면서, 모든 관심영역을 한꺼번에 격자로 생성하기에는 주어진 영역이 너무 복잡한 경우가 많아졌다. 이러한 복잡한 소재로 이루어진 건축물의 수치해석을 위해서는 전체 관심영역을 비교적 쉽게 격자를 생성할 수 있는 단순한 영역 혹은 재질별로 나누어 각 영역에 대하여 격자를 생성하는 방법을 많이 사용한다. 이러한 격자 생성은 빠른 시간에 전체 영역에 대한 격자를 손쉽게 생성할 수 있으나, 서로 인접하거나 만나는 독립된 영역들 간에 격자점들이 서로 일치하지 않는 비등각격자(non-conformal mesh)가 생성되며 이러한 면에서는 열에너지의 수수를 직접적으로 계산할 수 없게 된다. 이러한 불일치점을 해결하기 위해서는 양쪽격자사이에 가상의 격자(ghost cell)[8]을 사용하거나 불일치면에 격자를 서로 만나는 더 작은 면으로 나누는 클리핑(clipping)[9, 10] 방법이 있다. 가상격자를 사용하는 방법은 추가적인 기학학적인 정보를 필요로 하지 않으나 내삽하는 방법에 따라 해의 정확도가 바뀌는 단점과 근본적으로 해가 보존되지 않는다는 단점을 가지고 있다. 클리핑(clipping)방법은 복합물질의 접촉면의 격자를 양쪽이 서로 일치하는 더 작은 격자로 나누어 추가적인 기하학적 정보를 필요로 한다는 단점을 가지고 있으나 해가 근본적으로 보존된다는 장점을 가지고 있다.

본 연구에서는 복합재질로 구성된 전열문제에 대하여 각기 독립으로 격자를 생성하고 수치적 전열해석을 수행할 수 있도록 클리핑 방법을 구현하였다. 수치적 방법의 정확성을 검증하기 위하여 단순하지만 두면이 일치하지 않는 면의 격자에 대하여 본 연구에서 구현한 클리핑 방법을 적용하고 전열해석을 수행하였다. 향후 본 연구에서 구현한 방법을 좀 더 실제적인 형상을 지닌 창호 또는 단열재의 해석에 적용할 예정이다.

더 나아가, 이번 연구에서는 기존의 전열해석 연구와는 다르게 수치해석적인 알고리즘의 구현에도 목적을 두었다. 특히, 많은 연구자들이 사용하는 tetrahedron격자에 대해서 클리핑 알고리즘을 구현하였다. 본 연구는 전열해석 분야에서 수치해법을 연구하게 될 후학연구자들을 위한 기초연구에 도움이 될 것으로 생각된다.

비등각격자를 가진 복합재질 전열해석을 위해서는 서로 이웃한 경계면간에 열유속을 계산할 수 있도록 기학학적 정보를 필요로 한다. 즉 격자점이 서로 일치하지 않기 때문에 모든 경계면에 있는 격자면 들을 잘게 쪼개어 서로 면들이 일치하도록 해야 한다. 이러한 작업을 클리핑이라고 한다. 계산에 소요되는 시간을 최소화하기 위하여 quick reject algorithm을 이용하여 후보군을 정하고 이 후보군들 중에서 클리핑을 수행하게 된다. 클리핑이 끝이 나면 비등각격자면은 추가적인 작은 격자면 들로 나누어지며 이면은 이웃한 두 격자가 공유하는 면이 됨으로 열유속을 이용하여 전열계산을 수행할 수 있게 된다. 실제적으로 이러한 면들을 가진 이 해석 과정은 마치 다면체 격자를 이용한 해석과 같은 과정이 된다.

2장과 3장에서 클리핑을 위한 알고리즘과 열전달해석 방법을 설명하고, 클리핑과 열전달해석과정을 검증하기 위해 비등각격자면을 가진 형상에 대하여 열전도해석을 통한 시뮬레이션을 수행하였다.

컴퓨터 그래픽 분야에서의 클리핑은 윈도우(window)의 외부에 표시되지 않는 화상(image)의 일부를 없애는 과정이다. 즉, 출력하고자 하는 화면 내에 있는 형상만을 화면다각형을 이용해 제거하는 과정을 의미한다[9, 10]. 이러한 과정은 전산유체역학(CFD, Computational Fluid Dynamics)에서 격자점이 일치하지 않는 면을 공통의 작은 면으로 나누는 작업에 응용된다.

본 연구의 클리핑 방법은 사면체 격자의 삼각형 면에 제한함으로 quick reject algorithm과 격자점의 내부여부를 확인하는 winding number등을 단순화시킬 수 있었으며 계산시간 또한 줄일 수 있었다.

Quick reject algorithm은 같은 평면상에 놓인 삼각형 중에서 겹쳐질 가능성이 없는 면을 미리 제외시키는 알고리즘이다. Beaudoin and Jasak[9]은 사각형에 대하여 Quick reject algorithm의 한 종류인 Axis Aligned Bounding Bos(AABB)알고리즘을 이용하여 quick reject를 수행하였다. 하지만 Beaudoin and Jasak[2]은 직사각형에 대한 알고리즘을 사용하였다. 본 연구에서 고려한 격자는 사면체로, 비등각격자면에서 모양은 삼각형이기 때문에 변형된 방식으로 quick reject algorithm을 수행하였다. 개략도는 Fig. 1에 나타내었다. 우선 Fig. 1에서 polygon_P는 polygon_Q에 의하여 클리핑 되는 다각이라고 하면 polygon_P는 “clipee”라하고, polygon_Q는 “clipper” 혹은 “window”라고 정의한다. 본 연구에서는 clipper 대신 window라고 하겠다.

Fig. 1. 
Criterion of quick reject

Quick reject algorithm을 수행하기 위해 두 삼각형은 서로 같은 평면에 있다고 가정하였다. polygon_P와 window가 서로 겹치는지 확인할 필요가 있다. 우리는 polygon_P의 면적중심을 기준으로 하는 원을 이용하여 영역을 설정하였다. 이 영역 안에 window의 면적센터가 들어오게 되면 다음 계산단계로 넘어가는 방법을 수행하였다. 이 영역의 반지름은 다음과 같이 정의하였다.

여기서 k는 임의의 실수 상수이다. d_r은 polygon_P 모서리들의 산술평균 값이다. 본 연구에서 k값은 2~3사이의 값을 사용하였으며 이 값을 증가시킴에 따라 계산시간은 늘어나게 된다. (k값에 대한 보강, 테스트를 통해 3을 사용하였다.)이 값은 삼각형들의 크기에 따라서 자유롭게 바꿀 수 있지만, 효율적인 값은 여러 번의 계산을 통해서 경험적으로 정할수있다. 또한, 식(1)로부터 정해진 영역 안에 window 다각형이 존재한다고 해도 두 다각형이 겹쳐(overlap)지지 않을 수 있다. k값이 커질수록 겹쳐지지 않는 무의미한 삼각형의 개수가 늘어나기 때문에, k값의 선택은 계산의 효율에 중요하다.

Quick reject algorithm을 수행하면 후보군이 결정되며 이 후보군을 대상으로 winding number(ω)를 계산하게 된다.

Winding number는 임의의 점 주위를 곡선이 몇 바퀴 도는 지 나타내는 정수로서, 다음과 같이 정의가 된다.

여기서, γ는 곡선의 방정식, A는 임의의 점의 좌표, r은 점 A에서 곡선까지의 거리, 그리고 dφ는 극 좌표계에서의 미소 각도를 나타낸다.

Winding number는 클리핑 방법에 적용하면 window 삼각형 안에 임의의 점(clipee의 한 꼭지점)이 존재하는지 여부를 알 수 있게 해준다.[3] 만약 winding number의 값이 0이라면 점은 다각형 안에 존재하지 않게 된다.

클리핑 방법을 컴퓨터 프로그램의 알고리즘으로 구현하기 위해 Fig. 2과 같은 방법으로 winding number에 대한 식을 이산화 하였다. 점으로부터 각 꼭지점까지 선을 만들어서 세 개의 삼각형으로 나누었다. 각 각의 각도의 합이 2π와 같다면 winding number는 1 또는 –1의 값을 갖게 될 것이다. 하지만 다각형 밖에 존재할 경우 winding number는 주위의 곡선이 폐곡선을 이루지 못하므로 0의 값을 가지게 된다.

Fig. 2. 
Calculation of winding number (a) interior(ω≠0). (b) exterior(ω=0).

Fig. 3와 같이 임의의 두 선분이 같은 평면 위에 놓여있고, 교차를 하게 되면 교차점(I)을 만들게 된다. 이 교차점은 항상 클리핑 과정에서 만들어진 새로운 다각형의 점이 된다. 두 선분에 의한 교차점의 좌표는 벡터로 표현된 선분에 대한 식을 이용하여 얻을 수 있다. 알고리즘은 Table 1과 같다[4]. 여기서, α_P와 α_Q는 교차점 I의 위치를 나타내어 주는 비율값이고, <|>는 두 벡터의 내적을 의미한다. WEC는 Window Edge Coordinate의 약자로, 본 연구의 가상코드(pseudo code)에서 쓰인 변수의 이름이다.

Fig. 3. 
Calculation of intersection point

앞의 과정들을 이용하여 두 삼각형(window와 clipee)이 서로 겹쳐지는 다각형의 점들을 구할 수 있다. 점들의 좌표를 알기 때문에 면적을 구할 수는 있지만, 점의 순서가 불규칙하여 컴퓨터 계산상의 불편함이 따른다. 그래서 Fig. 4에서 볼 수 있듯이, 다각형을 이루는 점들에 대해, 다음 식과 같은 규칙성을 가지고 재배열 하였다.

Fig. 4. 
Rearrangement of the clipping points.

여기서 n→은 면의 법선벡터이고, V1→와 V2→는 점 C₁으로부터 다음 점들까지의 변위 벡터를 나타낸다. Fig. 5(a)에서는 C₁에서 점 C₂까지의 벡터를 V1→, 점 C₃까지의 벡터를 V2→라고 하였다. 그 다음 식(3)을 적용 시켰을 때 만족을 시키면 두 벡터는 변하지 않는다. 하지만 식(3)을 만족시키지 못할 경우에는 두 벡터를 바꿔서 점의 위치를 바꾸어 준다. 점의 위치를 바꾸는 예는 Fig. 5(b)와 같다. 같은 방법으로 모든 점들에 대해 적용을 시키게 되면 점의 원하는 순서로 나열 할 수 있다.

Fig. 5. 
Algorithm of the rearrangement.

재배열 된 점들은 삼각형의 면적을 구하는 식을 이용하여 다음과 같이 구할 수 있다.

여기서 N은 클리핑방법으로 얻어진 점의 개수이다. 컴퓨터 프로그램 상의 가상코드(pseudo code)는 Table 2와 같다.

앞에 설명한 전체 과정을 이용하여 임의의 모양을 가진 두 면에 대하여 검증을 진행하였다. 검증은 Fig. 6에 나와 있는 격자들을 이용하였는데, 임의의 법선 벡터를 가지고 있는 면에 대해 서로 다른 크기의 격자들로 나눈 것을 나타낸다. Fig. 6(a)와 (b)는 같은 면에 대한 격자들을 나타내고 있으며, 다른 크기의 삼각형들로 이루어져 있다. 식(1)에서의 k값은 3.0을 사용하였으며 격자들의 면의 개수와 node의 개수에 관한 자세한 내용은 Table 3에 나타내었다.

Fig. 6. 
Different cells on the same surface

클리핑을 수행하면 window도형들에 의해서 polygon_P는 여러 개의 조각으로 나뉘게 된다. 그 조각들의 면적을 구하여 합하게 되면 원래의 polygon_P의 도형의 면적과 같아진다는 조건을 이용하여 나누어진 방법에 대한 검증을 수행하였다. 상대오차는 다음의 식으로 계산하였다.

여기서A_exact는 원래그대로의 면적을 나타내고, A_clipping은 클리핑 방법을 이용하여 구한 면적들의 합을 나타낸다. 상대오차의 결과는 Fig. 7의 그래프로 나타내었다. Fig.7은 전체 계산영역에서 셀들이 가지고 있는 오차의 비율을 나타내고 있다.

Fig. 7. 
Relative error of the validation for the clipping method

Fig. 8(a)는 클리핑 되는 면의 격자 모양을 나타낸 것으로, polygon_P와 window들을 겹쳐서 보여주고 있다. 계산 영역의 모서리는 일치하지만 가운데 영역의 cell은 모서리가 서로 교차하면서 서로 다른 점들을 많이 만들어 내는 것을 Fig.8(b)와 (c)를 보면 알 수가 있다.

Fig. 8. 
(a) Windows and polygons_P on the test surface. (b) Windows and polygons_P in the black circle on the test surface. (c) Subfaces generated by the clipping method.

Fig. 9는 상대오차를 나타내었다. Fig.7은 셀들이 가지고 있는 오차의 비율을 나타내고 있는 반면, Fig.9에서는 각 각의 셀들이 가지고 있는 상대오차의 분포를 나타내고 있다. Fig.9를 통해 우리는 주요 오차발생위치와 원인을 파악하고자 하였고, 상대오차는 벽면 쪽이 상대적으로 높게 나오는 것을 확인할 수가 있었다. 벽면에서는 다각형 면들이 모서리에서 선분에 일치하게 되고 이는 winding number을 이용하여 교차점을 찾는데 어려움이 있으며 추가적인 작업을 필요로 하게 된다. 그러므로 이러한 영역에서 상대적으로 큰 오차를 만들어낸다는 것을 알 수 있다. 전체 계산영역에서의 계산된 상대오차의 최대값은 약 10^-4이고 최소값은 10^-5이다. 그러나 최대 오차가 O(10^-4)의 크기를 가지므로 전열해석에는 해의 정밀도를 떨어뜨리지 않을 것으로 예상된다.

Fig. 9. 
Contour of the relative error.

유한체적법(FVM)은 열전도 문제에서 흔히 쓰이는 수치기법이다. 유한제적법의 기본개념은 유한의 부피를 가지고 있는 체적으로 계산영역을 나누고 경계를 통과하는 열유속를 계산하는 것이다. 그러므로 충분히 수렴했다면 계산영역의 보존법칙은 만족되며 격자가 비교적 큰 경우 각 위치별로 정밀하게 표시하지는 못하지만 격자의 중심부에 해당하는 값은 정확하게 예측할 수 있수 있는 특징을 가지고 있다. 그러므로 비교적 작은 수의 격자를 사용해도 비교적 높은 정밀도를 유지하는 장점을 가지고 있다.

일반적으로 적분형태의 수송방정식은 종속변수를 φ라고 할 때, 다음과 같은 형태를 가진다.[11]

여기서 ρ는 밀도, φ는 어떤 값을 나타내는 종속변수, V는 검사체적의 부피, S는 검사체적의 면을 나타내는 검사표면, υ는 관측좌표계에서의 절대속도, υ_s는 검사표면의 속도 벡터,Γ_φ는 확산계수, Q는 원천(source)을 나타낸다. 좌변에서의 첫 번째 항은 시간에 의한 비정상항을 나타낸다.

만약, V와 S가 임의의 유한한 체적(V₀)과 그 체적의 표면 (S_j, j = 1, 2, …, n-1, n)을 나타낸다면 식(6)은 다음과 같이 쓸 수 있다.

위 식을 이산화하면 해석이 가능한 연립방정식을 얻을 수 있다.

FVM에서 서로 다른 물질이 만나는 interface에서는 열전도 계수와 비열을 정의할 수 없으며 단지 온도가 서로 일치해야 한다는 조건만이 있다. 일반격자에서는 격자의 비직교에 의한 교차확산을 고려해 주어야 한다.

여기서 q_b는 벽면에서의 열 플럭스를 나타낸다. S_b는 벽면에 있는 셀의 면적, n^은 면적의 단위법선벡터, h는 열전달계수와 기하학적 계수를 곱한 값 kSb/ds0⋅n^, T_b와 T_P는 각 각 셀의 경계면, 중심에서의 온도를 나타낸다. 식(8)와 (9)에서의 열 플럭스는 서로 반대이므로, 다음과 같이 나타낼 수 있다.

식 (11)를 이용하여 interface에서의 온도(T_b)를 구할 수 있다.

여기서 S₀과 S₁은 각 각 경계면에서 만나는 면적들을 나타내고 있다. 면의 온도 T_b를 온도 경계조건의 온도로 적용하였다.