Second law analysis of forced convection in a circular duct for non-Newtonian fluids

Second law analysis of forced convection in a circular duct for non-Newtonian fluids

ARTICLE IN PRESS Energy 31 (2006) 2226–2244 www.elsevier.com/locate/energy Second law analysis of forced convection in a circular duct for non-Newto...

337KB Sizes 6 Downloads 83 Views

ARTICLE IN PRESS

Energy 31 (2006) 2226–2244 www.elsevier.com/locate/energy

Second law analysis of forced convection in a circular duct for non-Newtonian fluids Shohel Mahmud, Roydon Andrew Fraser Department of Mechanical Engineering, University of Waterloo, 200 University Avenue West, Waterloo, Ont., Canada N2L3G1 Received 12 May 2003

Abstract The second law characteristics of fluid flow and heat transfer inside a circular duct under fully developed forced convection for non-Newtonian fluids are presented. Heat flux is kept constant at the duct wall. Analytical expressions for dimensionless entropy generation number (N S ), irreversibility distribution ratio (F), and Bejan number (Be) are obtained as functions of dimensionless radius (R), Peclet number (Pe), modified Eckert number (Ec), Prandtl number (Pr), dimensionless temperature difference (O), and fluid index (m or n). Spatial distributions of local and average entropy generation number, irreversibility ratio, and Bejan number are presented graphically. For a particular value of fluid index, n ¼ 1 (or m ¼ 2), the general entropy generation number expression for a non-Newtonian power-law fluid reduces to the expression for Newtonian fluid as expected. Furthermore, entropy generation minimization is applied to calculate an optimum fluid index (nEGM ). A correlation is proposed that calculates nEGM as a function of group parameter (Ec  Pr=O) and Peclet number (Pe) within 75% accuracy. Finally, for some selected fluid indices, the governing equations are solved numerically in order to obtain Nusselt number. It is observed that the numerically obtained asymptotic Nusselt number shows excellent agreement with the analytically obtained Nusselt number. r 2005 Published by Elsevier Ltd.

1. Introduction Entropy generation is associated with thermodynamic irreversibility, which is present in all types of heat transfer processes. Different sources of irreversibility are responsible for heat transfer’s generation of entropy, such as heat transfer across a finite temperature gradient, free expansion convective flows, and viscous friction for a simple compressible substance. Bejan [1,2] has focused on the different reasons behind entropy generation in applied thermal engineering. Generation of entropy destroys the available work of a system. Therefore, it makes good engineering sense to focus on the irreversibility [2] of heat transfer and fluid flow processes and to try to understand the function of entropy generation mechanisms. A second law (of thermodynamics) analysis [3] focusing on entropy generation and its minimization has recently been playing a dominant role in understanding the irreversibility in applied engineering and transport processes.

Corresponding author. Tel.: +1 519 888 4567x4764(W); fax: +1 519 888 6197.

E-mail address: [email protected] (R.A. Fraser). 0360-5442/$ - see front matter r 2005 Published by Elsevier Ltd. doi:10.1016/j.energy.2005.09.003

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

Nomenclature A Be C1 Cp D Ec GFFI GHTI k K L L0 m n NC NF NH NR NS p Pe Pr q r R Re SG T Tr T0 u u0 v um z Z

Cross-sectional area of the pipe (m2) Bejan number (see Eq. (40)) Integration constant (see Eq. (7)) Specific heat (J kg1 1C1) Diameter of the pipe (m) nþ1 Modified Eckert number ¼ um =ðrn1 0 C p DTÞ Global fluid friction irreversibility Global heat transfer irreversibility Thermal conductivity of the fluid (W m–1 1C–1) A constant (see Eq. (11)) Length of the pipe (m) An arbitrary axial length (m) Modified fluid index ¼ 1+1/n Fluid index Entropy generation number for axial conduction Entropy generation number due to fluid friction Local heat transfer irreversibility ¼ NC+NR Entropy generation number for radial heat transfer Total entropy generation number Pressure (N m–2) Peclet number ¼ Re  Pr Prandtl number ¼ m CP/k Constant heat flux at the wall (W m–2) Radial distance (m) Dimensionless radial distance ¼ r/ro Reynolds number ¼ rum D/m Entropy generation rate (W m–3 K–1) Temperature (1C) Reference temperature (1C) Reference temperature (1C) Axial velocity component (m s–1) Velocity at inlet of the pipe (m s–1) Radial velocity component (m s–1) Maximum velocity (m s–1) Axial distance (m) Dimensionless axial distance ¼ za/(r2oum)

Greek symbol a G m r t F Y O C

Thermal diffusivity of the fluid (m2 s–1) A constant (see Eq. (26)) Dynamic viscosity of the fluid (N s m–2) Density of the fluid (kg m–3) Shear stress (N m–2) Irreversibility ratio ¼ NF/(NR+NC) Dimensionless temperature ¼ (T–Tr)/DT Dimensionless temperature difference ¼ DT/T0 Convenient representation of group parameter ¼ Ec Pr/O

2227

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2228

Fluid flow and heat transfer inside a circular duct for different boundary conditions are a fundamental area of research in engineering [4,5]. Analyses of simple systems are often useful to understand important features of complex pattern forming processes in various fields of science and technology. Circular ducts appear in many engineering applications as a single unit or in combination. Bejan [1] presented a simplified analytical expression for entropy generation rate in a circular duct with constant heat flux at the wall for a Newtonian fluid. This analysis is then extended by calculating the optimum Reynolds number [2] as a function of Prandtl number and duty parameter. Further extension is done by Sahin [6] who introduced a second law analysis of viscous fluids in a circular duct under an isothermal boundary condition. In a more recent paper, Sahin [7] presented the effect of variable viscosity on entropy generation rate for a constant heat flux boundary condition for a circular duct. For a non-circular duct, Narusawa [8] gives a theoretical and numerical analysis of the second law (of thermodynamics) for flow and heat transfer inside a rectangular duct. Sahin [9] compared the entropy generation rate inside ducts with different geometric shapes (circular, triangular, square, etc.) and calculated the optimum duct shape subject to an isothermal boundary condition restricted to laminar flow. Nag and Kumar [10] present second law optimization techniques for convective heat transfer through a duct at a constant heat flux boundary condition. For other flow configurations, like periodic flow, Rocha and Bejan [11] present second law analyses of flow and heat transfer through a single tube and bundle of tubes. For other geometries, second law analyses as well as entropy generation profiles are available in the references by Bejan [1,2] and Drost and Zaworski [12]. In each reference mentioned in the previous paragraphs a Newtonian fluid is considered. In many practical situations a non-Newtonian fluid appears. Both the modeling and analysis involve additional complexity for non-Newtonian fluids. A second law analysis extends this complexity further. In this paper, the governing equations in cylindrical coordinates are simplified and solved analytically for a circular tube with constant heat flux at the wall for a non-Newtonian power law fluid. Velocity and temperature profiles are derived from the momentum and energy equations. Entropy generation number, irreversibility ratio, and Bejan number are also determined from a second law analysis.

