Significance of Rarefaction, Streamwise Conduction, and Viscous Dissipation on the ExtendedGraetz–Nusselt Problem: The Case of Finite-Length Microchannels with Prescribed Wall Heat Flux

%e 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. %e analytical solution proposed, valid for low-intermediate Peclet values, takes into account the presence of the thermal development region. %e influence of all transport parameters (Peclet Pe, Knudsen Kn, and Brinkman Br) 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 Br � 0, the best heating protocol is a decreasing wall heat flux function. In the presence of dissipation Br> 0, the best heating protocol is a uniform wall heat flux.

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][23][24][25] and temperature jump may occur at the wall. e 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 (z ≥ 0) as a fully developed isothermal flow, that is, T � T 0 at z � 0. However, this boundary condition at z � 0 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 z � 0 into the hydrodynamic development region z < 0. erefore, a temperature distribution is built up for z < 0 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. e presence of an upstream insulated region in microchannels has been also accounted for also by Knupp et al. [33][34][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. e 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. e solution, taking into account the presence of the thermal development region, is valid for low-intermediate Peclet values Pe 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 L − , wall heat flux Q w � 0), followed by a heat transfer region with length L and prescribed wall heat flux Q w (z), followed by a third thermally insulated region (length L + , wall heat flux Q w � 0). e 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 Pe and for different wall heat flux profiles. e 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 Pe, the Knudsen number Kn, and the Brinkman number Br. e influence of geometrical parameters, that is, the aspect ratio diameter-over-length D/L and the length of the upstream section L − /L, is also addressed. e solution, valid for small-intermediate values of Pe, is presented for cylindrical and flat rectangular microchannels (see Appendix A). A list of symbols is reported in Appendix B.

Statement of the Problem and Analytical Solution
Let us consider a Newtonian fluid with constant thermodynamic properties entering, at z � − L − , into a finite-length microtube of radius R, diameter D, and total axial length (L + + L + L − ), that is, L − ≤ z ≤ L + L + . (πR 2 )(ρ f c p U)T 0 is the inlet convective flux at z � L − , U being the average cross section velocity. e microtube is thermally insulated for L − ≤ z < 0 and for L < z ≤ L + L + . Let us indicate with Q w (z) the prescribed wall heat flux for the heat transfer region 0 ≤ z ≤ L and with Q av � (1/L) L 0 Q w (z)dz ≠ 0 its average value. Let V z (r) be the fully developed first-order velocity profile in slip-flow regime. In terms of the dimensionless radial coordinate, ρ � r/R, V z (r) attains the form where Kn is the Knudsen number.
A first-order model for rarefaction effects is adopted, being valid for 0.001 ≤ Kn ≤ 0.1 [36] with a tangential momentum accommodation coefficient F � 1 [30]. e velocity profile reduces to the classical no-slip parabolic velocity profile for Kn � 0.
By introducing the dimensionless space variables ρ � r/R, ζ � z/L, temperature θ � (T − T 0 )k/(Q av R), and velocity v(ρ) � V z /U, 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 Pe and the Brinkman number Br appear as together with the geometric dimensionless parameters α, β + , By further introducing the dimensionless axial convective flux e outlet boundary condition equation (6) at ζ � 1 + β + 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.

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 θ b (ζ) is the dimensionless bulk temperature (mixing cup temperature) and ϕ(ρ, ζ) is an auxiliary function satisfying the integral constrain: e auxiliary function ϕ(ρ, ζ) accounts for the temperature dependence on the radial coordinate, that is, (zθ/zρ) � (zϕ/zρ), 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 θ b and ϕ still appear. However, an independent equation for the bulk temperature θ b (ζ) 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 (z 2 ϕ/zζ 2 ) ≪ (d 2 θ b /dζ 2 ), one arrives at the following equation for the dimensionless bulk temperature θ b (ζ):

