Analytical solution to heat transfer in compressible laminar flow in a flat minichannel

International Journal of Heat and Mass Transfer 127 (2018) 975–988

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Analytical solution to heat transfer in compressible laminar flow in a flat minichannel Cheng Bao a,⇑, Zeyi Jiang a, Xinxin Zhang a, John T.S. Irvine b a b

School of Energy and Environmental Engineering, University of Science & Technology Beijing, Beijing 100083, PR China School of Chemistry, University of St Andrews, St Andrews, Fife KY16 9ST, Scotland, UK

a r t i c l e

i n f o

Article history: Received 1 June 2018 Received in revised form 6 July 2018 Accepted 20 August 2018

Keywords: Minichannel Heat transfer Compressible laminar flow Analytical solution Whittaker function

a b s t r a c t Heat transfer in compressible laminar flow in mini-/micro-channels, a classical and general topic in fields of fuel cells, electronics, micro heat exchanger, etc., is revisited. Based on a two-dimensional continuum flow model, analytical solutions of the dimensionless model are achieved in closed-form symbolic algebras of Whittaker eigenfunctions, corresponding to two kinds of boundary conditions with arbitrarily prescribed wall temperature or wall heat flux. As the eigenvalues and eigenfunctions are independent on the dimensionless quantities, which influence the along-the-channel behaviors, the algorithm reveals the common features of compressible laminar thermal flows. The algorithms do not require the assumption of a linear pressure distribution, which is proved to be untenable in some cases (e.g. constant wall heat flux). The algorithms are validated well by the exact (numerical) computations in exemplary cases of both small and moderate Reynolds number, Mach number and Eckert number of air. Although expressed in a series of eigenfunctions, only several terms (sometimes one or two terms) of solutions are required for a practical computation. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction With hydraulic diameters ranging from microns to millimeters, and correspondingly outstanding heat transfer performance due to the large specific surface area, minichannels and microchannels have been being increasingly used in many different applications, e.g. gas distributors or flow field plates in fuel cells, compact heat exchangers, heat sinks for electronics cooling, temperature control in injection molding of plastic or composite materials, micro heat pipes in micro cooling system, and on-chip laboratories, etc. [1]. Fluid flows and heat transfer in mini- and microchannels have been investigated a lot [1–5]. Besides large amounts of experimental studies and numerical calculations, analytical studies were performed to understand the fundamentals, including prediction of thermo-hydraulic performance, analysis of stability, parameter estimation in heat-transfer inverse problems, and also help in reducing the numerical computational cost. Brinkman [6] early calculated the temperature distribution in a capillary accounting for the energy dissipation of viscous flow under the conditions of a uniform wall heat flux or a homogeneous wall temperature. With an analytical solution in form of a series expansion, Siegel et al. [7] analyzed the heat transfer characteristics for a laminar forced ⇑ Corresponding author. E-mail address: [email protected] (C. Bao). https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.084 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.

convection flow in a circular tube with prescribed, i.e. uniform and arbitrary longitudinal variation of wall heat flux. Barron et al. [8] extended the Graetz problem, i.e. thermally developing heat transfer in laminar flow through a circular tube, to slip-flow regime accounting for the velocity slip and the temperature jump boundary conditions, and solved the eigenfunctions by the power-series method. With the similar mathematical method, for a gaseous flow in two-dimensional micro and nano-channels, Hadjiconstantinou and Simek [9] investigated the constant wall-temperature convective heat-transfer characteristics in the slip-flow regime 0
Kn
0.2, where Kn the Knudsen number, defined as the ratio between the molecular mean free path and the channel height. They found that it is necessary to include the effects of axial heat conduction in the continuum model due to the finite Peclet numbers. Rouizi et al. [10] presented an analytical solution of conjugate heat transfer in flat mini-channel based on Fourier transform of the temperature and normal flux and then used it inversely to estimate the fluid bulk temperature distribution from the external surface heat sources and the corresponding noised temperature profiles. Among these studies, incompressible fluids were considered, also for gaseous flow because of the small pressure drop and the small Mach number of flow in mini- and microchannels. For compressible flows, Prud’homme et al. [11] early presented an approximate solution for the expansion cooling of a compressible ideal-gas fluid, with a double perturbation

976

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

Nomenclature a0 coefficient defined in Eq. (63) a1n  a8n coefficients corresponding to nth eigenvalue defined in Eq. (31), (44), (44), (55), (56), (69), (75) and (75), respectively b coefficient defined in Eq. (25) bn coefficient corresponding to nth eigenvalue defined in Eq. (28) or (29) cp constant-pressure specific heat (J kg1 K1) de hydraulic diameter of the duct (m) D channel depth (m) E(n) E(n) = 1  4n2 Ec Eckert number L channel length (m) Mj,l(z) Whittaker’s function Nu Nusselt number Re Reynolds number p pressure (Pa) Pe Peclet number Pr Prandtl number q wall heat flux (J m2) Rg gas constant (J kg1 K1) T temperature (K) u x-direction velocity (m s1) W(n) W(n) = 3/2E(n) x x coordinate, i.e. along-the-channel direction y x coordinate, i.e. transverse direction

expansion of the radius to length ratio and the relative pressure drop. Vandenberg et al. [12] analytically investigated the velocity distribution of a steady, isothermal, compressible and laminar flow in a capillary, with a perturbation expansion of compressibility of the fluid. By setting the fluid inlet temperature equal to the homogenous wall temperature, Vandenberg et al. [13] then analytically investigated the thermal effects (expansion cooling and dissipation heating) of a compressible viscous flow. In this work, the eigenfunctions of the temperature distribution were expressed in form of power series, and the polynomial coefficients were numerically obtained from a set of linear equations. Harley et al. [14] made a theoretical investigation of compressible ideal-gas flow with low Reynolds number and high Mach number in microchannels and presented a one-dimensional analytical solution of the axial distribution of the gas density for an isothermal problem. Schlichting and Gersten [15] summarized the boundary-layer theory of incompressible and compressible flow, among the algorithms only some specific cases have the analytical solutions in closed form for compressible thermal flow. More physical applications of gaseous thermal flows in mini-/micro-channels refer to the recent review work [3,4]. In this article, we revisit the compressible gaseous thermal flow in a flat minichannel. The contents of the article are organized to first introduce the hydrodynamic equations, which are subsequently reduced and normalized based on scale analysis and dimensionless quantities. Then analytical solutions of the twodimensional model are performed in closed-form symbolic algebras of eigenfunctions. Both the cases of boundary conditions of arbitrarily prescribed wall temperature and wall heat flux are included. 2. Mathematical model Although theoretically, a fully three-dimension model is needed for the gas flow in a rectangle minichannel, we use a parallel-plate

Yn (DT)0 Greek

ap cdD cLD j

kn

l q

h ht

v

n

nth eigenfunction reference temperature rise, cf. Eq. (7)

thermal expansion coefficient, ap = 1/T for ideal gas cdD = de/D cLD = L/D thermal conductivity (J K1 m1) nth eigenvalue dynamic viscosity (Pa s) density (kg m3) dimensionless temperature dimensionless stagnant temperature or total temperature dimensionless x coordinate dimensionless y coordinate

Subscripts and superscripts m average or maximum (at centerline) n parameters corresponding to nth eigenvalue and eigenfunction 1 fluid properties of uniform flow from infinity w wall  dimensionless governing variables of u, p, q 0 zero-order solution 1 first-order solution

