Abstract
The article addresses the extended Graetz–Nusselt problem in finite-length microchannels for prescribed wall heat flux boundary conditions, including the effects of rarefaction, streamwise conduction, and viscous dissipation. The analytical solution proposed, valid for low-intermediate Peclet values, takes into account the presence of the thermal development region. The influence of all transport parameters (Peclet , Knudsen , and Brinkman ) and geometrical parameters (entry length and microchannel aspect ratio) is investigated. Performances of different wall heat flux functions have been analyzed in terms of the averaged Nusselt number. In the absence of viscous dissipation , the best heating protocol is a decreasing wall heat flux function. In the presence of dissipation , the best heating protocol is a uniform wall heat flux.
1. Introduction
A correct estimation of heat and mass transfer coefficients is a powerful tool in the design of heat exchangers [1], mass transfer equipment and reactors [2], and microdevices [3, 4] for chemical [5], biomedical [6, 7], and pharmaceutical applications [8, 9].
Focusing on laminar forced convection of an incompressible fluid in a duct, the estimation of transport coefficients requires the solution of the classical Graetz–Nusselt problem [2, 10]. Originally proposed for a sudden step change of the wall temperature at some positions along the duct and no axial diffusion, the Graetz–Nusselt problem is valid for both heat and mass transfer. It has been solved in transient and steady-state [11], for Dirichlet and Neumann boundary conditions [12], for different wall shape and curvature [9, 13, 14], for non-Newtonian fluids [15], and for counterflow streams [16], in the presence of high viscous dissipation [17], axial diffusion [18, 19], and simultaneous heat and mass transfer [20, 21].
As the size of the channel is reduced, the no-slip boundary condition needs to be modified because velocity slippage [22–25] and temperature jump may occur at the wall. The so-called “extended” Graetz problem in microtubes, accounting for rarefaction and viscous dissipation, has been recently investigated by Cetin et al. [26, 27] and by Jeong and Jeong [28, 29] by means of eigenfunction expansion (including streamwise conduction) and by Tunc and Bayazitoglu [30] using an integral transform technique (neglecting streamwise conduction).
In all the above-mentioned papers, it is assumed that the fluid enters the semi-infinite microchannel as a fully developed isothermal flow, that is, at . However, this boundary condition at may be extremely restrictive, especially in the case of laminar conditions and low Peclet values [31]. If axial conduction is important, then a sizeable amount of heat is conducted upstream the entrance cross section into the hydrodynamic development region . Therefore, a temperature distribution is built up for and this may significantly affect the temperature downstream. For this reason, Barletta and Magyari [32] addressed the problem of the thermal entrance forced convection in a circular duct with a prescribed wall heat flux distribution, including the effect of viscous dissipation but neglecting heat axial conduction. The presence of an upstream insulated region in microchannels has been also accounted for also by Knupp et al. [33–35] by combining the Generalized Integral Transform Technique (GITT) and a single domain reformulation strategy, aimed at providing hybrid numerical/analytical solutions to convection/diffusion problems. No viscous dissipation effects or slip boundary conditions are accounted for while wall conjugation effects are taken into account.
The present paper presents an analytical solution of the extended Graetz problem in finite-length microtubes including the effects of rarefaction, streamwise conduction, and viscous dissipation. The solution, taking into account the presence of the thermal development region, is valid for low-intermediate Peclet values and for the prescribed heat flux boundary conditions (no wall conjugation effects).
In dealing with finite-length channels, in the presence of axial dispersion and wall heat flux, one issue to be addressed is the proper boundary conditions at the inlet and outlet sections. For this reason, in the problem formulation, three distinct regions along the microchannel have been considered (see Figure 1): a thermally insulated region (length , wall heat flux ), followed by a heat transfer region with length and prescribed wall heat flux , followed by a third thermally insulated region (length , wall heat flux ).
The analytical solution is compared with the numerical results obtained by means of the finite elements method (FEM Comsol 3.5) in a wide range of and for different wall heat flux profiles. The range of validity of the analytical solution is investigated in detail.
From the temperature field, the local Nusselt axial profile and the average Nusselt number are obtained as a function of the transport parameters, that is, the Peclet number , the Knudsen number , and the Brinkman number . The influence of geometrical parameters, that is, the aspect ratio diameter-over-length and the length of the upstream section , is also addressed.
The solution, valid for small-intermediate values of , is presented for cylindrical and flat rectangular microchannels (see Appendix A). A list of symbols is reported in Appendix B.
2. Statement of the Problem and Analytical Solution
Let us consider a Newtonian fluid with constant thermodynamic properties entering, at , into a finite-length microtube of radius , diameter , and total axial length , that is, . is the inlet convective flux at , being the average cross section velocity. The microtube is thermally insulated for and for . Let us indicate with the prescribed wall heat flux for the heat transfer region and with its average value.
Let be the fully developed first-order velocity profile in slip-flow regime. In terms of the dimensionless radial coordinate, , attains the formwhere is the Knudsen number.
A first-order model for rarefaction effects is adopted, being valid for [36] with a tangential momentum accommodation coefficient [30]. The velocity profile reduces to the classical no-slip parabolic velocity profile for .
By introducing the dimensionless space variables , , temperature , and velocity , the starting point of the subsequent analysis is the steady state heat-balance equation, accounting for radial conduction, axial conduction, and convection and viscous dissipation:where the Peclet number and the Brinkman number appear astogether with the geometric dimensionless parameters , , and
By further introducing the dimensionless axial convective fluxthe following boundary conditions apply
The outlet boundary condition equation (6) at is an integral version of the Danckwerts outlet boundary condition, usually adopted for finite-length channels, and implies zero conductive axial heat flux at the outlet section. Its integral version [37] (integral over the radial cross-section) allows us to take into account the heat generated by viscous dissipation and is in agreement with the asymptotic condition usually adopted for infinitely long channels.
2.1. Analytical Solution
At low-intermediate values of the Peclet number, the temperature profile exhibits a weak dependence on the radial coordinate so that it can reasonably assumed that the dimensionless temperature can be written as the linear combination of two functions:where is the dimensionless bulk temperature (mixing cup temperature) and is an auxiliary function satisfying the integral constrain:
The auxiliary function accounts for the temperature dependence on the radial coordinate, that is, , and is a slowly varying function of the axial coordinate , so that
By substituting (9) into the balance equation (2) and into the boundary conditions equations (6)–(8) and making use of the simplifying assumption equation (11), one obtains
In (12), both and still appear. However, an independent equation for the bulk temperature can be obtained as follows. By integrating the balance equation (2) over the entire radial cross section, one obtains
By further enforcing solely the second-order simplifying assumption , one arrives at the following equation for the dimensionless bulk temperature :
Equations (16) and (17) can be analytically solved, thus obtaining
It can be observed that, by neglecting the viscous dissipation, that is, by setting , the bulk temperature profile is independent of the velocity profile, resulting in the same for slip and for no-slip flows .
By substituting the expression for , (18) into the transport equation (12), one arrives at the following equation for the auxiliary function :
By solving (24) along the radial coordinate, one arrives at the following expression for satisfying the boundary conditions equation (14):
The function can be determined by enforcing the integral constraint equation (10):
It can be observed that, in the case of very long-thin channels, that is, for , the analytical solution, equations (18)–(23),(25)–(26), strongly simplifies becauseand the dimensionless temperature profile attains the formrepresenting the asymptotic (infinite length of the heating section) fully developed temperature profile.
For (i.e., constant wall heat flux) and for (i.e., no upstream section) one recovers the fully developed temperature profile for a constant wall heat flux [26, 27], usually written in terms of the dimensionless axial coordinate :
3. Comparison with Numerical Results
In order to verify the correctness of the analytical solution equations (18)–(23),(25)–(26) and to identify the limits of validity in terms of and , the transport problem equations (2)–(8) have been solved with finite elements method (FEM, Comsol 3.5 a).
The convection-diffusion package in stationary conditions has been used. Lagrangian quadratic elements are chosen. The linear solver adopted is UMFPACK, with relative tolerance .
The number of finite elements is with a nonuniform mesh. The maximum element size in all the three subdomains (, , and ) is . Smaller elements have been located close to the boundaries. Specifically, a maximum element size of is chosen at the external boundary and at the internal cross sections (, ), that is, at the internal boundaries between the insulated and the heated regions.
Figure 2 shows the computational domain and the mesh element size adopted for the entire set of simulations with . The mesh adopted guarantees an accurate description of temperature gradients for all the wall heat flux functions adopted (continuous or discontinuous at ) and for all the values of analyzed.
As test cases, three different expressions for the wall heat flux function have been considered.(1)A triangular function: for which the corresponding integral functions and can be easily computed (nonreported here for the sake of brevity).(2)An exponential function [19], decreasing along for and increasing along for : The corresponding integral functions and attain the form(3)The uniform wall heat flux function that can be obtained as a limiting case of the exponential function for .
3.1. Temperature Profiles
Figure 3 shows the excellent agreement between numerical and analytical results for the temperature profiles in both the upstream and downstream sections for the triangular wall heat flux function. In this case, the function is continuous at and so that the simplifying assumption equation (11)fd11 works well in the entire transport domain. The results are shown for a high value at the limit of validity of the assumption equation (9)fd9.
(a)
(b)
It can be observed that there is a significant effect of the upstream thermally insulated section because the bulk temperature is order of unity, for , when the fluid enters the heated section. Even in the case of low dissipation, that is, , the temperature profile is significantly affected by the presence of the upstream section. The behaviour in the thermal development region (close to ) is well described by the analytical solution.
Figures 4(a)and 4(b) show the comparison between analytical and numerical results for increasing values of in terms of the dimensionless bulk temperature and wall temperature for the constant wall heat flux function.
(a)
(b)
The wall temperature , accounting for the temperature jump induced by the rarefaction effect, is given bywhere (thermal accommodation coefficient), , and are assumed to be the typical values for air [26–28].
It can be observed that model predictions for are actually very accurate in the whole range of Peclet values because (1), at low values, the simplifying assumption holds true and (2), at high values, the entire axial conduction term becomes negligible with respect to the axial convective contribution [38–40].
Model predictions for the wall temperature are accurate for small-intermediate values of , . Indeed, the analytical solution for follows quite closely the numerical solution also for , but small errors (related to the simplifying assumption equation (11)fd11 that fails at and for the uniform wall heat flux) are not negligible and are amplified when focusing on the spatial behaviour of local Nusselt number close to .
3.2. Local Nusselt Number
The local Nusselt number in the heating section () attains the formand is a function solely of the ratio and not of and separately. Moreover, the local Nusselt number is independent of the lengths and of the upstream and downstream thermally insulated regions.
The two limiting cases for (no axial convection) and for (Nusselt based on the fully developed temperature profile ) can be easily recovered from (34)fd35 by considering thatthus obtaining
Figures 5(a)and 5(b) show the comparison between numerical results (curves with dots) and analytical results for the local Nusselt number for low-intermediate values of for the triangular wall heat flux (continuous at ) and for the uniform wall heat flux (discontinuous at ). The analytical solution is very accurate for the triangular function in the whole range while deviations from the numerical solution can be observed for and for the uniform wall heat flux due to the simplifying assumption (11)fd11 that fails at discontinuity points .
(a)
(b)
However, the low accuracy of the analytical solution close to for a discontinuous wall heat flux has a very small impact of the average Nusselt number:
Figures 6(a)and 6(b) show the excellent agreement between numerical and analytical results for as a function of for the exponential wall heat flux function (discontinuous at ) in the presence of dissipation , for constant , axially increasing () and decreasing () wall heat flux. All data refers to an aspect ratio (long-thin channel and long preheating section) and show an excellent agreement with numerical data up to .
(a)
(b)
4. Limits of Validity of the Analytical Solution
The analytical expression proposed is reliable also for larger values of the aspect ratio (finite-length channel). Figure 7 shows the comparison between numerical and analytical results for and .
It can be observed that different curves, corresponding to different values of the aspect ratio, saturate towards the same limiting value, corresponding to the average Nusselt number evaluated on the basis of the fully developed temperature profile . However, the larger is the value of , the smaller is the range of validity, in terms of values, of the fully developed profile, and the influence of the thermal developing region must be necessarily accounted for (like in the analytical solution proposed) in order to have an accurate estimate of the average Nusselt number.
This observation becomes more evident by plotting the same data as in Figure 7 as a function of instead of (see Figure 8(a)).
(a)
(b)
Indeed, numerical (and analytical) results for collapse onto an invariant curve when plotted as function of for low-intermediate Peclet values ; see Figure 8 (A), and the asymptotic behaviour sets for .
On the other hand, for high Peclet values, streamwise conduction becomes negligible and numerical results for , for different aspect ratios, collapse onto a unique invariant curve for when plotted as a function of (see Figure 8(b)).
From these observations, it follows that the analytical solution proposed is actually reliable for that implies for (see Figure 8(b)). On the other hand, the solution based on the fully developed temperature profile is valid only in the range and that implies .
Figure 9 shows the range of validity, in terms of and values, of the analytical solution proposed in (34)fd35 (region with vertical bars) and of that based on the fully developed temperature profile (triangular region indicated by arrows), the latter reducing to an empty set for .
Data reported in Figures 8(a), 8(b), and 9 refer to a case in which the rarefaction effect is not present since . The same analysis can be performed by including the effect of rarefaction (). The region of stability of the analytical solution for almost coincides with that for .
The next section analyzes the influence of different parameters and of different wall heat flux functions on the average Nusselt number , focusing exclusively on the range of validity of the analytical expression equation (34)fd35.
5. Influence of Wall Heat Flux Function and Transport Parameters and
From a preliminary analysis of data reported in Figures 6(a)-6(b), it can be observed that, in presence of dissipation, that is, , the uniform wall heat flux function represents the best heating protocol (larger value of for low-intermediate values of ) with respect to both monotonically increasing or monotonically decreasing (exponential) wall heat flux functions.
For the axially increasing wall heat flux, the average Nusselt number exhibits a minimum for (see Figure 6(a)) and the minimum is more pronounced for increasing values of , corresponding to an increasing amount of energy furnished close to the channel outlet . For , no minimum is observed for decreasing wall heat flux functions (see Figure 6(b)).
For a constant wall heat function , is a monotonically decreasing function of , attaining values in the range .
Figures 10(a)and 10(b) show the upper bound and the lower bound as a function of for . The effect of dissipation is to lower by decreasing both the upper and lower bound. The effect of is to monotonically lower both the upper and the lower bounds in the absence of dissipation while, for , the upper and lower bounds exhibit a maximum as a function of .
(a)
(b)
In the absence of dissipation, that is, , both and become constant and independent of the wall heat flux function (see (37) and (38) for ). For this reason, the values of and reported in Figures 10 A-B for represent the asymptotic limits valid for any other wall heat flux function. However, may not be the minimum value attained by because it may exhibit a minimum for intermediate value of , depending on the wall heat flux function.
In order to investigate the role of the wall heat flux function, the following analysis focuses on the exponential function equation (31)fd31 that represents a decreasing function of for and an increasing function of for . For , it reduces to the constant function . For this exponential heat function, can be evaluated analytically for , thus obtaining
Figures 11(a)and 11(b) show the behaviour of as a function of for and for . The arrow indicates increasing values of , in the range . The larger is , the larger is the value of for low-intermediate values of . For , exhibits a neat minimum.
(a)
(b)
Therefore, in the absence of dissipation, a decreasing wall heat flux function () has to be preferred to a constant heat flux function and the larger part of energy must be provided at the entrance of the heating section. On the contrary, in the presence of viscous dissipation, a constant wall heat flux represents the best heating protocol, as shown in Figures 6(a)-6(b) and Figures 12(a) and 12(b) for . High values of amplify differences between increasing and decreasing wall heat flux functions.
(a)
(b)
6. Conclusions
The paper proposes an analytical solution for the extended Graetz–Nusselt problem in finite-length microchannels, including the effects of rarefaction, streamwise conduction, and viscous dissipation. The solution takes into account the presence of a thermally insulated upstream section. Different wall heat functions are analyzed.
A classical approach to the problem (see Jeong and Jeong [28, 29] and Tunc and Bayazitoglu [30]) suggested solving the nonhomogeneous energy balance equation for the dimensionless temperature by setting , where is the fully developed temperature profile and satisfies an homogeneous partial differential equation that can be solved by the method of separation of variables or by integral transform.
In the present case, the analysis focuses on finite-length channels. Therefore, a different “decoupling” of the problem is proposed leading to a more convenient way to express the temperature profile as the linear combination of the bulk temperature and of the auxiliary function , accounting for the dependence of the temperature on the radial coordinate and slowly varying with the axial coordinate . This represents the peculiarity and the strength of our approach to the problem, valid for low-intermediate Peclet values. A similar approach can be adopted to deal with the same transport problem with prescribed wall temperature condition (instead of a prescribed wall heat flux). A Dirichlet boundary condition would require a different nondimensionalization of the transport equation and the solution of an inverse problem; that is, given the analytic expression of the wall temperature for prescribed wall heat flux, find the wall heat flux function that guarantees the required wall temperature profile. This problem is already under investigation.
The range of validity of the analytical solution is investigated by comparing analytical and numerical results obtained with finite elements method. By analyzing different wall heating functions, different values of , and different channel aspect ratios , it is possible to conclude that the analytical solution is reliable for .
The influence of various transport parameters, that is, and , is analyzed in detail.
The solution proposed is valid in the whole range of Brinkman values considered in the microfluidic literature, that is, , and therefore valid also in the absence of dissipation effects. Moreover the solution is valid for slip () as well as for nonslip flows () and therefore also for viscous fluids for which dissipation effects may not be negligible.
Actually, a value of Brinkman less than is rather realistic for microfluidic applications and compatible with Peclet numbers less than 100.
Performances of different wall heat flux functions have been analyzed in terms of the averaged Nusselt number . In the absence of viscous dissipation, the best heating protocol is a decreasing wall heat flux function, where the larger part of energy is furnished at the entrance of the downstream-heated section.
In the presence of dissipation, that is, , the best heating protocol is a uniform wall heat flux and decreasing wall heat flux functions have to be preferred to increasing wall heat flux functions.
In Appendix A the analytical solution proposed for circular microchannels is extended to the case of flat 2-d rectangular microchannels. The influence of the cross section aspect ratio for rectangular microchannels (three-dimensional problem) on slip-flow heat transfer [41, 42] is not investigated.
Appendix
A. Flat Microchannels
The analytical solution proposed for circular microchannels (microtubes) can be extended to flat rectangular microchannels. Let be channel height (the distance between the two flat plates) and the channel width with , such that the hydraulic diameter .
The laminar velocity profile, depending exclusively on the cross section coordinate and satisfying the slip boundary condition at attains the form (for )
The transport equation and boundary conditions for the dimensionless temperature read aswhereand all the dimensionless parameters are the same as defined for microtubes by simply replacing the diameter with the channel height
The analytical solution can be sought as the linear combination of the bulk temperature and of an auxiliary function satisfying the conditions
By assumingthe transport equation (A.2) and the boundary conditions equations (A.3) and (A.4) can be rewritten in terms of and as follows:
By integrating (A.2) over the entire cross section and by enforcing the boundary conditions equation (A.4) at , one obtains the following equation for the bulk temperature :which can be solved analytically, thus obtaining
By substituting the expression for the bulk temperature into the balance equation (A.9), one obtains the following PDE for :which can be solved analytically for , thus obtaining
The analysis of the exact limits of validity of the analytical solution for a flat microchannel will be developed elsewhere, by following the same approach presented for microtubes.
B. List of Symbols
: Brinkman number, : fluid specific heat, : channel diameter, : tangential momentum and thermal accommodation coefficients, : integral function, (19)–(21), : heat transfer coefficient, : integral wall heat flux function, (22)-(23), : thermal conductivity, : Knudsen number, : length of the channel heated section, : upstream and downstream insulated channel lengths, : Nusselt number, : average Nusselt number, : Nusselt number for , : Nusselt number for , : cross-sectional Peclet number, : dimensionless wall heat flux function, : average wall heat flux, : radial and axial coordinates, : channel radius, : temperature, : inlet temperature at , : average axial velocity, : dimensional and dimensionless axial velocity profile, equation (1).
Greek symbols
: channel aspect ratio, , : dimensionless insulated upstream and downstream section lengths, : dimensionless temperature, : dimensionless bulk temperature, : molecular mean free path, : dimensionless axial coordinate entering equation (29), : auxiliary function, equation (9), : dimensionless radial and axial coordinates, : fluid density.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.