2. Physical model and derivation 2.1. First-law analysis Consider a circular tube of diameter D (radius ro ) and length L as shown in Fig. 1 carrying a nonNewtonian fluid. The flow is considered laminar and fully developed. We consider a situation which is axisymmetric in nature and where the swirling component of velocity (peripheral velocity component, uy ) is zero. The governing equations are: continuity equation :

1 qðrur Þ quz þ ¼ 0, r qr qz

(1)

    qur qur qp 1q qtrz r  momentum equation : r ur ðrtrr Þ þ þ uz ¼ þ , qr r qr qr qz qz

(2)

q r

ro

dz Fig. 1. Schematic view of the circular duct.

z

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244



  qp 1q qtzz ðrtrz Þ þ ¼ þ , qz r qr qz

(3)

      qT qT 1q qT q2 T þ uz r ¼k þ 2 þ mF. and energy equation : rC p ur qr qz r qr qr qz

(4)

quz quz þ uz z  momentum equation : r ur qr qz



2229

In Eqs. (1)–(4), ur and uz are the radial and axial velocity components, r, p, t, k, C p , and F are the fluid density, pressure, shear stress, fluid thermal conductivity, fluid specific heat, and viscous dissipation function, respectively. As we are considering the flow of a non-Newtonian fluid, we include the viscous stress in the diffusion terms of the radial and axial momentum equations (Eqs. (2) and (3)) instead of velocity gradient. Simplifications of Eqs. (1)–(4) are required in order to obtain analytical solutions to them. We assume a nearly unidirection flow through the pipe. At any axial location zL0 in the fully developed region, we have rD and uz u0 . Therefore, using Eq. (1), the radial velocity (ur ) in the fully developed region scales as ur Du0 =L0 . Assuming that the fully developed region is situated far away from the entrance, i.e., L0 bD, one can neglect ur (Du0/L051). For fully developed flow conditions, quz =qz must be zero, and given the resulting continuity equation one obtains qur =qr ¼ 0. The radial momentum equation reduces to qp=qr ¼ 0 which means papðrÞ; therefore, p ¼ pðzÞ. Now, by replacing qp=qz by dp=dz one obtains the simplified axial momentum equation as qðrtrz Þ dp ¼r qr dz and energy equation as   qT a q qT trz quz ¼ r uz . þ qz r qr qr rC p qr

(5)

(6)

From here, the symbol u is used instead of uz for the axial velocity component. Integrating Eq. (5) with respect to r gives the following relation:   r dp C1 . (7) t ¼ trz ¼ þ 2 dz r At r ¼ 0, the second term of the right-hand side of Eq. (7) becomes singular. But shear stress always has a finite value whatever the flow condition. Therefore, the second term of the right-hand side of Eq. (7) should be omitted to yield   r qp t¼ . (8) 2 qz According to the rheological classification of fluids [13], shear stress (t) of a non-Newtonian fluid can be expressed with a power law model by the following expression:  n qu . (9) t¼m qr The above expression covers a wide range of fluids. For pseudoplastic fluids like fine particle suspensions (no1), for dilatant fluids like ultrafine irregular particle suspensions (n41), and for Newtonian fluid like air (n ¼ 1). Replacing the shear stress (t) term from Eq. (8) with the expression for shear stress given in Eq. (9), and after rearrangement, the following expression is obtained:   qu 1 qp 1=n 1=n ¼ r . (10) qr 2m qz Integrating Eq. (10) with respect to r and applying the no-slip boundary condition at the wall (u ¼ 0 at r ¼ ro ) gives the expression for velocity as   m  r u¼K 1 . (11) ro

ARTICLE IN PRESS 2230

S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

In Eq. (11), the modified exponent m is equal to ðn þ 1Þ=n. For a constant pressure gradient along the axial direction, the coefficient K of Eq. (11) can be expressed as   1 qp n n nþ1=n r K¼ . (12) 2m qz n þ 1 o Alternatively, the velocity expression in Eq. (11) can also be expressed as   m  r u¼ 1 um ro

(13)

by using the maximum velocity um , where um is the centerline velocity which equals K for qp/qz ¼ constant. Using the definition of shear stress (t) from Eq. (9), the energy equation (Eq. (6)) can be expressed as     qT a q qT m qu nþ1 ¼ r . (14) u þ qz r qr qr rC p qr To get an expression for the temperature, Eq. (14) needs to be solved. However, it is difficult to get an analytical expression for temperature T from Eq. (14) due to its non-linear nature. To get an approximate analytical solution for T, the equation is put first into a dimensionless form. The velocity u is scaled with um , radial distance r is scaled with ro , axial distance z is scaled with r2o um =a, and the dimensionless temperature Y can be expressed as ðT2T r Þ=DT where T r is a reference temperature and DT is a reference temperature difference which is equal to qro =k for the present problem. Inserting the expression for of u from Eq. (13) into Eq. (14), followed by non-dimensional scaling, one obtains the following modified energy equation:     3n þ 1 qY 1 q qY 3n þ 1 nþ1 nþ1=n ½1  Rnþ1=n  ¼ R R , (15) þ PrEc nþ1 qZ R qR qR n where Pr is the Prandtl number, Ec is the Eckert number, and R is the dimensionless radial distance (r=ro ), respectively. The left-hand side of Eq. (15) represents convection, while the first term of the right-hand side represents thermal diffusion and the second term viscous dissipation. Neglecting the effect of viscous dissipation, Eq. (15) can be further modified into the following form:   3n þ 1 qY 1 q qY ½1  Rnþ1=n  ¼ R , (16) nþ1 qZ R qR qR subject to the following boundary conditions: at Z ¼ 0 and 0pRp1 : Y ¼ 0