channel flow (i.e. two-dimension model) as shown in Fig. 1 with a physical approximation by accounting for the negligible velocity and heat transfer in the direction normal to the x-y plane. For our problem, a uniform flow from infinity is first isothermally fully developed to obtain a parabola velocity profile, then within the range of x = 0 and x = L the fluid involves in heat exchange with the walls with prescribed wall temperature or wall heat flux. The general assumptions are: steady, continuum laminar flow with negligible volume forces and no external heat sources, ideal gas, Newtonian fluid with constant properties (dynamic viscosity, specific heat, and thermal conductivity), and no-slip flow at the wall. Some assumptions are further made as: (a) y-coordinate velocity (v) is safely discarded, i.e.

v = 0.

The hydrodynamic equations for conversation of mass, momentum, and energy in the two-dimensional Cartesian coordinate system are therefore

@ðquÞ ¼0 @x 0¼

qu

ð1Þ

@p l @ 2 u þ @y 3 @x@y

ð2Þ

@u @p 4l @ 2 u @2u ¼ þ þl 2 @x @x @y 3 @x2

qcp u

@T @p @2T @2T  ap Tu ¼j þ @x @x @x2 @y2

ð3Þ ! þU

ð4Þ

where the governing variables of q, T, p and u are density, static temperature, static pressure, and x-coordinate velocity of the fluid,

977

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

Fig. 1. Geometry of a flat minichannel.

respectively, ap the thermal expansion coefficient, and U the dissipation related to shear viscosity

ap ¼ 

 2  2   1 @q @u 4 @u ;U ¼ l þ l @y 3 @x q @T p

ð5Þ

where ap = 1/T accounting for the equation of the state of ideal gases, p = qRgT. (b) The gas pressure varies only along the Poiseuille flow direction, i.e. p = p(x). The safety of this assumption for mini-/ microchannel flow was validated even in the critical conditions of very compressible fluids and high Mach number [12–14]. Eq. (2) is therefore omitted to avoid its incompatibility with p = p(x), and the safety of this treatment was also validated by Vandenberg et al. [12,13]. Note that we do not ask for the linear pressure profile as that used in [12,13], and we will prove the assumption of linear pressure profile hereafter. Define the following dimensionless quantities

Re ¼

q1 u1 de de L lcp ; cdD ¼ ; cLD ¼ ; Pr ¼ ; Pe ¼ Re Pr; D l D j

Ec ¼

u21 cp ðDTÞ0

ð6Þ

where Re is the Reynolds number, de the hydraulic diameter and de = 2D for an infinite parallel-plate channel, Pr the Prandtl number, Pe the Peclet number, Ec is the Eckert number in which the reference temperature rise (DT)0 is defined as

8 R < 1 L T w ðxÞdx  T 1 ðWTÞ L 0 ðDTÞ0 ¼ RL : 2 q u1 cp D 0 qðxÞdx ðWHFÞ

ð7Þ

1

0¼

q
u



@p cdD @ 2 u
þ @n 3RecLD @ v@n



cLD cdD
@2u
@u dp 4 @2u ¼ þ þ 2 @v dv Re 3cLD @ v2 @n2

ð10Þ !

!
@h cLD cdD @ 2 h 1 @2h dp
þ Ecu q
u
¼ þ 2 þ Ec 2 2 @v @ d Pe v v c @n LD "   2 #
2
c c @u 4 @u þ 2  LD dD @n Re 3cLD @ v

ð11Þ

ð12Þ

For the geometry of a minichannel, the depth of the channel (D) is much smaller than its length (L), i.e. 1/cLD = D/L is normally in the order of 103–102. Therefore, the items of o2/ox2 and (ou/ox)2 can be omitted and the system holds the framework as the classical boundary-layer equations [15], based on the scale analysis, as below. For more accurate solutions of Eqs. (1)-(4), one should take into account the effect of the terms in series of 1/cLD that are neglected in the present analysis, i.e. f(v,n) = f0(v,n) + (1/cLD)f1(v, n) + (1/cLD)2f2(v,n) + O((1/cLD)3) [11,12], where f denotes the dimensionless velocity, pressure, density, and temperature. On the other hand, the dimensionless form of the y-direction momentum equation (Eq. (10)) can help to check the assumption of p = p (x) accounting for the magnitude of 1/RecLD.

q
u




cLD cdD @ 2 u
@u dp ¼ þ @v dv Re @n2

ð13Þ

q
u


 
@h cLD cdD @ 2 h dp c c @ u
2
¼ þ Ecu þ Ec LD dD 2 @v dv @n Pe @n Re

ð14Þ

Here, ‘WT’ and ‘WHF’ represent the two cases of boundary conditions as shown in Fig. 1, i.e. prescribed wall temperature (WT) or prescribed wall heat flux (WHF), in which Tw(x) is the wall temperature in the WT case, q(x) the wall heat flux in the WHF case. When Tw(x) = T1 or q(x) = 0, (DT)0 can be arbitrarily taken as a non-zero value, like (DT)0 = T1, to avoid the singularity in mathematics, but in this case the Eckert number does not denote the ratio of twice of the adiabatic temperature rise to the actual temperature rise [15]. Further, define the dimensionless variables

Correspondingly, the original elliptic system becomes a parabolic one, as a result, only the boundary condition at the channel inlet (x = 0) is needed in the x direction. We will mainly solve the profiles based on the simplified set of equations, cf. Eqs. (13) and (14), and make improvements with some minor corrections accounting for the axial thermal conduction in Eq. (12). The corresponding boundary conditions are as follows accounting for no-slip flow at the wall and symmetry at the centerline, the dimensionless boundary conditions are

x y q u p  p1 T  T1 v ¼ ; n ¼ ; q
¼ ; u
¼ ; p
¼ ; h¼ q1 L D u1 q1 u21 ðDTÞ0


ð0; nÞ ¼ u

ð8Þ

where q1, p1, T1 and u1 are the gas density, pressure, temperature and mean velocity at the entrance of the channel, respectively, we then get the dimensionless form of Eqs. (1)(4)


Þ
u @ðq ¼0 @v

ð9Þ

 
ðv; 0Þ 3 @u 1
v; ¼ 0; ð1  4n2 Þ; ¼ 0; u 2 @n 2 @hðv; 0Þ hð0; nÞ ¼ 0; ¼ 0; @n 8   < h v; 12 ¼ hw ðvÞ ¼ ½T w ðvÞ  T 1 =DT 0 ðWTÞ 1 : @hðv;2Þ ¼  D qðxÞ ðWHFÞ @n jðDTÞ 0

ð15Þ

978

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988 1 1 Y n ðnÞ ¼ pffiffiffi M kn ;1 ð2kn n2 Þ; Y n ðn ! 0Þ ¼ ð2kn Þ4 n 8 4

3. Analytical approximation To get analytical solutions, the main assumption is further made as (c) The mass flow is locally fully developed, i.e. the profile of the mass flux is the same as the one obtained for a fully developed incompressible flow. Vandenberg et al. [12,13] used and validated such an assumption in isothermal compressible capillary flow, Guo and Wu [16] also presented that the locally fully developed approximation is valid even if the flow Mach number is up to be moderate. In mathematical form, considering the continuity equation (cf. Eq. (9)) and the parabolic velocity profile of an incompressible parallel flow in a straight channel, that is

3 EðnÞ; EðnÞ ¼ 1  4n2 2

q
ðv; nÞu
ðv; nÞ ¼ WðnÞ ¼

Eq. (14) then becomes

3 @h cLD cdD ð1  4n2 Þ ¼ 2 @v Pe

 2 #
cLD cdD @ u
@ h dp
þ Ec u þ dv @n Re @n2

ð16Þ

"

2

ð17Þ

The two items in the right hand side of the above equation denote heat conduction and dissipation heat, therefore in physics, the excess temperature is contributed by two effects: the effect of gas/wall heat exchange without dissipation, and the effect of frictional heat with respect to the adiabatic wall [15]. In mathematics, we take a linear superposition of zero-order and first-order parts with respect to Ec, which is normally small enough to be as the perturbation variable (e.g. u1 = 10 m s1, cp = 1006 J kg1 K1 for air, (DT)0 = 10 K, Ec = 0.01)

