Applied Thermal Engineering 31 (2011) 921e931
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
A numerical study on heat transfer performance of microchannels with different surface microstructures Y. Liu a, *, J. Cui a, Y.X. Jiang b, W.Z. Li a a b
Key Laboratory of Ocean Energy Utilization and Energy Conservation of Ministry of Education, Dalian University of Technology, Dalian 116024, Liaoning Province, China Haier Refrigerator R&D Department, Qingdao Haier Joint Stock Co., Ltd., 266103, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 July 2010 Accepted 10 November 2010 Available online 18 November 2010
Forced convection heat transfer occurring in microchannel is numerically studied in this paper using the CFD (computational fluid dynamics) and LB (lattice Botlzmann) approaches. Simulation results of these two methods are compared and tested against available experimental correlation, and a good agreement is achieved. It suggests that both methods are suitable to describe the liquid flow in microchannels. The influence of microchannel geometric shape on heat transfer performance is investigated by evaluating fluid thermophysical parameter and Nusselt number of the high-temperature surface. It is found that the inflow liquid temperature raises intensively along the flow direction at the entry region, and it develops gradually into equilibrium stage approaching the exit for each microstructure studied in this study. The heat exchange efficiency increases with inlet Reynolds number. The results also imply that the shieldshaped groove microchannel possesses the highest heat exchange performance. Compared with the lowest heat transfer efficiency of the plain surface structure, the averaged Nusselt number can be increased by about 1.3 times. Through the field synergy principle analysis, we find that it is the synergy between temperature gradient and velocity that results in different heat transfer performance for different microstructures. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Microchannel Heat transfer Numerical simulation CFD LBM
1. Introduction The application of microchannel heat sinks is drawing more and more attention as one of the most promising high efficiency heat exchange technologies in, e.g., the cooling of electronic devices, automotive heat exchangers, laser process equipment and aerospace technology, etc. The MEMS (Micro-Electro-Mechanical Systems) microchannel heat sinks, which can be directly packaged with the electronic chip, have also been found to be very promising for situations that require high heat flux in relatively small spaces, such as electronic cooling applications. Recognizing the significant potential of microchannel heat sinks in removing heat, there has been a great deal of interest in the study of fluid flow and heat transfer in microchannels. The pioneering numerical work can be found in Ref. [1], where the concept of microchannel heat sink was first proposed. Many research works have been done afterwards [2e8], as reviewed by [9]. There are presently three main approaches to simulate the fluid flow and heat
* Corresponding author. Tel.: þ86 13234060684. E-mail address:
[email protected] (Y. Liu). 1359-4311/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2010.11.015
transfer in microchannels. The first one is the LBM (Lattice Boltzmann method), the second method is the MD (Molecular dynamics) approach and the third one is the conventional CFD method. Lattice Boltzmann (LB) model [10] is a class of CFD methods for flow simulation. It is well known that the NaviereStokes equation cannot adequately describe gas flows in the transition and freemolecular regimes. In these regimes, the Boltzmann equation (BE) of kinetic theory is invoked to govern the flows. Instead of solving the NaviereStokes equation, the discrete BE is solved to simulate the Newtonian fluid flow with collision models such as BhatnagareGrosseKrook (BGK) [11]. By simulating streaming and collision processes across a limited number of particles, the intrinsic particle interactions evince a microcosm of viscous flow behavior applicable across the greater mass. Regarding to heat transfer process, the TLBM (Thermal Lattice Boltzmann methods) proposed by [12], is developing very fast in the recent five years. This approach is quite appropriate to describe microscale flow and the involved heat transfer problems, because the arithmetic itself is developed from the microscale flow mechanism. This method has been extensively studied in many microscopic flow systems and it has been introduced into the calculations of microchannel flow as well. The main
922
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
Fig. 1. Schematic of the microstructure and coordinate system (structure A: Ridge-shaped groove; structure B: V-shaped groove; structure C: shield-shaped groove; structure D: straight slot groove; structure E: plain surface).
focus was cast onto the rarefied and compressible gas flow [13e16], because this method can well describe gaseous flows in microchannels, which is beyond the describing capability of traditional CFD approaches. Some MD attempts were also made, such as the works of [17e19]. For liquid flow in microchannels the NaviereStokes and energy equations and the no-slip boundary condition based on the continuum assumption are still valid and can be used to predict the flow and heat transfer characteristics with reasonable accuracy [20]. The heat transfer process in two microchannels under laminar flow condition was modeled in [20]. The results showed that different micro-geometries can have great effect on heat transfer performance for the microchannel system. Toh et al. [21] performed a numerical study on flow behavior and heat transfer characteristics in a silicon microchannel heat sink where the microstructure was a wide rib. The CFD method was employed to account for the problem. It was found that the thermal resistances were well predicted compared with experimental data, and specific heat and thermal conductivity of the fluid have insignificant effect on flow parameters. Ji et al. [22] investigated a pressure-driven gas flow system in microchannel between two long parallel flat plates via CFD approach. Wall roughness effect on the rarefied and compressible gas flows was evaluated. It was found that rarefied flow parameters are sensitive to the shape of roughness elements. An increase of the relative roughness height leads to a pronounced decrease in the local heat flux for both rarefied and compressible flows. For rarefied flows, the roughness elements do not improve the average Nu number, i.e., the efficiency of heat transfer. Chen et al. [23] investigated the flow parameters of three different microstructures by CFD method and it mainly focused on the pressure drop and overall heat power, the detailed heat exchange parameters were not well presented or analyzed. However, this paper as well as the work of [20] pointed out that the microstructures of the microchannel can have great effect on the flow characteristic, and alerting the microstructures can significantly change the flow traits and heat exchange performance. Thus, to make an intensive study regarding the microchannel structure’s effect on heat transfer performance is indeed very important. This point is further extended in this paper to numerically explore which kind of structure would possess the best heat exchange efficiency.
2. Model development 2.1. Problem statement We consider cold water with a constant density of 1000 kg/m3, flowing over five different microstructure planes with a constant high temperature. The flow behavior and heat transfer performance are investigated using numerical approaches. Fig. 1 illustrates the corresponding cross-sectional geometries and dimensions of the microstructure, above which the cool water flow through a 60 mm length. As can be seen, the microstructure appears periodically, and in most previous simulation works, only the flow at one characteristic structure was modeled. To avoid the geometry unsymmetry effect on the flow field and to eliminate potential numerical errors, 10 columns (flow passage width ¼ 10 s) of the periodic microstructure flow passages are modeled and simulated in the present research. In the present investigation, the CFD method and LBM are employed to account for the above physical problem. The mathematical formulations are going to be described in the next sections. Primarily, the aim of this comparative study is to test the reliability of the two approaches in the microscale liquid flow. Then, the heat exchange performance of different structures is tested and a relatively better structure is expected to be found.
2.2. CFD method Flow in microchannel is a microscopic problem, and temperature jump and velocity slip processes may take place herein, which will invalidate the traditional CFD method. However, as pointed out by [20,22], the rarefaction, temperature jump and velocity slip’s effect on flow characteristics and heat transfer behavior only occur for gaseous flow. These effects can be omitted for liquid flow in microchannels with minimum dimension of about one tenth millimeter. Thus the CFD method is still appropriate and justified to describe the liquid microscopic flow considered in this study and it is adopted to investigate the heat transfer performance of different microstructures herein. Present model is derived from continuum-based steady-stage conservation equations of mass, momentum and energy, with the following several basic assumptions:
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
(1) Continuous Newtonian fluid, Reynolds with steady laminar flow and heat transfer. (2) Specific heat, thermal conductivity and viscosity are single variable functions of temperature. (3) Negligible gravity, radiation heat transfer and natural convection heat transfer.
923
Table 1 Constants in properties vs. temperature correlations. x
bx
cx
dx
0 1 2 3
8958.9 40,535 0.11243 1.014 104
0.58166 6.3556 103 7.964 106
2.414 105 247.8 413.0
The resulting continuity, momentum and energy governing equations are:
! V$ u ¼ 0
(1)
Cp ¼ b0 þ b1 ðT 273Þ þ b2 ðT 273Þ2 þb3 ðT 273Þ3
(5)
! rð! u $V u Þ ¼ Vp þ V$s
(2)
k ¼ c0 þ c1 ðT 273Þ þ c2 ðT 273Þ2
(6)
rCp ! u $VT ¼ V$ðkVTÞ
(3)
m ¼ d0 10d1 =ðTd2 Þ
(7)
! where r is water density, p is pressure, u is the velocity, k is the water thermal conductivity and s is the stressestrain tensor expressed in equation (4).
2 !T ! s ¼ m V! V$ u I u þVu 3
(4)
The above equations are subject to boundary conditions and the detailed boundary conditions are listed as follows: (1) Velocity-inlet boundary: Seven inlet flow velocities uv,in are selected in the present study and the corresponding Re is 50, 100, 150, 200, 300, 400 and 500, respectively.
uu ¼ 0; uv ¼ uv;in ; uw ¼ 0; T ¼ Tinlet (2) Outflow boundary: The diffusion fluxes for all flow variables are considered to be zero which means that the values at the outflow plane are extrapolated from within the domain and have no impact on the upstream flow.
vuu vuv vuw vT ¼ 0 ¼ 0; ¼ 0; ¼ 0; vy vy vy vy (3) Wall with constant temperature boundary:
uu ¼ 0; uv ¼ 0; uw ¼ 0; T ¼ Twall (4) Symmetry wall boundary: Zero normal velocity and zero normal gradient of temperature at the symmetry plane.
u! ¼ 0; n
vT ! ¼ 0 vn
In the present study, the properties of water are allowed to vary with temperature, and the functional relations between specific heat, thermal conductivity and viscosity with temperature are given as:
In equations (5)e(7), thermal conductivity k is in W/m K, viscosity m is in N s/m2 and temperature T is in K, respectively. The various constants are listed in Table 1. The Prandtl number, defined as Pr ¼ y/ k, is also changeable with temperature, where y and k are kinematical viscosity and thermal diffusivity, respectively. The mass, momentum and temperature governing equations are discretized using the second order upwind scheme. The SIMPLE algorithm is applied for the pressureevelocity coupling. The governing equations are solved iteratively until they are convergent and the convergence criterion is 103 for all variables. To evaluate the grid size effect on the accuracy of numerical solutions, gridindependence tests are performed for the five microstructures as the grid size is refined until acceptable difference between the last two grid sizes is found. It is found that 2,254,200, 711,600, 1,001,700, 738,600 and 242,400 computational cells, for structures AeE respectively, are sufficient to produce results independent of the grid size. They are thus employed in this study as displayed in Fig. 2. Based on the above meshes, the governing equations are solved on the Fluent 6.3.26 platform [24], and an 88-node Intel Xeon E5420 Linux cluster is employed to accomplish the solution iteration processes. Take the structure A for example. It takes about 10e13 h to get convergence using the Fluent 8-node parallel calculation mode. 2.3. LB model This physical problem is also accounted for by using a DDF (double distribution function) LBM [25,26] in this section. The D3Q15 model is used to describe the flow filed, while the D3Q6 model is adopted to simulate the temperature field. The EDF (equilibrium distribution function) of the temperature field is calculated from the macroscopic velocity, which is directly decided from velocity field. Consequently, the flow field determines temperature field and it is not affected by the temperature field. The configurations of lattice velocity for the two models ! [27,28] are given in Fig. 3, where the velocity vectors e i for the D3Q15 lattice are (0, 0, 0), (1, 1, 1), (1, 0, 0), (0, 1, 0), (0, 0, * 1), and the velocity vectors vi for the D3Q6 lattice are (1, 0, 0), (0, 1, 0), (0, 0, 1). The LBM dynamics for the velocity field are given by:
i 1h ! ! ! ! ! fi ð x ; tÞ fieq ð x ; tÞ fi ð x þ e i dt ; t þ dt Þ ¼ fi ð x ; tÞ
sn
(8)
! where x is the position on the lattice, t is the time, fi represents the PDF (particle density function) in the ith velocity direction of the lattice, and feq i is the EDF described as:
924
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
Fig. 2. Computational mesh.
" fieq ¼ wi r 1 þ
# ! ! ! !2 ! ! e i$ u e i$ u Þ u$u þ ð 2c4s c2s 2c2s
(9)
r¼
i¼0
fi ;
! u ¼
14 X
! fi e i =r
(10)
i¼0
The evolution equation for temperature is given as follow:
i 1h ! ! ! ! gi ð x ; tÞ gieq ð x ; tÞ gi ð x þ vi dt ; t þ dt Þ ¼ gi ð x ; tÞ
sT
(11)
The LBGK model for PDFs in the temperature field is similar to that in velocity field but a simplified EDF can be used: eq
gi
¼
! ! T v $u 1 þ i2 6 cs
(12)
Temperature can be calculated by:
a
6 X
gi
(13)
i¼1
with the weights wp 0¼ ffiffiffi 2/9, wi ¼ 1/9 for i ¼ 1e6, wi ¼ 1/71 for i ¼ 7e14 and cs ¼ c= 3 is defined as the lattice speed of sound, here c ¼ dx/dt is the lattice speed. The density and velocity can be obtained through moment summations in the velocity space as: 14 X
T ¼
The ChapmaneEnskog expansion for the temperature distribution function recovers the macroscopic energy equation. This gives the kinematic viscosity y and thermal diffusivity k in terms of the single relaxation:
y ¼ c2s dt ðsy 0:5Þ;
k ¼ c2s dt ðsT 0:5Þ
(14)
Here, the above LBM governing equations are employed to investigate the flow behavior and heat transfer characteristics. The uniform 900 300 60 lattice is found to be sufficient to produce results independent of lattice size in the grid-independence tests, and thus used in our LBM simulations. The solid obstacles with five microstructures are placed on the bottom walls to form the corresponding groove microchannels, as shown in Fig. 4, where the red and blue regions represent the solid obstacles and the liquid fluid, respectively (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.). As to the boundary conditions, the velocity-inlet and the constant-low-temperature conditions are employed for the inlet, while the outflow condition is imposed at the exit. The top surface is treated as symmetry wall. For the above three surfaces, the
b
Fig. 3. Lattice geometry and velocity vectors of (a) D3Q15 model and (b) D3Q6 model.
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
925
Fig. 4. Solid obstacles in the lattice system for the five microstructure grooves.
macroscopic variables are determined firstly, as described in Section 2.2. Then, the nonequilibrium extrapolation rule proposed by Guo et al. [29], which is of second order and has good numerical stability, is used to obtain the distribution functions. The standard bounce-back rule is employed for the bottom surface wall and the solid obstacles’ borders, and their temperature keeps at a constant high value. The periodic boundary condition is adopted for both the left and right hand surfaces. The working fluid is water, and its properties are set the same as in the CFD approach, see in Section 2.2 for detail. Accordingly, the dimensionless relaxation times, for both flow field sy and temperature field sT, are also variable. They can be determined by equation (14). The accuracy and stability of numerical simulation can be guaranteed by adjusting the time step dt and lattice speed c, and the details can be found in [30], where it is suggested that the relaxation times should fall into the range between 0.5 and 2. The developed serial LBM code for the problem considered in this paper is optimized specifically for parallel calculation environment considering the computational cost [31]. A parallel computational code [32] is eventually developed and compiled to solve the above equations using 32 nodes for each run. The 88-node Intel Xeon E5420 Linux cluster is also employed to accomplish the LBM solution iteration process. It takes about 32e38 h to get convergence for each case. 3. Results and discussion 3.1. Velocity distribution The simulation result of flow in the grooves reveals the presence of secondary flows consisting of vortices. The flow fields in the grooved structures are shown in Fig. 5 and the results are from the CFD calculations. The upward motion above the peak and the downward flow towards the valley can be found in structures A and B, which is evident and in agreement with the experimental findings of [33] and the DNS study of [34].
The dimensionless temperature distributions at four xez crosssections along the channel are shown in Fig. 6 for the Re ¼ 500 case, where the four cross-sections (y0 ¼ 0.167, y0 ¼ 0.33, y0 ¼ 0.5, y0 ¼ 1) are all placed normal to the flow direction. It can be seen that the isotherms are closest at the entrance region and then spread out gradually along the channel for each microstructure. It means the temperature gradient, e.g. heat transfer, is highest at the entrance and it decreases along the flow direction. It can be also noticed that the temperature field tends to be more “flat” when approaching the exit. This implies the heated upstream water is flushed downstream and the convection effect overtakes that of temperature diffusion. It can be also seen that the temperature field almost reaches developed flow after the y0 ¼ 0.5 cross-section and there is no distinct difference between y0 ¼ 0.5 and y0 ¼ 1 slices. Besides, it can be also found that the temperature of most computational cells in the groove C and groove D nearly reaches the wall temperature, while only some cells of the grooves A and B are filled with the near-wall temperature fluid. This implies that microstructures C and D have better local heat exchange performance, and may have better overall heat exchange performance than A and B, which will be quantitatively introduced later. It should be noted that the temperature contour results shown in Fig. 6, are derived from the LBM solution and those from the CFD model have quite similar trend, thus not qualitatively presented here. The heat transfer results of both models are also going to be presented quantitatively next. Apart from the above local temperature distribution analysis, the outlet temperature is a more direct way for one to evaluate the overall heat exchange intensity. Fig. 7 shows the dimensionless temperature hoist at the exit, as defined by equation (16). As the equation suggests, the higher outlet temperature is, the higher dimensionless temperature hoist is, which symbolizes a better heat exchange performance. It can be seen from Fig. 7 that both the CFD and LBM produce agreed and reasonable results. The structure C possesses the highest dimensionless temperature hoist for all Reynolds numbers, while the flat plane structure E has the lowest heat exchange intensity. The microstructures D, A and B take the second, third and fourth position.
3.2. Temperature distribution The temperature results are presented and analyzed in this section. To facilitate the comparison between the CFD approach and LBM results, the simulation results obtained are nondimensionalized as follows:
T Tinlet T0 ¼ ; Twall Tinlet
y y0 ¼ L
(15)
where T0 and y0 are dimensionless temperature and dimensionless distance, respectively.
0 Toutlet ¼
Toutlet Tinlet Twall Tinlet
(16)
3.3. Heat exchange characteristics In order to evaluate the local heat transfer characteristics at the microstructure surface along the flowing direction, the local Nusselt numbers must be defined. The definition of local Nusselt number for strong conjugate heat transfer problem suggested in [7]
926
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
Fig. 5. Flow field results (y ¼ L/2 slice).
Fig. 6. Dimensionless temperature contour profile for different longitudinal cross-sections (Re ¼ 500).
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
927
Fig. 7. Dimensionless temperature hoist at the exit. Fig. 8. Comparison of predicted Nu with experimental correlation of [35] for microstructure D.
is adopted in this study for the CFD model as shown in equation (17):
qwall $L Nu ¼ ðTwall Tbulk Þk
(17)
where Twall is the wall temperature, Tbulk is the fluid bulk temperature, L is the flowing length, and qwall is the local heat flux of the microstructure channel surface given by:
vTðx; y; zÞ qwall ¼ k ! vn
(18)
For the results obtained by LBM, the Nusselt number can be calculated by the following expression:
Nu ¼
Twall T 0 L dzðTwall Tbulk Þk
(19)
where T0 is the temperature of the lattice located nearest to the wall surface. Single-phase forced convective heat transfer and flow characteristics of water in microchannel with small rectangular channels were studied experimentally by [35], where correlations for Nusselt number under laminar and turbulent flow were proposed. It is expressed in equation (20) for laminar flow condition, where s, s2 and h are the geometrical parameters indicated in Fig. 1.
Nu ¼ 0:1165
2s2 $h sðs2 þ hÞ
0:81 0:79 h Re0:62 Pr 1=3 s2
(20)
Predicted Nusselt number in the thermal developed region under different Reynolds number is compared with this correlation to validate the present numerical models for microstructure D in Fig. 8. It can be seen that a reasonable agreement is achieved between the model predicted value and experimental correlation value for both the two approaches. It is worth noting that the predicted values of both approaches slightly deviate from correlation value for high Re cases. This can be explained as: the hydraulic diameter of microchannel D in the present study is 0.67 mm, which is larger than the original experimental hydraulic diameter range of correlation 20, i.e. 0.133e0.367 mm. However, this is the only available and the closest to our structure dimension correlation in the open literature. Fig. 8 implies that the present models can be used to describe the flow and heat transfer in microchannel and this experimental correlation can actually be extended to larger
hydraulic diameter microchannels under low Re number cases. The comparisons of other microstructures are not preformed because there is presently no experiment investigating the heat transfer efficiency for other microstructures at the best of the authors’ knowledge. However, it is believed that with the justified predicted heat transfer results for microstructure D, the present models can be evaluated to be accurate and reliable and it can be extended to other micro-geometries, and this is just what the CFD tool was intended for. Fig. 9 shows the Nusselt number variation of the five microstructures along the flowing direction, with the Reynolds number ranging from 50 to 500. It can be clearly observed that the local heat transfer intensity is very high at the inlet, and it decreases with the fluid going downstream and gradually approaches constant, just as in a channel with a conventional dimension. It can be also deduced that the Nusselt number of microstructure C is about 1.3 times more than that of microstructure E from an averaged sense, which symbolizes a much better heat transfer performance. From the numerical results it can be found that after the fluid flowing past a certain length, the local Nusselt number becomes constant, indicating the occurrence of fully developed regime. The thermal entry length of a circular tube Sthermal can be estimated by equation (21), as described in [36]. For the microchannel studied in this paper, the thermal entry lengths can be estimated by replacing the tube diameter d with the characteristic length: groove spacing s indicated in Fig. 1.
Sthermal ¼ 0:05RePrd
(21)
For the five microchannel configurations investigated, assuming Pr ¼ 7.0, Reynolds number ranging from 50 to 500, the corresponding dimensionless thermal entry lengths (S0 ¼ Sthermal/L) are 0.06, 0.12, 0.175, 0.233, 0.35, 0.466 and 0.583, respectively. Comparing the above results with Fig. 9(a)e(e), it can be seen that both the CFD and LBM approaches have quite similar results, and they both overestimated the thermal entry length for low Reynolds number case, while a quite reasonable agreement is achieved for cases with Reynolds number greater than 200. Thus, the entrance length of microchannel fluid flow and heat transfer is still a problem open to discussion. It can be also observed from Fig. 9 that the heat transfer performance for high Reynolds number cases is better than that for low Reynolds number cases for each microstructure flow. This is because a higher fluid velocity will
928
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
Fig. 9. Longitudinal Nusselt numbers variation for various Reynolds numbers. a. Structure A; b. structure B; c. structure C; d. structure D; e. structure E.
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931 90 89 88 87
β, degree
86 85
Structure A Structure B Structure C Structure D Structure E
84 83 82 81 80 0
50
100
150
200
250
300
350
400
450
500
550
Re Fig. 10. Synergy angle b varying with Reynolds number.
produce a higher near-wall temperature gradient and eventually result in a better heat exchange performance.
3.4. Further discussion To gain a better understanding of the previous results and particularly concerning why the five microstructures have different heat exchange performance, the field synergy principle is introduced, which is a novel concept of heat transfer enhancement [37e39]. The general form of 3D energy equation can be expressed as:
rcp uu
vT vT vT þ uv þ uw vx vy vz
¼
v vT v vT k þ k vx vx vx vy
v vT _ þ k þF vx vz
(22)
It can be rewritten by integration over the whole computational domain into the following form:
Z
rcp uu U
þ
vT vT vT v vT v vT k þ k þuv þuw vx vy vz vx vx vx vy
v vT _ dU ¼ k vT k F ! vx vz v n wall
(23)
where U is the computational domain and F stands for internal heat source. For liquid flow with Peclet number greater than 100, the heat conduction in the liquids can be neglected compared with ! the convection term, and the term kðvT=v n Þjwall represents fluid thermal conduction effect at the surface wall. Thus, we can finally write:
Z U
vT ! ð u $VTÞdU ¼ k ! v n wall
(24)
! The vector dot product, u $VT in the dimensional integration in equation (24) can be expressed as:
! ! u $VT ¼ j u jjVTjcos b
(25)
where b is the intersection angle, defined as synergy angle, between the velocity and the temperature gradient (heat flow vector). Equation (25) shows that in the convection domain there
929
! ! are two vector fields, u and VT, or three scalar fields, j u j, jVTj and cos b. The heat exchange strength of the convection heat transfer not only depends on the velocity field, the temperature field, but also on their synergy level. In order to have a insight view of the above temperature distribution results and to find what it is that makes the microstructure C have the best heat exchange performance, the volume averaged synergy angle for the five microchannels flow is calculated. Fig. 10 shows the results of synergy angle b varying with Reynolds number. It can be seen that the structure C possesses the highest synergy angle covering all Reynolds numbers, meaning that the synergy degree is best for this microstructure. The plain microstructure E has the lowest synergy angle, meaning that the synergy between velocity and temperature gradient is worst in this configuration and consequently the heat exchange performance is worst. The microstructures D, A and B take the second, third and fourth place regarding to the synergy angle, i.e. heat exchange intensity, which is the same conclusion as drawn by the above temperature analysis. This insight analysis reveals that it is the synergy level that lies behind and decides the heat exchange performance for different microstructure microchannels. To enhance the heat exchange intensity in microchannels, the synergy between velocity field and temperature gradient must be improved. The various surface microstructures can all be treated as methods intended to increase the synergy level. The microstructure C can produce best synergy between velocity and temperature gradient and this configuration is thus recommended for enhanced heat transfer purpose among the five microstructures studied in this research. 4. Conclusions In this paper, CFD and LBM are employed to simulate the flow and heat transfer processes of liquid flow in microchannels, where five microstructures are modeled and analyzed. The two approaches have quite similar and reasonable results compared with experimental data. The distributions of the temperature field and longitudinal Nusselt number in microchannels are presented for heat transfer performance. In addition, the comparison of the heat exchange efficiencies is conducted among the five microchannels using outlet dimensionless temperature hoist and longitudinal Nusselt number. Besides, the field synergy principle is introduced to have an insight view of the heat exchange process, and what leads to different microstructure’s different overall heat exchange performance is found. The conclusions of this paper can be summarized as: (1) Both CFD and LBM can produce reasonable accurate results compared to experimental data. (2) The temperature of liquid fluid is increased along the flow direction because the heating effect of the high-temperature surface. The water temperature field reaches thermal developed stage just past a short flowing distance. (3) The heat exchange is very intensive at the inlet for each microchannel, and it decreases gradually along the flowing direction towards the exit. (4) Microstructure C has the highest overall Nusselt number and heat exchange performance, the microstructure E has the lowest heat exchange efficiency. Quantitatively, the Nusselt number of microstructure C is about 1.3 times more than that of microstructure E from an averaged sense. (5) It is the synergy level between velocity and temperature gradient that decides the heat exchange performance. Various surface microstructures can all be treated as methods intended to increase the synergy level. The microstructure C has the highest synergy level covering all Reynolds numbers and the
930
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931
best heat transfer performance is achieved for this microchannel. Microstructure E has the lowest heat transfer efficiency, which can be also reflected in the field synergy analysis. Nomenclature
Cp c cs ! e f g h k L Nu ! n Pr p q Re S s s1 s2 t T ! u ! v w x y z
specific heat, kJ/kg K lattice speed, e lattice speed of sound, e velocity vectors for the D3Q15 lattice, e distribution function for flow field, e distribution function for temperature field, e geometrical parameter, mm thermal conductivity, W/m K flowing length, mm Nusselt number, dimensionless normal direction vector, e Prandtl number, dimensionless pressure, Pa heat flux, W/m2 Reynolds number, dimensionless thermal entry length, mm geometrical parameter, mm geometrical parameter, mm geometrical parameter, mm discrete time, e temperature, K velocity, m/s velocity vectors for the D3Q6 lattice, e weight factor, e coordinate, e coordinate, e coordinate, e
Greek letters synergy angle, lattice time step, e lattice spacing step, e internal heat source, W/m3 geometrical angle, density, kg/m3 stressestrain tensor, Pa relaxation time for flow field, e relaxation time for temperature field, e
b dt dx F q r s sv st
Subscripts i index of discrete lattice velocity u x direction in Cartesian coordinate v y direction in Cartesian coordinate w z direction in Cartesian coordinate Superscripts dimensionless variable eq equilibrium state 0
References [1] D.B. Tuckerman, R.F.W. Pease, High-performance heat sinking for VLSI, IEEE Electron Device Letters 2 (5) (1981) 126e129. [2] L.A. Jiang, M. Wong, Y. Zohar, Forced convection boiling in a microchannel heat sink, Journal of Microelectromechanical Systems 10 (1) (2001) 80e87. [3] S.G. Kandlikar, W.J. Grande, Evolution of microchannel flow passages e thermohydraulic performance and fabrication technology, Heat Transfer Engineering 24 (1) (2003) 3e17.
[4] A.G. Fedorov, R. Viskanta, Three-dimensional conjugate heat transfer in the microchannel heat sink for electronic packaging, International Journal of Heat and Mass Transfer 43 (3) (2000) 399e415. [5] P.S. Lee, S.V. Garimella, D. Liu, Investigation of heat transfer in rectangular microchannels, International Journal of Heat and Mass Transfer 48 (9) (2005) 1688e1704. [6] S.J. Kim, D. Kim, Forced convection in microstructures for electronic equipment cooling, Journal of Heat Transfer e Transactions of the ASME 121 (3) (1999) 639e645. [7] J. Li, G.P. Peterson, P. Cheng, Three-dimensional analysis of heat transfer in a micro-heat sink with single phase flow, International Journal of Heat and Mass Transfer 47 (19e20) (2004) 4215e4231. [8] T.S. Zhao, Q. Liao, Thermal effects on electro-osmotic pumping of liquids in microchannels, Journal of Micromechanics and Microengineering 12 (2002) 962e970. [9] C.B. Sobhan, S.V. Garimella, A comparative analysis of studies on heat transfer and fluid flow in microchannels, Nanoscale and Microscale Thermophysical Engineering 5 (4) (2001) 293e311. [10] Z. Guo, T.S. Zhao, Lattice Boltzmann model for incompressible flows through porous media, Physical Review E 66 (3) (2002) 36304. [11] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical Review 94 (3) (1954) 511e525. [12] X. He, S. Chen, G.D. Doolen, A novel thermal model for the lattice Boltzmann method in incompressible limit* 1, Journal of Computational Physics 146 (1) (1998) 282e300. [13] X.D. Niu, S.A. Hyodo, T. Munekata, K. Suga, Kinetic lattice Boltzmann method for microscale gas flows: issues on boundary condition, relaxation time, and regularization, Physical Review E 76 (3) (2007) 036711. [14] Z.W. Tian, C. Zou, H.J. Liu, Z.L. Guo, Z.H. Liu, C.G. Zheng, Lattice Boltzmann scheme for simulating thermal micro-flow, Physica A e Statistical Mechanics and Its Applications 385 (1) (2007) 59e68. [15] Y.T. Chew, X.D. Niu, C. Shu, Three-dimensional lattice Boltzmann BGK model and its application to flows with heat transfer in a rectangular microchannel, International Journal for Numerical Methods in Fluids 50 (11) (2006) 1321e1334. [16] N. Verma, D. Mewes, A. Luke, Lattice Boltzmann study of velocity, temperature, and concentration in micro-reactors, International Journal of Heat and Mass Transfer 53 (15e16) (2010) 3175e3185. [17] C.Y. Ji, Y.Y. Yan, A molecular dynamics simulation of liquidevapouresolid system near triple-phase contact line of flow boiling in a microchannel, Applied Thermal Engineering 28 (2e3) (2008) 195e202. [18] S.V. Nedea, A.J. Markvoort, A.A. van Steenhoven, P.A.J. Hilbers, Heat transfer predictions for micro-/nanochannels at the atomistic level using combined molecular dynamics and Monte Carlo techniques, Journal of Heat Transfer 131 (3) (2009) 033104. [19] G. Nagayama, M. Kawagoe, A. Tokunaga, T. Tsuruta, On the evaporation rate of ultra-thin liquid film at the nanostructured surface: a molecular dynamics study, International Journal of Thermal Sciences 49(1) 59e66. [20] Z. Li, W.Q. Tao, Y.L. He, A numerical study of laminar convective heat transfer in microchannel with non-circular cross-section, International Journal of Thermal Sciences 45 (12) (2006) 1140e1148. [21] K.C. Toh, X.Y. Chen, J.C. Chai, Numerical computation of fluid flow and heat transfer in microchannels, International Journal of Heat and Mass Transfer 45 (26) (2002) 5133e5141. [22] Y. Ji, K. Yuan, J.N. Chung, Numerical simulation of wall roughness on gaseous flow and heat transfer in a microchannel, International Journal of Heat and Mass Transfer 49 (7e8) (2006) 1329e1339. [23] Y. Chen, C. Zhang, M. Shi, J. Wu, Three-dimensional numerical simulation of heat and fluid flow in noncircular microchannel heat sinks, International Communications in Heat and Mass Transfer 36 (9) (2009) 917e920. [24] Fluent, FLUENT 6.3 User’s Guide. [25] Q. Li, Y.L. He, Y. Wang, W.Q. Tao, Coupled double-distribution-function lattice Boltzmann method for the compressible NaviereStokes equations, Physical Review E 76 (5) (2007) 56705. [26] G.H. Tang, W.Q. Tao, Y.L. He, Thermal boundary condition for the thermal lattice Boltzmann equation, Physical Review E 72 (1) (2005) 16703. [27] A. Dupuis, J.M. Yeomans, Lattice Boltzmann modelling of droplets on chemically heterogeneous surfaces, Future Generation Computer Systems 20 (6) (2004) 993e1001. [28] C.S. Azwadi, S. Syahrullail, A three-dimension double-population thermal lattice BGK model for simulation of natural convection heat transfer in a cubic cavity, WSEAS Transactions on Mathematics 8 (9) (2009) 561e571. [29] Z. Guo, B. Shi, C. Zheng, A coupled lattice BGK model for the Boussinesq equations, International Journal for Numerical Methods in Fluids 39 (4) (2002) 325e342. [30] M. Wang, N. Pan, Modeling and prediction of the effective thermal conductivity of random open-cell porous foams, International Journal of Heat and Mass Transfer 51 (5e6) (2008) 1325e1331. [31] G. Wellein, T. Zeiser, G. Hager, S. Donath, On the single processor performance of simple lattice Boltzmann kernels, Computers & Fluids 35 (8e9) (2006) 910e919. [32] C. Körner, T. Pohl, U. Rüde, N. Thürey, T. Zeiser, Parallel lattice Boltzmann methods for CFD applications, Numerical Solution of Partial Differential Equations on Parallel Computers (2006) 439e466.
Y. Liu et al. / Applied Thermal Engineering 31 (2011) 921e931 [33] Y. Suzuki, N. Kasagi, Turbulent drag reduction mechanism above a riblet surface, AIAA Journal 32 (1994) 1781e1790. [34] H. Choi, P. Moin, J. Kim, Direct numerical simulation of turbulent flow over riblets, Journal of Fluid Mechanics 255 (1993) 503e539. [35] X.F. Peng, G.P. Peterson, Convective heat transfer and flow friction for water flow in microchannel structures, International Journal of Heat and Mass Transfer 39 (12) (1996) 2599e2608. [36] F.P. Incropera, D.P. DeWitt, T.L. Bergman, A.S. Lavine, Fundamentals of heat and mass transfer.
931
[37] W.Q. Tao, Z.Y. Guo, B.X. Wang, Field synergy principle for enhancing convective heat transfer e its extension and numerical verifications, International Journal of Heat and Mass Transfer 45 (18) (2002) 3849e3856. [38] Z.Y. Guo, W.Q. Tao, R.K. Shah, The field synergy (coordination) principle and its applications in enhancing single phase convective heat transfer, International Journal of Heat and Mass Transfer 48 (9) (2005) 1797e1807. [39] Y.P. Cheng, Z.G. Qu, W.Q. Tao, Y.L. He, Numerical design of efficient slotted fin surface based on the field synergy principle, Numerical Heat Transfer, Part A: Applications 45 (6) (2004) 517e538.