(17a)

at ZX0 and R ¼ 0 : Y is finite; and

(17b)

at ZX0 and R ¼ 1 : qY=qR ¼ 1:

(17c)

In order to obtain a solution to Eq. (16) a separation of variables solution is assumed in the following form: YðR; ZÞ ¼ R1 ðRÞZ1 ðZÞ þ R1 ðRÞ þ Z 1 ðZÞ.

(18)

The first term on the right-hand side of Eq. (18) is significant for a decaying initial transient and entrance effect, the second term is significant for an axial temperature rise due to accumulated wall heat flux, and the third term is significant for radial temperature variations that drive wall heat flux into the fluid. Difficulties are imminent in the determination of the decaying initial transient. Neglecting entrance effects and assuming that the system has already passed the decaying initial transient, that is assuming initial transients eventually decay to insignificance, then only Y ¼ R1 þ Z1 is left. Inserting Y ¼ R1 þ Z 1 into Eq. (18) leaves two separate

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2231

ordinary differential equations connected by a constant l: qZ 1 ¼ l, qZ

(19)

  m 1 q qR1 R ¼ l. m þ 2 Rð1  Rm Þ qR qR

(20)

Integrating Eqs. (19) and (20) and applying the boundary conditions described in Eqs. (17a)–(17c), an expression for dimensionless temperature (Y) is obtained of the following form: Y ¼ YðR; ZÞ ¼ G þ

2ðm þ 2Þ mþ2 2 2 Zþ R  Rmþ2 . m 2m mðm þ 2Þ

(21)

In the above equation G is an integration constant. To evaluate G the mixing-cup (bulk mean) temperature [14] defined by the equation Z 1 uT dAc (22) Tm ¼ Ac uav Ac is used, where Ac and uav are the cross-sectional area of the pipe and average velocity, respectively. The quantity uav is evaluated from the following equation [14]: Z 2 ro m um : uav ¼ 2 ur dr ¼ (23) ro 0 mþ2 The dimensionless form of Eq. (22) is given by Z 2ðm þ 2Þ 1 Ym ¼ uYR dR, mum 0

(24)

where uav is replaced by um by using Eq. (23). Inserting Eqs. (13) and (21) into Eq. (24) yields Ym ¼ 4G þ

2ðm þ 2Þ m2 þ 6m þ 12 Zþ . m ðm þ 2Þðm þ 4Þ

(25)

Since Eq. (17a) requires Ym ¼ 0 at Z ¼ 0, the unknown constant G becomes G¼

1 m2 þ 6m þ 12 . 4 ðm þ 2Þðm þ 4Þ

(26)

Inserting G into Eq. (25), the asymptotic mixing-cup temperature becomes Ym ¼

2ðm þ 2Þ Z m

(27)

and the asymptotic difference between wall and mixing-cup temperatures is Yw  Ym ¼

1 m2 þ 10m þ 20 . 4 ðm þ 2Þðm þ 4Þ

(28)

In Eq. (28), Yw is the dimensionless wall temperature and is calculated from Eq. (21) by substituting 1 for R. The dimensionless temperature difference in Eq. (28) is then used to calculate an asymptotic Nusselt number (Nu1 ) as follows:  2hro 2 qY ¼ . (29) Nu1 ¼ Yw  Ym qR R¼1 k

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2232

Table 1 Limiting Nusselt numbers for different type of flow situations Type of flow

n

m

NuN, analytical (Eq. (30))

NuN, numerical (Section 3.5)

Fully developed slug flow [15] No specific name No specific name No specific name Fully developed Newtonian flow [15] No specific name No specific name No specific name Fully developed flow with linear velocity profile [14]

0 0.1 0.25 0.5 1 2 3 5 N

N 11 5 3 2 1.5 1.33 1.2 1

8.0 6.22 5.30 4.75 4.36 4.14 4.05 3.98 3.87

Not possible 6.1515 5.2800 4.7359 4.3589 4.1324 4.0498 Diverge Not possible

Use of Eqs. (21) and (28) in Eq. (29) gives the analytic asymptotic Nusselt number to be Nu1 ¼ 8

m2 þ 6m þ 8 15n2 þ 8n þ 1 ¼8 . 2 m þ 10m þ 20 31n2 þ 12n þ 1

(30)

In the heat transfer literature, three types of flow situations (for example, flow with flat velocity profile (slug flow), parabolic velocity profile (Newtonian flow), and linear velocity profile) are usually discussed in connection with fully developed heat transfer in duct flow [5,14,15]. For each case the Nusselt number’s asymptotic value is available, for example, see Arpaci and Larsen [14] and Bird et al. [15]. In this work, a single equation (Eq. (30)) is able to predict Nusselt number values for all cases without performing a case specific analysis. The above three cases are summarized in Table 1. The results in Table 1 may be used for accuracy and validation of present modeling. The dimensional temperature T can be derived from the definition of Y and Eq. (21). Therefore, finally, the equation for temperature T is "      mþ2 # qro 2ðm þ 2Þ a mþ2 r 2 2 r T¼ Gþ  (31) zþ þ T r. m r2o um 2m ro mðm þ 2Þ ro k

2.2. Second-law analysis Based on the second law of thermodynamics and assumptions already made the local volumetric rate of entropy generation, S G (W m 3 K1), in cylindrical coordinates, is given in Eq. (32). For a detailed derivation see Bejan [2] "   2 #   k qT 2 qT m qu m=m1 SG ¼ 2 þ . (32) þ qz qr T qr T The first square bracketed term on the right-hand side is due to heat transfer, and the second term is due to viscous dissipation. Differentiating Eq. (13) with respect to r and Eq. (31) with respect to r and z, one obtains the following expressions: qu mum ¼  m rm1 , qr ro

(33)

"     # qT q m þ 2 r 2 r mþ1 ¼  qr k m ro m ro

(34)

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2233

and qT 2ðm þ 2Þ qa ¼ . qz m kro um

(35)