h ¼ h0 þ Ech1

ð18Þ

The Eckert number can be further related to the Mach number of the infinity flow (Ma1) with respect to the infinity velocity (u1) and the sound velocity at T1, i.e. Ec = (c1)(Ma1)2T1/(DT)0, in which c is the specific heat ratio, c = 1.4 for diatomic gas molecules. At the low and moderate Mach number, or even at the high Mach number but meanwhile the high temperature rise (DT)0, Ec is small and the gas/solid heat exchange overwhelms the dissipation heat to dominate the gas temperature rise. 3.1. Zero-order solutions 3.1.1. Prescribed wall-temperature (WT) case The zero-order system is

8 < 3 ð1  4n2 Þ @h0 ¼ cLD cdD 2 @v Pe : BC : h ð0; nÞ ¼ 0; 0

@ 2 h0 @n2

@h0 ðv;0Þ @n

¼ 0; h0





v; 12 ¼ hw ðvÞ

ð24Þ

Table 1 lists the values of the first fifteen non-zero kn (n = 1–15). It presents a rapidly equally-spaced feature between two neighbor eigenvalues (knkn1  8), corresponding to the periodical oscillation of Mx/8,1/4(x/2), as shown in Fig. 2. Expanding the inhomogeneous item (last term) of Eq. (21) with P the basis of the eigenfunctions, i.e. f(v,n) = fn(v)Yn(n), it yields

dxn ¼ bk2n xn þ f n dv

  2c c b ¼ LD dD 3Pe

ð25Þ

of which the solution is in the form

xn ðvÞ ¼ C n ebn v þ ebn v

Z v 0

ebn s f n ðsÞds

ð26Þ

Here bn denotes the damping coefficient of xn. To get a more accurate prediction of the axial profiles (i.e. relaxation behaviors due to system damping), we determine the value of bn from the characteristic equation corresponding to the second-order differential equation (i.e. the dimensionless form of Eq. (12), which keeps the terms of (1/cLD)2 d2xn/dv2)



1

c2LD

2

bn þ

Pe

cLD cdD

2 bn þ k2n ¼ 0 3

ð27Þ

Table 1 The first 15 eigenvalues and coefficients for WT case. n

kn

knkn1

a1n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

3.363190644477972 11.339714691790149 19.336484925020809 27.335322885215088 35.334747130698553 43.334410648649573 51.334192972667623 59.334042089371408 67.333932137332994 75.333848912529149 83.3337840124531 91.3337321717271 99.3336899352377 107.3336549485964 115.3336255561839

 7.976524047312177 7.996770233230659 7.998837960194280 7.999424245483464 7.999663517951021 7.999782324018050 7.999849116703786 7.999890047961586 7.999916775196155 7.999935099923960 7.999948159274027 7.999957763510523 7.999965013358704 7.999970607587585

0.745652186583203 0.137087067881769 0.064491968114589 0.039510622244401 0.027469857769128 0.020574322246476 0.016184902033439 0.013182960291618 0.011020446519237 0.009400044885553 0.008147738226243 0.007155497289913 0.006353025717266 0.005692784230661 0.005141573998427

ð19Þ

Defining

h0 ðv; nÞ ¼ hw ðvÞ þ xðv; nÞ

ð20Þ

it yields a system with homogeneous boundary conditions

8 < 3 ð1  4n2 Þ @ x ¼ cLD cdD

@2 x @n2

 32 dhdvw ð1  4n2 Þ : BC : xð0; nÞ ¼ hw ð0Þ; @ xðv;0Þ ¼ 0; xv; 1 ¼ 0 2 @n @v

2

Pe

ð21Þ

Based on expansion in series of the eigenfunctions, cf. Appendix B, we construct the solution as

xðv; nÞ ¼

1 X

xn ðvÞY n ðnÞ

ð22Þ

n¼1

where the eigenfunctions and eigenvalues are

  kn ¼0 4 2

M kn ;1 8

ð23Þ Fig. 2. Zeros and periodically oscillation of Whittaker function Mx/8,1/4(x/2).

979

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

the zero-order system with homogeneous boundary conditions is obtained by

So it yields

0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  2 cLD @ Pe 8 2 Pe A bn ¼ þ kn  3 2 cdD cdD

ð28Þ

For a limiting case, it reduces to the characteristic constant corresponding to the first-order differential equation (cf. Eq. (25))

bn ¼

ðwhen cLD ! 1Þ

bk2n

ð29Þ

According to Sturm-Liouville theory, the coefficients fn(v) and Cn are determined by weighted integral of the inhomogeneous item of Eq. (21) and the boundary condition at the inlet (x = 0), respectively, using the orthogonality of the eigenfunctions

f n ðvÞ ¼ a1n

dhw ; C n ¼ hw ð0Þa1n dv

Z

1 2

kY n k Z

kY n k2 ¼

1 2

ð1  4n2 ÞY n ðnÞdn

ð1  4n2 ÞY 2n ðnÞdn

0

h0 ðv; nÞ ¼ hw

!

n¼1

¼ a1n bn ebn v

Z v 0

þ

hm ðvÞ ¼ 2

1 2

0

un Y n ; un ðvÞ ð33Þ

P

ð34Þ

n¼1

The local Nusselt number and the mean Nusselt number can then be obtained by

  cdD @h0 NuðvÞ ¼ hw ðvÞ  hm ðvÞ @n n¼1 2 " #1 N N X X cdD ffiffiffi ¼ p a1n kY n k2 ðun  a1n hw Þ ðun  a1n hw Þ 3 2 n¼1 n¼1   kn ð35Þ  ðkn þ 2ÞM kn þ1;1 4 8 2 Z Num ¼ 0

1

NuðvÞdv

ð36Þ

D

jðDTÞ0

qðvÞn2 þ xðv; nÞ

ð40Þ

        kn 1 kn kn 1 kn Mkn þ1;1 M kn ;1 þ þ  ¼ 0 ðn P 1Þ 8 8 4 4 4 2 2 4 2 2

k0 ¼ 0;

The function xn(v) holds the same framework as Eqs. (25) and (26). According to Sturm-Liouville theory, the coefficients fn(v) and Cn are determined by the weighted orthogonal eigenfunctions

f n ðvÞ ¼ Cn ¼

  D dq 2ba2n qðvÞ þ a3n dv jðDTÞ0

Dqð0Þ

ð37Þ

ð42Þ

ð43Þ

a3n

jðDTÞ0

where

a2n ¼

Z

1 kY n k2

1 2

Y n ðnÞdn; a3n ¼

0

1 kY n k2

Z

1 2

n2 ð1  4n2 ÞY n ðnÞdn

ð44Þ

0

Table 2 lists the first 16 eigenvalues and coefficients a2n and a3n. The final solution of zero-order dimensionless total temperature for WHT case is

h0 ðv; nÞ ¼ 

N X D qðvÞ n2  a3n Y n jðDTÞ0 n¼0 N X

þ

!

un ðvÞY n ðnÞ; un ðvÞ

n¼0

¼

  Z v D ð2ba2n þ a3n bn Þebn v ebn s qðsÞds jðDTÞ0 0

ð45Þ

P in which the first item vanishes as a3nYn(n) = n2 when N ? 1, anyway, it is still kept accounting for the truncation error as limited items should be taken for a practical numerical calculation. The average dimensionless temperature, local Nusselt number, and the mean Nusselt number can then be obtained by

hm ðvÞ ¼ 2

Z

NuðvÞ ¼ 

3.1.2. Prescribed wall heat flux (WHF) case Defining

h0 ðv; nÞ ¼ 

ð39Þ

ð41Þ

n¼1

N X a1n kY n k2 ðun  a1n hw Þ

