Contaminant spreading by natural convection in a box

Vietnam Journal of Mechanics, VAST, Vol. 38, No. 2 (2016), pp. 141 – 152 DOI:10.15625/0866-7136/38/2/7521 CONTAMINANT SPREADING BY NATURAL CONVECTION IN A BOX Tran Van Tran1,∗, Nguyen Ngoc Thang2, Nguyen Thi Thuy1 1VNU University of Science, Hanoi, Vietnam 2University of Fire Fighting and Prevention, Hanoi, Vietnam ∗E-mail: trantv@vnu.edu.vn Received December 11, 2015 Abstract. In this paper the spreading of a contaminant accompanied with natural con- vection in a box is numerically simulated. The box may be considered as a cooking room or a working place where some sources of heat and contaminant are in the simultane- ous action. The box floor is supposed to be divided into several domains with different boundary conditions for temperature or heat flux. Here the purpose of the simulation is to understand the contaminant spreading process in the box under the influence of a convective motion. The model can be also applied for an enclosure with separated parts differentially heated by the sunlight on its boundaries. A good knowledge of this process is very useful for setting an efficient ventilation scheme. In this paper the finite difference method based on the Samarski scheme with ADI technique is applied for numerical sim- ulation. Here the box floor is divided into two domains of equal sizes but with different temperature or heat flux. The contaminant source locates in the middle of the box bottom. The simulation shows that over the part of the floor where temperature or heat flux is greater the contaminant concentration is lager. That result is in the accordance with the experiment done in the framework of this investigation. Keywords: Box, three-dimensional, natural convection, contaminant spreading, numerical simulation, finite difference method, Samarski scheme, ADI, multigrid. 1. INTRODUCTION The spreading process of a contaminant accompanied by a convective motion in the air is a much known phenomenon that often occurs almost everywhere. The con- taminant source may be “independent” on the heat source or they are merged together. Combustion reactions are often of the second case. In an industrial enclosure heat sources usually also issue one or several contaminant matters simultaneously. The heat convec- tive motion as expected always makes the contaminant spreading more quick. But it is not its unique effect. In the case when the temperature or heat flux is not homogeneous c© 2016 Vietnam Academy of Science and Technology 142 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy on boundaries, in our case on the floor of the box, the contaminant distribution in the enclosure is strongly influenced too. This issue is the primary interest of this paper. Natural convection in an enclosure caused by heating from below or the difference in temperature of the side walls has been theoretically and experimentally investigated intensively from several decades ago. This problem has been attractive for the theoreti- cal investigation as well as application. One of the earliest studies of three-dimensional natural convection in a box with differential side heating by numerical simulation was carried out by Millinson and De Vahl Davis [1]. They revealed the steady air motion for moderate Rayleigh number (Ra). And the motion is essentially three-dimensional. In [2] the same problem was considered for Ra ranged from 103 to 2.1016. The laminar flow was observed again at not very large Ra. This problem was also solved experimentally in [3] for Ra from 104 to 2.107, and numerically by finite difference method in [4] for Ra not exceeded 106. It is interesting to note that in [5] the transition from the steady flow to the time-periodical natural convective motion in a box was observed. The natural convection considered in [6] is different from that of [1–5] by the heating way. Namely, in [6]the box is heated from below and the Rayleigh number is from 3500 to 104. Four different stable convective structures were recognized. Orhan Aydin and Wen-jei Yang [7] studied nat- ural convection in a two-dimensional rectangular enclosure with localized heating from below and symmetrical cooling from the sides. Four dimensionless heat source length of 1/5, 2/5, 3/5, 4/5 were taken for numerical simulation at Ra from 103 to 106. Recently natural convection of nanofluids has been investigated [8]. In this paper we consider the natural convection in a box caused by the non- homogeneity of the temperature or heat flux applied to different parts of the box bottom. And the way a contaminant spreads in the box in the presence of this convective motion is our main interest. 2. THE PROBLEM FORMULATION The numerical simulation in this investigation will be carried out on the base of the Boussinesq approximation of the Navier-Stokes equations. The reason for this choice is that, first we want to understand the influence of the temperature non-homogeneity at the floor on the contaminant spreading process at moderate Rayleigh number. Sec- ond, our computational facilities are beneath computing the solution of the full three- dimensional turbulent problem. So with this assumption the system of equations of our problem consists of the following [9, 10] ∂~V ∂t + ( ~V.∇ ) ~V = − 1 ρ0 ∇p + ν∆~V + (gβT′ + gβcC′)~k div~V = 0 ∂T′ ∂t + ( ~V.∇ ) T′ = χ∆T′ ∂C′ ∂t + ( ~V.∇ ) C′ = χc∆C′ where ρ0, ν are the air density and viscosity respectively at a referenced temperature T0. T′, C′ are the deflections of the temperature and contaminant from their referenced value Contaminant spreading by natural convection in a box 143 respectively, and β, βc are the temperature and contaminant coefficients of the expan- sion whilst χ,χc are their diffusivity coefficients respectively. Making the usual non- dimensionalizing procedure for the natural convection problem [10] we have ∂u ∂t + u ∂u ∂x + v ∂u ∂y + w ∂u ∂z = −∂p ∂x + ∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 , (1) ∂v ∂t + u ∂v ∂x + v ∂v ∂y + w ∂v ∂z = −∂p ∂y + ∂2v ∂x2 + ∂2v ∂y2 + ∂2v ∂z2 , (2) ∂w ∂t + u ∂w ∂x + v ∂w ∂y + w ∂w ∂z = −∂p ∂z + ∂2w ∂x2 + ∂2w ∂y2 + ∂2w ∂z2 + Gr T + Grc C, (3) ∂u ∂x + ∂v ∂y + ∂w ∂z = 0, (4) ∂T ∂t + u ∂T ∂x + v ∂T ∂y + w ∂T ∂z = 1 Pr ( ∂2T ∂x2 + ∂2T ∂y2 + ∂2T ∂z2 ) , (5) ∂C ∂t + u ∂C ∂x + v ∂C ∂y + w ∂C ∂z = 1 Sc ( ∂2C ∂x2 + ∂2C ∂y2 + ∂2C ∂z2 ) , (6) where Gr, Pr are the Grashof and Prandtl number for the temperature effect respectively whilst Grc, Sc are the similar numbers for the contaminant. The boundary conditions for Eqs. (1)-(6) will be discussed later. Here for saving the computer memory as well as reducing the computation time the finite difference method is applied. To do this, first we rewrite the governing equations (1), (2) and (3) in the vorticity variable introducing the vorticity vector by the form ~Ω (Ωx,Ωy,Ωz) = rot~V, Ωx = ∂w ∂y − ∂v ∂z , Ωy = ∂u ∂z − ∂w ∂x , Ωz = ∂v ∂x − ∂u ∂y . (7) By differentiating (1), (2) and (3) in respect to appropriate spatial variables and taking account of (4) we have Ωxt +uΩ x x+vΩ x y+wΩ x z−Ωxux−Ωyuy−Ωzuz− ( Ωxxx+Ω x yy+Ω x zz ) =GrTy + GrcCy, (8) Ωyt +uΩ y x+vΩ y y+wΩ y z−Ωxvx−Ωyvy−Ωzvz− ( Ωyxx+Ω y yy+Ω y zz ) =−GrTx−GrcCx, (9) Ωzt + uΩ z x + vΩ z y + wΩ z z −Ωxwx −Ωywy −Ωzwz − ( Ωzxx +Ω z yy +Ω z zz ) = 0. (10) Next by the differentiation and transformation of (4) one can get the following equations for the velocity components uxx + uyy + uzz = Ω y z −Ωzy, (11) vxx + vyy + vzz = Ωzx −Ωxz , (12) wxx + wyy + wzz = Ωxy −Ωyx. (13) Note that in Eqs. (8)-(13) the subscript letter means the derivative in respect to the corresponding variable. These equations will be integrated for determining the vector of vorticity and velocity of our problem. In Fig. 1 shown the enclosure considered in this 144 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy paper. It is a box of length L, so after the non-dimensionlizing procedure it becomes the unit box. We denote the sides of the box by ST, SB, SL, SR, SF, SK for top, bottom, left, right, front and back respectively. Now the boundary conditions for the system of Eqs. (5)-(13) are set as follows. x y z S1 S2 S3 S4 S4 Fig. 1. The box and its base divided into domains with different boundary conditions for temperature and contaminant u = v = w|all sides = 0, (14) Ωx|Sl ,SR = 0, Ωx|SF ,SB = wy, Ωx|SB,ST = −vz, (15) Ωy|SF,SK = 0, Ωy|SL,SR = −wx, Ωy|SB,ST = uz, (16) Ωz|SB,ST = 0, Ωz|SL,SR = vx, Ωz|SF,SK = −uy, (17) Tn = Cn|SF,SK ,SL,SR,ST = 0. (18) To study the influence of the heat non-homogeneity at the bottom side on the con- taminant spreading we divide the box base into several parts along the axis 0x (Fig. 1). The area ratio of S1, S2, S3, and S4 to the base are 1/16, 3/8, 3/8 and 3/16, respectively. The contaminant source occupies whole central part S1. On the other parts we impose the following conditions Tn|S2 = −1, T|S3,S4 = 0, (19) Cn|S2,S3,S4 = 0, Cn|S1 = −1, (20) T|S1 = 1, (21) T|S1 = 0. (22) Now we denote the problem that consists of (5)-(20) and (21) by A-problem whilst the problem consists of (5)-(20) and (22) by B-problem. Condition (19) shows that S2 is the hot domain while S3, and S4 are relatively cool ones. Condition (20) means that the contaminant is issued only from domain S1. Condition (21) indicates that S1 is a dual source i.e. heat and contaminant simultaneously originate from this domain while (22) shows that S1 is a monotonous source of contaminant only. Contaminant spreading by natural convection in a box 145 3. THE NUMERICAL METHOD To integrate the transport equations for vortices (9)-(11) here we apply the ADI and time splitting technique for the finite difference method based on the Samarski scheme. We describe this numerical procedure in detail for Eq. (9) as follows. We split every time step of the integration (τ) into three substeps. At the first substep we integrate Eq. (9) with all derivatives of Ωx in respect to variable x in the left part (x-direction) (Ωx)n+1/3−(Ωx)n τ/3 +u (Ωx)n+1/3x˜ −0.5 |u| k1{Ωxx −Ωxx¯}n+1/3−a (Ωx)n+1/3xx¯ −ux(Ωx)n+1/3 = = −v (Ωx)ny˜ − w (Ωx)nz˜ + uy (Ωy)n + uz (Ωz)n + { Ωxyy¯ +Ω x zz¯ }n + GrTny + GrcC n y , (23) At the second substep the equation in y-direction is integrated (Ωx)n+2/3−(Ωx)n+1/3 τ/3 +v (Ωx)n+2/3y˜ −0.5 |v| k1 { Ωxy−Ωxy¯ }n+2/3−b(Ωx)n+2/3yy¯ −ux(Ωx)n+2/3 = =−u (Ωx)n+1/3x˜ −w (Ωx)n+1/3z˜ +uy (Ωy)n+uz(Ωz)n+{Ωxxx¯+Ωxzz¯}n+1/3+GrTny +GrcCny , (24) And finally, at the third substep the equation in z-direction is solved (Ωx)n+1−(Ωx)n+2/3 τ/3 +w (Ωx)n+1z˜ −0.5 |w| k1{Ωxz−Ωxz¯}n+1−c (Ωx)n+1zz¯ −ux(Ωx)n+1 = =−u (Ωx)n+2/3x˜ −v (Ωx)n+2/3y˜ +uy (Ωy)n+uz(Ωz)n+ { Ωxxx¯+Ω x yy¯ }n+2/3 +GrTny +GrcC n y , (25) In (23)-(25) we denote (ux)ijk = ( ui+1,jk − uijk ) /hi, (ux¯)ijk = ( uijk − ui−1,jk ) /hi−1, (ux˜)ijk = 0.5 (ux + ux¯) , (uxx¯)ijk = 2 (ux − ux¯) / (hi−1 + hi) , hi = xi+1 − xx, a = 1/ ( 1 + ∣∣uijk∣∣ k2hiRe) , b = 1/ ( 1 + ∣∣vijk∣∣ k2hiRe) , 1/ (1 + ∣∣wijk∣∣ k2hiRe) . (26) The analogous procedure is applied for integrating equations (10) and (11). It is obviously from (26) that the central scheme for both the first and second derivatives re- sults in the case when k1 = k2 = 0 while in the case with k1 = 1, k2 = 0 one has the so called upwind scheme. Finally we have the Samarski scheme [10, 11] taking k1 = k2 = 1. In this paper the last scheme is applied. To calculate the solution of (11)-(13) the second order central finite difference scheme is used for the Laplace operator and the multigrid method [12–14] is applied to solve the system of finite difference equations. The multigrid method extremely reduces the com- putational time of our numerical simulation. 146 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy 4. NUMERICAL RESULTS AND DISCUSSION The numerical simulation in this paper is carried out for the case when the contaminant-carbon dioxide spreads by the natural convection in a box filled with the air. The Grashof number (Gr) is from 104 to 5.105, and Grc is fixed at value 105. The common fixture that our simulation shares with the other studies mentioned in Section 1 of this paper is the existence of the steady air flow at moderate Rayleigh number de- spite the fact that the boundary condition for heat in our case is much different from that imposed in the mentioned works. Moreover the natural convective motion in our study is also interacted with the contaminant spreading process. The recording flow param- eters at a series of time moments in three points P1(0.5, 0.5, 0.175), P2(0.5, 0.5, 0.5) and P3(0.5, 0.5, 0.875) helps to determine the kind of the motion: steady or unsteady. In Fig. 2 shown the change in time of the velocity, temperature and contaminant concentration in these points. It is obviously that the flow in both cases becomes stationary after a relatively short time interval. -15.0 -10.0 -5.0 0.0 5.0 10.0 15.0 20.0 0.0 0.5 1.0 1.5 2.0 u1 u2 u3 v1 v2 v3 w1 w2 w3 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.5 1.0 1.5 2.0 T1 T2 T3 10C1 10C3 10C2 (a) A-problem -10.0 -8.0 -6.0 -4.0 -2.0 0.0 2.0 4.0 6.0 8.0 0.0 0.5 1.0 1.5 2.0 V1 w3 w2 v2 w1 v3 u1 u2 u3 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.0 0.5 1.0 1.5 2.0 100C1 T1 T2 T3 100C2 100C3 (b) B-problem Fig. 2. The evolution of velocity components u, v, w, temperature T and contaminant C in points P1, P2 and P3 at Gr = 2. 104, Grc = 105 The effect of the boundary condition (21) and (22) on the amount of the contami- nant released from the source S1 are clearly shown in Fig. 2. This amount in case of the Contaminant spreading by natural convection in a box 147 A-problem (S1 is the source of heat and contaminant simultaneously) is nearly ten times bigger than that of the B-case when no heat issued from S1. This is naturally reasonable because heat always assists the emission and spreading contaminant. Next, in the A-case on the axis of the symmetry of the box the temperature decreases from point P1 to point P3 while in the B-case we have T2 > T1 > T3. This result is correct too because in the first case the hot air rises directly from S1 domain while in the second the air layer adjacent to S1 is always relatively cool due to the effect of the Rayleigh-Taylor stability. The bound- ary condition (21) and (22) also strongly effect on the structure of the air flow in the box. This can be seen in Fig. 3 where shown the velocity field on the middle horizontal section of the cube. (a) A-problem (b) B-problem Fig. 3. The velocity field on section z = 0.5 at Gr = 2.104, Grc = 105 The flow type changes when Gr increases to 2.105. Fig. 4 shows that at Gr = 2.105 and Grc = 105 the flow in the A-case is clearly unsteady while a time-periodic motion is resulted in the B-case. The existence of the time-periodic flow in the B-case is more clearly indicated in Fig. 5 where the complex structure of this motion is shown also. As seen in Fig. 5 at every point of P1, P2 and P3 each of the flow variables such as u, v, w, T and C has its characteristic phase that is different from that of the other parameters. This phase also changes from point to point. When the Grashof number increases the unsteady air flow appears. In Fig. 6 shown a flow of that kind for the B-case at Gr = 3.105 and Grc = 105 . The simulations indicate that when the Grashof number increases the frequency of the oscillation increases too. The complex structure of the natural convection in the box caused by a complicated boundary condition for heat on the cube base can be demonstrated by Fig. 7. As shown in [6] for the natural convection in a box heated from below there are four different stable structures. In our case the convective structure is more complex because of the non- homogeneity of heat boundary condition on the box bottom. 148 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy -80 -60 -40 -20 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 u w v -0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.2 0.4 0.6 0.8 1.0 T 10*C (a) A-problem -40 -30 -20 -10 0 10 20 30 40 0 2 4 6 8 10 w v u 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0 2 4 6 8 10 T 10*C (b) B-problem Fig. 4. The evolution of velocity components, temperature and contaminant concentration in P1 for Gr = 2.105, Grc = 105 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 0 2 4 6 8 10 u1 u2 u3 -60 -40 -20 0 20 40 60 0.0 2.0 4.0 6.0 8.0 10.0 v1 v2 v3 -60 -40 -20 0 20 40 60 80 0 2 4 6 8 10w1 w2 w3 Fig. 5. The evolution of velocity components in P1, P2, P3 for B-problem. Gr = 2.105, Grc = 105 -80 -60 -40 -20 0 20 40 60 80 0 1 2 3 4 5 6 7 8 9 10 v w 0.00 0.05 0.10 0.15 0.20 0.25 0 2 4 6 8 10 T 50C Fig. 6. The velocity, temperature and contaminant evolution for the B-case at Gr = 3.105, Grc = 105 Contaminant spreading by natural convection in a box 149 (a) (b) (c) Fig. 7. The vector field on: (a) x-sections, x = 0.125, x = 0.5, x = 0.875; (b) y-section, y = 0.125, y = 0.5, y = 0.875; (c) z-section, z = 0.125, z = 0.5, z = 0.875 for the B-case, Gr = 3.105, Grc = 105 Our next concern is the distribution of the contaminant over the box base in the presence of the non-homogeneous heating from below expressed by conditions (19), (21) and (22). The simulation here indicates that if the contaminant issue from the central source S1 is continuous then at any moment the amount of the contaminant over a high- temperature region of the bottom is lager than that over a relatively lower-temperature one. In Fig. 7 presented the distribution of temperature and contaminant on three sections x = 0.875, z = 0.125 and z = 0.875. It is necessary to remind that the hot region is 0 ≤ x ≤ 1, 0.625 ≤ y ≤ 1 while the cooler domain is 0 ≤ x ≤ 1, 0 ≤ y ≤ 0.625. The temperature distribution on any horizontal section, as seen in Fig. 8, is very reasonable namely temperature of the air over the hot region is higher than that over the cooler one. 150 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy And the temperature gradient in y-direction decreases with z-increasing owing to the convective motion in the box. Fig. 8 also shows that the contaminant distribution has the same characteristics as the temperature. Fig. 8. A-case: isolines of temperature (three first pictures) and contaminant on sections: x = 0.875, z = 0.125 and z = 0.875 respectively for Gr = 2.105, Grc = 105 (a) (b) (c) Fig. 9. B-case: isolines of temperature (two first pictures of each line) and contaminant on sections: z = 0.125 and z = 0.875 respectively for (a) steady, (b) periodic, (c) unsteady flow Contaminant spreading by natural convection in a box 151 For the B-problem these characteristics of the temperature and contaminant remain the same. And this is correct for all three types of the flow: steady, periodic and unsteady. This conclusion is demonstrated in Fig. 9. Fig. 10. The experimental illustration of contaminant spreading in a natural convection in a box Fig. 10 shows pictures that illustrates the above conclusion on the way of contam- inant spreading in a box in the presence of a natural convective flow. In this experiment the box bottom is divided into two equal parts. The left part surface is kept cool by ice below while the right part one is heated by a lamp. In the middle picture the contami- nant source locates at the center of the box base (close to our model). In the left and right picture the source is lain in the cool or hot region respectively. As seen in the pictures the smoke always tends to spread more into the hot region of the air. 5. CONCLUSION At moderate Grashof number the three-dimensional natural convection with spread- ing contaminant in an enclosure can be modeled by Navier-Stokes equation in the Boussi- nesq approximation. Using the finite difference method based on the Samarski scheme to calculate the solution of the transport-type equations of vorticity, energy and contami- nant and the multigrid method to compute the velocity shows to be efficient for numeri- cal simulation. The simulation indicates that the mentioned above motion can be steady, time-periodic or unsteady. The transition from one type of the motion to the other de- pends on the Grashof number value as well as the boundary condition on the bottom of the box. The method applied in this study may be applicable to problems with very difficult setting the boundary condition. REFERENCES [1] G. D. Mallinson and G. D. V. Davis. Three-dimensional natural convection in a box: a nu- merical study. Journal of Fluid Mechanics, 83, (01), (1977), pp. 1–31. [2] N. C. Markatos and K. A. Pericleous. Laminar and turbulent natural convection in an en- closed cavity. International Journal of Heat and Mass Transfer, 27, (5), (1984), pp. 755–772. 152 Tran Van Tran, Nguyen Ngoc Thang, Nguyen Thi Thuy [3] W. J. Hiller, S. Koch, and T. A. Kowalewski. Three-dimensional structures in laminar natural convection in a cubic enclosure. Experimental Thermal and Fluid Science, 2, (1), (1989), pp. 34– 44. [4] T. Fusegi, J. Hyun, K. Kuwahara, and B. Farouk. A numerical study of three-dimensional natural convection in a differentially heated cubical enclosure. International Journal of Heat and Mass Transfer, 34, (6), (1991), pp. 1543–1557. [5] R. J. A. Janssen, R. A. W. M. Henkes, and C. J. Hoogendoorn. Transition to time-periodicity of a natural-convection flow in a 3D differentially heated cavity. International Journal of Heat and Mass Transfer, 36, (11), (1993), pp. 2927–2940. [6] J. Pallares, I. Cuesta, F. X. Grau, and F. Giralt. Natural convection in a cubical cavity heated from below at low Rayleigh numbers. International journal of heat and mass transfer, 39, (15), (1996), pp. 3233–3247. [7] O. Aydin and W.-J. Yang. Natural convection in enclosures with localized heating from below and symmetrical cooling from sides. International Journal of Numerical Methods for Heat & Fluid Flow, 10, (5), (2000), pp. 518–529. [8] H. F. Oztop and E. Abu-Nada. Numerical study of natural convection in partially heated rectangular enclosures filled with nanofluids. International Journal of Heat and Fluid Flow, 29, (5), (2008), pp. 1326–1336. [9] L. D. Landau and E. M. Lifshitz. Mechanics of continous mediums. Nauka, Moscow, (1986). (in Russian). [10] V. I. Polezaev, A. V. Bune, N. A. Verezub, G. S. Glusco, V. I. Griaznov, K. G. Dubovic, S. A. Nikitin, A. I. Prostomolotov, A. I. Fedoseev, and S. G. Cherkasov. Mathematical modeling of convective heat and mass transfer on the base of Navier-Stokes equation. Nauka, Moscow, (1987). (in Russian). [11] A. A. Samarski. Introduction to the theory of finite difference method. Nauka, Moscow, (1971). (in Russian). [12] W. Hackbusch. Multi-grid methods and applications. Springer-Verlag, Berlin/New York, (1985). [13] V. John and L. Tobiska. Numerical performance of smoothers in coupled multigrid methods for the parallel solution of the incompressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 33, (4), (2000), pp. 453–473. [14] U. K. N. G. Ghia, K. N. Ghia, and C. T. Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48, (3), (1982), pp. 387–411.