Viscous dissipation
Slip/no−slip laminar flow Equations (16) and (17) can be analytically solved, thus obtaining It can be observed that, by neglecting the viscous dissipation, that is, by setting Br � 0, the bulk temperature profile is independent of the velocity profile, resulting in the same for slip (Kn ≠ 0, χ ≠ 1) and for no-slip flows (Kn � 0, χ � 1).
By substituting the expression for θ b , (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): e function C(ζ) 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 α ⟶ 0, the analytical solution, equations (18)-(23),(25)-(26), strongly simplifies because and the dimensionless temperature profile attains the form representing the asymptotic (infinite length of the heating section) fully developed temperature profile.

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 Pe and α, the transport problem equations (2)- (8) have been solved with finite elements method (FEM, Comsol 3.5 a). e convection-diffusion package in stationary conditions has been used. Lagrangian quadratic elements are chosen. e linear solver adopted is UMFPACK, with relative tolerance 10 − 12 .
e number of finite elements is 3 × 10 5 with a nonuniform mesh. e maximum element size in all the three subdomains (ζ < 0, 0 ≤ ζ ≤ 1, and ζ > 1) is 1 × 10 − 2 . Smaller elements have been located close to the boundaries. Specifically, a maximum element size of 2 × 10 − 4 is chosen at the external boundary (0 ≤ ζ ≤ 1, ρ � 1) and at the internal cross sections (ζ � 0, 1, 0 ≤ ρ ≤ 1), 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 e mesh adopted guarantees an accurate description of temperature gradients for all the wall heat flux functions adopted (continuous or discontinuous at ζ � 0, 1) and for all the values of Pe analyzed.
As test cases, three different expressions for the wall heat flux function q w (ζ) have been considered.
(1) A triangular function: for which the corresponding integral functions g(ζ) and IF w (ζ) can be easily computed (nonreported here for the sake of brevity). (2) An exponential function [19], decreasing along ζ for n > 0 and increasing along ζ for n < 0: e corresponding integral functions g(ζ) and IF w (ζ) attain the form (3) e uniform wall heat flux function q w (ζ) � 1 that can be obtained as a limiting case of the exponential function for n ⟶ 0. 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 F w (ζ) is continuous at ζ � 0 and ζ � 1 so that the simplifying assumption equation (11)fd11 works well in the entire transport domain. e results are shown for a high Pe value Pe � 100 at the limit of validity of the assumption equation (9)fd9. 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 Br � 0.1, when the fluid enters the heated section. Even in the case of low dissipation, that is, Br � 0.01, the temperature profile is significantly affected by the presence of the upstream section. e behaviour in the thermal development region (close to ζ � 0) 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 Pe in terms of the dimensionless bulk temperature θ b (ζ) and wall temperature θ w (ζ) for the constant wall heat flux function.

Temperature Profiles.
e wall temperature θ w (ζ), accounting for the temperature jump induced by the rarefaction effect, is given by International Journal of Chemical Engineering where F T � 1 (thermal accommodation coefficient), c � 1.4, and Pr � 0.7 are assumed to be the typical values for air [26][27][28].
Model predictions for the wall temperature θ w are accurate for small-intermediate values of Pe, Pe ≤ 100. Indeed, the analytical solution for θ w follows quite closely the numerical solution also for Pe � 10 3 , but small errors (related to the simplifying assumption equation (11)fd11 that fails at ζ � 0 and ζ � 1 for the uniform wall heat flux) are not negligible and are amplified when focusing on the spatial behaviour of local Nusselt number close to ζ � 0, 1.
However, the low accuracy of the analytical solution close to ζ � 0, 1 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 〈Nu〉 as a function of Pe for the exponential wall heat flux function (discontinuous at ζ � 0, 1) in the presence of dissipation Br > 0, for constant n � 0, axially increasing (n < 0) and decreasing (n > 0) wall heat flux. All data refers to an aspect ratio α � 0.02 (long-thin channel and long preheating section) and show an excellent agreement with numerical data up to Pe≃20 − 40.

