Rayleigh-B□nard 대류 발생에 관한 수치모사

조홍제, 신치범1, 최창균*

서울대학교 응용화학부, 1아주대학교 화학공학과

*Corresponding Author 이메일:ckchoi@snu.ac.kr

Numerical simulation of the onset of Rayleigh-B□nard convection

H. J. Cho, C. B. Shin1, and C. K. Choi*

School of Chemical Engineering, Seoul National University

1Department of Chemical Engineering, Ajou University

*Corresponding Author e-mail:ckchoi@snu.ac.kr

서론

초기에 정지되어 있는 수평 유체층의 밑면을 가열하면 부력에 의한 자연대류가 발생한다. 자연대류 발생 시점은 Mahler 등[1]이 제시한 것과 같이 계가 불안정해지는 시간과 성장률을 이용해 해석할 수 있다. Choi 등[2, 3]이 제안한 전파이론(propagation theory)은 자연대류 발생에 관한 연구에 있어서 많은 계에 성공적으로 적용되어 왔다. 본 연구에서는 Prandtl수가 매우 큰 경우에 대하여 기존의 전파이론을 통해 얻어진 교란 분포를 초기조건으로 하여 주어진 계를 유한요소법(finite element method; FEM)으로 수치모사를 하였다. 이로부터 자연대류 발생 시점을 나타낼 수 있는 성장률을 새롭게 정의하여 본 수치모사의 타당성을 보이고, 실험결과와 비교하여 보았다.

이론

본 연구는 수직 방향의 거리가 L인 두 평판 사이에 밑면이 가열되는 비압축성 뉴튼유체를 고려하였다. 밑면의 온도는[eqn: { T}_{b }] , 윗면의 온도는[eqn: { T}_{i }] 라 하고, 점성소산에 의한 에너지 손실을 무시하였다. Boussinesq 가정을 도입하여 2차원 연속방정식,  운동량식과 에너지식을 무차원화하면 다음과 같다.

[eqn: {PARTIAL v} over {partial y }+{partial w} over{partial z}=0]                                                   (1)

[eqn: {partial v } over {partial tau }+ v{partial v } over {partial y}+ w{partial v } over {partial z }=-{partial p } over {partial y}+P r left ({{partial}^{2}v } over { {partial y}^{2}}+ {{partial}^{2}v } over {{partial z}^{2} } right)]                 (2)

[eqn: {partial w} over {partial tau }+ v{partial w } over {partial y}+ w{partial w } over {partial z }=-{partial p } over {partial z}+ P r left ({{partial}^{2}w } over { {partial y}^{2}}+ {{partial}^{2}w } over {partial z}^{2} } right) + P rRa theta]      (3)

[eqn: {partial theta } over {partial tau }+ v{partial theta } over {partial y}+ w{partial theta } over {partial z }= {{partial}^{2}theta } over {{partial y}^{2}}+ {{partial}^{2}theta } over {{partial z}^{2} }]                             (4)

위 식에서 v, w, [eqn: theta] 는 무차원화된 y방향 속도, z방향 속도, 온도이다. 길이는 유체층 깊이 L, 시간은 [eqn: { L}^{2 }/ alpha], 속도는 [eqn: alpha /L], 압력은 [eqn: {alpha}^{2}{rho}_{0}/{L}^{2}] 그리고 온도는 [eqn: DELTA T] (=[eqn: { T}_{b }-{ T}_{i}])으로 무차원화하였다. 중요한 무차원군은 Prandtl수 Pr(=[eqn: nu /alpha]), Rayleigh수 Ra(=[eqn: g beta DELTA T {L}^{3}/(alpha nu)])이다. 여기서 g, [eqn: alpha], [eqn: beta], [eqn: nu], [eqn: {rho}_{0}]는 중력가속도, 열확산계수, 부피팽창계수, 동점도, 밀도를 나타낸다.

주어진 계에 롤(roll) 형태의 셀(cell)이 생긴다는 가정 아래 유한요소법으로 수치모사를 하였다. 무차원화된 z방향의 길이는 1이고, y방향의 길이는 전파이론에서 구한 Ra에 따른 무차원 파수 ac를 이용하여 하나의 셀이 생기도록 [eqn: pi/]ac로 정하였다. 그리고 y방향으로 주기적인 셀의 모습이 생기므로 식 (1)-(4)를 풀기 위한 초기조건은  다음과 같이 정의하였다.

[eqn: ( theta(0,y,z),w(0,y,z) ) = ({theta}_{*}f(z),{w}_{*}g(z) )] cos(ac y)                     (5)

