Indian J Phys https://doi.org/10.1007/s12648-017-1154-4

ORIGINAL PAPER

Convective heat transfer for a gaseous slip flow in micropipe and parallel-plate microchannel with uniform wall heat flux: effect of axial heat conduction Y Haddout*, E Essaghir, A Oubarra and J Lahjomri Laboratory of Mechanics, Faculty of Science Ain Chock, University Hassan II, B.P. 5366, 20100 Maarif, Casablanca, Morocco Received: 20 July 2017 / Accepted: 18 October 2017

Abstract: Thermally developing laminar slip flow through a micropipe and a parallel plate microchannel, with axial heat conduction and uniform wall heat flux, is studied analytically by using a powerful method of self-adjoint formalism. This method results from a decomposition of the elliptic energy equation into a system of two first-order partial differential equations. The advantage of this method over other methods, resides in the fact that the decomposition procedure leads to a selfadjoint problem although the initial problem is apparently not a self-adjoint one. The solution is an extension of prior studies and considers a first order slip model boundary conditions at the fluid-wall interface. The analytical expressions for the developing temperature and local Nusselt number in the thermal entrance region are obtained in the general case. Therefore, the solution obtained could be extended easily to any hydrodynamically developed flow and arbitrary heat flux distribution. The analytical results obtained are compared for select simplified cases with available numerical calculations and they both agree. The results show that the heat transfer characteristics of flow in the thermal entrance region are strongly influenced by the axial heat conduction and rarefaction effects which are respectively characterized by Pe´clet and Knudsen numbers. Keywords: Heat transfer; Forced convection; Slip flow; Axial heat conduction; Self-adjoint formalism PACS Nos.: 44.10. ? i; 44.27. ? g; 44.05. ? e; 44.90. ? c

1. Introduction The analysis of heat transfer problems at the microscale is among the major themes in the Micro-Electro-Mechanical Systems (MEMS) and plays central role in transport of heat in fluid which flows through a microfluidic thermal system [1]. Microfluidic is referred to devices that have length typically between 1 lm and 1 mm such as micropipe, parallel-plate microchannel, micro-exchangers etc. This study is devoted to the classical nontrivial problem in convection–diffusion: The Graetz problem extended to microscale. Special consideration is given to flow of gas in a micropipe and in a parallel plate microchannel for slip flow regime. The basic parameter characterizing the slip flow in microchannels is the Knudsen number Kn = k/Dh, which is the ratio of the molecular mean free path k to the

hydraulic diameter Dh. When Kn is considered as the basic parameter, the interval of slip flow regime may be estimated as [1]: (10-3 B Kn B 0.1) and the gas becomes slightly rarefied. In this case, it can be shown that the flow and the heat transfer phenomena can be safely modeled using continuum conservation equations together with slip boundary conditions. These boundary conditions are essentially characterized by discontinuities of the velocity and the temperature fields at the fluid-wall interface due to rarefaction. The fundamental useful data concerning heat transfer rates in microchannels may be obtained from the determination of the Nusselt number. Therefore, the need for accurate evaluation of this number is essential in the design and manufacturing of microchannel. One of the objectives of this work is to determine thermal characteristics for a micropipe and a parallel plate microchannel applicable for a fully developed gas flow. These two configurations are basic geometries which represent the limiting cases of annuli and rectangular channels. The thermal

*Corresponding author, E-mail: [email protected]

Ó 2017 IACS

Y Haddout et al.

characteristics results for micropipe and parallel plate microchannel are of theoretical interest and practical importance in engineering and may be used, for example, in cooling of microelectronic components [2]. Although these are the simplest configurations, they could nevertheless give insight on the heat transfer in the thermal entrance region of various microchannels for more complex geometries. Although, the effect of axial heat conduction could be often ignored in the classical convection of liquids and gas in the channels of macroscopic dimension, this simplification is not justifiable in the slip flow regime. This is because, for microscale geometries, the Pe´clet number is small and the effect of axial heat conduction plays a significant role. The effects of axial heat conduction in the conventional forced convection inside channels for a fully developed flow, known as the extended Graetz problem, with different thermal boundary conditions, have been extensively investigated in the past. For more detail, relative of the effect of axial heat conduction, see the books of Weigand [3]. Previous investigations in macroscale using classical methods, revealed the existence of some difficulties: the associated eigenvalue problem for the elliptic energy equation is non-selfadjoint and yields an incomplete set of eigenvalues. Consequently, the eigenfunctions appearing in the series representation of the solution are not orthogonal and the determination of the related expansion coefficients by a classical method becomes difficult. For example, by using a separation of variables method and considering an infinite axial domain, Hsu [4] solved the problem of thermal entry region heat transfer for laminar pipe and parallel-plate channel flow with axial heat conduction and uniform wall heat flux. He used a cumbersome Gram–Schmidt procedure to determine the expansion coefficients appearing in the series representation of the solution. However, in series of innovative papers in macroscale [5–7], Papoutsakis et al., using a self-adjoint formalism, successfully solved this problem by decomposing the energy equation into a pair of first order partial differential equations. They showed that the analytical solution obtained is computationally simple and efficient as the solution for the parabolic problem. This method has been used by Weigand et al. [8–10] and Weigand [3] for the extended Graetz problem with constant heat flux wall boundary conditions to solve a forced convection of laminar and turbulent flow in a pipe and a parallel plate channel. Lahjomri et al. [11] also used a method of selfadjoint formalism to solve the extended Graetz problem in magnetohydrodynamic channel of Hartmann flow heated by uniform wall heat flux. The problem was treated by considering a semi-infinite axial domain (0 B x \ ? ?) and assuming uniform temperature at the inlet (x = 0).

