International Communications in Heat and Mass Transfer 33 (2006) 908 – 916 www.elsevier.com/locate/ichmt
Numerical simulation of conjugate, turbulent mixed convection heat transfer in a vertical channel with discrete heat sources☆ Roy N. Mathews, C. Balaji ⁎ Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600 036, India Available online 27 March 2006
Abstract This paper presents the results of a comprehensive numerical study to analyze turbulent mixed convection in a vertical channel with a flush-mounted discrete heat source in each channel wall. The conjugate heat transfer problem is solved to study the effect of various parameters like the thermal conductivity of the wall material (ks), the thermal conductivity of the flush-mounted discrete heat source (kc), Reynolds number (Res), modified Richardson number (Ri⁎) and the aspect ratio of the channel (AR). The standard k–ε turbulence model, modified by including buoyancy effects, without wall functions, has been used for the analysis. The twodimensional governing equations are discretised on a semi-staggered, non-uniform grid, using the finite volume method. The asymptotic computational fluid dynamics (ACFD) technique has been then applied to obtain a correlation for the non-dimensional ¯ max, which can be used for a wide range of parameters. maximum temperature θ © 2006 Elsevier Ltd. All rights reserved. Keywords: Turbulent flow; Vertical channel; k–ε model; Mixed convection; Conjugate; Numerical; ACFD
1. Introduction Electronic equipment have become an integral part of almost every phase of modern living. The reliability of these systems is of great importance to the electronic industry. The principal objective of the thermal control system of any electronic component is to maintain a relatively constant temperature, which should be below the maximum service temperature specified by the manufacturer. As the heat fluxes from the components increase, the fluid velocities are to be increased accordingly to maintain the temperature within the prescribed limits. When the heat fluxes become very high, consequent upon the miniaturization of electronic components, values of the Rayleigh number and the Reynolds number may cross the laminar limits. It is in this context that studies on turbulent mixed convection heat transfer are becoming increasingly relevant. The aim of the present study is to analyze turbulent mixed convection
☆
Communicated by A.R. Balakrishnan and S. Jayanti. ⁎ Corresponding author. E-mail address:
[email protected] (C. Balaji).
0735-1933/$ - see front matter © 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2006.02.013
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
909
Nomenclature AR Gr⁎ H k kc kf kn ks ¯ P p ¯ q˙ ‴ Res Ri Ri⁎ S T T∞ ΔTq ¯ U u ¯ ¯ V v ¯ V∞ X Y
Aspect ratio of the channel (H / S) Modified Grashof number, based on volumetric heat generation (βgS3ΔTq / υ2) Height of the channel (m) Turbulent kinetic energy (m2/s2) Thermal conductivity of the heat source (W/m-K) Thermal conductivity of the fluid (W/m-K) Non-dimensional turbulent kinetic energy (k/(V∞)2) Thermal conductivity of the channel wall (W/m-K) 2 Non-dimensional pressure (p¯ / ρV∞ ) 2 Time averaged pressure (N/m ) Rate of volumetric heat generation from the heat source (W/m3) Reynolds number based on S (V∞S / υ) Richardson number (Gr / Re2) Modified Richardson number (Gr⁎s / Res2) Spacing between the two channel walls (m), width of the square enclosure (m) Temperature at any point (K) Uniform inlet temperature (K) Reference temperature difference (q˙ ‴Hδ / kf) Non-dimensional horizontal velocity (u¯ / V∞) Time averaged x-component of velocity (m/s) Non-dimensional vertical velocity (v¯ / V∞) Time averaged y-component of velocity (m/s) Uniform inlet velocity (m/s) Non-dimensional horizontal coordinate (x / S) Non-dimensional vertical coordinate (y / S)
Greek β δ ε εn θ¯ υ νt τ
Coefficient of volume expansion (K− 1) Thickness of channel wall (m) Dissipation rate of k (m2/s3) Non-dimensional ε (εS / (V∞)3) Non-dimensional temperature (T − T∞ / ΔTq) Kinematic viscosity (m2/s) Turbulent eddy viscosity (m2/s) Non-dimensional time (tv∞ / S)
s
heat transfer in a vertical channel for a wide range of parameters and to obtain a correlation for θ¯ max, using the ACFD approach. Gururaja Rao et al. [1] carried out a study on conjugate mixed convection heat transfer for a similar geometry limited to the laminar regime. Premachandran and Balaji [2] studied conjugate mixed convection heat transfer from a horizontal channel with four protruding heat sources and obtained a correlation for θ¯ max, using the ACFD approach, for the laminar regime. They concluded that the thermal conductivity of the protruding heat source affects the maximum temperature to the extent of 40%. Kuyper et al. [3] and Velusamy et al. [4] studied buoyant turbulent flows in square and rectangular enclosures. Many investigators have presented experimental results for mixed convection turbulent flows in various geometries (see, for example, Carr et al. [5], Huang et al. [6]and Zang and Dutta [7]). Nakajima et al. [8] presented experimental and theoretical results for combined free and forced convection between vertical parallel plates. Agrawal et al. [9] and Cotton and Jackson [10], numerically, studied turbulent mixed convection flows in a driven cavity and vertical tube respectively. Balaji and Herwig [11] and Premachandran and Balaji [2] successfully
910
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
employed the ACFD approach to multimode heat transfer problems and concluded that it is an attractive option in the analysis of such problems in heat transfer. The studies of Gururaja Rao et al. [1] and Premachandran and Balaji [2] clearly demonstrated the effect of conjugate heat transfer from a vertical channel and a horizontal channel, respectively, in the laminar regime. The present study intends to explore the effect of conjugate mixed convection heat transfer from a vertical channel with flush-mounted discrete heat sources, in the turbulent regime, and to obtain a correlation for θ¯ max, for a wide range of parameters. The parameters considered are Res, Ri*, Pr, ratio of the thermal conductivities (kc/kf and ks/kf) and the aspect ratio (AR). The number of parameters is large and a full-scale parametric study to develop a correlation for the maximum temperature of the heat sources will be extremely time-consuming, for the turbulent flow regime. Hence, in this study we use the technique of asymptotic expansions to develop a correlation for θ¯ max, with a limited parametric study. 2. Mathematical formulation The geometry for the problem under consideration, along with the system of coordinates, is shown in Fig. 1. There are two flush mounted discrete heat sources, one on each wall. They are of the same size and are centrally located. The outer boundaries of the channel walls are assumed to be perfectly insulated. The flow is assumed to be turbulent and incompressible with constant properties, outside of density, for which the Boussinesq approximation is assumed to hold good. The turbulent flow in the channel is described by the time averaged Reynolds equations. The standard k–ε model is adopted for the prediction of turbulent fluxes due to its simplicity and successful application in turbulent buoyant flows [3,4].
δ
S
δ
g
H1
Substrate having thermal conductivity k s W/m-K
H
Y, V X, U
Flush mounted, uniform heat generating electronic chip, having thermal conductivity k c W/m-K
Fluid entering with uniform velocity V∞ (m/s), temperature T∞ (K) and thermal conductivity kf (W/m-K) Fig. 1. Schematic of the vertical channel with a flush-mounted discrete heat source in each wall.
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
911
2.1. Governing equations The Reynolds averaged Navier–Stokes (N–S) equations along with the energy equation in non-dimensional form are given by BU¯ BV¯ þ ¼0 BX BY
ð1Þ
BU¯ BU¯ BU¯ 1 BP¯ 1 B BU¯ 1 B BU¯ BV¯ þ U¯ þ V¯ ¼− þ 2ð1 þ ttn Þ ð1 þ ttn Þ þ þ Bs BX BY q BX Res BX BX Res BY BY BX
ð2Þ
BV¯ BV¯ BV¯ 1 BP¯ 1 B BU¯ BV¯ 1 B BV¯ þ U¯ þ V¯ ¼− þ ð1 þ ttn Þ þ 2ð1 þ ttn Þ þ þ Ri* h¯ Bs BX BY q BY Res BX BY BX Res BY BY
ð3Þ
Bh¯ Bh¯ Bh¯ 1 B þ U¯ þ V¯ ¼ Bs BX BY Res BX
1 ttn þ Pr rT
Bh¯ BX
1 B þ Res BY
1 ttn þ Pr rT
Bh¯ BY
ð4Þ
The turbulent viscosity υt is computed from a velocity scale (k1 / 2) and a length scale (k3 / 2 / ε) which are predicted at each point in the flow, via a solution of the following transport equations for the turbulent kinetic energy (k) and its dissipation rate (ε), in non-dimensional form, Bkn Bkn Bkn 1 B þ U¯ þ V¯ ¼ Res BX Bs BX BY
1þ
ttn rk
Bkn 1 B ttn Bkn 1þ þ þ Pkn þ Gkn −en Res BY BX rk BY
ttn Ben 1 B ttn Ben 1þ 1þ þ Res BY re BX re BY en þ ½Ce1 ðPkn þ Ce3 Gkn Þ−Ce2 en kn
ð5Þ
Ben Ben Ben 1 B þ U¯ þ V¯ ¼ Res BX Bs BX BY
ð6Þ
Where " 2 2 # 2 ttn BU¯ BV¯ BU¯ BV¯ ttn Ri* Bh¯ Bh¯ k2 þ þ Pkn ¼ 2 þ2 þ ; Gkn ¼ − ; ttn ¼ Res Cl n ; BX BY BY BX Res rT Res BX BY en ¯ Ce3 ¼ tanh V U¯ The X-derivative of temperature appearing in Gkn term has been neglected by many researchers. However, this term has been taken into account in the present study. The constants appearing in the turbulence model are assigned the following values [12] Cμ = 0.09, Cε1 = 1.44, Cε2 = 1.92, σT = 1.0, σk = 1.0, σε = 1.314. 2.2. Boundary conditions 2.2.1. Inlet boundary conditions The fluid is assumed to be entering the channel with uniform velocity and temperature profiles and the appropriate boundary conditions are U¯ ¼ 0; V¯ ¼ 1; h¯ ¼ 0 The k and ε at the inlet are calculated from a prescribed turbulence intensity of 5%.
912
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
2.2.2. Outlet boundary conditions At the exit of the computational domain, the following boundary conditions are used for any variable ϕ B2 / ¼0 BY 2 The above is the familiar smoothing condition and with the use of this, one can dispense with the extension of the computational domain beyond the channel. 2.2.3. Wall boundary conditions Adiabatic boundary conditions are imposed at the outer surfaces of the walls. By doing an energy balance, assuming the fin approximation to be valid (neglecting temperature gradient across the wall), the interface temperature of the wall and fluid can be computed from the equations given below. For the heat generating part of the wall B2 h¯ kf S Bh¯ kf S 2 þ ¼0 þ BY 2 ks d BX ks Hd
ð7Þ
For that part of the wall that does not generate heat B2 h¯ kf S Bh¯ ¼0 þ BY 2 ks d BX
ð8Þ
The boundary conditions for all other variables are U¯ ¼ 0; V¯ ¼ 0; kn ¼ 0 εn = ∞ (For numerical implementation it has been set as 1015) 3. Solution procedure The two dimensional governing equations are discretised on a semi-staggered, non-uniform, grid using a finite volume method. The velocities and pressures are iteratively calculated based on the SIMPLE algorithm. The discretised equations are then solved using an explicit pseudo-time marching technique. Wall functions are not used in the study and sufficiently fine grids have been employed in the inner layer to enable integration up to the walls [3]. The grids in the X direction are generated using a geometric progression. The smallest non-dimensional grid size used close to the wall is 5 × 10− 6. Separate cosine grids are employed along the Y direction for the bottom and top non-heat generating portion and the middle heat generating portion. Comparatively more number of grids is used along the heat generating part. Based on a detailed grid sensitivity test, a grid size of 91 × 181 is used with 41 grids along the heat source. 3.1. Validation The turbulent code has been validated for (i) Turbulent natural convection in a square enclosure in the laminar and turbulent regimes, by comparing the values of average Nusselt number with the published results of Velusamy et al. [4], Markatos and Pericleous [12] and De Vahl Davis [13]. The results were found to be in good agreement. ii) Turbulent mixed convection between vertical parallel plates by comparing the vertical velocity and temperature distribution at the channel exit with the results of Nakajima et al. [8] for Re = 5870 and Ri = 2.2 × 10− 3, and again the results were found to be in good agreement.
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
913
4. Results and discussion Numerical results have been obtained for a wide range of Res, Ri*, kc/kf, ks/kf and AR. For the present study, air is considered as the coolant and the heat generating electronic component is modeled as a discrete heat source flush mounted to the PCB. Two such PCBs vertically mounted constitute the geometry under consideration, as shown in Fig. 1. The geometric parameters are H = 0.3 m, δ = 1.5 × 10− 3 m, H1 = 0.03 m and S that varies according to AR. Calculations have been performed for the following range of parameters: 3000 ≤ Res ≤ 20,000, 0.00375 ≤ Ri* ≤ 0.1125, 41 ≤ kc/kf ≤ 2066, 4 ≤ ks/kf ≤ 413 and 8 ≤ AR ≤ 15. For the entire study, only one parameter is allowed to change at a time and all other parameters are kept at the reference values. The values of the parameters chosen for the base line case are Res = 12,000, Ri* = 0.03825, AR = 15, kc/kf = 325 and ks/kf = 41. 4.1. Effect of the thermal conductivity ratio kc/kf Fig. 2 shows the temperature distribution along the wall for three values of the thermal conductivity ratio of the wall material. An increase of thermal conductivity of the material from 1 to 25 W/m-K decreases the maximum temperature from 329 to 325 K. This means a decrease in the temperature rise (Tmax − T∞) by about 14%. When the thermal conductivity increases further, the decrease in ¯ max is marginal. θ 4.2. Effect of the thermal conductivity ratio ks/kf ¯ max is much more significant than the conductivity ratio It is observed that the effect of thermal conductivity ratio ks/kf on θ kc/kf. By increasing the conductivity ks from 1 to 10 W/m-K, the maximum temperature rise reduces by as much as 20%. Since there are practical limitations for raising the conductivity of the substrate to very high values, in the present study the maximum limit for ks is set as 10 W/m-K. As ks increases, the temperature gradient at the interface shows a decreasing trend, as expected. 4.3. Effect of the aspect ratio, AR The aspect ratio is varied by changing the width of the channel. It was observed that as the AR increases the θ¯ max decreases almost linearly. For higher aspect ratios, however, the dimensional temperature increases due to a decrease in spacing and poorer heat transfer.
335
k c=1 (W/m-K)
330
k c=10 (W/m-K) k c=25 (W/m-K)
Temperature, T (K)
325 320 315 310 305 300 295 0.00
0.05
0.10
0.15
0.20
0.25
0.30
Vertical distance, y (m) Fig. 2. Distribution of temperature along the channel wall for different thermal conductivities.
0.35
914
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
380
Temperature, Τ (Κ)
360
Ri* =0.00375 Ri* =0.03825 Ri* =0.1125
340
320
300
280 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
Vertical distance, y (m) Fig. 3. Distribution of temperature along the channel wall for different Ri*.
4.4. Effect of Reynolds number, Res As the Reynolds number increases, the θ¯ max decreases non-linearly, when all the other parameters are constant at their reference values. 4.5. Effect of modified Richardson number, Ri⁎ The value of Ri* is varied by varying q˙ ‴, while Res is kept constant. The θ¯ max increases non-linearly with Ri*. From Fig. 3, one can clearly see that as the Ri* increases the portion of channel wall above the heat source is being heated up but there is not much change in the temperature of the bottom portion. The non-linear variation of the maximum temperature with Ri* is also seen in the figure. 4.6. Correlation for the maximum temperature The principal objective of the cooling system for any electronic equipment is to maintain the temperature below the maximum allowable limit set for the safe operation of the equipment. Hence, a correlation for predicting the maximum temperature for a wide range of parameters that influence temperature is highly useful for the designers of such systems. Asymptotic expansion is a powerful tool for obtaining a good correlation with a limited parametric study. Practically, a majority of the present day cooling systems operate in the turbulent regime. So, a correlation obtained with a limited parametric study will reduce the computational effort considerably for such problems. The independent variables, that influence the maximum temperature, identified for the present study are AR, Res, Ri*, kc/kf and ks/kf. The reference values selected for these variables are AR = 15, Res = 12,000, Ri⁎ = 0.03825, kc/kf = 325 and ks/kf = 41. The θ¯ max,ref obtained with these values assigned to independent variables is 0.005515 and this corresponds to a Tmax of 325 K. A preliminary ¯ max can be expanded linearly with respect to a new set of variables. A correlation can now be expressed in the study suggested that θ ¯ max,ref as, vicinity of θ h¯ max ¼ h¯ max;ref þ
Bh¯ max Bh¯ max Bh¯ max Bh¯ max Bh¯ max ðe1 −e1;ref Þ þ ðe2 −e2;ref Þ þ ðe3 −e3;ref Þ þ ðe4 −e4;ref Þ þ ðe5 −e5;ref Þ Be1 Be2 Be3 Be4 Be5
ð9Þ
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
915
where 0:6 * Þ0:32 ; e e1 ¼ ðRe=Reref Þ0:35 ; e1;ref ¼ 1:0; e2 ¼ ðRi* =Riref 2;ref ¼ 1:0; e3 ¼ ðAR=ARref Þ ; e3;ref ¼ 1:0; e4 ¼ ðkc =kf Þ=ðkc =kf Þref Þ0:1 ; e4;ref ¼ 1:0; e5 ¼ ðks =kf Þ=ðks =kf Þref Þ0:33 ; e5;ref ¼ 1:0
The coefficients of the above equation are evaluated one at a time by obtaining solutions to the original governing equation with ¯ max with different ¯ max/∂ε1 is evaluated by solutions obtained for θ all but one parameter fixed at the reference values. For example, ∂θ values of ε1, while ε2 to ε5 be kept at their reference values. The evaluated coefficients are now substituted and the correlation is derived as, h¯ max ¼ 0:005515 þ 9:065 10−3 ð1−e1 Þ−4:195 10−5 ð1−e2 Þ þ 1:283 10−2 ð1−e3 Þ þ 2:888 10−3 ð1−e4 Þ þ 1:015 10−3 ð1−e5 Þ
ð10Þ
This equation is valid for the range of parameters given earlier, when air is used as the cooling medium. A parity plot showed that ¯ max obtained from the numerical solution and θ ¯ max obtained from the above correlation compare very the maximum temperature, θ well (error band ±5% ).
5. Conclusions In the present work, conjugate turbulent mixed convection heat transfer in a vertical channel has been analyzed numerically. The procedure for combining the method of asymptotic expansions with CFD, called ACFD, has been applied to the above-mentioned problem. The main conclusions of the study are: 1. A reasonably fine grid coupled with the boundary conditions k → 0, ε → ∞ and the use of a semi-staggered grid gives convergent solutions and through this approach the use of wall functions can be avoided. 2. The maximum temperature of electronic components can be reduced significantly by increasing the effective thermal conductivity of the PCBs. 3. For the range of parameters tested, the maximum temperature can be reduced by increasing the effective thermal conductivity of the components up to around 25 W/m-K. Increasing the conductivity beyond this does not significantly affect the maximum temperature. 4. The correlation developed using the ACFD approach with a limited parametric study gives fairly accurate predictions for the maximum temperature in the channel over a wide range of parameters that are of engineering interest. References [1] C. Gururaja Rao, C. Balaji, S.P. Venkateshan, Effect of surface radiation on conjugate mixed convection in a vertical channel with a discrete heat source in each wall, International Journal of Heat and Mass Transfer 45 (2002) 3331–3347. [2] B. Premachandran, C. Balaji, Mixed convection Heat transfer from a horizontal channel with protruding heat sources, Heat and Mass Transfer 41 (2005) 510–518. [3] R.A. Kuyper, Th.H. van der Meer, C.J Hoogendoorn, R.A.W.M. Henkes, Numerical study of laminar and turbulent natural convection in an inclined square cavity, International Journal of Heat and Mass Transfer 36 (1993) 2899–2911. [4] K. Velusamy, T. Sundararajan, K.N. Seetharamu, Interaction effects between surface radiation and turbulent natural convection in square and rectangular enclosures, ASME Journal of Heat Transfer 123 (2001) 1062–1070. [5] A.D. Carr, M.A. Connor, H.O. Buhr, Velocity, temperature and turbulence measurements in air for pipe flow with combined free and forced convection, ASME Journal of Heat Transfer 95 (1973) 445–452. [6] T.M. Huang, C. Gau, W. Aung, Mixed convection flow and heat transfer in a heated vertical convergent channel, International Journal of Heat and Mass Transfer 38 (1995) 2445–2456. [7] X. Zhang, S. Dutta, Heat transfer analysis of buoyancy-assisted mixed convection with asymmetric heating conditions, International Journal of Heat and Mass Transfer 41 (1998) 3255–3264. [8] M. Nakajima, K. Fukui, H. Ueda, T. Mizushina, Buoyancy effects on turbulent transport in combined free and forced convection between parallel plates, International Journal of Heat and Mass Transfer 23 (1980) 1325–1336. [9] L. Agrawal, J.C. Mandal, A.G. Maathe, Computations of laminar and turbulent mixed convection in a driven cavity using pseudocompressibility approach, Computers and Fluids 30 (2001) 607–620. [10] M.A. Cotton, J.D. Jackson, Vertical tube air flow in the turbulent mixed convection regime calculated using a low-Reynolds-number k–ε model, International Journal of Heat and Mass Transfer 33 (1990) 275–286.
916
R.N. Mathews, C. Balaji / International Communications in Heat and Mass Transfer 33 (2006) 908–916
[11] C. Balaji, H. Herwig, The use of ACFD approach in problems involving surface radiation and free convection, International Communications in Heat Mass Transfer 30 (2003) 251–259. [12] N.C. Markatos, K.A. Pericleous, Laminar and turbulent natural convection in an enclosed cavity, International Journal of Heat and Mass Transfer 27 (1984) 755–772. [13] G. De Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, International Journal for Numerical Methods in Fluids 3 (1983) 249–264.