Substituting Eqs. (33), (34), and (35) into Eq. (32) and rearranging, the volumetric entropy generation rate becomes 2(   ) (     )2  2    q2 4 4ðm þ 2Þ2 a 2 mþ2 r 2 r mþ1 m mT 0 r m SG ¼ þ  , (36) þ r o um m ro m ro m2 kT 20 kðDTÞ2 ro where we made the additional assumption that temperature variations over the pipe’s cross-section are negligible [1] compared to the absolute temperature, hence T  T 0 where T 0 is a characteristic (reference) absolute temperature. It is convenient to non-dimensionalize Eq. (36) and define the entropy generation number N S as       4ðm þ 2Þ2 1 mþ2 2 mþ1 2 2 Ec Pr m R  R R NS ¼ þ m þ , (37) Pe2 m m O m2 where Pe is the Peclet number (a product of Reynolds number (Re) and Prandtl number (Pr)), Ec is the Eckert number, and O is the dimensionless temperature difference equal to DT=T 0 . On the right-hand side of Eq. (37) the first square bracketed term represents entropy generation by heat transfer due to axial conduction (N C ), the second square bracketed term accounts for entropy generation due to heat transfer in the radial direction (N R ), and the last term is the fluid friction contribution to entropy generation (N F ). 2.3. Special case: fully developed Newtonian flow (n ¼ 1 or m ¼ 2) The entropy generation profile expressed by Eq. (37) is valid for a wide range of fluid indices mð¼ 1 þ 1=nÞ. Starting from a flat velocity profile (slug flow) at n ¼ 0 (m ¼ 1) to a linear profile at n ¼ 1 (m ¼ 1). For validation purposes one may consider the special case of a Newtonian fluid (n ¼ 1 or m ¼ 2) with parabolic flow profile. In this case Eq. (37) reduces to     16 4Ec Pr 2 3 2 NS ¼ R  þ þ ½2R  R (38) Pe2 O which is the same solution given by Bejan [1]. 2.4. Fluid friction vs. heat transfer irreversibility Entropy is generated in a process or system due to the presence of an irreversibility [1,2]. In a convection problem both fluid friction and heat transfer contribute to the rate of entropy generation. Eq. (37) generates the spatial entropy profile, but fails to give any idea whether fluid friction or heat transfer dominates. According to Bejan [1], the irreversibility distribution ratio (F) takes care the above problem and is, by definition, equal to the ratio of entropy generation due to fluid friction (N F ) to that of heat transfer (N C þ N R ). For the present problem, the irreversibility distribution ratio can be expressed as follows: F¼

½m4 E Pr Rm =O . ½4ðm þ 2Þ2 =Pe2  þ ½ðm þ 2ÞR  2Rmþ1 2

(39)

Heat transfer dominates for 0pFo1 and fluid friction dominates when F41. For F ¼ 1, both heat transfer and fluid friction contribute the same to entropy generation. In many engineering design and optimization problems [16], the contribution of heat transfer entropy generation (N C þ N R ) to overall entropy generation rate (N S ) is required. Therefore, as an alternative irreversibility distribution parameter, Paoletti et al. [17] define the Bejan number (Be) as the ratio of heat transfer entropy generation to total entropy generation.

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2234

Mathematically the Bejan number is Be ¼

NC þ NR 1 . ¼ 1þF NS

(40)

For the present problem, an expression for Bejan number can be derived using Eqs. (37), (39), and (40) to yield Be ¼

½4ðm þ 2Þ2 =Pe2  þ ½ðm þ 2ÞR  2Rmþ1 2 . ½4ðm þ 2Þ2 =Pe2  þ ½ðm þ 2ÞR  2Rmþ1 2 þ ½m4 E Pr Rm =O

(41)

Bejan number inherently ranges from 0 to 1. Accordingly, Be ! 1 is the limiting case where heat transfer irreversibility dominates, Be ! 0 is the opposite limit where the irreversibility is dominated by fluid friction effects, and Be ¼ 1=2 is the case where the heat transfer and fluid friction entropy generation rates are equal. 3. Results and discussion 3.1. Flow and thermal fields Although the main focus of this paper is entropy generation, graphical presentations of velocity and temperature profiles are required in order to understand the mechanisms of entropy generation inside the pipe. The dimensionless velocity (u=um ) and temperature (Y) are plotted as a function of dimensionless radial distance (R) in Figs. 2 and 3(a), respectively, for different fluid indices. The parabolic velocity profile (n ¼ 1 in Fig. 2) for Newtonian fluid flow in a pipe is a well-discussed case presented in all fluid mechanics textbooks (for example, see White [4]). The nature of pipe and channel flows for a non-Newtonian power law fluid (na1) and several other non-Newtonian models (for example, truncated power law model, modified power law model, Carreau model, etc.) are discussed in several textbooks (for example, Bird et al. [15]). A reduction in the shear layer (when compared with Newtonian fluid flow) is a characteristic feature of non-Newtonian fluids when no1. A non-Newtonian fluid having an index no1 is sometime referred as a shear-thinning fluid. At the 1

0.8

n= 0.6

0

u/um

0.25 0.5 0.75

0.4

1 3 5 0.2

10 ∞

0 0

0.2

0.4

0.6

0.8

r/ro Fig. 2. Dimensionless velocity profiles at different n.

1

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2235

1.2

1.2

1

0.8

10

∂Θ/∂R

1

3

0.6

5 0.5

Θ

0.7

0.4 0.2

0.6

n= 0 0.25 0.5 0.75

5

n=0

0.2

n= 1 3 5 10 ∞

0.9

∞ 5

0.3 pipe flow channel flow (n=1)

0

-0.2 (a)

0

0.2

0.4

0.6 r/ro

0.8

0

1

0

0.2

(b)

0.4

0.6

0.8

1

r/ro

Fig. 3. Dimensionless temperature profiles at different n.