This method has also been used for dipolar fluids flow by Akyildiz and Bellout [12]. The forced convection problem of slip flow with the effect of axial heat conduction in a microchannel has also been extensively investigated in the literature, by using either analytical methods [13–19] or numerical ones [20–27]. These studies show that the axial heat conduction plays significant role in the thermal entrance region of microchannels. A good review of experimental, numerical and analytical studies, concerning convective heat transfer of slip flow in microchannels, with various geometries and different heat boundary conditions, can be found in the work of Colin [28]. In this review, the author critically analyzed and compared the previous studies. Additionally, the recent theoretical works in the heat transfer for gaseous slip flow through microchannels in thermally and hydrodynamically fully developed region, with heat flux boundary condition, were undertaken in [29–31]. Recently, Haddout and Lahjomri [32] extended a self-adjoint formalism method to solve the forced convection of gaseous slip flow in microchannels with an isothermal heating section of finite length, including the effects of the axial heat conduction, viscous dissipation, and pressure work. This extension was done using a new three dimensions matrix operator and a suitable scalar product between two vectors in the Hilbert space. To our knowledge, the existing analytical solutions treating the problem of forced convective heat transfer in microchannels with axial heat conduction and a prescribed heat flux are often based on the variables separation method and all past works [13–16] considered a semi-infinite axial domain (0 B x \ ? ?), or others were only concerned with the thermally fully developed solution. In this configuration, the flow has been assumed, immediately in the beginning, hydrodynamically fully developed and enters with uniform temperature distribution T0 at x = 0. However, when the effect axial conduction is significant, this thermal inlet boundary condition and the hypothesis of the development of the flow may be unrealistic and unjustifiable, this results in an erroneously infinite value of Nusselt number at x = 0 [14–16, 23, 25]. Indeed, for small Pe´clet number, the heat is transferred upstream and exceeds the inlet cross section leading to a transverse variation of the fluid temperature in the region (x \ 0), which certainly will affect the heat transfer in the downstream heat-exchange section (x C 0). Consequently, the solution of the problem in downstream region considering an infinite axial domain (-? \ x \ ? ?) is different from that obtained with the assumption of uniform temperature at the inlet when only semi-infinite axial domain (0 B x \ ? ?) is considered, except in the thermally fully developed region for which the solutions for the two problems are identical. On the other hand, recall that the classical separation of

Effect of axial heat conduction Fig. 1 Geometrical configuration and coordinate system

variables method, used in microscale flow, leads to an eigenvalues problem with a non-orthogonal eigenfunctions. So, the associated eigenvalues could be, in principle, complex. This method is subject to uncertainties and therefore questionable, since it could not justify that the eigenvalues are real. The objective of the present contribution is to provide an analytical closed form expressions of the thermally developing of temperature field and Nusselt number for slip flow in the thermal entrance region of a micropipe as well as of a parallel plate microchannel, including axial heat conduction with prescribed uniform heat flux at the wall. In this study, the first-order slip model boundary conditions for velocity and temperature fields at the wall is employed. To avoid the difficulty of the formulation of the thermal inlet boundary condition, we consider an infinite axial domain (-? \ x \ ? ?) for which physically, the walls of microchannels are considered as adiabatic for x \ 0, while for x C 0 are uniformly heated. Furthermore, the fluid enters the microchannels with a flat temperature profile at x = - ?, since there is no axial conduction and heat transfer at this location. Beside this the hydrodynamical hypothesis of fully developed flow at the thermal entrance section x = 0 is entirely justifiable for the problem considered. The analytical solution is based on a selfadjoint formalism employed in macroscale by Papoutsakis et al. [6], Weigand et al. [8] and Weigand [3] and extended here for the first time to the slip flow regime with Neumann boundary condition. This powerful method which is built from purely physical arguments has the advantage to lead to a selfadjoint problem, even if the initial problem is apparently not a self-adjoint one. Therefore, it allows, on the one hand, to rigorously prove that the eigenvalues are real even though the eigenfunctions are not orthogonal and on the other hand, to determine more easily (than the separation of variables method [4, 33]) the explicit expressions of the related expansion coefficients. The effects of different dimensionless parameters, involved in the problem of the heat transfer characteristics will be quantified and studied in detail.

2. Statement of the problem We consider an hydrodynamically fully developed laminar flow in a long micropipe and a parallel plate microchannel with some specific velocity profile ux(r). The fluid enters with a uniform temperature distribution T0 at x = - ?. Figure 1 shows the geometrical configuration and the coordinate system. The variable n and the characteristic length L denote respectively both the radial coordinate r and the radius R for the flow in a micropipe, and the transversal coordinate y and half of the channel height h for the flow in a parallel plate microchannel, i.e., (n = r or y) and (L = R or h). The microchannels are thermally insulated for x \ 0 and are uniformly heated for x C 0 with heat flux qw (see Fig. 1). So, the density distribution of heat flux along the wall of microchannels has the following form:  qw if x  0 qw ¼ ð1Þ 0 if x\0 It is assumed that the fluid is flowing in the slip-flow regime, corresponding to Knudsen numbers in the range 10-3 B Kn B 10-1. The boundary conditions, representing two major characteristics of slip flow, are velocity-slip and temperature-jump conditions at the fluid-wall interface, which are important in microscale at ordinary pressure and in rarefied gases at low pressure. If we consider the firstorder slip model these two boundary conditions can be expressed respectively as [1]:  oux  ux;s ¼ bv k  ð2Þ on n¼L  oT  Ts ¼ Tw  bt k  ð3Þ on n¼L where ux,s, Ts and Tw are, respectively, the slip velocity, the fluid temperature at the wall and the unknown wall temperature, with bv ¼ ð2  Fv Þ=Fv and bt ¼ ð2  Ft Þ=Ft ð2c=Pr ðc þ 1ÞÞ such that Fv and Ft are the tangential momentum and thermal accommodation coefficients. c is the ratio of specific heats and Pr is the Prandtl

Y Haddout et al.

number. Experimentally, Fv and Ft have been found to vary between 0.79 and 1 with most of the values are close 1 [1]. Our model assumes, as first approximation, that the variations of temperature, pressure and density are not too large (DT  T0, DP  P0 and Dq  q0) relative to the reference state of the fluid (q0, P0, T0) evaluated at large axial distance upstream at the inlet (x ? - ?). Therefore, these variations have negligible effects on the hydrodynamical flow and all physical properties can be assumed as constant in the intervals of the temperature and pressure considered. Although, the viscosity could be varied through the channel and because the assumption of not too large variations of the temperature DT  T0 as discussed above, this variation will be small. Indeed, experience shows, however, that good results are often obtained, for example, by assuming the viscosity as constant and by assigning to it the value corresponding to an average temperature of the fluid. For example, for air at standard pressure the viscosity is l = 1.85710-5 Pa s at 300 K, and l = 1.93510-5 Pa s at 320 K given a relative change of about 4%. So, for an average temperature of the fluid at 310 K, the viscosity can be evaluated as: l = 1.88910-5 Pa s. In addition, the velocities are generally small enough compared to the speed of sound, to neglect the effects of compressibility in microchannel. It should be noted that, if the value of the Mach number along the channel is less than 0.2 and for relatively low Reynolds number and small hydraulic diameter, and when the pressure drop is less than 10% of the outlet pressure, the inertial effects can be neglected [1, 34, 35]. Therefore, an hydrodynamical fully developed flow is possible if the above conditions are satisfied. However, as pointed out by Duan [35], for practical engineering applications, the fluid can be considered as incompressible for small pressure drops, which is approximately consistent with experimental results of Jang and Werley [36]. To have an idea of the orders of magnitude of the physical parameters in the engineering applications at microscale for gaseous slip flows in microtube, for example for air, nitrogen, helium and argon, at standard atmospheric and experiment conditions [37], we estimate the physical properties and parameters for small velocities as: q & 1 kg/m3, Cp & 103 J/(kg K), l & 210-5 Pa s, k & 310-2 W/(m K), Pr & 0.7, k * 85 nm, Dh * 1–100 lm, Um * 1–10 m/s, qw * 105–106 W/m2, ReDh & 0.05–50, PeDh & 0.033–33, Kn & 0.001–0.085, Br & 410-7–410-4. It is clear from these values that the Brinkman number Br is very small, so the effects of viscous dissipation and pressure work can be neglected. Therefore, the theoretical assumptions of our analysis will be used in accordance with the estimate of the conditions of the experiment; and in the Sect. 4 below, we simulate our

