Introduction

본 문제에서 다루고자 하는 파라미터 영역은 R_e≫1, P_r∼O(1), δ ≺ O(1), ϵ ≺ O(1) 이다. 즉, 유동에 의해 발생하는 경계층의 두께가 용기의 길이척도에 비해 매우 작은 경우이다. 따라서 용기 내부 대부분의 유동영역(bulk zone)을 비점성 영역으로 가정할 수 있다. 이론해석은 이 영역에 대해 과도상태에서 발생하는 온도변화를 예측하는데 초점을 맞출 것이다.

물리적인 관점에서 살펴보면 본 문제에서 다루는 유동 및 열전달 현상이, R_e≫1인 경우를 고려하면, 에너지방정식[식 (5)]으로 부터 구할 수 있는 열확산 시간척도, 즉, 긴시간변수(long time variable) t~=t/R_e 와 운동량방정식[식 (3)]의 우변에 나타나는 강제항 (forcing term) -2πδ(1-ϵT)cos(2πt)으로 부터 구할 수 있는 짧은 시간변수(short time variable) t 에 관계하는 이중시간척도 문제인 것을 알 수 있다.

따라서 지배방정식의 좌변에 나타나는 시간미분항은 다음과 같이 2-변수 전개법(two variable expansion method)으로 처리할 수 있다.

따라서 근사해를 구하기 위해 앞서 파라미터 ϵ(ϵ≪1)을 점근해석을 위한 매우 작은 섭동파라미터(perturbation parameter)로 이용하고 지배방정식 (2)~(5)와 초기 및 경계조건 (6), (7)을 만족하는 다음의 급수해를 생각하자.

식 (9a)~(9d)를 지배방정식 (2)~(7)에 대입하고, 각 차수의 해를 구하면 다음과 같다.

3.1.1 0-차해

경계조건(7a)와 식 (2) & (4)로부터

U0=V0=0,

(10)

을 얻을 수 있다. 사실 속도장에 대한 0-차해는 밀폐된 용기를 채우고 있는 밀도 일정인 비압축성 유체의 운동이 고체의 운동과 같다는 물리적인 사실로부터도 쉽게 유도할 수 있다.

압력에 대한 해를 구하기 위해 식 (10)을 식 (3)에 대입하면 P₀에 대한 방정식

∂P0∂x=-2πδcos⁡(2πt)

(11)

이 구해진다. 따라서 P₀에 대한

P0=-2πδcos⁡(2πt)x

(12)

의 해를 구할 수 있다.

온도에 대한 0-차 방정식과 초기 및 경계조건은 식 (5), (6) & (7b)로부터

∂T0∂t~=1Pr∇2T0,

(13)

T0(t~=0,x,y)=0,T0t~,x2+y2=1=-1,

(14)

와 같이 얻을 수 있고, 식 (13)으로 나타나는 열전도방정식에 대한 해를 구하면

T0=-1+2∑n=1∞ exp⁡-βn2Prt~J0x2+y2βnβnJ1βn

(15)

의 고전적인 해를 얻는다.⁽⁷⁾ 여기서 β_n은 0-차 Bessel 함수의 n번째 0-차 해이다.

이후의 그림에 나오는 모든 급수해 계산은 Mathematica를 이용해 수행했고 급수의 30번째 항까지 덧셈을 수행한 결과임을 밝혀둔다. 그리고 Prandtl 수는 15℃, 1기압 물을 기준으로 P_r = 8.0으로 고정했다.

식 (15)에 주어진 0-차 해(T₀) 는 예상한 것처럼 긴시간변수 t~에만 의존하는 정지한 유동에 대한 단순 열전도 방정식의 해이다. 시간이 지남에 따라 축대칭 형태의 냉각이 내부로 서서히 진행되어가는 것을 Fig. 3에서 확인할 수 있다.

Fig. 3 
Behavior of the 0-th order temperature T₀. (P_r = 8.0)

