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

Authors: Antonio Brasiello; Alessandra Adrover
Journal: International Journal of Chemical Engineering
 (2020)

---

## 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 (PecletPe, 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.


---

## Body

## 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 microchannelz≥0 as a fully developed isothermal flow, that is, T=T0 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. Therefore, 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. 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 valuesPe 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 Figure1): a thermally insulated region (length L−, wall heat flux Qw=0), followed by a heat transfer region with length L and prescribed wall heat flux Qwz, followed by a third thermally insulated region (length L+, wall heat flux Qw=0).Figure 1
Schematic representation of the physical domain and boundary conditions.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 ofPe 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 numberPe, the Knudsen number Kn, and the Brinkman number Br. The 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.The solution, valid for small-intermediate values ofPe, 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, atz=−L−, into a finite-length microtube of radius R, diameter D, and total axial length L++L+L−, that is, L−≤z≤L+L+. πR2ρfcpUT0 is the inlet convective flux at z=L−, U being the average cross section velocity. The microtube is thermally insulated for L−≤z<0 and for L<z≤L+L+ . Let us indicate with Qwz the prescribed wall heat flux for the heat transfer region 0≤z≤L and with Qav=1/L∫0LQwzdz≠0 its average value.LetVzr be the fully developed first-order velocity profile in slip-flow regime. In terms of the dimensionless radial coordinate, ρ=r/R, Vzr attains the form(1)Vzρ=2U1−ρ2+42−F/FKn1+82−F/FKn=2U1−ρ2+4Kn1+8Kn=U21−ρ2χ+1−χ,χ=11+8Kn,Kn=λD,where Kn is the Knudsen number.A first-order model for rarefaction effects is adopted, being valid for0.001≤Kn≤0.1 [36] with a tangential momentum accommodation coefficient F=1 [30]. The 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−T0k/QavR, and velocity vρ=Vz/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:(2)−Peα4vρ∂θ∂ζ+1ρ∂∂ρρ∂θ∂ρ+α24∂2θ∂ζ2+32Brχ2ρ2=0,−β−<z<1+β+,where the Peclet number Pe and the Brinkman number Br appear as(3)Pe=UDk/ρfcp,Br=μU2QavD,together with the geometric dimensionless parameters α, β+, and β−(4)α=DL,β−=L−L,β+=L+L.By further introducing the dimensionless axial convective flux(5)Jcρ,ζ=Peα4vρθρ,ζ,the following boundary conditions apply(6)Jcζ=−β−=0,∂∂ζ∫01Jc2πρdρ=16πBrχ2,atζ=1+β+,(7)∂θ∂ρρ=0=0,∂θ∂ρρ=1=Fwζ,(8)Fwζ=0,for−β−≤ζ<0,qwζ=Qwζ/Qav,for0≤ζ≤1,0,for1<ζ≤1+β+.The 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.