analytical solution by choosing the numerical values of the modeling parameters belonging the range given above. Under previous assumptions of incompressible flow, the steady-state energy equation, including axial heat conduction, neglecting viscous dissipation and pressure work and considering of the symmetry of the flow about the x-axis, for 0 \ n \ L and all x, is given by:   oT o2 T k o p oT ¼k 2 þ p n q Cp ux ðnÞ ð4Þ ox ox n on on The index p appearing in Eq. (4) is equal to 0 for a parallel plate microchannel and equal to 1 for a micropipe. The boundary conditions associated to Eq. (4) are oT ¼ 0 for n ¼ L and x\0 on oT qw ¼ for n ¼ L and x  0 on k oT ¼ 0 for n ¼ 0 and all x on lim T ¼ T0

ð5Þ ð6Þ ð7Þ ð8Þ

x!1

The velocity distribution ux(n) in Eq. (4) could be calculated from the momentum equation by neglecting inertial effects and assuming hydrodynamically fully developed flow   op l d p dux ¼ p n ð9Þ ox n dn dn with the boundary conditions dux ¼0 dn

for n ¼ 0

ux ¼ bv k

dux dn

for n ¼ L

ð10Þ ð11Þ

Additionally, the conservation of mass in integral form Z L Um Lpþ1 ¼ ðp þ 1Þ ux np dn ð12Þ 0

must be satisfied. Using the Eqs. (10)–(12), by eliminating the pressure gradient, the fully developed velocity distribution in the micropipe and parallel plate microchannel can be derived in term of the mean velocity Um as   Um n2 i 1  2 þ 23p bv Kn ð13Þ ux ðnÞ ¼ h 2 L 3p b Kn þ 2 v pþ3 where Kn = k/Dh is the Knudsen number, with Dh = 22-pL denoting the hydraulic diameter. By introducing the dimensionless slip velocity us = ux(L)/ Um, the dimensionless velocity profiles can be written:

Effect of axial heat conduction

uðgÞ ¼

 ux ð1  us Þ  1  g2 ¼ us þ ðp þ 3Þ 2 Um

ð14Þ

n b Kn ðp þ 3Þ and us ¼ uð1Þ ¼ p2 v g¼ L ð2 þ bv Knðp þ 3ÞÞ

with

ð15Þ In the limiting case of no slip, Kn ? 0 and us = 0, we find, for the two configurations considered p = 0 and p = 1, the classical expression of the Hagen–Poiseuille flow. From Eq. (14), we can deduce the expression of the friction coefficient Cf and the maximum value of u(g) at g = 0. Cf ¼

2sw 23p ðp þ 3Þð1  us Þ ¼ 2 ReDh qUm

ð16Þ

ð1  us Þ ð17Þ Umax ¼ us þ ðp þ 3Þ 2   where sw ¼ ldu dn n¼L denotes the wall shear stress and 2p ReDh ¼ 2 ReL is the Reynolds number based on the hydraulic diameter (ReL ¼ qUm L=l). We can also define the Poiseuille number Po for the hydrodynamically developed flow as the product of the friction coefficient and the Reynolds number: Po ¼ Cf ReDh ¼ 23p ðp þ 3Þð1  us Þ

ð18Þ

The expression of Cf agrees well with that obtained by Shah and London [38] in the case of no slip flow (us = 0) for a parallel-plate channel, for which Cf = 24/ReDh and for a macro-pipe, Cf = 16/ReDh. By introducing the following dimensionless variables and parameters: n g¼ L

x 1 n¼ L PeL qCp Um L PeL ¼ k

ux uðgÞ ¼ Um

T  T0 H¼ L qw =k

ð19Þ

into Eqs. (4) and (5)–(8), the energy equation and the boundary conditions can be cast into the following form:  oH 1 o oH 1 o2 H ¼ p gp uðgÞ þ 2 2 ð20Þ on g og og PeL on ðfor 0\g\1 and all nÞ oH ¼0 og

for g ¼ 1; n\0

ð21Þ

oH ¼1 og

for g ¼ 1 and n  0

ð22Þ

oH ¼0 og

for g ¼ 0 and all n

ð23Þ

lim H ¼ 0

n!1

3. Analytical solution Among the different methods which may be used for solving the elliptic energy Eq. (20), there exists a powerful method based on the self-adjoint formalism of functional analysis theory, which consists to decompose this equation into two first order partial differential equations by defining the dimensionless axial energy flow as an auxiliary function following a purely physical argument as used earlier in macroscale, see Refs. [3, 6, 8]:  Z g 1 oH 0p 0 0 uðg ÞH  2 ð25Þ g dg /ðn; gÞ ¼ PeL on 0 In this study, we extend this method in the slip flow regime to derive an analytical solution of the energy Eq. (20) with boundary conditions (21)–(24). For a description and details on the principle of this method, see [3, 6, 8]. Here, we only describe the outline of the method. Differentiation of Eq. (25) with respect to g and n leads to the following system oH Pe2 o/ ¼ Pe2L uðgÞH  pL on g og

ð26Þ

o/ oH ¼ gp on og

ð27Þ

Equations (26) and (27) now may be written in matrix notation as !

! oF ¼ LF on

ð28Þ !