xn ðvÞY n ðnÞ

n¼0

ð32Þ

WðnÞh0 ðv; nÞdn

¼ hw þ 3

1 X

xðv; nÞ ¼

The eigenvalues when n > 0 are determined by the boundary condition at the wall, cf. Appendix B, and there is also a trivial eigenvalue when n = 0, i.e.

in which the first item vanishes as a1nYn(n) = 1 when N ? 1, anyway, it is still kept accounting for the truncation error when limited items (i.e. N < 1) are taken for the numerical calculation. Table 1 also lists the values of the first 15 coefficients, a1n, which show a feature of positive and negative alternation and a gradual decrease in magnitude. The average dimensionless temperature weighted by the mass flow rate is

Z

Also expanding in a series of eigenfunctions, we construct

1 Y 0 ðnÞ ¼ 1; Y n ðnÞ ¼ pffiffiffi Mkn ;1 ð2kn n2 Þ ðn P 1Þ n 8 4

N X

ebn s hw ðsÞds

ð38Þ

ð31Þ

Therefore, the zero-order dimensionless total temperature for WT case is obtained by N X 1 a1n Y n

D

jðDTÞ0

ð30Þ

0

1 2

> : BC : xð0; nÞ ¼

h i þ jðDDTÞ  2cLDPecdD qðvÞ þ 32 ddqv ð1  4n2 Þn2 0 @ xðv;1Þ ðv;0Þ qð0Þn2 ; @x@n ¼ 0; @n 2 ¼ 0 @2 x @n2

Note that, different from Eq. (22), the index of n in Eq. (39) starts from n = 0. This difference is due to both the Newman-type boundary conditions at the two ends (n = 0 and 1/2) in the WHF case, which allows the existence of a nontrivial direct-current item corresponding to n = 0, i.e.

where

a1n ¼

8 > < 32 ð1  4n2 Þ @@xv ¼ cLDPecdD

Z Num ¼ 0

1

0

1 2

WðnÞh0 ðv; nÞdn ¼ 

3Db jðDTÞ0

Z v 0

qðsÞds

ð46Þ

qðvÞd  e

jðDTÞ0 h0 v; 12  hm ðvÞ

ð47Þ

NuðvÞdv

ð48Þ

980

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

Table 2 The first 16 eigenvalues and coefficients for WHF case. n

kn

knkn1

a2n

a3n

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0 8.574449891262043 16.607448955054519 24.621212125443339 32.629043392171262 40.634194492735453 48.637883309996965 56.640677408765349 64.642879897453426 72.644668596367097 80.646155277311152 88.647414016896690 96.648496014183735 104.6494378846336 112.6502665659481 120.6510023546735

 8.574449891262043 8.032999063792476 8.013763170388820 8.007831266727923 8.005151100564191 8.003688817261512 8.002794098768383 8.002202488688077 8.001788698913671 8.001486680944055 8.001258739585538 8.001081997287045 8.000941870449893 8.000828681314474 8.000735788725436

1.5 0.415571171308799 0.259767859813622 0.194925373536150 0.158440349195970 0.134716049103952 0.117905609077432 0.105294085812198 0.095438852691162 0.087498182354645 0.080945877577779 0.075435054203879 0.070727183475642 0.066652512961491 0.063086739382730 0.059936621273379

0.05 0.031699493520337 0.008889805049549 0.004085635816366 0.002327688456663 0.001496984555122 0.001040849607436 0.000764220478390 0.000584142193557 0.000460521842797 0.000372077359867 0.000306669699167 0.000256968853360 0.000218338785437 0.000187732406313 0.000163080971660

and

3.2. Profiles of pressure, density and velocity

" Going back to the momentum equation with the dimensional governing variable of u and p, it is

qu

@u dp q1 u1 cLD cdD @ 2 u ¼ þ @v dv Re @n2

ð49Þ

Considering the equation of state of ideal gas and Eq. (16), the gas velocity is related to the local temperature as



qRg T ¼ pðvÞ q u1 Rg WðnÞ )u¼ 1 Tðv; nÞ pðvÞ qu ¼ q1 u1 WðnÞ

ð50Þ

where Rg is the gas constant of the fluid, Rg = 287 J kg1 K1 for air. Substituting Eq. (50), Eq. (49) becomes

! 1 W 2T dp @T p  ¼ W 2 p2 dv @v q21 u21 Rg þ

cLD cdD d2 W

dW @T @2T þW 2 T þ2 2 dn @n dn @n

Re

!

ð51Þ Considering

1 q ðDTÞ0 þ 12 ½ðDTÞ0 2 2T 1 2T 1

ð52Þ

Further, taking an approximation of the gas temperature,

Tðv; nÞ ¼ T 1 þ ðDTÞ0 h  T 1 þ ðDTÞ0 h0

ð53Þ

substituting it into Eq. (51) and integrating in the n direction with a weight of function W(n), we got

2 !2 3 1 d 4 54 q1 u21 5 p 2q21 u21 Rg dv 35 bq ¼

 X 12cLD cdD T 1 dun cLD cdD a4n þ þ FðvÞ þ ðDTÞ0 a5n un Re dv Re n ð54Þ

where

Z

1=2

a4n ¼ 2

W 3 Y n ðnÞdn ¼

0

Z a5n ¼ 2

W 0

27 4

Z

1=2

d W dn2

3

ð1  4n2 Þ Y n ðnÞdn

ð55Þ

0

Z

2

1=2

ð57Þ for the WT case, or

FðvÞ ¼

D

j

! N 3 X dq a3n a4n  70 n¼0 dv

ð58Þ

for the WHF case, in which a5n = 0 (n > 1) has been used accounting for the orthogonality between Y0(n) = 1 and Yn(n) (n > 1). The funcP tion F(v) is denoted as the truncation error function, as a1na4P P a1na5n = 12 in the WT case, and a3na4n = 3/70 n = 54/35 and in the WHF case, when N ? 1. Eq. (54) principally gives us a nonlinear profile of the gas pressure and correspondingly a linear profile in the cases of moderate temperature rises. 3.3. First-order total temperature The first-order system for the total temperature is

T 1 1  ¼ ; bq p2 pqavg Rg pq1 bq Rg ¼1

FðvÞ ¼ ðDTÞ0

# ! ! N N X 54 X dhw cLD cdD 12 þ  þ a1n a4n  a1n a5n hw ðvÞ 35 n¼1 dv Re n¼1

1=2

ð1  4n2 ÞY n ðnÞdn

Y n ðnÞdn ¼ 36 0

ð56Þ

8 2 c c @2 h 2 @h > >
ddvp
þ cLDRecdD @@nu
> 32 ð1  4n Þ @ v1 ¼ LDPedD @n21 þ u > < ðv;0Þ BC : h1 ð0; nÞ ¼ 0; @h1@n ¼ 0; > > >   > : h v; 1 ¼ 0 ðWTÞ or @h1 ðv;12Þ ¼ 0 ðWHFÞ 1 2 @n

ð59Þ

The system with the homogeneous boundary conditions still holds the Whittaker’s feature as the same as the zero-order system, i.e. Eqs. (23) and (24) or and Eqs. (40) and (41) hold for the eigenvalues and eigenfunctions for the WT and WHF case, respectively. Substituting Eqs. (50) in (59), it yields


cLD cdD 3 @h1 cLD cdD @ 2 h1 q1 Rg WT dp ð1  4n2 Þ ¼ þ þ 2 p dv @v Pe @n2 Re 2 2 2 q 1 Rg @T 12nT þ W  @n p2

ð60Þ

Using the temperature profile (cf. Eq. (53)) and the pressure profile (cf. Eq. (54)), expanding all the items in Eq. (59) in the basis of the eigenfunctions, following the Sturm-Liouville theory for weighted orthogonality of the eigenfunctions, we can in principle solve the first-order solution. Considering the tedious mathematical deviation, we instead present a simplified algorithm introducing below.

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