[eqn: v(0,y,z) = - {{w}_{*}}over{{a}_{c}}{d g(z)} over{dz}] sin(ac y)                         (6)

Fig. 1. Finite element meshes

and boundary conditions.

여기서 f(z)와 g(z)는 전파이론에서 얻은 대류 불안정성이 시작되는 시점에서의 z방향에 대한 온도, 속도 교란분포를 각각의 최대값으로 규격화한 것으로 z만의 함수이다. [eqn: {theta}_{*}],[eqn: {w}_{*}]는 온도, 속도 교란의 초기 크기이다. 주어진 문제를 풀기 위한 경계조건과 z방향으로는 36, y방향으로는 40인 [eqn: 36 times 40]인 메쉬(mesh)는 Figure 1에 나타나 있다. 식 (1)-(4)를 Hughes 등[4]이 제안한 predictor/multicorrector 알고리즘으로 풀었다.

결과 및 토론

위의 계를 해석하기 위해 평균장 근사를 도입하여 [eqn: theta=<theta>+ { theta}^{' }]으로 두었다. 여기서 [eqn: <theta>]은 [eqn: { INT _{ 0}^{pi/{a}_{c} }theta d y } / {INT _{ 0}^{pi/{a}_{c}} d y }] 으로 수평 방향의 평균온도를 의미한다. 따라서 [eqn: <theta>]은 [eqn: tau]와 z의 함수이며,[eqn: { theta}^{' }]은 이에 대한 국지적인 차이를 나타낸다. 또한 온도에 대한 제곱근평균제곱(rms=root mean square)값을 기준으로 성장률[eqn: { r}_{0 }],[eqn: { r}_{1 }]을 다음과 같이 정의한다.

[eqn: { r}_{0 }={1}over{<theta> }_{r ms } { d{<theta>}_{r ms}} over {d tau }], [eqn: { r}_{1 }={1}over{{theta}{'} }_{r ms } { d{theta {'}}_{r ms}} over {d tau }]                (7)

Figure 2에서 보면, 시간이 작을 때는 기본온도에 의한 성장률 [eqn: { r}_{0 }]이 교란온도에 의한 성장률 [eqn: { r}_{1 }]보다 크다가 [eqn: { r}_{0 } = { r}_{1 }]이 되는 시간이 존재한다. 즉 교란온도의 성장률이 기본온도의 성장률과 같게 되는 가장 빠른 시간을 [eqn: { tau}_{c }]로, 그리고 교란온도의 성장률이 최대가 되는 시간을[eqn: { tau}_{r1,max }]라 정의한다.[eqn: { tau}_{c }]이후의 시간에는 교란온도의 성장률이 기본온도에 의한 성장률보다 커지는 시간이므로 대류 불안정성이 발생하여 교란이 성장하고 있다고 하겠다. 따라서 [eqn: { tau}_{c }] 는 계의 대류 불안정성이 시작되는 시간이라 말할 수 있다. 이 그림은 초기온도 크기의 효과도 보이고 있는데, [eqn: Ra={ 10}^{5 }],[eqn: {w}_{*}]= 5일 때,[eqn: {theta}_{*}]=[eqn: { 10}^{-3 }] , [eqn: { 10}^{-4 }], [eqn: { 10}^{-5 }]로 가정하고 얻은 성장률 곡선을 나타낸다. 세 경우에 대하여 [eqn: { r}_{0 }]는 모두 일치하였고, [eqn: { r}_{1 }]은 [eqn: { tau}_{r1,max }]이전까지는 동일한 선상에 있었다. 이는 대류 불안정성이 시작되는 시간은 초기온도 크기에 관계없이 나타나는 불변의 시간이라는 것을 의미한다. 그리고 이 값을 Patrick과 Wragg[5], Inoue 등[6]의 실험에서 얻어진 Nu수가 전도 상태에서 벗어나기 시작하는 undershoot시간, 즉[eqn: {tau}_{u}]와 비교하였다. 그 결과 교란온도의 성장률이 최대가 되는 시간[eqn: { tau}_{r1,max }]는[eqn: {theta}_{*}]=[eqn: { 10}^{-3 }]일 때[eqn: {tau}_{u}]와 거의 일치하였다. 따라서[eqn: { tau}_{r1,max }]는 명백한 자연대류가 발생하는 시간이라고 생각된다. 시간이 무한히 흐르면 [eqn: { r}_{0 }], [eqn: { r}_{1 }] 모두 0에 가까이 간다는 것을 볼 수 있는데, 이는 시간이 많이 흐르면 계가 안정해져서 변화가 없는 쪽으로 진행된다는 것을 보여준다고 하겠다.