with the two-component vector F ðn; gÞ and the operator L given by 2 3 " # Hðn; gÞ Pe2 o ! Pe2L uðgÞ  gpL og 6 7 ð29Þ 5; L ¼ F ¼4 o gp og 0 /ðn; g The most remarkable aspect of L is that it gives rise to a selfadjoint operator even if the original equation is nonselfadjoint. This fact is of course dependent on the sort of inner product between two vectors, which will be used. Indeed, if we define the inner product between two vectors ! ~ / ¼ ½/1 ðgÞ; /2 ðgÞt and w ¼ ½w1 ðgÞ; w2 ðgÞt as 

 Z 1 p ! ! g 1 / ðgÞw1 ðgÞ þ p /2 ðgÞw2 ðgÞ dg /; w ¼ g Pe2L 1 0 ð30Þ and the following domain in the Hilbert spaces for L,

ð24Þ

Y Haddout et al.

n o ~ 2 H; / ð0Þ ¼ / ð1Þ ¼ 0 DðLÞ ¼ ~ / 2 H : L/ 2 2

ð31Þ

Therefore, the solution can be obtained, using an expansion theorem for any vector ~ F 2 H: 2D D E E D E 3 þ  ~ ~ ~ ~ 1 1 F; ~ /n F; ~ /n F ; / X X n þ 7 6 ~ F¼ /n ¼ /n þ 2 ~ /n 5 4 þ 2 ~ 2 ~  n¼1 ~ n¼1 /n /n /n ~ ~ ð32Þ where ~ /n is the norm of ~ /n . Then, for a given fully developed velocity profile u(g), we can show that the solution can be represented as an  infinite series of eigenfunctions /þ n1 and /n1 (i.e. the first  components of vectors ~ /n ) by applying the theorem (32) 1 X An expðk2n nÞ/þ ð33Þ Hðn; gÞ ¼ n1 ðgÞ n\0 n¼1

Hðn; gÞ ¼ QðgÞ þ ðp þ 1Þn þ

1 X

  Bn exp b2n n / n1 ðgÞ

n¼1

n0 ð34Þ where Z g Z g0 ðpþ1Þ 1 QðgÞ ¼ þ g00p uðg00 Þdg00 dg0 2 0p PeL 0 g 0 "Z # Z 0 Z 1 g 1 g 00p 00 00 0 p ðpþ1Þ uðgÞg g uðg Þdg dg dg 0p 0 0 g 0 and the expansion coefficients An and Bn can be obtained explicitly as: " !#1 /þ0 3 d n1 ð1Þ An ¼ 2 kn ; dkn k2n ð35Þ " !#1 0 d / ð1Þ n1 Bn ¼ 2 b3n dbn b2n The kn and bn designate the real eigenvalues associated  with the eigenfunctions /þ n1 and /n1 respectively. The þ  eigenfunctions /n1 and /n1 are then the solutions of the following differential equations:      d p d/n1 p  ln g  uðgÞ / ð36Þ þ g ln n1 ðgÞ ¼ 0; dg dg Pe2L subject to the boundary conditions /0 n1 ð0Þ

¼0

/0 n1 ð1Þ ¼ 0

ð38Þ

þ 2 þ The explicit forms for / n ðgÞ and /n ðgÞ, with ln ¼ kn and 2  ln ¼ bn can be obtained as:

  bn 2 pþ1 uð0Þ b3 p þ 1  bn  n2 ; ; bn tg2 tg M / n1 ðgÞ ¼ exp  4 4t 2 4PeL t 2

ð39Þ 

 ikn 2 pþ1 uð0Þ ik3n p þ 1 2  ik þ ; ik M tg /þ ðgÞ ¼ exp  tg ; n n n1 4 4t 2 4Pe2L t 2

ð40Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffi pffiffiffiffiffiffiffi where t ¼ ðp þ 3Þð1  us Þ 2, i ¼ 1, u(0) denotes the value of dimensionless velocity at the centerline, and M(a, b; z) is the Kummer confluent hypergeometric function of the first kind [39]. Now that the temperature field is known, we can compute two quantities of practical importance, the bulk temperature: RL R1 p Hugp dg 0 Tux n dn Tb ¼ R L ð41Þ Hb ¼ R0 1 p p 0 ux n dn 0 ug dg and the local Nusselt number: qw ¼ 22p Nu ¼ kðT T Þ w

 ¼

b

Dh 2p

2 0

1 Hw Hb

qw qw ðHw  Hb Þ

Introducing the predicted temperature field in the flow, we obtain the following expressions for the bulk temperature and the wall temperature which this latter is calculated from Eq. (3) in all regions of the microchannels: n\0: Hb ¼ ðp þ 1Þ

1 X



An expðk2n nÞ

n¼1

Hw ¼

1 X

k2n Pe2L

Z

1

g 0

p

/þ n1 dg



An expðk2n nÞ/þ n1 ð1Þ

ð43Þ

ð44Þ

n¼1

n  0: Hb ¼

ðp þ 1Þ þ ðp þ 1Þn  ðp Pe2L  2 Z 1  1 X b þ 1Þ Bn expðb2n nÞ n2 gp / dg n1 PeL 0 n¼1

Hw ¼ Qð1Þ þ ðp þ 1Þn þ ð37Þ

ð42Þ

if n  0 if n\0

1 X

ð45Þ

Bn expðb2n nÞ/ n1 ð1Þ

n¼1

þ 22P bt Kn

ð46Þ

The Nusselt number is zero for n\0 because there is zero

Effect of axial heat conduction Table 1 The first ten eigenvalues kn, bn and the corresponding expansion coefficients An and Bn for a micropipe (p = 1) obtained from analytical method using Maple and numerical method for PeDh = 1 and Kn = 0.1 bn

n

kn

Bn

An

Analytical method (using Maple) 1

1.332055821

0.186040308

2

1.833875310

- 0.071453865

0.4995995387

3

2.222821646

0.040174912

1.912966328

0.063955669

4

2.552544506

- 0.026551884

2.288504236

- 0.037218024

5

2.844047822

0.019206339

2.609940728

0.025046300

6 7

3.108213262 3.351557087

- 0.014719661 0.011744193

2.895671133 3.155517529

- 0.018320466 0.014147443

8

3.578351860

- 0.009652016

3.395472695

- 0.011349642

9

3.791580291

0.008115030

3.619516697

0.009366516

10

3.993422977

- 0.006946783

3.830454230

- 0.007900621

1.439081021

7.966181274 - 0.151804479

Numerical method 1

1.332055737

0.186040338

0.499599530

7.966181853

2

1.833874852