(a) Isothermal contours at T₀ = 0.0, 0.1, 0.2, ∙∙∙, 0.9 . (b) Temperature distributions along radial direction (r=x2+y2) when t~ = 0.05, 0.1, 0.3, 0.6, 0.9 , 1.2.

3.1.2 1-차해

R_e≫1인 경우 식 (3) & (4)의 우변의 두 번째 항은 벽 근처의 경계층을 제외한 대부분의 유동영역에서 무시해도 좋을 만큼 작다. 그러나 이 항은 벽면근처의 O(R_e^-1/2)-두께의 경계층에서는 그 크기가 다른 항과 비슷해 무시할 수 없고, 수학적인 관점에서 볼 때 이 항이 미분의 최고차항이기 때문에 벽면경계조건을 위해 필요하다. 그러나 본 연구에서는 유동영역의 대부분을 차지하는 내부 비점성영역에서의 온도장에 대해 관심이 있기 때문에 이후의 해석에서 점성항을 생략하고, 경계층 영역에 대한 자세한 해석은 생략한다.

식 (8) & (9)를 식 (3) & (4)에 대입하여 1-차 방정식을 구하면

∂U1∂t=-∂P1∂x+2πδT0cos⁡2πt,

(16)

∂V1∂t=-∂P1∂y,

(17)

을 얻을 수 있다.

식 (16) & (17)과 식 (2)로부터 P₁에 대한 방정식을 구하면,

∇2P1=2πδ∂T0∂xcos⁡(2πt)

(18)

을 얻을 수 있다. 식 (18)의 해 P₁의 해석해를 직접 구하기는 힘들고 수치 해석적인 방법이 필요하다. 그러나 여기서는 식 (18)의 해를 직접적으로 구하는 대신 다른 이론적인 접근을 통해 Park⁽²⁾이 발견한 열전달 촉진 메카니즘을 설명할 수 있는 논리적 근거를 찾는 형태로 우회한다.

x-축(y=0)과 y-축(x=0)을 따라 식 (18)의 해를 구하면 다음과 같다.

(1) x-축(y =0) 상의 P₁ 해

∂P1(x,0)∂x=2πδT0(x,0)cos⁡(2πt),

(19a)

(2) y-축(x =0) 상의 P₁ 해

P10,y=0.

(19b)

이제 식 (19a)와 (19b)를 식 (16) 우변의 첫째항에 넣고 U₁을 계산하면

(1) x-축(y =0) 위의 해

U10,y=0

(20a)

(2) y-축(x =0) 위의 해

U10,y=δT00,ysin⁡2πt,

(20b)

을 얻을 수 있다.

식 (20a)와 (20b)이 주는 물리적 의미는 다음과 같다. 무한 수직벽 사이에 갖힌 유체인 경우[상기 (1)의 경우]는 용기의 좌우교반과 x-축방향 온도 기울기가 있어도 가속도가 만드는 체적력과 온도변화에 의해 만들어지는 압력구배가 상쇄되어 유속이 발생하지 않는다. 반면, 무한 수평벽 사이에 갖힌 유체의 경우[상기 (2)의 경우]는 온도 기울기가 미미해 압력구배가 0에 가깝고 결과적으로 식 (16)의 우변 두 번째 항, 즉, 강제진동에 의한 체적력만 살아남아 식 (20b)와 같은 유동을 발생시킨다.

위 결과와 x-축에서 멀어질수록 ∂T₀/∂x 값이 급격히 작아지는 T₀에 대한 등온선 그림[Fig. 3(a)]으로부터 x-축에서 먼 유동영역에서는 식 (16)의 우변 두 항의 크기가 ∂P1∂x≪2πδcos⁡(2πt) 인 것을 알 수 있다.

위 사실로부터 식 (16) 우변의 첫 번째 항을 무시하면 x-축에서 먼 점(0,1) 근처의 내부 비점성영역에 대한 수평방향 속도성분 U₁ 구하면 다음과 같다.