3.4. Simplified algorithm

where

In the preceding context, the fully two-dimensional profile of the gas density, i.e. q = q(x,y) was considered for strict mathematical deviations. As a matter of fact, due to the large length-depth ratio of a minichannel, the gas density in the y direction develops rapidly approaching to be homogenous except for the region very close to the channel inlet. Therefore, just in this section, the gas density is assumed to be just function of x coordinate, i.e. q = qm(x). Although it is incompatible with the equation of the state, i.e. q = p(x)/RgT(x,y) for ideal gases, the weakness of this assumption is considered to be insignificant accounting for the thermal expansion coefficient, ap = 1/T, which means a relative variation of the gas density in y direction in order of 103.

a6n ¼

3.4.1. Profiles of velocity and density Still with the parabolic velocity profile, i.e. u = um(x)E(y) where um(x) denotes the centerline (maximum in the y direction) velocity, the continuity equation leads to

q
m ðvÞu
ðv; nÞ ¼ q
m ðvÞu
m ðvÞEðnÞ ¼ WðnÞ ) q
m u
m ¼

3 2

ð61Þ

The zero-order temperature solution in the preceding mathematics is not affected. Revisiting the momentum equation, multiplying with q(x) at the both sides of Eq. (13), we got

Z

1 kY n k

2

981

1=2

n2 Y n ðnÞdn

ð69Þ

0

The index n starts from n = 1 or 0 in the WT and WHF cases, respectively, to N1, the number of terms. In the WHF case, a1n = 0 (n > 1) holds due to the orthogonality between Y0(n) = 1 and Yn(n) (n > 1). Note that again, we take the same symbol sets for the WT and the WHF cases, but the values of the variables are different because of the different eigenvalues {kn} and eigenfunctions {Yn}. The above algorithm can be further simplified. Multiplying Eq.
and summing it with Eq. (14), we got (13) with u

q
u


   
@ht cLD cdD @ 2 ht cLD cdD 1 @ @u
u ¼ þ Ec 1  Pr @n @n @v Pe @n2 Re

ð70Þ

where ht denotes the dimensionless total (or stagnant) temperature

Tt ¼ T þ

u2 Tt  T1 1
2 ; ht ¼ ¼ h þ Ecu 2 2cp ðDTÞ0

ð71Þ

With the linear perturbation expansion of ht

1 2
ht ¼ ht;0 þ Echt;1 ; ht;0 ¼ h0 ; ht;1 ¼ h1 þ u 2 0

ð72Þ

ð62Þ

in which the zero-order velocity u0 = um(x)E(y), the differential system of the first-order total temperature is therefore obtained by

Via integrating the above equation in the y direction weighted by the mass flow rate, it yields

ð73Þ


mu
m Þ2 E2 ðq


m
cLD cdD @ ln q dp @2E
m ¼ q þ q
m u
m 2 dv @v Re @n

 
m
du 35 dp 70cLD cdD
m  ¼ a0 u a0 ¼ 36 dv dv 9Re

ð63Þ

Accounting for the pressure profile, which is generally a combination of power functions and exponent functions (cf. Eq. (54) and see also the Section 4), i.e. K 2 X
X dp ¼ b vj þ am e
hm v dv j¼0 j j¼0 K

ð64Þ

" # j K 3 a0 v 35 X ð1Þjþ1 j! a0 v X ð1Þk ðjÞk jk
m ðvÞ ¼ e  b e þ v u kþ1 2 36 j¼0 j ða0 Þjþ1 k¼0 ða0 Þ 

36

j¼0

aj

a0 
hj

ðe
hj v  ea0 v Þ

Compared to Eq. (67), the differential system is simplified without the pressure gradient. By using the same mathematical framework, it yields the corresponding solution as

ut;n ðvÞY n ðnÞ; ut;n ðvÞ

n¼1;or0

  Z 16cLD cdD 1 bn v v
2m ðsÞds e a8n ebn s u 1 Pr 3Re 0

ð74Þ

where

ð65Þ

a7n ¼ ¼

q
m ðvÞ ¼ 3=2u
m

N1 X

ht;1 ¼

¼ a7n ebn v 

we can get the velocity and density profiles as

K2 35 X

8 @h @2 h 3 >
2m ð1  12n2 Þ > ð1  4n2 Þ @ vt;1 ¼ cLDPecdD @n2t;1  8cLDRecdD ð1  Pr1 Þu > 2 > < 2 ðv;0Þ BC : ht;1 ð0; nÞ ¼ 98 ð1  4n2 Þ ; @h1@n ¼ 0; > > > 1   > @h v ; ð Þ : h v; 1 ¼ 0ðWTÞ or t;1 2 ¼ 0ðWHFÞ t;1 2 @n

ð66Þ

Z

1 kY n k

2

kY n k

0

Z

1 2

1=2

3 9 ð1  4n2 Þ Y n ðnÞdn; a8n 8

1=2

ð1  12n2 ÞY n ðnÞdn

ð75Þ

0

The first-order temperature is finally simplified to 3.4.2. First-order total temperature When q = qm(x), Eq. (59) becomes


3 @h1 cLD cdD @ 2 h1 dp
m ð1  4n2 Þ ð1  4n2 Þ ¼ þu 2 dv @v Pe @n2 64cLD cdD 2 2
m n þ u Re

2 1 2
ð1  4n2 Þ h1 ¼ ht;1  u 2 m

4. Applications

ð67Þ

Following the preceding mathematical framework, the firstorder temperature profile is easily obtained by expanding all the items in the basis of the eigenfunctions as

h1 ¼

N1 X

u1n ðvÞY n ðnÞ; u1n ðvÞ

n¼1;or0

¼ ebn v

Z v 0

ebn s




2 dp 128cLD cdD 2
m a1n þ
m a6n ds u u 3 ds 3Re

ð76Þ

ð68Þ

In this section, we apply the preceding mathematical framework to two specific cases, with the prescribed wall temperature or wall heat flux in power-law profile, which is without losing of generality accounting for Taylor-series expansion of continuous rational functions. Then the analytical solutions are compared with the exact (numerical) computations, in the cases of constant wall temperature or constant heat flux. In the text below, only the zero-order solutions are included for a compact expression, as the first-order solutions contribute very little in cases of small values of Ec. Finally, a calculation is employed to validate the

982

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

algorithms at a moderate Reynolds number and Mach number, in which situation the friction heat (i.e. the first-order solution) could dominate the gas temperature profile. Case 1: the excess temperature of the wall presents a powerlaw distribution, in this case, similar solutions exist in the classical boundary-layer problems.

T w ðvÞ  T 1 ¼

M X m¼0

cm vm ; hw ðvÞ ¼

M X

c m vm

m¼0

, ðDTÞ0 ; ðDTÞ0 ¼

M X

cm m þ1 m¼0 ð77Þ

Following Eq. (33) and Eq. (B3),

  X M N N X 1 X c m vm 1  a1n Y n þ un Y n ðnÞun ðvÞ ðDTÞ0 m¼0 n¼1 n¼1 " # M m X a1n X ð1Þk ðmÞk mk ð1Þmþ1 m! bn v ¼ cm v þ e m k ðDTÞ0 m¼0 ðbn Þ ðbn Þ k¼0

h0 ðv; nÞ ¼

ð78Þ Using Eqs. (35) and (36), the Nusselt number is

"

N cdD X

NuðvÞ ¼  pffiffiffi 3 2 

N X n¼1

a1n kY n k

n¼1

un 

2

M X a un  1n c m vm ðDTÞ0 m¼0 !