- 0.071453916

1.439080939

- 0.151804504

3

2.222820451

0.040174975

1.912965887

0.063955714

4

2.552542138

- 0.026551956

2.288503073

- 0.037218082

5

2.844043782

0.019206419

2.609938410

0.025046367

6

3.108206996

- 0.014719748

2.895667162

- 0.018320542

7

3.351547990

0.011744288

3.155511354

0.014147527

8

3.578339282

- 0.009652117

3.395463714

- 0.011349733

9

3.791563537

0.008115137

3.619504260

0.009366613

10

3.993401313

- 0.006946895

3.830437645

- 0.007900724

heat flux at the wall. For n  0 the Nusselt number is given by 8 ðpþ1Þ 2P > < Qð1Þ þ 2 bt Kn  Pe2L 2p h  2 R i Nu ¼ 2 1 P bn 1 p  > :þ Bn expðb2n nÞ / n1 ð1Þ þ ðp þ 1Þ Pe2 0 g /n1 dg L

n¼1

91 > = > ;

ð47Þ The asymptotic value of the local Nusselt number Nuas, corresponding to the thermally fully developed region can be derived from Eq. (47) by taking the limit n ! þ1 Nuas ¼

22p Qð1Þ þ 22P bt Kn  ðpþ1Þ Pe2

ð48Þ

L

Nuas is dependent of Kn, bv and bt.

4. Results and discussion 4.1. Calculation procedure, accuracy and validation For this analytical study, we calculate the eigenvalues and the eigenfunctions by two methods: first, using numerical

method and specifically high-accuracy numerical program written in FORTRAN to solve a differential Eq. (36) with boundary conditions (37) and (38). Second, to validate and examine the accuracy of this numerical code we compare the results obtained with the analytical method by using the explicit forms of the eigenfunctions [Eqs. (39) and (40)] with the help of the symbolic algebra software MAPLE. One advantage of the numerical method is that it can deal with any arbitrary velocity profile and the code can be used for the resolution of others microchannel problems such as: slug flow, Hartmann flow, flow in porous media, electroosmotic flow, etc. The numerical method is based on a Newton–Raphson algorithm coupled with a bisection method [40]. Simultaneously the eigenfunctions were determined by solving the differential Eqs. (36)–(38) by using a finite difference method with standard second-order difference for approximating the derivatives. Then a Newton’s method is used to solve the resulting system of linear equations which is a tridiagonal matrix. Each eigenfunctions may be normalized by the requirement that they satisfy /n1 ð0Þ ¼ 1 and /0n1 ð0Þ ¼ 0. These two conditions are used to construct the initial approximation of the solution vector which is obtained from Taylor expansion

Y Haddout et al.

around a point g ¼ 0 for the determination of the first vector component and iteratively by using finite difference equation to obtain the other components of the solution vector to attain a point g ¼ 1. This initial approximation of the solution vector is necessary for starting Newton’s iterations and then constructs a sequence of the solution vector by solving multiple times a tridiagonal matrix equation with the help of Thomas algorithm to converge efficiently with high accuracy to the final solution of the boundary value problem. Thus, starting with an arbitrary value of the nth eigenvalue one can step forward from g ¼ 0 to g ¼ 1 and vary the value of kn or bn to satisfy the required fundamental boundary condition: /0n1 ð1Þ ¼ 0. The results for the eigenvalues kn, bn and the corresponding expansion coefficients An and Bn, for a micropipe (p = 1), obtained from analytical method (using Maple) and numerical method for Pe´clet number, PeDh ¼ 1, (PeDh ¼ 22p PeL ), Kn = 0.1 are compared in Table 1. The results found for the eigenvalues and the corresponding expansion coefficients agree well within a relative error less than 10-7 for the eigenvalues, and less than 10-6 for the expansion coefficients. Figure 2 shows the effect of different truncation orders (N = 10, 20, 100, and 1000) of the series on the accuracy of the distribution of local Nusselt number Nu, for Kn = 0.1, PeDh = 1 and PeDh = 10. This result shows that, for example, for PeDh = 1, by increasing N the accuracy of the calculation increases, and only N = 100 terms are required to attain the convergence of Nu around

n ¼ 103 . Clearly, beyond this value of N, the distribution of Nu did not change. The validity of our analytical solution could be checked, in the simplified case of continuum flow Kn = 0, by comparing the predicted variation of the temperature profiles for various axial distance n in a parallel plate channel (p = 0) with analytical results of Weigand et al. [8]. This comparison is illustrated in Fig. 3 for Pe´clet number, PeDh = 5. Clearly, that the agreement between both calculations is excellent. Table 2 depicts the variation of fully developed Nusselt number, Nuas which is independent of the magnitude of Pe´clet number as predicted by Eq. (48), as a function of Knudsen number for a micropipe and parallel-plate microchannel. This table shows also a comparison of our analytical results, in the range of applicability of the slip flow regime, with numerical results of Aydın and Avcı [21] who have neglected the effect of axial heat conduction in a micropipe. It can be seen again that both solutions are in good agreement and therefore supports the validity of our analytical solution, within a relative error less than 10-5. Furthermore, the results of Table 2 show that in the slip flow regime, the asymptotic Nusselt number decreases monotonically with increasing Knudsen number. It is also noted that for Kn = 0, which represents the classical noslip continuum flow in which both the velocity-slip and temperature-jump conditions are absent, one recovers the classical values of Nuas for macro-pipe flow (e.g., 4.3636) and for parallel plate channel flow (e.g., 8.2353).

Fig. 2 Effect of different truncation orders N on the accuracy of the distribution of local Nusselt number in micropipe for PeDh = 1, PeDh = 10 and Kn = 0.1

Fig. 3 Comparison of temperature profiles for a parallel plate channel with analytical results of Weigand et al. [8] in continuum flow

Effect of axial heat conduction Table 2 Variation of asymptotic Nusselt number with Knudsen number for a micropipe and parallel-plate microchannel and comparison with numerical results of Aydın and Avcı [21] Kn

Nuas for micropipe (p = 1)

Nuas for parallel-plate microchannel (p = 0)

This work

Aydın and Avcı [21]

This work

0

4.363636363

4.364

8.235294118

0.02

4.070750059

4.071

6.819172479

0.04 0.06

3.748762889 3.438582539

3.749 3.439

5.724212768 4.894216406

0.08

3.155646984

3.156

4.256894521

0.1

2.903659447

2.904

3.757486137