위의 결과를 근거로,[eqn: {theta}_{*}]=[eqn: { 10}^{-3 }]으로 하여 초기속도에 대한 영향을 살펴보았다. 그 영 향은 Figure 3에 나타나 있다.[eqn: {w}_{*}]=0일지라도 성장률[eqn: { r}_{0 }]과[eqn: { r}_{1 }]은 모두 일치하였다. 따라서 초기속도 크기에 대한 영향은 거의 없다고 하겠다. 본 연구에서는 [eqn: {w}_{*}]=5로 두고, Ra=1707.762에서[eqn: { 10}^{7 }]까지 계산하여 보았다.

일반적으로 Ra[eqn: >=]1707.762이면 자연대류가 발생한다고 알려져 있다.[eqn: { r}_{1}={r}_{0}]이 되는 시간은 Figure 4에서 구할 수 없다. 그리고 시간이 많이 경과한 후에는 두 곡선이 일치한다. 일치하는[eqn: { r}_{0 }] 과[eqn: { r}_{1 }]의 값이 거의 0이라는 사실을 통해 1707.762가 자연대류가 발생하는 임계 Ra수이며, 계에는 불안정한 시기가 나타나지 않는다고 할 수 있다.

Figure 5는 Patrick과 Wragg[5], Inoue 등[6]의 실험결과와 이를 비교한 그림이다. 이 그림에서 보면[eqn: { tau}_{r1,max }]는 실험에서 얻은 [eqn: { tau}_{u}]와 거의 일치한다. 이는 본 연구의 결과가 실제 현상을 잘 설명해준다는 것을 뒷받침한다. 그리고 위의 결과들을 요약하면 Table 1과 같다. Table 1에서는 Ra가 증가함에 따라[eqn: { tau}_{r1,max}]/[eqn: { tau}_{c }]값이 어느 일정한 값으로 수렴한다는 것을 알 수 있다. 즉, Ra가 큰 경우 [eqn: { tau}_{r1,max}]/[eqn: { tau}_{c }][eqn: SIMEQ]5.7이 된다. 이는 곧 [eqn: { tau}_{c }]에서 발생한 대류 불안정성이 가시화 되기까지 일정한 성장기간이 필요하다는 것을 의미한다. 좀 더 구체적으로 말하면, Figure 5의[eqn: { tau}_{c }] 곡선을 이 수렴값 만큼 곱하여 도시하면 [eqn: { tau}_{r1,max }]의 곡선을 얻을 수 있다.

이와 같이 본 연구에서 이용된 척도, 즉[eqn: {r}_{0}={r}_{1}]이라는 조건은 주어진 계의 자연대류 발생 시점에 관한 실험결과를 잘 설명해준다고 하겠다.

참고문헌

1. Mahler, E. G., Schechter, R. S. and Wisseler, E. H.: Phys. Fluids, 11, 1901 (1968).

2. Kim, M. C., Park, H. K., and Choi, C. K.: Theoret. Comput. Fluid Dynamics, 16, 49(2002).

3. Yang, D. J. and Choi, C. K.: Phys. Fluids, 14, 930(2002).

4. Hughes, T. J. R., Pister. K. S. and Taylor, R. L.: Comput. Meth. Appl. Mech. Engng., 17/18, 159(1979).

5. Patrick, M.A. and Wragg, A.A.: Int. J. Heat Mass Transfer, 18, 1397(1975).

6. Inoue, Y., Akutagawa, S., Saeki, S., and Ito, R.: Kagagu Kagagu Rhobunshu,      9, 359(1983).

Fig. 3. Effect of amplitude of initial velocity.

Fig. 2. Effect of amplitude of initial temperature.

Fig. 4. Growth rates at Ra=1707.762.

Fig. 5. Comparison of present results

with available experiments.

Table 1. Results of numerical simulation.

Ra	[eqn: {10 }^{4}]	[eqn: {10 }^{5}]	[eqn: {10 }^{6}]	[eqn: {10 }^{7}]
[eqn: { tau }_{c }] (A)	1.20[eqn: times{10 }^{-2}]	2.50[eqn: times{10 }^{-3}]	0.52[eqn: times{10 }^{-4}]	1.10[eqn: times{10 }^{-4}]
[eqn: { tau }_{r1,max }] (B)	6.80[eqn: times{10 }^{-2}]	1.37[eqn: times{10 }^{-2}]	2.97[eqn: times{10 }^{-3}]	6.30[eqn: times{10 }^{-4}]
B/A	5.67	5.48	5.71	5.73