!#1

  M a1n X kn cm vm ðkn þ 2ÞM kn þ1;1 8 4 ðDTÞ0 m¼0 2

ð79Þ

Fig. 3. Computational results for the constant WT case: (a) analytical and exact profiles of dimensionless temperature, density and velocity along the channel, (b) analytical and exact profiles of dimensionless pressure gradient (upper) and Nusselt number (lower).

983

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

Following Eq. (54), the pressure profile is

2

u21 d 4 p 54
þ 1  p 2Rg dv q1 u21 35bq ¼

5

N

X 12cLD cdD T 1 c c a4n /n þ LD dD a5n un þ FðvÞ þ ðDTÞ0 Re Re n¼1

where

c

dD ffiffiffi NuðvÞ ¼  p 3 2

!2 3

ð80Þ

!



ð81Þ Specially, for a constant wall temperature, i.e. Tw = Const and hw = 1 (m = 0), there are

h0 ðv; nÞ ¼ 1 

N X

a1n ebn v Y n ðnÞ

n¼1

!1

  N X kn bn v a1n ðkn þ 2ÞMkn þ1;1 e 4 8 2 n¼1 ð83Þ

N M X 54 X þ a1n a4n mcm vm1 35 n¼1 m¼1 X  N M X c c  LD dD 12 þ a1n a5n cm vm /n ðvÞ Re n¼1 m¼0 " # M m1 X ð1Þk ðmÞ dun a1n X ð1Þm m! bn v kþ1 mk1 ¼ ¼ cm v þ e m1 k dv ðDTÞ0 m¼0 ðbn Þ ðbn Þ k¼0

FðvÞ ¼

N X a21n kY n k2 ebn v

2 !2 3 u21 d 4 p1 54 5
þ  p 2Rg dv q1 u21 35bq ¼

1

X 12cLD cdD T w c c  ðDTÞ0 a1n a4n bn þ LD dD a5n ebn v Re Re n¼1

ð84Þ

In the case of the moderate temperature difference, (DT)0, the second term in the right hand side of Eq. (84) is normally small, a linear pressure profile can then be obtained

 
dp 12cLD cdD ðDTÞ0

1þ ¼ dv Re 1  54 cMa2 T1 35bq

ð85Þ

1

which can be further simplified at the moderate Mach number as

 
12cLD cdD dp ðDTÞ0 ¼ 1þ dv Re T1

ð86Þ

ð82Þ

The velocity and density profile can then be calculated by Eq. (50) and (16). The centerline velocity can be obtained based on the simplified algorithm, cf. Eq. (65)

Fig. 4. Computational results for the constant WT case: profiles of (a) dimensionless mass flux and (b) dimensionless temperature, in the y direction at different locations along the channel.

Fig. 5. Computational results for the constant WT case: profiles of (a) dimensionless average temperature and (b) dimensionless centerline temperature at different numbers of items.

n¼1

984

3 a0 v 35 e  b ð1  ea0 v Þ; 2 36a0 0   12cLD cdD ðDTÞ0

1þ b0 ¼ 54 T1 Re 1  35b cMa21 q

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988


m ðvÞ ¼ u

ð87Þ

which can be used to calculate the first-order temperature, cf. Eq. (68). Fig. 3 shows the analytical and exact along-the-channel profiles of the dimensionless gas temperature at the centerline, the average density, centerline velocity, pressure gradient and the local Nusselt number in the constant WT case. The parameters are: D = 103 m, L = 0.05 m, u1 = 0.5 m s1, p1 = 105 Pa, T1 = 300 K, Tw = 400 K, (DT)0 = 100 K, which correspond to the dimensionless parameters

as cLD = 50, Pr = 0.744, Re = 64.9, Ma1 = 0.0014, Ec = 2.48  106 for air. The exact solutions were obtained by using the commercial package, COMSOL [17], and the number of items N = 5 was taken in default for the analytical computations. The algorithms, including the simplified algorithm, were validated well. As shown in Fig. 3b, the Nusselt number approaches to Nu = 7.54, which agrees very well with the theoretical value for the fully-developed thermal flow in a rectangular straight channel with an infinite width/ height ratio of the duct. As shown in Fig. 4a, the profile of mass flux holds consistent with a parabolic one, which demonstrates the validity of the main assumption, i.e. Eq. (16) in this work. At different sections along the channel, even in the region very close to the channel inlet

Fig. 6. Computational results for the constant WHF case: (a) analytical and exact profiles of dimensionless temperature, density and velocity along the channel, (b) analytical and exact profiles of dimensionless pressure gradient (upper) and Nusselt number (lower).

985

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

(v = 0.01), the analytical temperature profile in the y direction matches well with the exact one, as shown in Fig. 4b. Fig. 5 further compares the profiles of the average and centerline temperatures with different numbers of eigenfunctions. An overshooting at the entrance region occurs at very few items and gradually disappears with more items, anyway, even the computation as N = 1 can predict well the main profiles due to the fast convergence of Whittaker functions. The overshooting at the entrance is because of the inconsistence between the inlet and wall boundary conditions, T(0,y) = T1 and T(x,D/2) = Tw at the top left corner (x = 0 and y = D/2). Such an inconsistence of the boundary conditions at the top left corner will propagate in the closed-form analytical solutions, while it doesn’t influence the numerical computation (e.g. the control-volume method does not use the info at the corner of meshes). This also explains the error of pressure gradient distribution, cf. Fig. 3b, in the region close to the channel inlet. Case 2: the distribution of the wall heat flux is in power-law form, i.e.

qðvÞ ¼

M X m¼0

cm vm ; ðDTÞ0 ¼

M 2cLD X cm q1 u1 cp m¼0 m þ 1

ð88Þ

Following Eq. (45) and Eq. (B3),

" # X N M M X X D cm 2 m mþ1 n  h0 ðv; nÞ ¼  a3n Y n cm v þ 3b v jðDTÞ0 mþ1 n¼0 m¼0 m¼0 M X D ð2ba2n þ a3n bn Þ cm jðDTÞ0 n¼1 m¼0 " # k mþ1 m X ð1Þ ðmÞ ð1Þ m! k mk ð89Þ  v þ ebn v mþ1 kþ1 ðbn Þ ðbn Þ k¼0

þ

N X

un Y n ðnÞun ðvÞ ¼ 

Using Eq. (47), the Nusselt number is

NuðvÞ ¼  ¼

M X d  e 1

c vm h ðvÞ jðDTÞ0 h0 v; 2  hm ðvÞ m¼0 m m M 3Db X cm vmþ1 jðDTÞ0 m¼0 m þ 1

ð90Þ

Following Eq. (54), the pressure profile is

2 !2 3 u21 d 4 p1 54 5
þ  p 2Rg dv q1 u21 35bq

12cLD cdD T 1 þ FðvÞ Re ! M M 3bD 54 X 12cLD cdD X cm m mþ1 þ c v þ v j 35 m¼0 m Re mþ1 m¼0

¼

 ðDTÞ0

N X a4n /n

ð91Þ

n¼1

where

FðvÞ ¼

D

j

¼

! N M X 3 X dun  a3n a4n mcm vm1 /n ðvÞ ¼ 70 n¼0 dv m¼1 D

jðDTÞ0

ð2ba2n þ a3n bn Þ

"



M m1 X X ð1Þk ðmÞ ð1Þm m! bn v kþ1 mk1 cm v þ e m kþ1 ðbn Þ ðbn Þ m¼0 k¼0

# ð92Þ

Specially, for a constant wall heat flux, i.e. q(v) = q0 (c = q0 and m = 0), there is

" ! N X Dq0 2 h0 ðv; nÞ ¼  n  a3n Y n jðDTÞ0 n¼0 #  N  X 2ba2n bn v þ3bv þ þ a3n ð1  e ÞY n ðnÞ bn n¼1