lower extreme of the fluid index (i.e., n ¼ 0), shear stress becomes constant causing a flat velocity profile similar to that seen in inviscid flow. Several explanations for this shear-thinning behavior at no1 have been advanced. One explanation is that asymmetric particles or molecules are progressively aligned with streamlines, an alignment that responds nearly instantaneously to changes in the imposed shear; after complete alignment at high shear the apparent viscosity becomes constant. In another explanation, a shear thinning fluid has a structure when undisturbed that breaks down as shear stresses are increased. Nevertheless, a large portion the fluid shows a small velocity gradient near the centerline of the pipe when n decreases and approaches zero. As later seen in this paper such a small velocity gradient results in very small fluid friction irreversibility. Correspondingly, one observes the opposite behavior, a thickening of shear layer, when n41. One physical explanation of shear thickening is that as the fluid particles come close to each other some sort of bonds and structure form. Hence increasing shear leads to a thicker fluid. Alternatively, it has been suggested that fluid particles rubbing against one another can acquire an electric charge with the resultant electric attraction leading to an increase in apparent viscosity. Whatever the situation, a non-zero velocity gradient always contributes to a fluid friction irreversibility for n41. At the upper extreme of the fluid index (n ! 1), the velocity profile is linear. Dimensionless temperature (Y) and radial temperature gradient (qY=qR) profiles are presented in Fig. 3(a) and (b), respectively. In Fig. 3(a) the dimensionless axial distance is Z ¼ 0:1 in Fig. 3(a). Since the radial temperature gradient is independent of Z, Fig. 3(b) holds for all Z. For a particular n, the shape of the temperature profiles in Fig. 3(a) remains the same for all Z but displays a shift to increasing temperature for increasing Z. Due to the imposed boundary conditions, fluid inside the pipe is heated by thermal energy coming from the pipe wall at all n as confirmed in Fig. 3(b) by the positive temperature gradient at all n. For all fluid indices, the temperature gradient at the pipe centerline is zero due to symmetry, and at the wall is 1 due to the isoflux boundary conditions. When n ¼ 0, the velocity profile is flat (see Fig. 2) and independent of radial distance. In this case the energy equation reduces to the conduction equation in cylindrical coordinate with a constant source term. The solution to this n ¼ 0 conduction equation is a parabolic temperature profile (see Fig. 3(a)) with a linear variation in the radial temperature gradient (see Fig. 3(b)). In order to understand the increasing (or decreasing) tendency of qY=qR with increasing (or decreasing) n at any radial location inside the fluid region, as shown in Fig. 3(b), consider a control volume of axial length dz. A first law of thermodynamics energy

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2236

balance (see Bejan [18]) over this control volume yields     qT 2 q 2 3n þ 1 q ¼ ¼ . qz ro rC p uav ro rC p n þ 1 um

(42)

Eq. (42) implies that the temperature at all radial distances varies linearly with z with the axial temperature gradient being proportional to q=uav (or q=um ). As previously discussed the shear layer becomes thinner (or thicker) with decreasing (or increasing) n (see Fig. 2), which in turn causes a higher (or lower) average velocity and lower (or higher) qT=qz with decreasing (or increasing) n. The radial variation of temperature must adjust itself to respond to changes in qT=qz, thus one expects a decreasing (or increasing) qT=qr with decreasing (or increasing) n in order to ensure conservation of energy is satisfied. 3.2. Local entropy generation The expressions for N S , F, and Be (Eqs. (37), (39), and (41)) can be further simplified by neglecting the conduction term (proportional to Pe2). Except for liquid metals (Proo1) and creeping flow (low Reynolds number), the magnitude of the Peclet number (Pe) is usually large in value (negligible Pe2). In Fig. 4 the dimensionless entropy generation number (N S ) is plotted as a function of dimensionless radial distance for different fluid indices n. In all cases, entropy generation is zero along the centerline of the pipe, as are the velocity and temperature gradients due to flow symmetry. Also fluid friction irreversibility is negligible because a zero group parameter situation is being considered and the temperature gradient at the wall is fixed by the imposed isoflux boundary condition that is independent of any fluid index. Therefore, at the wall, the entropy generation rate is the same in magnitude for all values of n. As for N S , its magnitude is higher for higher values of n. Inside the fluid region the finite temperature gradient dominates entropy generation when the group parameter is negligible and it has already been seen (Fig. 3(b)) that the radial temperature gradient increases with increases in the fluid index at any radial location except at the wall and centerline of the pipe. For n ¼ 1, the maximum N S occurs inside the fluid at R ¼ ð2=3Þ1=2 where, due to the effect of wall curvature, the

Pe-1=0, Ec.Pr/Ω=0

1.2

1

NS

0.8

0.6 n=1 n=2

0.4

n=3 n=5 0.2

n=∞

0 0

0.2

0.4

0.6

0.8

1

r/ro Fig. 4. Entropy generation number (N S ) as a function of dimensionless radial distance (r=ro ) at different fluid index (n).

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2237

maximum radial temperature gradient occurs [1]. Temperature gradient for a Newtonian fluid flow through a parallel-plate channel is additionally plotted in Fig. 3(b) in order to understand the wall curvature effect. The expression for the temperature distribution in a channel flow with an isoflux boundary condition is available in advanced heat transfer textbooks (for example, see Arpaci and Larsen [14]) and is not derived here. In a parallel plate channel the theoretical radius of the channel wall is infinity, or the wall is free of any curvature. In this case, the maximum temperature gradient occurs at the parallel plate channel wall as shown in Fig. 3(b). In the presence of wall curvature the maximum temperature gradient occurs in the fluid region near the wall except for n ¼ 0 where the gradient is constant throughout the fluid. The maximum value of N S shifts towards the center of the duct with increasing n. For n ! 1 the maximum entropy generation occurs at R ¼ 3=4. The radial location (R ) where N S is maximum for large Pe can be estimated from the following equation:    1 mþ2  ln R ¼ exp . (43) m 2ðm þ 1Þ Note that the above discussion is restricted to large Peclet number and negligible group parameter (Ec Pr=O). The group parameter determines the importance of viscous effects [2] and cannot be neglected in real flow situations. Fig. 5 shows the distribution of N S as a function of radial distance for different values of group parameter ranging from 0 to 1 in the limit where the irreversibility due to axial conduction is negligible. Entropy generation is zero at the centerline of the duct for all values of group parameter due to flow symmetry (i.e., both velocity and temperature gradients are zero at the pipe centerline). In all cases, the pipe wall region acts as a strong concentrator of irreversibility. The temperature gradient at the pipe walls is fixed due to the imposed boundary condition and is independent of n’s variation. Although the velocity gradient varies with n (see Fig. 2); for the situation shown in Fig. 5 the velocity gradient at the pipe wall is fixed (n is fixed in Fig. 5). Any change in fluid index or group parameter must change the entropy generation rate at the pipe wall as well as in the fluid region. As viscous effects become important (i.e., as Ec Pr=O increases), the radial location where maximum entropy generates shifts closer to the wall due to the strong influence of the near wall fluid