### 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:(9)θρ,ζ=θbζ+ϕρ,ζ,θbζ=1π∫01θvρ2πρdρ,where θbζ is the dimensionless bulk temperature (mixing cup temperature) and ϕρ,ζ is an auxiliary function satisfying the integral constrain:(10)1π∫01ϕρ,ζvρ2πρdρ=0.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(11)∂ϕ∂ζ≪dθbdζ,∂2ϕ∂ζ2≪d2θbdζ2.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(12)∂ϕ∂ζ≪dθbdζ,∂2ϕ∂ζ2≪d2θbdζ2,(13)θbζ=−β=0,Peα4dθbdζζ=1+β+=16Brχ2,(14)∂ϕ∂ρρ=0=0,∂ϕ∂ρρ=1=Fwζ.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(15)−Peα4dθbdζ+α24d2θbdζ2+∂2ϕ∂ζ2+16Brχ2+2Fwζ=0.By further enforcing solely the second-order simplifying assumption∂2ϕ/∂ζ2≪d2θb/dζ2, one arrives at the following equation for the dimensionless bulk temperature θbζ:(16)−Peα4dθbdζ+α24d2θbdζ2+16Brχ2+2Fwζ=0,(17)θb−β=0,Peα4dθbdζζ=1+β+=16Brχ2.Equations (16) and (17) can be analytically solved, thus obtaining(18)θbζ,(19)gζ<0=PeαePeζ/α∫01e−Peζ′/αqwζ′dζ′,(20)g0≤ζ≤1=PeαePeζ/α∫ζ1e−Peζ′/αqwζ′dζ′,(21)gζ>1=0,(22)IFwζ<0=0,IFwζ>1=1,(23)IFw0≤ζ≤1=∫0ζqwζ′dζ′.It can be observed that, by neglecting the viscous dissipation, that is, by settingBr=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 ϕρ,ζ:(24)1ρ∂∂ρρ∂ϕ∂ρ=2gζ+8Brχ2vρ−1+2Fwζ+8Brχ2−32Brχ2ρ2.By solving (24) along the radial coordinate, one arrives at the following expression for ϕρ,ζ satisfying the boundary conditions equation (14):(25)ϕρ,ζ=Fwζ+8Brχ21+χ+gζχρ22+−8Brχ21+χ+gζχρ44+Cζ.The functionCζ can be determined by enforcing the integral constraint equation (10):(26)Cζ=gζχχ24−16+Fwζχ12−14+Brχ2χ23−χ−43.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(27)limPe/α⟶∞gζ,Peα=Fwζ,and the dimensionless temperature profile attains the form(28)θ∞ρ,ζ=8Pe2Fwζ+PeαIFwζ+8Brχ2ζ+β−︸θbζ++1+χFwζ+8Brχ2ρ22−8Brχ21+χ+Fwζχρ44︸ϕρ,ζ++Fwζ−14−χ12+χ224+Brχ2−43−χ+χ23︸Cζ,representing the asymptotic (infinite length of the heating section) fully developed temperature profile.Forqwζ=1 (i.e., constant wall heat flux) and for β−=0 (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 ξ=z/PeR:(29)θ∞ρ,ζ=41+8Brχ2ξ+1+8Brχ21+χρ22+−8Brχ21+χ+χρ44+8Pe2++−14−χ12+χ224+Brχ2−43−χ+χ23.


## 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:(9)θρ,ζ=θbζ+ϕρ,ζ,θbζ=1π∫01θvρ2πρdρ,where θbζ is the dimensionless bulk temperature (mixing cup temperature) and ϕρ,ζ is an auxiliary function satisfying the integral constrain:(10)1π∫01ϕρ,ζvρ2πρdρ=0.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(11)∂ϕ∂ζ≪dθbdζ,∂2ϕ∂ζ2≪d2θbdζ2.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(12)∂ϕ∂ζ≪dθbdζ,∂2ϕ∂ζ2≪d2θbdζ2,(13)θbζ=−β=0,Peα4dθbdζζ=1+β+=16Brχ2,(14)∂ϕ∂ρρ=0=0,∂ϕ∂ρρ=1=Fwζ.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(15)−Peα4dθbdζ+α24d2θbdζ2+∂2ϕ∂ζ2+16Brχ2+2Fwζ=0.By further enforcing solely the second-order simplifying assumption∂2ϕ/∂ζ2≪d2θb/dζ2, one arrives at the following equation for the dimensionless bulk temperature θbζ:(16)−Peα4dθbdζ+α24d2θbdζ2+16Brχ2+2Fwζ=0,(17)θb−β=0,Peα4dθbdζζ=1+β+=16Brχ2.Equations (16) and (17) can be analytically solved, thus obtaining(18)θbζ,(19)gζ<0=PeαePeζ/α∫01e−Peζ′/αqwζ′dζ′,(20)g0≤ζ≤1=PeαePeζ/α∫ζ1e−Peζ′/αqwζ′dζ′,(21)gζ>1=0,(22)IFwζ<0=0,IFwζ>1=1,(23)IFw0≤ζ≤1=∫0ζqwζ′dζ′.It can be observed that, by neglecting the viscous dissipation, that is, by settingBr=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 ϕρ,ζ:(24)1ρ∂∂ρρ∂ϕ∂ρ=2gζ+8Brχ2vρ−1+2Fwζ+8Brχ2−32Brχ2ρ2.By solving (24) along the radial coordinate, one arrives at the following expression for ϕρ,ζ satisfying the boundary conditions equation (14):(25)ϕρ,ζ=Fwζ+8Brχ21+χ+gζχρ22+−8Brχ21+χ+gζχρ44+Cζ.The functionCζ can be determined by enforcing the integral constraint equation (10):(26)Cζ=gζχχ24−16+Fwζχ12−14+Brχ2χ23−χ−43.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(27)limPe/α⟶∞gζ,Peα=Fwζ,and the dimensionless temperature profile attains the form(28)θ∞ρ,ζ=8Pe2Fwζ+PeαIFwζ+8Brχ2ζ+β−︸θbζ++1+χFwζ+8Brχ2ρ22−8Brχ21+χ+Fwζχρ44︸ϕρ,ζ++Fwζ−14−χ12+χ224+Brχ2−43−χ+χ23︸Cζ,representing the asymptotic (infinite length of the heating section) fully developed temperature profile.Forqwζ=1 (i.e., constant wall heat flux) and for β−=0 (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 ξ=z/PeR:(29)θ∞ρ,ζ=41+8Brχ2ξ+1+8Brχ21+χρ22+−8Brχ21+χ+χρ44+8Pe2++−14−χ12+χ224+Brχ2−43−χ+χ23.