Fig. 7. Computational results for the constant WHF case: profiles of (a) dimensionless mass flux and (b) dimensionless temperature, in the y direction at different locations along the channel.

ð93Þ

Fig. 8. Profiles of the dimensionless centerline temperature at different numbers of items in the constant WHF case.

986

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

NuðvÞ ¼ 

qðvÞd 3Dbq0  e

; hm ðvÞ ¼  v jðDTÞ0 jðDTÞ0 h0 v; 12  hm ðvÞ

ð94Þ

!2 3 u21 d 4 p1 54 5
þ  p 2Rg dv q1 u21 35bq   12cLD cdD T 1 3bDq0 54 12cLD cdD þ v þ ¼ 35 Re j Re D

j

q0

N X

a4n ð2ba2n þ a3n bn Þebn v



1 54 35bq

cMa21

  12cLD cdD 3bDq0 54 12cLD cdD þ v þ Re jT 1 35 Re ð96Þ

2

þ


dp ¼ dv 1 

Note that, the constant wall heat flux case shows a parabolic pressure profile, instead of a linear one. As a result, the centerline velocity is obtained by

35 35b1 ðb a0  b1 Þ  v 36a0 36a20 0   3 35 þ þ ðb a0  b1 Þ ea0 v 2 36a20 0


m ðvÞ ¼  u ð95Þ

n¼1

The pressure profile can be further simplified because of the rapid decay of the exponent items,

ð97Þ

Fig. 6 validated the analytical profiles of the dimensionless gas temperature at the centerline, the average density, centerline

Fig. 9. Comparison between the analytical and exact solutions at moderate Reynolds number and Mach number with along-the-channel profiles of (a) zero- and first-order centerline gas temperature, (b) average gas temperature and wall temperature, (c) gas centerline velocity and (d) pressure gradient.

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

987

conditions with arbitrarily prescribed wall temperature or wall heat flux. As the eigenvalues and eigenfunctions are independent on the dimensionless quantities, which influence only the along-the-channel behaviors, our algorithm reveals the common features of all kinds of compressible laminar thermal flows. (3) The algorithm is robust without assuming a linear distribution of pressure in literature. The linear pressure profile was proved to be invalid in this work in the case of constant wall heat flux. (4) Then the algorithm was applied to two specific cases with the prescribed wall temperature or wall heat flux in power-law profile. The analytical solutions, both the zeroorder and first-order parts corresponding to the perturbation variable of Ec, were exemplarily validated well via comparison with the exact (numerical) solutions, in the cases of small Ma1 or Ec with boundary conditions of constant wall temperature or constant wall heat flux, and in the case of moderate Re and Ma1 (or Ec) with boundary condition of parabolic wall heat flux. (5) In the form of a series of eigenfunctions and eigenvalues, the analytical solutions theoretically involve many terms. At this point, the analytical solutions in this work are not marginal better the solutions in literature. Anyway, it provides the ability for arbitrarily prescribed wall temperature or wall heat flux. On the other hand, validated in the cases of constant wall temperature or constant heat flux, only several terms are required except for the region very close to the channel inlet.

velocity, and pressure gradient along the channel in the constant WHF case. The parameters are employed as: D = 103 m, L = 0.05 m, u1 = 0.5 m s1, p1 = 105 Pa, q0 =  600 W m2, (DT)0 = 102.67 K, which correspond to the dimensionless parameters as cLD = 50, Pr = 0.744, Re = 64.9, Ma1 = 0.0014, Ec = 2.419  106 for air. The algorithms work well again with good consistence between the analytical, simplified and exact solutions. The result in Fig. 6b presents a linear profile of pressure gradient, i.e. a parabolic pressure profile in the constant WHF case. This denotes that the assumption of a linear pressure profile is invalid in the constant WHF case. Fig. 7 further validates the profiles of the dimensionless mass flux and temperature. Combination of results in Fig. 6a and Fig. 7b (also Fig. 3a and Fig. 4b in the WT case), in x- and y-direction (at different values of x coordinate) respectively, gives a clear picture of the 2D temperature profile. For this reason, the 2D parameter distribution, e.g. contour distributions were not provided further. Without inconsistence between the boundary conditions, the close-form expansion with only two eigenfunctions (N = 1) works very well, and even only the directcurrent items (N = 0) can predict well the main profiles, as shown in Fig. 8. Case 3: The algorithms are finally validated at moderate Re and Ma1 (i.e. Ec). The parameters are designed as D = 104 m, L = 5  103 m, u1 = 100 m s1, p1 = 1.146  105 Pa (correspondingly the outlet pressure of 105 Pa), q(v) = 1.4026  105(v20.5v) W m2, (DT)0 = 8.726 K, which correspond to the dimensionless parameters as cLD = 50, Re = 1487.5, Ma1 = 0.288, Ec = 1.139 for air. Fig. 9a shows the zero-order and first-order solutions of the dimensionless centerline temperature along the channel. Corresponding to the small value of the overall temperature rise, (DT)0, the heat transfer between gas and wall is overwhelmed by the dissipation effect due to the high gas velocity. Therefore, the first-order solution dominates the temperature profile over than the zero-order solution. Here q(v) > 0 means heat flows from the gas to the wall and vice versa. It is natural that the gas is cooled at the first half of the channel as q(v) > 0. In the region of v = 0.5–0.9, the gas is heated by the wall as q(v) < 0, however, the centerline temperature gradually decreases because of the stronger thermal effect of expansion. Fig. 9a and b further show the analytical and exact along-the-channel profiles of the average gas temperature, wall temperature, gas velocity and pressure gradient. The algorithms work well again in this case of moderate Reynolds number and Mach number. However, Because the first-order solution was obtained based on the simplified algorithm, i.e. q = qm(x) instead of q = q(x,y), such a simplification leads to relatively larger errors compared to those in case 1 and 2 at small Mach numbers.

This work was supported by the Beijing Science and Technology Project [grant number Z181100004518004] and the Fundamental Research Funds for the Central Universities [grant number FRF-GF-17-B31]. CB thanks for the China Scholarship Council (CSC) fellowship support. CB also thanks Dr. Ruifeng Dou and Prof. Xiangjun Liu for helpful discussions.

5. Conclusions

Appendix A. Whittaker function and its properties

The classical problem of compressible laminar thermal flow in a flat mini- or micro-channel was revisited in this article, based on a two-dimensional continuum flow model expressed in the form of dimensionless quantities of length/depth ratio (cLD), Reynolds number (Re), Prandtl number (Pr), Mach number (Ma1) and Eckert number (Ec). The main conclusions are summarized as follows.

The standard solutions of Whittaker’s equation, Mj,l(z) and Wj,l(z) are [18]

(1) The main assumption of a parabolic profile of the mass flux, similar to the velocity profile of an incompressible parallel flow in a straight channel, was thought to be reasonable in the exact (numerical) computations. (2) Then analytical solutions of the dimensionless model were achieved in closed-form symbolic algebras of Whittaker eigenfunctions, corresponding to two kinds of boundary

This work is helpful to understand heat transfer in compressible laminar flow in mini- or micro-channels, a general problem in fields of fuel cells, electronics, micro heat exchangers, etc. Analytical solutions including the effects of gas injection/suction will be discussed in the next.

Conflict of interest The authors declared that there is no conflict of interest. Acknowledgments

12z

Mj;l ðzÞ ¼ e

z

1þ 2

l

1 X s¼0

1

 þlj s s z ð1 þ 2lÞs s! 2

ðA1Þ

When 2l is not an integer,

Cð2lÞ Cð2lÞ  M j;l ðzÞ þ 1  Mj;l ðzÞ C 2lj C 2þlj

W j;l ðzÞ ¼ 1

ðA2Þ

where C(∙) is the gamma function. When z ? 0, there is limiting form