5

Pe-1=0, n=1

Ec .P 1.0 r/Ω=

4 0.7 0.6

3 NS

0.5

2

0.3

0.4

0.2

0.1

0 1

0

0

0.2

0.4

0.6

0.8

1

r/ro Fig. 5. Entropy generation number (N S ) as a function of dimensionless radial distance (r=ro ) at different group parameter.

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2238

n=0

1 0.6

1

0.4

1.2 0.8

Be

1.5 2

0.6

4

3 6



0.4

Pe-1=0, Ec.Pr/Ω = 0.1 0.2

0

0

0.2

0.4

0.6

0.8

1

r/ro Fig. 6. Bejan number (Be) as a function of dimensionless radial distance (r=ro ) at different fluid index (n).

friction irreversibility. Maximum entropy generation occurs at the wall except for a group parameter equal to zero. Bejan number, a measure of heat transfer irreversibility, is plotted as a function of radial distance in Fig. 6 for a group parameter equal to 0.1 and for fluid indices (n) ranging from 0 to N. The velocity profile is flat for n ¼ 0 leaving no contribution of fluid friction to the entropy generation rate. Bejan number is constant and equal to its maximum value ( ¼ 1) at this n. For pseudoplastic fluids (0pno1), the velocity profile near the center region of the duct is flat, a flatness that reduces with increasing n. This flat velocity profile yields a near zero contribution from fluid friction to entropy generation in this region, and thus the Bejan number is a maximum ( ¼ 1) near the center region of the pipe. A near wall nonzero velocity gradient contributes to the overall entropy generation rate and Bejan number is no longer equal to unity near the pipe wall for no1. For a Newtonian fluid (n ¼ 1), the Bejan number is a maximum adjacent to the center; it then decreases up to the wall of the pipe. Dilatant fluids (n41) show an interesting Bejan number distribution. Whatever the value of n, at the center ‘qu=qy’ is equal to zero necessitating a maximum ( ¼ 1) Bejan number at the center; in Fig. 6 this maximum value at r=ro ¼ 0 is not shown for clarity reasons as there is a rapid fall in the Bejan number from r=ro ¼ 0 to very small radial distances from the centerline, this fall depends on n. The higher the value of n the lower the value of Be. For n ¼ 1, Be drops to 0. After this drop, Be increases with increasing radial distance and approaches a maximum before it decreases again as the wall is approached. Proceeding from the pipe centerline towards the pipe wall, both the velocity and temperature gradients increase. Up to a certain radial distance increases in the velocity gradient is not significant when compared to the temperature gradient. Therefore, the Bejan number increases with increasing R and shows a maximum near the peak temperature gradient. The magnitude of the velocity gradient is largest near the pipe wall which in turn causes a decrease in the Bejan number’s variation. For n ¼ 1:2 the maximum Be occurs at R ¼ 0:367, for n ¼ 1:5 the maximum Be occurs at R ¼ 0:49, for n ¼ 2 the maximum Be occurs at R ¼ 0:57, and so on. When n approaches N the maximum Be occurs at R ¼ 0:75. Irreversibility ratio (F) and Bejan number (Be) are plotted as functions of radial distance in Fig. 7 for different values of group parameter. A larger value of group parameter indicates a greater contribution from

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2239

4 Pe-1=0, n=1 1

Ec×Pr/Ω=0.0 0.1 3 0.2

0.8 0.3

0.5

2

0.6

0.7

Be

Φ

0.4

1.0 1.0

1

0.4 0.7

0.5 0.3 0.2

0.1 0

0 0

0.2

0.4

0.6

0.8

1

r/ro Fig. 7. Irreversibility distribution ratio (F, solid lines) and Bejan number (Be, dashed lines) as a function of radial distance at different group parameter.

fluid friction irreversibility. A lower magnitude of Bejan number is observed for higher values of group parameter at a given radial distance. For a particular value of group parameter, the Bejan number falls gradually along the radial distance except for Ec Pr=O ¼ 0 where Be is constant ( ¼ 1) due to the absence of fluid friction irreversibility. In Fig. 7 the single ‘dash-dot’ line corresponds to F ¼ 1. Above this line, entropy generation is dominated by fluid friction irreversibility while heat transfer irreversibility dominates entropy generation below this line. 3.3. Global entropy generation The dimensionless form of volumetric averaged entropy generation rate ([NS]av) can be evaluated using the following equation: Z 1 N S d8, (44) ½N S av ¼ 8 8 where 8 represents the volume of the channel. Inserting Eq. (37) into Eq. (44), and after some arithmetic operations, the global or averaged entropy generation rate becomes       mþ2 2 mþ2 mþ2 4 2 Ec Pr m2  þ ½N S av ¼ 2 þ2 þ 2 (45) mPe m2 4 m þ 4 ðm þ 2Þ2 O mþ2 in terms of modified fluid index m or       2 3n þ 1 2 nð3n þ 1Þ 3n þ 1 4n 2n2 Ec Pr ðn þ 2Þ2  þ þ2 þ 2 ½N S av ¼ Pe n þ 1 4n 5n þ 1 ð3n þ 1Þ2 O nð3n þ 1Þ ðn þ 1Þ2

(46)

in terms of fluid index n. The first two terms in Eqs. (45) and (46) represent the global heat transfer irreversibility (GHTI) and third term is the global fluid friction irreversibility (GFFI), respectively. In order to understand the effect of GHTI and GFFI on the global entropy generation rate, GHTI, GFFI, and [NS]av are

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2240

103

103 Pe=1,

Pe=0.5, EcPrΩ-1=1

[Ns]av GHTI GFFI

[Ns]av GHTI GFFI

102 GHTI, GFFI, [NS]av

102 GHTI, GFFI, [NS]av

EcPrΩ-1=1

101

100

101

100

10-1 10-3

10-2

10-1

100

(a)

101

10-1 10-3

102

n

10-2

(b)

10-1

100

101

102

n

Fig. 8. (a). Average entropy generation as a function of fluid index. (b). Average entropy generation as a function of fluid index.

103 Pseudoplastic

Dilatant