## 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 Pe 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 tolerance10−12.The number of finite elements is3×105 with a nonuniform mesh. The 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.Figure2 shows the computational domain and the mesh element size adopted for the entire set of simulations with β−=β+=1. The 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.Figure 2
Color plot of the mesh element size [dimensionless units] for the computational domain−β−≤ζ≤1+β+,0≤ρ≤1 for β−=β+=1.As test cases, three different expressions for the wall heat flux functionqwζ have been considered.(1)
A triangular function:(30)qwζ=4ζ,for≤ζ<1/2,qwζ=4−4ζ,for1/2≤ζ≤1,for which the corresponding integral functionsgζ and IFwζ 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:(31)qwζ=ne−nζ1−e−n.The corresponding integral functionsgζ and IFwζ attain the form(32)gζ<0=nePeζ/α1−e−nPe/αPe/α+n1−e−Pe/α+n,gζ>1=0,g0≤ζ≤1=ne−nζ1−e−nPe/αPe/α+n1−e−Pe/α+n1−ζ,IFwζ<0=0,IFwζ>1=1,IFw0≤ζ≤1=1−e−nζ1−e−n.(3)
The uniform wall heat flux functionqwζ=1 that can be obtained as a limiting case of the exponential function for n⟶0.


### 3.1. Temperature Profiles