U1=δT0sin⁡2πt.

(21)

그리고 식 (17)로부터

V1=0.

(22)

를 얻을 수 있다.

결론적으로 점(0,1) 근처의 내부 비점성영역에서의 수평방향 속도성분 U₁은 외부에서 가해지는 가진가속도와 90⁰의 위상차를 갖고 수평방향으로 주기적인 운동을 한다는 것을 알 수 있다. 물론 이와 같은 상대운동이 가능한 이유는 유체의 온도변화에 의한 지엽적인 밀도변화 (local density variation)가 발생했기 때문이다 (이 효과는 식 (21)의 해에서 T₀에 반영되어 있다).

내부 비점성영역에서 온도장에 대한 1-차 방정식은 식 (9)를 식 (5)에 대입하여 아래와 같이 얻을 수 있다.

∂T1∂t+U1∂T0∂x=0.

(23)

식 (21)를 식 (23)에 대입하고 해를 구하면

T1=δ2πcos⁡2πtT0∂T0∂x,

(24)

을 얻을 수 있다.

앞서 구한 1-차해를 살펴보면 벽면 온도변화가 전도 열전달을 통해 (0-차 해(T₀)를 통해) 내부유체에서 축대칭 원형 형태의 온도변화를 일으키고[Fig. 3(a)], 이것이 다시 가진방향(x-축방향)으로 밀도구배(density gradient)를 일으켜 동시에 가진에 의해 발생하는 체적력의 지엽적인 변화(local variation)를 발생시키고, 궁극적으로 유동(U₁)을 발생시킨다는 것을 알 수 있다.

내부 영역에서는 발생하는 유동(U₁)의 크기는 식(21)에서 볼 수 있듯이 0-차 온도(T₀)에 비례한다. T₀ 는 변수 x 에 대해 우함수이고(y-축을 기준으로 좌우의 유체가 동일한 방향으로 움직인다는 의미) 시간에 따라 싸인함수로 형태의 운동을 한다는 점을 고려하면, 시간평균 관점에서 1-차 해 U₁에 의해 발생하는 유동이 열전달을 직접 촉진하지는 않는다는 것을 알 수 있다.

위 사실은 온도에 대한 1-차해[식 (24)]에서 다시 확인할 수 있다. 식 (24) 우변에 나타나는 합수 값(T₀∙∂T₀/∂x)의 분포를 살펴보면 더욱 분명해진다. Fig. 4에서 알 수 있듯이 식 (24)에서 구한 T₁은 시간평균 뿐만 아니라 공간평균 값도 0 인 함수이다. 결론적으로 주기운동에 의해 발생하는 1-차 해는 그 자체로는 열전달 촉진에 기여하지 않는다.

Fig. 4 
Distribution of the function, ∂∂xT02∂T0∂x, value near the point (0,1) at t~ = 0.05. P_r = 8.0.

앞서 구한 이론해를 이용하여 1-cycle 가진에 대한 온도 평균값을 구하자.

외부가진 1-cycle 동안에 대한 온도의 시간평균값은 다음과 같이 주어진다.

여기서 시간평균값 정의는 ϕ¯=∫tt+1 ϕ(t,t~,x,y)dt이고, ϕ는 온도에 대한 각 차수해를 의미한다.

각 차수의 온도에 대한 평균값을 구하면 식 (15)로부터 T0¯=T₀ 이고, 식 (24)로부터 T1¯=0 그리고 식 (28)로부터 T2¯=∫01 T2(t,t~,x,y)dt=-δ216π2∂∂xT02∂T0∂x 이다.

따라서

식 (30) 우변의 첫 번째 항은 단순 열전도방정식의 해이고, 두 번째 항은 수평가진에 의해 발생한 온도변화 효과이다. Fig. 5로 부터 우변 두 번째 항(ϵ2T2¯)의 값이 +값이라는 것을알 수 있다.