Ec.Pr/Ω= 1.0 Newtonian

0.8 102

0.6

[NS]av

0.4 0.3 0.2 0.1 0.0

Pe=1.0 101

10-3

10-2

10-1

100

101

102

n Fig. 9. Average entropy generation number at different group parameters.

plotted as a function of fluid index in Figs. 8(a) and (b) for two different Peclet numbers and a constant group parameter. It has already been shown that the local heat transfer irreversibility (N C þ N R ¼ N H ) increases with increases in the axial and radial temperature gradients and fluid index at any radial location inside the fluid. The GHTI, a spatial average of N H , is a function of fluid index and Peclet number. Note that a higher

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2241

Peclet number signifies a reduction in axial conduction as revealed by the conduction irreversibility (N C ). Therefore, at a constant Peclet number, GHTI increases with increases in fluid index as observed in Figs. 8(a) and (b). A high near wall velocity gradient contributes a major portion to fluid friction entropy generation. As previously observed that the near wall velocity gradient is higher in magnitude when n is lower. Therefore, GFFI decreases with increases in fluid index as observed in Figs. 8(a) and (b). These increasing and decreasing trends of GHTI and GFFI with increasing n are responsible for the [NS]avn profile which leads one to consider entropy generation minimization as discussed later. Fig. 9 shows the variation of [NS]av with exponent n for selected values of group parameter and constant Peclet number ( ¼ 1.0). Depending on the magnitude of n, regions for pseudoplastic, dilatant, and Newtonian fluids are identified in the figure. For a particular value of n, the magnitude of [NS]av is larger for larger values of group parameter as expected. For pseudoplastic fluids [NS]av is larger in magnitude at lower values of n. [NS]av, falls rapidly with increasing n, shows a minimum, then increases but at lower rate than the fall. The magnitude of n where the minimum value of [NS]av occurs is larger for larger group parameter values. For two consecutive average entropy generation profiles the magnitude of the difference between [NS]av’s decreases as n ¼ 1 (Newtonian fluid) is approached. For dilatant fluids, a very small increase in [NS]av occurs for 1pnp10. For larger values of n, [NS]av profiles show asymptotic behavior. The effect of Peclet number on [NS]av is shown in Fig. 10. Usually a larger value of Pe results in a negligible incremental contribution to the entropy generation rate. For a circular pipe with constant heat flux boundary condition Bejan [1] neglects the effect of Pe for PeX4. Similarly, for the present study, in most of the results, larger Pe is taken as negligible in its effect. However, for this particular case, small values of Pe are selected to check its effect on average entropy generation rate. The magnitude of [NS]av is higher for lower values of Pe due to the axial conduction irreversibility. For lower n, profiles of [NS]av merge with each other in the pseudoplastic fluids region. In the dilatant fluid region [NS]av profiles show asymptotic behavior for nX10. The difference between the magnitude of [NS]av between two consecutive profiles increases with decreases in Pe.

104 Ec.Pr/Ω=0.1

Pe=0.1

103 0.2

0.3

0.4

102 [NS]av

0.7

0.5 1.0 2.0

101



100 Pseudoplastic Dialatant 10-1 10-1

10-2

10-1

100

101

102

n Fig. 10. Average entropy generation number at different Peclet numbers.

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2242

101

100

Pipe Flow

4.0

2.0

1.0

10-1

nEGM

0.5

10-2

0.1 0.05 0.01

10-3 0.005

10-4

10-5 -3 10

.001

Pe=0

10-2

10-1

100

Ψ Fig. 11. Optimum fluid index (nEGM ) as a function of group parameter at different Peclet number.

3.4. Entropy generation minimization For a particular combination of Peclet number and group parameter, the [NS]avn profile shows a minimum (see Figs. (8)–(10)). The corresponding value of the fluid index (n) where [NS]av is a minimum can be termed the optimum fluid index (nEGM ). To determine nEGM , Eq. (46) is differentiated with respect to n and set equal to zero. After simplification, an equation for nEGM is obtained in the following form: ð125CPe2  17Pe2  5400Þn7EGM þ ð575CPe2  59Pe2  7560Þn6EGM þ ð1065CPe2  76Pe2  4176Þn5EGM þ ð1011CPe2  44Pe2  1136Þn4EGM þ ð519CPe2  11Pe2  152Þn3EGM þ ð141CPe2  Pe2  8Þn2EGM þ 19CPe2 nEGM  CPe2 ¼ 0.

ð47Þ

The left-hand side of Eq. (47) is a seventh-order polynomial of nEGM where the variable C, for convenience, represents the group parameter (Ec Pr/O). Commercial software MAPLE7 [19] is used to solve Eq. (47). For a wide range of group parameters and Peclet numbers, the distribution of nEGM is shown in Fig. 11. The linear nature of the nEGM 2Ec Pr=O profile on a log–log plot facilitates the proposal of a correlation of the following form:   Ec  Pr 0:5 nEGM ¼ FðPeÞ  , O   Ec  Pr p1 and 0:001pPeo4:0 . ð48Þ FðPeÞ ¼ 0:3945  Pe1:0249 103 p O For a given value of Peclet number and group parameter the correlation of Eq. (48) can be used to calculate the value of nEGM to within a bias of 75%. Eq. (48) is a very important result obtained from the present analysis. The group parameter and Peclet number in Eq. (48) are functions of fluid properties (for example, viscosity, thermal conductivity, etc.), flow properties (for example, fluid velocity, pressure gradient, etc.), and

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2243

Table 2 Power law parameters for aqueous solutions [15] Fluid

Temperature (K)

m (Pa sn)

n (dimensionless)

20% hydroxyethylcellulose

293 313 333

93.5 59.7 38.5

0.189 0.223 0.254

0.5% hydroxyethylcellulose

293 313 333

0.84 0.30 0.136

0.509 0.595 0.645

1.0% polyethylene oxide

293 313 333

0.994 0.706 0.486

0.532 0.544 0.599