4.2. Main hydrodynamic results Figure 4 represents the dimensionless velocity profiles calculated for different values of Knudsen number Kn. The increase in Knudsen number Kn results in an increase of slip velocity at the wall, while the fluid velocity decreases in the central core of the microchannels. This behavior can be explained by the conservation of mass flow rate. The velocity profile corresponding to the zero value of Kn represents the usual parabolic profile of Poiseuille flow in a macrochannels. Figure 5a, b show the variation of the friction coefficient and the slip velocity, respectively, as a function of the Knudsen number Kn in a parallel-plate microchannel and a

(a)

micropipe. Figure 5a clearly indicates that the friction coefficient decreases monotonically with increasing Knudsen number. For example, by increasing Kn from the value 0–0.1 in the case of a parallel plate microchannel, we can note a relative decrease of about 54.58%, whereas in the case of a micropipe we have a relative decrease of 44.5% compared to the case without slip velocity. Figure 5b shows that the slip velocity increases with increasing the rarefaction effect, which gives a relative increase of 54% for parallel plate microchannel and 44% for micropipe, compared to the no-slip case. Consequently, this increase in the slip velocity leads to a decrease in the velocity gradient at the wall, which is the reason for the reduction in the friction coefficient in microchannels.

(b)

Fig. 4 Dimensionless velocity profiles for various values of Knudsen number Kn (a) for a micropipe and (b) for a parallel-plate microchannel

Y Haddout et al.

(a)

(b)

Fig. 5 Effect of Knudsen number (a) on the friction coefficient and (b) on the slip velocity for a microchannels

4.3. Main results of heat transfer Typical local temperature profiles H ¼ ðT  T0 Þ=ðqw L=kÞ for thermally developing flow at various axial distances n ¼ x=ðL PeL Þ in a micropipe, for different values of Knudsen number, with Pe´clet number PeDh = 1; 10 and also in the case of negligible axial heat conduction (parabolic calculation PeDh = ?) are presented in Fig. 6a–c, respectively. These figures show that the temperature profiles retain globally the same shape as in the absence of rarefaction effect. We observe that for large values of Pe´clet number (Fig. 6b, c), near the entrance region, precisely for axial distance n B 0.01 and for PeDh C 10, the fluid temperature in the central core of the micropipe is close to T0 and is not influenced by Knudsen number except near the wall, for which the thermal boundary layer dth develops increasingly until a value of n 0:1. It can be also seen that this thickness decreases as Kn increases. This effect leads to a decrease of the fluid temperature at the wall when Kn increases. For small Pe´clet number (PeDh = 1), the effect of axial heat conduction becomes pronounced near the entrance region. Therefore, the heat which is coming from the wall to the fluid in the downstream region, a part of this heat is transferred towards the upstream region in which the walls of the microchannels are insulated. The results also show that at any cross-section of the flow, the Knudsen number does not affect the fluid temperature near the central core of the flow, while its effect becomes important near the wall. This can be explained by the effect of slip velocity at the wall (see

Fig. 5b). Indeed, as Kn increases, the slip velocity increases, therefore, the fluid particles near the wall are less heated than in the case of continuum flow. The effect of Knudsen number on the distribution of the bulk temperature Hb and on the wall temperature Hw in a micropipe for PeDh = 1, PeDh = 10 and PeDh = ? are further presented in Fig. 7a–c, respectively. It is clear that for the small value of Pe´clet number (PeDh = 1) and for any fixed value of Kn, both Hb and Hw tend to zero for a larger negative axial distance (i.e. n ¼ 25). At this location, there is no axial conduction and heat transfer. By increasing the value of Pe´clet number to 10, the effect of axial heat conduction becomes small and the temperatures on the upstream region decreases. The dimensionless temperatures profiles start closely to n ¼ 0 with an approximate value of zero and vanish completely in the extreme case of parabolic calculation. Figure 7a, b corresponding to PeDh = 1 and PeDh = 10 respectively, show that the Knudsen number has no large influence on the bulk temperature, so that for the overall heat balance, the linear increase of Hb must follow in the downstream region. This behavior is often encountered for prescribed uniform wall heat flux boundary condition. For PeDh ? ?, Hb increases also linearly with n and is not affected by Knudsen number. These results could be understood by considering the expression of Hb given by Eq. (45), where for large PeDh, the term containing the effect of rarefaction disappears. The results of Fig. 7 show that, in the fully developed region, because of the boundary condition of uniform heat

Effect of axial heat conduction

(b)

(a)

(c) Fig. 6 Effect of Knudsen number on the temperature profiles (a) for PeDh = 1, (b) for PeDh = 10 and (c) for PeDh = ? (Parabolic calculation) with various axial distance n