Limits of Validity of the Analytical Solution
e 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 α � 0.02 and α � 0.1.
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 〈Nu ∞ 〉 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 Pe 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.
is observation becomes more evident by plotting the same data as in Figure 7 as a function of Pe/α instead of Pe (see Figure 8 On the other hand, for high Peclet values, streamwise conduction becomes negligible and numerical results for 〈Nu〉, for different aspect ratios, collapse onto a unique invariant curve for (Peα/4) > 10 − 1 when plotted as a function of (Peα/4) (see Figure 8(b)).    From these observations, it follows that the analytical solution proposed is actually reliable for 0 < (Peα/4) ≤ 10 − 1 that implies for 0 < Pe ≤ (4 · 10 − 1 /α) (see Figure 8(b)). On the other hand, the solution 〈Nu ∞ 〉 based on the fully developed temperature profile θ ∞ is valid only in the range (Pe/α) ≥ 10 2 and (Peα/4) ≤ 10 − 1 that implies 10 2 α ≤ Pe ≤ (4 · 10 − 1 /α). Figure 9 shows the range of validity, in terms of Pe and α values, of the analytical solution proposed in (34)fd35 (region with vertical bars) and of that based on the fully developed temperature profile Nu ∞ (triangular region indicated by arrows), the latter reducing to an empty set for α ⟶ 0.063.
Data reported in Figures 8(a), 8(b), and 9 refer to a case in which the rarefaction effect is not present since Kn � 0. e same analysis can be performed by including the effect of rarefaction (Kn ≠ 0). e region of stability of the analytical solution for Kn ≠ 0 almost coincides with that for Kn � 0. e next section analyzes the influence of different parameters and of different wall heat flux functions on the average Nusselt number 〈Nu ∞ 〉, focusing exclusively on the range of validity of the analytical expression equation (34)fd35.

Influence of Wall Heat Flux Function and Transport Parameters Br and Kn
From a preliminary analysis of data reported in Figures 6(a)-6(b), it can be observed that, in presence of dissipation, that is, Br > 0, the uniform wall heat flux function represents the best heating protocol (larger value of 〈Nu ∞ 〉 for low-intermediate values of Pe ) 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 Pe≃5α (see Figure 6(a)) and the minimum is more pronounced for increasing values of n, corresponding to an increasing amount of energy furnished close to the channel outlet ζ � 1. For Br > 0, no minimum is observed for decreasing wall heat flux functions (see Figure 6(b)).
valid for any other wall heat flux function. However, 〈Nu ∞ 〉 may not be the minimum value attained by 〈Nu〉 because it may exhibit a minimum for intermediate value of Pe, 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 n > 0 and an increasing function of ζ for n < 0. For n � 0, it reduces to the constant function q w � 1. For this exponential heat function, 〈Nu〉 can be evaluated analytically for Br � 0, thus obtaining 〈Nu〉 � 2 (n +(Pe/α))A(Kn) +(Pe/α)B(Kn) n + Pe α + log 1 + 1 − e − (n − (Pe/α)) (Pe/α)B(Kn) (n +(Pe/α))A(Kn) ,

Conclusions
e 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. e 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 θ � θ 1 (ρ, ζ) + θ ∞ (ρ), where θ ∞ is the fully developed temperature profile and θ 1 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. erefore, 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 θ b (ζ) and of the auxiliary function ϕ(ρ, ζ), accounting for the dependence of the temperature on the radial coordinate ρ and slowly varying with the axial coordinate ζ. is 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 q w (ζ) that guarantees the required wall temperature profile. is problem is already under investigation. e 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 Pe, and different channel aspect ratios α, it is possible to conclude that the analytical solution is reliable for Pe < (1/α). e influence of various transport parameters, that is, Br and Kn, is analyzed in detail. e solution proposed is valid in the whole range of Brinkman values considered in the microfluidic literature, that is, Br ≤ 0.1, and therefore valid also in the absence of dissipation effects. Moreover the solution is valid for slip (Kn > 0) as well as for nonslip flows (Kn � 0) and therefore also for viscous fluids for which dissipation effects may not be negligible.
Actually, a value of Brinkman less than 0.1 − 0.2 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 〈Nu〉. 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, Br > 0, 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. e influence of the cross section aspect ratio for rectangular microchannels (three-dimensional problem) on slip-flow heat transfer [41,42] is not investigated.