geometric parameters (for example, length and diameter of the channel). For a given fluid, the fluid index (n) is fixed. See Table 2 for the power law parameters (viscosity and fluid index) of some non-Newtonian fluids [15]. Once the fluid is selected, one can vary the flow properties and geometric parameters in order to generate different combinations of group parameter and Peclet number. For a given fluid and a particular combination of group parameter and Peclet number which yields nEGM corresponds to fluid index for minimum irreversibility. 3.5. Numerical calculation and validation In order to validate some of our results Eqs. (1)–(4) are numerically solved using finite element method. Details regarding the implementation of the finite element approximation to solve partial differential equations are available in standard textbooks [20]. The entire solution methodology is developed using the commercially available software FEMLAB2.3 [21], which can be run either as a programmable toolbox for the development of finite element solutions using MATLAB6 [22], or as a simple graphical user interface (GUI)based integrated environment for the solution of partial differential equations using the finite element technique. The Chemical Engineering Toolbox of FEMLAB2.3 has the option to solve non-Newtonian power law fluid flow and heat transfer problems. The computational domain is initially discretized into a triangular mesh using Lagrangian elements of second order (quadratic elements). The density of these elements is higher near the wall and gradually decreases towards the pipe centerline with an expansion factor of 1.1. A comparatively finer mesh is used near the wall when the fluid index is smaller. All the calculations are performed on a Pentium IV (2.8 GHz) personal computer. For seven selected fluid indices (n ¼ 0:1, 0.25, 0.5, 1.0, 2, 3, and 5), the asymptotic Nusselt numbers (Nu1 ) are calculated from the converged solutions using MATLAB6 and reported in Table 1 for comparison with the analytically obtained results. The current version of FEMLAB is incapable of handling the zero and infinite fluid indices. Furthermore, poor convergence rate of the FEMLAB2.3 solver prevented simulation at n45. Except n ¼ 0:1, the difference between NuN,analytical and NuN,numerical is less than 1%. 4. Conclusions The second law of thermodynamics is applied to forced convection in a circular duct with constant heat flux at the wall for non-Newtonian fluids in order to obtain general expressions for entropy generation number, irreversibility distribution ratio, and Bejan number. Each derived expression has three main parts. The first part is related to axial conduction which is inversely proportional to Peclet number. The second part is related to radial heat transfer which is proportional to the radial distance. The third part is related to fluid friction which is proportional to radial distance, and a group parameter defined by the product of Eckert number, Prandtl number, and the inverse of a dimensionless temperature difference. In many practical situations the axial conduction term has negligible effect on the total entropy generation rate. Group parameter (Ec Pr O–1)

ARTICLE IN PRESS S. Mahmud, R.A. Fraser / Energy 31 (2006) 2226–2244

2244

affects significantly the entropy generation rate. For a particular fluid index n the entropy generation rate is larger for larger values of group parameter at any radial distance except at the center of the duct. For a fixed group parameter the entropy generation rate is larger for larger fluid indices n at any radial position. Differences between the rate of entropy generation for two consecutive n’s approaches zero when n approaches N. Bejan number plays a significant role for characterizing heat transfer entropy production. Its value is unity at the center of the duct for all n’s. Maximum Bejan number occurs at the center of the duct and then decreases for fluids having no1. For fluids having nX1, a second maximum in Bejan number is observed inside the fluid region, which strongly depends on n. The radial position of the second-maximum in Bejan number ranges from 0pRp0:75 for 1pnp1. For pseudoplastic fluids the average entropy generation rate strongly depends on group parameter; the larger the value of group parameter, the larger the magnitude of the average entropy generation rate as long as the fluid index (n) is constant. The average entropy generation number is independent of group parameter for dilatant fluids. Finally, an optimum fluid index (nEGM ) is calculated by applying the technique or entropy generation minimization, and a proposed correlation for nEGM is provided that works well within a bias of 75%. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

Bejan A. A study of entropy generation in fundamental convective heat transfer. J Heat Transfer 1979;101:718–25. Bejan A. Entropy generation minimization. New York: CRC Press; 1996. Bejan A. Second law analysis in heat transfer. Energy 1980;5:721–32. White FM. Viscous fluid flow. New York: McGraw-Hill; 1974. Shah RK, London AL. Laminar flow forced convection in ducts. Advances in heat transfer. New York: Academic Press; 1978. Sahin AZ. Second law analysis of laminar viscous flow through a duct subjected to constant wall temperature. J Heat Transfer 1998;120:76–83. Sahin AZ. Effect of variable viscosity on the entropy generation and pumping power in a laminar fluid flow through a duct subjected to constant heat flux. Heat Mass Transfer 1999;35:499–506. Narusawa U. The second-law analysis of mixed convection in rectangular ducts. Heat Mass Transfer 2001;37:197–203. Sahin AZ. A second law comparison for optimum shape of duct subjected to constant wall temperature and laminar flow. Heat Mass Transfer 1998;33:425–30. Nag PIK, Kumar N. Second law optimization of convective heat transfer through a duct with constant heat flux. Int J Energy Res 1989;13:537–43. Rocha LAO, Bejan A. Geometric optimization of periodic flow and heat transfer in a volume cooled by parallel tubes. J Heat Transfer 2001;123:233–9. Drost MK, Zaworski JR. A review of second law analysis techniques applicable to basic thermal science research. ASME AES 1988;4:7–12. Kumar KL. Engineering fluid mechanics. New Delhi: Eurasia Publishing House; 1992. Arpaci VS, Larsen PS. Convection heat transfer. Englewood Cliffs, NJ: Prentice-Hall, Inc.; 1984. Bird RB, Stewart WE, Lightfoot EN. Transport phenomena. New York: Wiley; 2002. Bejan A, Tsatsaronis G, Moran M. Thermal design and optimization. New York: Wiley; 1996. Paoletti S, Rispoli F, Sciubba E. Calculation of exergetic losses in compact heat exchanger passages. ASME AES 1989;10:21–9. Bejan A. Convection heat trasnfer. New York: Wiley; 1984. MAPLE7. Maplesoft, 615 Kumpf Drive, Waterloo, Ont., Canada N2V 1K8. See also http://www.maplesoft.com/products/maple. Zienkiewicz OC, Taylor RL. The finite element method. New York: McGraw-Hill; 1989. FEMLAB2.3. Cosmol, Inc., 8 New England Executive Park Suite, 310 Burlington, MA 01803, USA. See also http:// www.comsol.com/products/femlab. MATLAB6. The MathWorks, 3 Apple Hill Drive Natick, MA 01760-2089, USA. See also http://www.mathworks.com/products/ matlab.