flux, the slopes of the two curves of Hw and Hb are identical (i.e. dHw =dn ¼ dHb =dn). Thus, the temperature profiles of Hw and Hb are both continuous, parallel, never intersect (i.e. Hw  Hb 6¼ 0Þ and are different for finite Pe´clet number at n ¼ 0 except for the case of infinite Pe´clet number for which Hw ¼ Hb at n ¼ 0. Hence, from Eq. (42

or 47), the curve of Nusselt number for finite PeDh must be discontinuous but finite at n ¼ 0 (see Fig. 8 below), and presents a singularity at the point n ¼ 0 in the particular case when PeDh = ?. This singularity does not exist for Nu at finite value of Pe´clet number. From Fig. 7, we can also distinguish two types of variations of the wall

Y Haddout et al.

(a)

(b)

(c)

Fig. 7 Bulk and wall temperatures with respect to axial distance n for various values of Knudsen number Kn and (a) for PeDh = 1, (b) for PeDh = 10 and (c) for PeDh = ? (Parabolic calculation) in micropipe

temperature in the direction of the flow. In the thermally developing region, Hw presents a non-linear variation with n, thus the axial temperature gradient at the wall will depend on the axial distance, while in the thermally fully developed region the variation of Hw is linear, therefore, the axial temperature gradient at the wall becomes constant in this region. Figure 7 indicates also that for a given Pe´clet number, the Knudsen number affects the wall temperature Hw . As

seen in Fig. 7, for PeDh B 10 an increase in Kn leads to an increase of Hw in downstream region, while a decrease with Kn is observed in upstream region. These results can be understood by considering the expression of Hw given by Eq. (3), Hw ¼ Hs þ Hj where Hs and Hj are the dimensionless fluid temperature at the wall and the temperature jump, respectively. Indeed, as we have seen previously in Fig. 6, the fluid temperature at the wall Hs decreases with increasing Kn. The temperature jump,

Effect of axial heat conduction

Fig. 8 Influence of Knudsen number Kn and Pe´clet number PeDh on the local Nusselt number in micropipe

which is defined as Hj ¼ 2 bt Kn in the downstream region (n  0), increases with Kn; while in the upstream region Hj is zero because the walls are considered as thermally insulated. Therefore, the wall temperature in the downstream region increases due to an increase of temperature jump with Kn. For a parabolic calculation, the effect of Knudsen number on Hw , in the downstream region, is like to that for finite Pe´clet number. Figure 8 shows the variation of the local Nusselt number Nu with the axial distance for different values of Pe´clet and Knudsen numbers in a micropipe. It is noted that, for a given value of Kn and PeDh, Nu decreases monotonically with n in the thermal entrance region until attaining its asymptotic value corresponding to the established thermal regime. These asymptotic values for different Kn are shown in the figure. We note also that for a fixed value of Kn and for finite value of PeDh, the different curves of Nu have an inflection point and intersect at n 0:03 for which Nu has the same value for each Kn. This behavior has been also observed by Hennecke [41] and our result agrees with this author for Kn = 0 in the continuum flow (n in Hennecke’s definition is equivalent to n=4 in our definition). We also note that the results found here, in the developing region, for Nusselt number disagree with those obtained in Refs. [14–16, 23, 25] who considered a semi-infinite axial domain, except for the results obtained in the thermally fully developed region for which the asymptotic values of Nu are identical. This disagreement is due to the restrictive assumption of uniform temperature boundary condition at the inlet at n ¼ 0 by these authors. It should also be noted

that the curve of the particularly value PeDh = 1 and for all values of Kn, presents another point of inflexion located approximately at n 0:15. As discussed above, for a fixed value of Kn, Nu has a finite value at n ! 0 and increases with PeDh. For a given PeDh, the effect of Knudsen number leads to decrease the local Nusselt number. Indeed, as Kn increases, the temperature jump increases which leads to increase the wall temperature (as discussed above in Fig. 7) and because the bulk temperature is insensible to Kn. Consequently, from the definition of Nu (i.e. Nu ¼ 2=ðHw  Hb Þ), Nu decreases with Kn. To quantify the importance of the effects of axial heat conduction, temperature jump and slip velocity, we present in Figs. 9 and 10 the percentages variations of Nusselt number with axial distance n in the presence and in absence of these three effects. For  this, we define  these percentages variations as: DEC ¼ Nuelli  Nupara  Nuelli which represent the relative deviation from an elliptic calculation (with axial conduction PeDh = 1) and a parabolic calculation (without axial conduction) to characterize the effect of axial conduction; Ps ¼ NuðKn 6¼ 0Þ=NuðKn ¼ 0Þ for bt = 0 to characterize the effect of the slip velocity in the absence of temperature jump; and Pjs ¼ NuðKn 6¼ 0Þ=NuðKn ¼ 0Þ for bt = 1.667 to characterize the effect of temperature jump in the presence of the effect of slip. The results of Fig. 9, clearly show that the deviation DEC is significant in the developing region, indicating the importance of axial conduction in this region. Thereafter, this deviation decreases until the point of inflexion n & 0.03 is reached and vanishes at this location, since the

Fig. 9 Relative deviation between the local Nusselt number for an elliptic (with axial heat conduction PeDh = 1) and parabolic calculation (without axial heat conduction) in micropipe

Y Haddout et al.

Fig. 10 Effect of Knudsen number on the percentage variations Ps and Pjs along the micropipe in the heating section, in the presence and in the absence of temperature jump for PeDh = 1

Nusselt number is independent of Pe´clet number there. Beyond this location, the deviation increases and attains a maximum approximately at n & 0.15, which corresponds to the point of inflexion for the curve of Nu for PeDh = 1. After this, DEC vanishes completely in the thermal fully developed region in which the effect of axial conduction becomes negligible. Figure 9, indicates also that effect of axial heat conduction diminishes by increasing the effect of the rarefaction. The results of Fig. 10 illustrate the percentage variations Ps and Pjs defined above, with different values of Knudsen number Kn in the range of the slip flow regime for the micropipe and for PeDh = 1. Where bt = 1.667 is a typical value for air corresponding to Ft = 1, c = 1.4 and Pr = 0.7. Clearly, the slip velocity and the temperature jump have opposite effects on the heat transfer. As Kn increases from the value 0–0.1, for example at n = 0.001, the slip velocity increases the heat transfer of about 40% of relative increase, while the temperature jump tends to decrease it by giving a relative decrease of about 47%. Therefore, we can conclude that the decrease of heat transfer in microchannel flow is mainly due, to the effect of temperature jump which attenuates its increase by the slip velocity.

5. Concluding remarks In this paper, the analytical solution of the forced convection heat transfer for gaseous slip flow regime through a micropipe and a parallel plate microchannel including axial

heat conduction with uniform wall heat flux is obtained by using a powerful method of a self-adjoint formalism. The analytical solution obtained is valid for all values of Pe´clet number and considers an infinite axial domain. This solution is different from those found earlier in microchannel flow research which consider a semi-infinite axial domain with non-justifiable assumptions of hydrodynamically developed flow immediately in the beginning x = 0 and uniform temperature boundary condition at the inlet. Our solution is computationally more efficient than numerical methods. It solves efficiently the singularity and the convergence difficulties that all numerical methods face near boundary discontinuities for the wall heat flux at x = 0. This is clearly demonstrated by the results obtained near the entrance heating section (see Fig. 2). It is clear, that the solution methodology could be extended to arbitrary velocity profile and for any prescribed heat flux distributions. Although the solution has some limitation, neglecting the effect of compressibility, it has the potential to provide qualitative information about the physical phenomena involved in the heat transfer with axial heat conduction and rarefaction and can serve as verification for numerical models in compressible flow. It should be noted that in this study, we are concerned only with the first order slip model boundary conditions given by Eqs. (2) and (3) in the slightly rarefied slip flow regime 0.001 \ Kn \ 0.1 for which the Navier–Stokes equations and these boundary conditions are both valid and are first order accurate in Knudsen number Kn (i.e. O(Kn)). In the case of isoflux boundary condition, the effect of thermal creep appears as a second order in Knudsen number (i.e. O(Kn2)) in the slip boundary conditions which is generally used for the second order slip model. This boundary condition can be written in dimensionless form as [1]:  ou 3 c  1 Kn2 dHw 2p 2  Fv us ¼ 2 Kn  þ 242p og g¼1 2p c Br dn Fv ð49Þ  2

where Br ¼ Ec Pr ¼ lUm qw L is the Brinkman number and Ec ¼ k Um2 Cp qw L is the Eckert number. For a given Br, the term due to the thermal creep can play significant role at large values of Kn especially in transitional flow regime (i.e. 0.1 \ Kn \ 10), where the Navier–Stokes equations are not valid there and in which second order or higher slip boundary conditions are applicable. Beside this, we recall that in this study, we have completely neglected the effect of viscous dissipation, pressure work and shear work at the wall (see for example Ref. [32]) which means that Brinkman number tends to zero Br ? 0. The term of thermal creep becomes however dominant or infinitely larger compared to the first term

Effect of axial heat conduction

which is physically non-acceptable. Thus, the use Navier– Stokes and energy equations with second-order term in boundary conditions is inconsistent with the basic assumption of hydrodynamically fully developed flow and the approximation of negligible viscous dissipation and pressure work which are both of order of Br [32]. Therefore, maintaining the second-order term in the boundary condition (49) does not seem justified with the hypothesis made in the present analysis. Indeed, near the thermal entrance region (i.e. thermally developing region), the tangential axial temperature gradient due to the thermal creep effect, is a priori unknown and non-constant and must vary with axial distance n (or x) of the direction of the flow. This variation leads necessarily to a variation of the slip velocity with x and consequently also for the velocity field of the fluid flow, which is also inconsistent with the starting hypothesis of hydrodynamically fully developed flow. In problem formulation, because of the above reasons, we did not consider the second-order term in Kn2 in the boundary conditions for characterizing the effect of slip. However, we note that an analytical solution using the second order slip model boundary conditions assuming both hydrodynamically and thermally fully developed flow with non-negligible viscous dissipation can be easily obtained [42], since the tangential axial temperature gradient at the wall is constant in thermal fully developed region. In contrast, a difficulty arises for the solution of the problem in the thermal entrance region with the effect of thermal creep and necessitates further research in the future.

References [1] G Karniadakis, A Beskok and N Aluru Microflows and Nanoflows: Fundamentals and Simulation (New York: Springer) Ch 2 p 60 Ch 10 p 400 (2005). [2] J M Ha and G P Peterson J. Heat. Transf. 120 1064 (1998). [3] B Weigand Analytical Methods for Heat Transfer and Fluid Flow Problems (Berlin Heidelberg: Springer) Ch 5 p 149 (2015). [4] C J Hsu AICHE J. 17 732 (1971). [5] E Papoutsakis and D Ramkrishna and H C Lim Appl. Sci. Res. 36 13 (1980). [6] E Papoutsakis, D Ramkrishna and H C Lim AICHE J. 26 779 (1980). [7] E Papoutsakis and D Ramkrishna Trans. ASME J. Heat Transfer 103 429 (1981).

[8] B Weigand, M Kanzamar and H Beer Int. J. Heat and Mass Transfer 44 3941 (2001). [9] B Weigand and F Wrona Heat Mass Transfer 39 313 (2003). [10] B Weigand and D Lauffer Int. J. Heat Mass Transfer 47 5303 (2004). [11] J Lahjomri, K Zniber, A Oubarra and A Alemany Energy Conversion and Management 44 11 (2003). [12] F T Akyildiz and H Bellout Int. J. Heat Mass Transfer 47 2747 (2004). [13] G Tunc and Y Bayazitoglu Int. J. Heat Mass Transfer 44 2395 (2001). [14] H E Jeong and J T Jeong Int. J. Heat Mass Transfer 49 2151 (2006). [15] H E Jeong and J T Jeong J. Mech. Sci. Technol. 20 158 (2006). [16] B Cetin, A G Yazicioglu and S Kakac Int. J. Therm. Sci. 48 1673 (2009). [17] A K Satapathy Int. J. Therm. Sci. 49 153 (2010). [18] B Cetin, and O Bayer Thermal Science 15 103 (2011). [19] B C¸etin and S Zeinali J. Eng. Math. 89 13 (2014). [20] G Tunc and Y Bayazitoglu Int. J. Heat Mass Transfer 45 765 (2002). [21] O Aydın and M Avcı Int. J. Heat Mass Transfer 49 1723 (2006). [22] O Aydın and M Avcı Int. J. Therm. Sci. 46 30 (2007). [23] W Sun, S Kakac and A G Yazicioglu Int. J. Therm. Sci. 46 1084 (2007). [24] B Cetin, A G Yazicioglu and S Kakac Int. Commun. Heat Mass Transfer 35 535 (2008). [25] J V Rij, T Ameel and T Harmann Int. J. Therm. Sci. 48 271 (2009). [26] A Aziz and N Niedbalski Int. J. Therm. Sci. 50 332 (2011). [27] Y Kabar, R Bessaih and M Rebay Superlattices and Microstructures 60 370 (2013). [28] S Colin Trans. ASME J. Heat Transfer 134 020908 (2012). [29] H M Kushwaha and S K Sahu Sa¯dhana¯ 41 653 (2016). [30] H M Kushwaha and S K Sahu Fluid Mechanics and Fluid Power–Contemporary Research (Springer India) 1331 (2017). [31] H M Kushwaha and S K Sahu Journal of The Institution of Engineers (India): Series C 98 553 (2017). [32] Y Haddout and J Lahjomri Int. J. Heat Mass Transfer 80 673 (2015). [33] J Lahjomri and A Oubarra Trans. ASME J. Heat Transfer 121 1078 (1999). [34] G L Morini, Y Yang, H Chalabi and M Lorenzini Exp. Therm. Fluid Sci. 35 849 (2011). [35] Z Duan J. Fluids Eng. 133 074501.1 (2011). [36] J Jang and S T Wereley Microfluid. Nanofluid. 1 41 (2004). [37] L Jiang, Y Wang, M Wong and Y Zohar J. Micromech. Microeng. 9 422 (1999). [38] R K Shah and A L London Laminar flow forced convection in ducts (New York: Academic Press) Ch 17 p 394 (1978). [39] M Abramowitz and I A Stegun Handbook of Mathematical functions (New York: Dover) (1972). [40] W H Press, S A Teukolsky, W T Vetterling and B P Flannery Numerical Recipes in Fortran 77 (Cambridge University Press) (2nd ed.) (1992). [41] D K Hennecke Heat Mass Transfer 1 177 (1968). [42] B C¸etin Trans. ASME J. Heat Transfer 135 101007 (2013).