Fig. 5 
Distribution of the function, ∂∂xT02∂T0∂x, value near the point (0,1) at t~ = 0.05. P_r = 8.0.

이것은 수평방향으로 교반가진을 가하는 경우 (0,1) 근처 비점성영역에 있는 내부유체는 단순 열전도문제에 비해 냉각이 느려진다는 것을 의미한다. 그 이유는 온도변화로 인해 발생하는 수평방향 속도성분 U₁, U₂ 때문이다. 경계층유동에서 벽면온도변화가 경계층을 가로질러 내부비점성영역 유체로 쉽게 전달되지 않는 경우와 같다. 벽면의 온도변화(냉각)에 의해 일어나는 유체와 벽사이의 열전달 통로가 수평방향 대류의 영향으로 벽에 평행한 방향을 따라 이루어지기 때문에 나타난 현상이다.

Figs. 4 & 5에서 비록 0.5≤y≤1.0 영역에 대한 함수값을 그렸지만 이전의 이론전개 과정을 고려하면 여기서 구한 이론해의 엄밀성은 (0,1) 근처 영역에서 잘 맞는다는 점을 다시 한 번 강조한다.

이제 지금까지 구한 이론해 U₁, U₂ 을 이용해 수치해석 결과(Fig. 2)에서 알 수 있는 점 (0,1)에서 먼 x-축에 가까운 내부 비점성영역 유체가 겪는 수평방향 열전달 촉진 현상을 설명하고자 한다.

수평가진이 가해진 경우 수평방향으로 열전달이 전도 열전달에 비해 더 빨리 일어나는 현상을 설명한다.

식 (24)와 식 (26)에 주어진 U₁, U₂를 다시 쓰면

이다. 외부가진이 1-cycle 일어나는 동안을 짧은 시간 변수 t 관점에서 살펴보면, 식(31)의 싸인함수 앞에 있는 계수함수 δT₀와 δ24πT0∂T0∂x는 준정상상태이다. 따라서 계수함수 δT₀와 δ24πT0∂T0∂x는 고정된 t~ 값에 대해 공간변수 x, y 만의 함수로 생각할 수 있다.

계수함수의 부호는 x-축을 기준으로 유동영역 좌, 우에서 아래와 같다:

(1) -1≤x <0 에서는 δT₀ <0, δ24πT0∂T0∂x<0

(2) 0 <x≤1에서는 δT₀ <0, δ24πT0∂T0∂x>0

위 결과는 x-축을 기준으로 좌측영역(-1≤x<0)에서는 속도성분 U₁, U₂가 동상(in-phase)으로 움직이고, 우측영역(0<x≤1)에서는 180도 위상차이가 나는 이상(out of phase)으로 움직인다는 것을 보여준다.

수평방향유속이 두 속도의 합산, 즉, ϵU₁+ϵ²U₂ 로 계산됨을 고려하고, 섭동파라미터의 크기가 적절히 작아 두 항의 크기가 비교될 정도일 때를 가정하면, 점(0,1) 근처 영역에서의 유동은 두 조화함수 sin(2πt)와 sin(4πt)의 (영역 -1≤x<0에서는) 합으로, (영역 -1≤x<0에서는) 차로 나타나는 특성을 갖는다는 것을 알 수 있다.

Fig. 6으로부터 두 함수 sin(2πt)와 sin(4πt)합과 차의 크기의 최대치가 다음 구간에서 일어난다는 것을 알 수 있다.

Fig. 6 
Time history of two function arithmatic

첫째는 동상인 경우이다. 이 경우 시간구간 0≤t≤0.25와 0.75≤t≤1.0에서 최대값이 나타난다. 이때 용기에 가해지는 가속도의 방향은 오른쪽이다. 식 (31)을 이용하여 수평유속 U₁, U₂을 살펴보면 유동영역 왼쪽 (-1≤x<0)에서는 서로 같은 방향(- 방향)이고, 유동영역 오른쪽 (0<x≤1)에서는 서로 반대방향(U₁은 - 방향, U₂는 + 방향)이다. 다시 말해 시간구간 0≤t≤0.25와 0.75≤t≤1.0에서는 유동이 좌측영역(-1≤x<0)에서 활성화된다는 것을 의미한다.