Figure3 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 Fwζ is continuous at ζ=0 and ζ=1 so that the simplifying assumption equation (11)fd11 works well in the entire transport domain. The results are shown for a high Pe value Pe=100 at the limit of validity of the assumption equation (9)fd9.Figure 3
Contour plot of the dimensionless temperatureθρ,ζ for the triangular wall heat flux for Pe=100, α=0.02, and β−=β+=1. Comparison between numerical results (black lines with dots) and analytical results (continuous red lines). (a) Br=0.1, Kn=0.0; (b) Br=0.01, Kn=0. Labels indicate some reference temperature values.
(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, forBr=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. The behaviour in the thermal development region (close to ζ=0) is well described by the analytical solution.Figures4(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.Figure 4
Comparison between numerical (blue dots) and analytical results (red continuous lines) for the constant wall heat fluxqwζ=1, Br=0.1, Kn=0.0, α=0.02, β−=β+=1. (a) θb versus ζ for Pe=100,101,102,103,104,105; (b) θw versus ζ for Pe=10−2,10−1,100,101,102,103. The arrows indicate increasing values of Pe.
(a)(b)The wall temperatureθwζ, accounting for the temperature jump induced by the rarefaction effect, is given by(33)θwζ=θ1,ζ+2bKn∂θ∂ρρ=1=θbζ+ϕ1,ζ+2bKnFwζ,b=2−FTFT2γγ+11Pr,where FT=1 (thermal accommodation coefficient), γ=1.4, and Pr=0.7 are assumed to be the typical values for air [26–28].It can be observed that model predictions forθbζ are actually very accurate in the whole range of Peclet values Pe∈10−2,105 because (1), at low Pe values, the simplifying assumption ∂2ϕ/∂ζ2≪d2θb/dζ2 holds true and (2), at high Pe values, the entire axial conduction term ∂2θ/∂ζ2 becomes negligible with respect to the axial convective contribution Peα/4vρ∂θ/∂ζ [38–40].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=103, 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.


### 3.2. Local Nusselt Number

The local Nusselt numberNuζ in the heating section (0≤ζ≤1) attains the form(34)Nu0≤ζ≤1=2qwθw−θb=2qwζϕ1,ζ+2bKnqwζ=2qwζqwζ1/4+χ/12+2bKn+8Brχ21+χ+gζ,Pe/αχ1/12+χ/24,and is a function solely of the ratio Pe/α and not of Pe and α separately. Moreover, the local Nusselt number is independent of the lengths β− and β+ of the upstream and downstream thermally insulated regions.The two limiting casesNu0 for Pe/α⟶0 (no axial convection) and Nu∞ for Pe/α⟶∞ (Nusselt based on the fully developed temperature profile θ∞) can be easily recovered from (34)fd35 by considering that(35)limPe/α⟶0gζ,Pe/α=0,limPe/α⟶∞gζ,Pe/α=qwζ,thus obtaining(36)Nu0=2qwζqwζ1/4+χ/12+2bKn+8Brχ21+χ1/12+χ/24,(37)Nu∞=2qwζqwζ1/4+χ/6+χ2/24+2bKn+8Brχ21+χ1/12+χ/24.Figures5(a)and 5(b) show the comparison between numerical results (curves with dots) and analytical results for the local Nusselt number Nuζ for low-intermediate values of Pe=10−2,10−1,100,101 for the triangular wall heat flux (continuous at ζ=0,1) and for the uniform wall heat flux (discontinuous at ζ=0,1). The analytical solution is very accurate for the triangular function in the whole ζ range while deviations from the numerical solution can be observed for ζ<0.015 and ζ>0.99 for the uniform wall heat flux due to the simplifying assumption (11)fd11 that fails at discontinuity points ζ=0,1.Figure 5
Comparison between numerical (blue curves with dots) and analytical results (red continuous lines) forNuζ versus ζ. The arrows indicate increasing values of Pe. Br=0.1, Kn=0.0, α=0.02, β+=β−=1, Pe=10−2,10−1,100,101. (a) Triangular wall heat flux. (b) Uniform wall heat flux.
(a)(b)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:(38)Nu=∫01Nuζ′dζ′.Figures6(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.Figure 6
Comparison between numerical (red curves with dots) and analytical results (blue continuous lines) for the average NusseltNu as a function of Pe for the exponential wall heat flux function qwζ=ne−nζ/1−e−n, Br=0.1, Kn=0.0, α=0.02, β+=β−=2. (a) Constant wall heat flux n=0 and exponentially increasing wall heat flux n=−1,−2,−3,−5. (b) Constant wall heat flux n=0 and exponentially decreasing wall heat flux n=1,2,3.
(a)(b)


## 3.1. Temperature Profiles

Figure3 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 Fwζ is continuous at ζ=0 and ζ=1 so that the simplifying assumption equation (11)fd11 works well in the entire transport domain. The results are shown for a high Pe value Pe=100 at the limit of validity of the assumption equation (9)fd9.Figure 3
Contour plot of the dimensionless temperatureθρ,ζ for the triangular wall heat flux for Pe=100, α=0.02, and β−=β+=1. Comparison between numerical results (black lines with dots) and analytical results (continuous red lines). (a) Br=0.1, Kn=0.0; (b) Br=0.01, Kn=0. Labels indicate some reference temperature values.
(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, forBr=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. The behaviour in the thermal development region (close to ζ=0) is well described by the analytical solution.Figures4(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.Figure 4
Comparison between numerical (blue dots) and analytical results (red continuous lines) for the constant wall heat fluxqwζ=1, Br=0.1, Kn=0.0, α=0.02, β−=β+=1. (a) θb versus ζ for Pe=100,101,102,103,104,105; (b) θw versus ζ for Pe=10−2,10−1,100,101,102,103. The arrows indicate increasing values of Pe.
(a)(b)The wall temperatureθwζ, accounting for the temperature jump induced by the rarefaction effect, is given by(33)θwζ=θ1,ζ+2bKn∂θ∂ρρ=1=θbζ+ϕ1,ζ+2bKnFwζ,b=2−FTFT2γγ+11Pr,where FT=1 (thermal accommodation coefficient), γ=1.4, and Pr=0.7 are assumed to be the typical values for air [26–28].It can be observed that model predictions forθbζ are actually very accurate in the whole range of Peclet values Pe∈10−2,105 because (1), at low Pe values, the simplifying assumption ∂2ϕ/∂ζ2≪d2θb/dζ2 holds true and (2), at high Pe values, the entire axial conduction term ∂2θ/∂ζ2 becomes negligible with respect to the axial convective contribution Peα/4vρ∂θ/∂ζ [38–40].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=103, 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.


## 3.2. Local Nusselt Number

The local Nusselt numberNuζ in the heating section (0≤ζ≤1) attains the form(34)Nu0≤ζ≤1=2qwθw−θb=2qwζϕ1,ζ+2bKnqwζ=2qwζqwζ1/4+χ/12+2bKn+8Brχ21+χ+gζ,Pe/αχ1/12+χ/24,and is a function solely of the ratio Pe/α and not of Pe and α separately. Moreover, the local Nusselt number is independent of the lengths β− and β+ of the upstream and downstream thermally insulated regions.The two limiting casesNu0 for Pe/α⟶0 (no axial convection) and Nu∞ for Pe/α⟶∞ (Nusselt based on the fully developed temperature profile θ∞) can be easily recovered from (34)fd35 by considering that(35)limPe/α⟶0gζ,Pe/α=0,limPe/α⟶∞gζ,Pe/α=qwζ,thus obtaining(36)Nu0=2qwζqwζ1/4+χ/12+2bKn+8Brχ21+χ1/12+χ/24,(37)Nu∞=2qwζqwζ1/4+χ/6+χ2/24+2bKn+8Brχ21+χ1/12+χ/24.Figures5(a)and 5(b) show the comparison between numerical results (curves with dots) and analytical results for the local Nusselt number Nuζ for low-intermediate values of Pe=10−2,10−1,100,101 for the triangular wall heat flux (continuous at ζ=0,1) and for the uniform wall heat flux (discontinuous at ζ=0,1). The analytical solution is very accurate for the triangular function in the whole ζ range while deviations from the numerical solution can be observed for ζ<0.015 and ζ>0.99 for the uniform wall heat flux due to the simplifying assumption (11)fd11 that fails at discontinuity points ζ=0,1.Figure 5
Comparison between numerical (blue curves with dots) and analytical results (red continuous lines) forNuζ versus ζ. The arrows indicate increasing values of Pe. Br=0.1, Kn=0.0, α=0.02, β+=β−=1, Pe=10−2,10−1,100,101. (a) Triangular wall heat flux. (b) Uniform wall heat flux.
(a)(b)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:(38)Nu=∫01Nuζ′dζ′.Figures6(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.Figure 6
Comparison between numerical (red curves with dots) and analytical results (blue continuous lines) for the average NusseltNu as a function of Pe for the exponential wall heat flux function qwζ=ne−nζ/1−e−n, Br=0.1, Kn=0.0, α=0.02, β+=β−=2. (a) Constant wall heat flux n=0 and exponentially increasing wall heat flux n=−1,−2,−3,−5. (b) Constant wall heat flux n=0 and exponentially decreasing wall heat flux n=1,2,3.
(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). Figure7 shows the comparison between numerical and analytical results for α=0.02 and α=0.1.Figure 7
Comparison between numerical (red dotted lines) and analytical results (blue continuous lines) for the average NusseltNu as a function of Pe for qwζ=1 with Br=0 and Br=0.1, Kn=0. Two different values of α are considered, α=0.02 and α=0.1. The arrows indicate increasing values of α.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 numberNu∞ 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.This observation becomes more evident by plotting the same data as in Figure7 as a function of Pe/α instead of Pe (see Figure 8(a)).Figure 8
Comparison between numerical (red dotted lines) and analytical results (blue continuous lines) for the average NusseltNu as a function of Pe/α (a) and Peα/4 (b) for qwζ=1. Two different values of α are considered, α=0.02 and α=0.1. Figures 8(a)-8(b) refer to Br=0 and Br=0.1, Kn=0. Reference values (vertical lines) are Pe/α=102 and Peα/4=10−1. The arrows indicate increasing values of α.
(a)(b)Indeed, numerical (and analytical) results forNu collapse onto an invariant curve when plotted as function of Pe/α for low-intermediate Peclet values Pe<100; see Figure 8 (A), and the asymptotic behaviour sets for Pe/α≃102.On the other hand, for high Peclet values, streamwise conduction becomes negligible and numerical results forNu, 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 for0<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/α≥102 and Peα/4≤10−1 that implies 102α≤Pe≤4⋅10−1/α.Figure9 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.Figure 9
Region of validity of the analytical solution equation (35) (region with vertical bars below the red curve Pe=4⋅10−1/α) and of Nu∞ (triangular region below the red curve Pe=4⋅10−1/α and above the blue curve Pe=102α).Data reported in Figures8(a), 8(b), and 9 refer to a case in which the rarefaction effect is not present since Kn=0. The same analysis can be performed by including the effect of rarefaction (Kn≠0). The region of stability of the analytical solution for Kn≠0 almost coincides with that for Kn=0.The next section analyzes the influence of different parameters and of different wall heat flux functions on the average Nusselt numberNu∞, focusing exclusively on the range of validity of the analytical expression equation (34)fd35.


## 5. Influence of Wall Heat Flux Function and Transport ParametersBr and Kn

From a preliminary analysis of data reported in Figures6(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 forPe≃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)).For a constant wall heat functionqw=1, Nu∞ is a monotonically decreasing function of Pe/α, attaining values in the range Nu∞,Nu0.Figures10(a)and 10(b) show the upper bound Nu0 and the lower bound Nu∞ as a function of Kn for Br=0,0.05,0.1. The effect of dissipation is to lower Nu by decreasing both the upper and lower bound. The effect of Kn is to monotonically lower both the upper and the lower bounds in the absence of dissipation Br=0 while, for Br>0, the upper and lower bounds exhibit a maximum as a function of Kn.Figure 10
Upper boundNu0 and lower bound Nu∞ for Nu as a function of Kn for Br=0,0.05,0.1. Constant wall heat flux qwζ=1.
(a)(b)In the absence of dissipation, that is,Br=0, both Nu0ζ and Nu∞ζ become constant and independent of the wall heat flux function (see (37) and (38) for Br=0). For this reason, the values of Nu0 and Nu∞ reported in Figures 10 A-B for Br=0 represent the asymptotic limits 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 qw=1. For this exponential heat function, Nu can be evaluated analytically for Br=0, thus obtaining(39)Nu=2n+Pe/αAKn+Pe/αBKnn+Peα+log1+1−e−n−Pe/αPe/αBKnn+Pe/αAKn,AKn=14+χ12+2bKn,BKn=χ12+χ224.Figures11(a)and 11(b) show the behaviour of Nu as a function of Pe/α for Br=0 and for Kn=0,0.1. The arrow indicates increasing values of n, in the range n∈−3,3. The larger is n, the larger is the value of Nu for low-intermediate values of Pe. For n<0, Nu exhibits a neat minimum.Figure 11
Nu versus Pe/α for qwζ=ne−nζ/1−e−n. No viscous dissipation Br=0. The arrow indicates increasing values of n. n=0 (black broken line) corresponds to a constant function qwζ=1. The blue curves correspond to increasing wall heat flux. The red curves correspond to decreasing wall heat flux. (a) Kn=0; (b) Kn=0.1.
(a)(b)Therefore, in the absence of dissipation, a decreasing wall heat flux function (n>0) 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 Br=0.1. High values of Kn amplify differences between increasing and decreasing wall heat flux functions.Figure 12
Nu versus Pe/α for qwζ=ne−nζ/1−e−n with viscous dissipation Br=0.1. n=0 (black broken line) corresponds to a constant function qwζ=1. Blue lines with dots correspond to axially increasing wall heat flux, n<0. Red continuous lines correspond to axially decreasing wall heat flux, n>0. (a) Kn=0; (b) Kn=0.1.
(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 θ=θ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. 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 θbζ 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 qwζ 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 ofPe, and different channel aspect ratios α, it is possible to conclude that the analytical solution is reliable for Pe<1/α.The influence of various transport parameters, that is,Br and Kn, is analyzed in detail.The 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 than0.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 numberNu. 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. The influence of the cross section aspect ratio for rectangular microchannels (three-dimensional problem) on slip-flow heat transfer [41, 42] is not investigated.



---

*Source: 8884393-2020-08-25.xml*