Mj;l ðzÞ ¼ z2þl ð1 þ OðzÞÞ; 2l–  1; 2; 3; ::: 1

ðA3Þ

988

C. Bao et al. / International Journal of Heat and Mass Transfer 127 (2018) 975–988

The differential of Mj,l(z) is



   dM j;l ðzÞ 1 1 1 z M jþ1;l ðzÞ  jz1  M j;l ðzÞ ¼ jþlþ 2 2 dz

ðA4Þ

Note that, unlike the WT case, there are 0th (n = 0) eigenvalue and eigenfunction in the WHF case, which is presented in the text and not expressed here. R Appendix C. Integral of ebx ebssmds

Appendix B. Determination of eigenfunctions and eigenvalues Making the integral Corresponding to the homogeneous part of the differential system Eq. (21) and (38)

x 0

3 @ x cLD cdD @ 2 x ¼ ð1  4n2 Þ 2 @v Pe @n2

ðB1Þ

using separation-of-variables technique, we construct

xðv; nÞ ¼

Z

F m ðxÞ ¼

1 X

xn ðvÞwn ðnÞ

(

ebs sm ds ðb–0Þ

Using integration by parts, it is easy to get

F m ðxÞ ¼ 1b ½ebx xm  mF m1 ðxÞ ðm P 1Þ; F 0 ðxÞ ¼ 1b ðebx  1Þ

ðB2Þ

ebx F m ðxÞ ¼ ð1Þmþ1

which satisfies 2

dn2

¼ k2n ð1  4n2 Þwn

ðB3Þ

ðB4Þ

Eq. (B4) gives an exact solution in the form of Whittaker function, a kind of confluent hypergeometric function [18],

i 1 h wn ðnÞ ¼ pffiffiffi C 1n M kn ;1 ð2kn n2 Þ þ C 2n Mkn ;1 ð2kn n2 Þ 8 4 8 4 n

ðB5Þ

The standard solution of Whittaker’s equation is generally expressed as a linear combination of the two linearly independent solutions Mj,l(z) and Wj,l(z), see also Appendix A. Some studies used such mathematics for incompressible thermal flow accounting for the parabolic profile of velocity [19,20]. Alternatively, we instead use a linear combination of Mj,l(z) and Mj,l(z), as they are linear independent with respect to l = 1/4, and it gives us a simpler way to determine the eigenvalues and correspondingly more compact form of eigenfunctions. Considering the limiting form of Mj,l(z), cf. Eq. (A3), there is 1

wn ðn ! 0Þ ¼ C 1n ð2kn Þ4 n þ C 2n ð2kn Þ4

ðB6Þ

As a result, the boundary condition at the centerline (n = 0), cf. Eq. (B4), gives us C1n = 0. Therefore, the eigenfunctions{Yn} are determined in the form 1 1 Y n ðnÞ ¼ pffiffiffi Mkn ;1 ð2kn n2 Þ; Y n ðn ! 0Þ ¼ ð2kn Þ4 ðn P 1Þ n 8 4

ðB7Þ

In the WT case, the boundary condition at the wall (n = 1/2) gives

  kn ¼0 4 2

M kn ;1 8

b

mþ1

ebx þ

m X ðmÞ ð1Þk kþ1k xmk b k¼0

ðC3Þ

where

ðC4Þ

1

  dwn 2 dwn ð0Þ 1 ¼ 0 ðWHFÞ ¼ 0; wn ¼ 0 ðWTÞ; or dn dn 2

3

m!

ðmÞk ¼ mðm  1Þ . . . ðm  k þ 1Þ; ðmÞ0 ¼ 1

The boundary conditions are

ðC2Þ

As a result of the above the recurrence relation,

n

d wn

ðC1Þ

ðB8Þ

which determines the eigenvalues {kn}. Similarly, in the WHF case, the eigenvalues are determined by dYn/dn|n=1/2 = 0 according to the differential property of the Whittaker function, cf. Eq. (A4)

        kn 1 kn kn 1 kn Mkn þ1;1 M kn ;1 þ ¼ 0 ðn P 1Þ þ  4 4 8 8 4 2 2 4 2 2

ðB9Þ

References [1] S.G. Kandlikar, S. Garimella, D. Li, S. Colin, M.R. King, Heat Transfer and Fluid Flow in Minichannels and Microchannels, Elsevier, 2014. [2] J. Sarkar, A critical review on convective heat transfer correlations of nanofluids, Renew. Sust. Energy Rev. 15 (6) (2011) 3271–3277. [3] T. Dixit, I. Ghosh, Review of micro- and mini-channel heat sinks and heat exchangers for single phase fluids, Renew. Sust. Energy Rev. 41 (2015) 1298– 1311. [4] A.A. Hussien, M.Z. Abdullah, M.A. Al-Nimr, Single-phase heat transfer enhancement in micro/minichannels using nanofluids: theory and applications, Appl. Energy 164 (2016) 733–755. [5] M.M. Awad, A.S. Dalkilic, S. Wongwises, A critical review on condensation heat transfer in microchannels and minichannels, J. Nanotechnol. Eng. Med. 5 (1) (2014), 010904-010904-25. [6] H.C. Brinkman, Heat effects in capillary flow-I, Appl. Sci. Res. 2 (2) (1950) 120– 124. [7] R. Siegel, E.M. Sparrow, T.M. Hallman, Steady laminar heat transfer in a circular tube with prescribed wall heat flux, Appl. Sci. Res. 7 (5) (1958) 386–392. [8] R.F. Barron, X.M. Wang, T.A. Ameel, R.O. Warrington, The Graetz problem extended to slip-flow, Int. J. Heat Mass Transfer 40 (8) (1997) 1817–1823. [9] N.G. Hadjiconstantinou, O. Simek, Constant-wall-temperature Nusselt number in micro and nano-channels, J. Heat Trans.-T Asme 124 (2) (2002) 356–364. [10] Y. Rouizi, D. Maillet, Y. Jannot, Fluid temperature distribution inside a flat mini-channel: semi-analytical wall transfer functions and estimation from temperatures of external faces, Int. J. Heat Mass Transfer 64 (2013) 331–342. [11] R.K. Prudhomme, T.W. Chapman, J.R. Bowen, Laminar compressible flow in a tube, Appl. Sci. Res. 43 (1) (1986) 67–74. [12] H.R. Vandenberg, C.A. Tenseldam, P.S. Vandergulik, Compressible laminar-flow in a capillary, J. Fluid. Mech. 246 (1993) 1–20. [13] H.R. Vandenberg, C.A. Tenseldam, P.S. Vandergulik, Thermal effects in compressible viscous-flow in a capillary, Int. J. Thermophys. 14 (4) (1993) 865–892. [14] J.C. Harley, Y.F. Huang, H.H. Bau, J.N. Zemel, Gas-flow in micro-channels, J. Fluid. Mech. 284 (1995) 257–274. [15] H. Schlichting, K. Gersten, Boundary-Layer Theory, ninth ed., Springer, Berlin, 2016. [16] Z.Y. Guo, X.B. Wu, Further study on compressibility effects on the gas flow and heat transfer in a microtube, Microscale Therm. Eng. 2 (2) (1998) 111–120. [17] COMSOL Multiphysics, User’s Guide Version 3.5a. [18] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark, NIST Handbook of Mathematical Functions, Cambridge, 2010. [19] A. Aubert, F. Candelier, C. Solliec, Semi-analytical solution for heat transfer in a water film flowing over a heated plane, J. Heat Trans.-T ASME 132 (6) (2010) 064501. [20] A.E. Quintero, M. Vera, Laminar counterflow parallel-plate heat exchangers: an exact solution including axial and transverse wall conduction effects, Int. J. Heat Mass Transfer 104 (2017) 1229–1245.