두 번째는 이상인 경우이다. 이때는 최대치가 시간구간 0.25≤t≤0.75, 즉, 용기에 가해지는 가속도의 방향이 왼쪽을 향하는 시간 동안에 일어난다. 앞서의 설명과 동일한 이유로 시간구간 0.25≤t≤0.75 동안에는 유동이 우측영역(0<x≤1)에서 활성화된다. 위 설명은 x-축에서 먼 점(0,1) 근처의 경계층 밖 비점성영역에 대한 것이라는 사실을 다시 한 번 강조한다.

이상의 이론적인 분석을 바탕으로 전 유동영역에서 일어나는 유동특성에 대해 지금부터 살펴보자.

가진에 주기운동을 하는 용기가 오른쪽으로 가속되는 시간은 0≤t≤0.25와 0.75≤t≤1.0 이다. 이때 벽근처의 경계층 내부에서는 벽면냉각에 의해 내부유체보다 밀도가 커진 유체가 더 큰 관성력을 받아 용기가 가속되는 반대방향(왼쪽, 즉, -방향)을 향하게 된다. 한편 경계층 밖 비점성영역에서는 유동이 시간구간 0≤t≤0.25에서는 왼쪽 방향(- 방향), 시간구간 0.75≤t≤1.0에서는 오른쪽 방향(+방향)을 향한다 (이것은 앞서 구한 이론해로부터 쉽게 확인할 수있다).

따라서 시간구간 0≤t≤0.25에서는 경계층 유동의 방향과 경계측 밖 비점성 영역 유동의 방향이 동일한 병류(cocurrent flow)가 형성되어, Fig. 7(a)에서와 같은 매우 큰 순환유동이 만들어진다. 상반부를 기준으로 설명하면 반시계 방향의 순환유동이 만들어지고 순환류의 중심을 좌측에 존재한다 (지금 논의하는 시간구간에서는 유동이 활성화되는 영역이 -1≤x<0이라는 것을 상기하자).

Fig. 7 
Flow patterns : Schematics (upper frame) and Streamline of numerical result when A^* = 0.025 , f = 0.1 & (a)t = 0.2, (b)t = 0.7(lower frame).

반면 시간구간 0.75≤t≤1.0일 때는 경계층 유동은 – 방향을 향하고, 경계층 밖 비점성 영역의 유동은 +방향을 향하는 유동방향이 서로 반대인 countcurrent flow가 형성되어 전체적인 유동강도를 감소시켜 내부 전영역에 걸친 순환류를 만들지 않게 된다.

용기가 왼쪽으로 가속되는 시간 0.25≤t≤0.75 동안에는 우측 유동영역 0<x≤1에서 유동이 활성화된다. 이때 우측 유동영역에서도 앞서와 같은 이유로 시간구간 0.25≤t≤0.5 에서는 경계층유동과 경계층 밖 비점성 유동이 서로 반대인 countcurrent flow가 형성되어 내부영역의 유동이 약화되고, 시간구간 0.5≤t≤0.75 에서는 병류(cocurrent flow)가 발생하여 Fig. 7(b)에서와 같은 상반부를 기준으로 매우 큰 시계방향 순환유동이 만들어지고 순환유동의 중심은 우측에 존재하게 된다[Fig. 6(b)].

결과적으로 1 주기 동안 번갈아 발생하는 두 순환유동에 의한 대류효과가 x-축 방향으로의 열전달을 촉진해서 Fig. 2에서와 같이 가진방향으로 냉각을 촉진하게 된다.

이상의 이론적인 설명은 수치계산 결과와 비교한 Fig. 6에서 명확히 확인할 수 있다.