- Email: [email protected]
Forced turbulent heat convection in a square duct with non-uniform wall temperature
International Communications in Heat and Mass Transfer 38 (2011) 844–851
Contents lists available at ScienceDirect
International Communications in Heat and Mass Transfer j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / i c h m t
Forced turbulent heat convection in a square duct with non-uniform wall temperature☆ G.A. Rivas ⁎, E.C. Garcia, M. Assato Divisão de Engenharia Mecânica, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil
a r t i c l e
i n f o
Available online 20 April 2011 Keywords: Heat transfer Numerical methods Forced turbulent heat convection Non-uniform wall temperature Square ducts
a b s t r a c t In the present work, the numerical simulation to calculate the problem of the turbulent convection with nonuniform wall temperature in a square cross-section duct was adopted. To solve this problem some assumptions for the flow, such as: the condition of fully developed turbulence and incompressible flow have been assumed. The methodology of the dimensionless energy equation was used to calculate the fluid temperature field in the square cross-section in function of the non-uniform wall temperatures prescribed. Numerical simulations were done using two different turbulent models to resolve the momentum equations and two more models to resolve the energy equation. The models of turbulence k–ε Nonlinear Eddy Viscosity Model (NLEVM) and the Reynolds Stress Model (RSM) were used to determine the turbulent intensities as well as the profiles of axial and secondary mean velocities. The turbulence model RSM was simulated using a commercial software. The thermal field was determined from other two models: Simple Eddy Diffusivity (SED), based in the hypothesis of the constant turbulent Prandtl number; and Generalized Gradient Diffusion Hypothesis (GGDH). In this last model, as the turbulent heat transfer depends on the shear tensions, the anisotropy is considered. These two last equation models of the energy equation of the fluid have been implemented in FORTRAN, a code of programming. The performances of the models were evaluated by validating them based in the experimental and numerical results published in the literature. Two important parameters of great interest in engineering are presented: the friction factor and the Nusselt number. The results of this investigation allow the evaluation of the behavior of the turbulent flow and convective heat fluxes for different square cross-sectional sections throughout the direction of the main flow, which is mainly influenced by the temperature distribution in the wall. Crown Copyright © 2011 Published by Elsevier Ltd. All rights reserved.
1. Introduction Square ducts are widely used in heat transfer devices. For instance in compact heat exchangers, gas turbine cooling systems, cooling channels in combustion chambers and nuclear reactors. The forced turbulent heat convection in a square duct is one of the fundamental problems in the thermal science and fluid mechanics. Recently, Qin and Plecther [1] showed that the Prandtl's secondary flow of the second kind has a significant effect on the transport of heat and momentum as revealed by the recent Large Eddy Simulation (LES) technique. Several experimental and numerical studies have been conducted on turbulent flow though a non-circular duct, specifically, (Nikuradse [2]; Gessner and Emery [3]; Gessner and Po [4]; Melling and Whitelaw [5]; Nakayama et al. [6]; Myon and Kobayashi [7]; Assato [8]; Assato and De Lemos [9]; Home et al. [10]; Luo et al. [11]; Ergin et al. [12] and others). Similarly important works in the turbulent heat convection ☆ Communicated by W.J. Minkowycz. ⁎ Corresponding author. E-mail addresses: [email protected] (G.A. Rivas), [email protected] (E.C. Garcia).
were developed (Launder and Ying [13]; Emery et al. [14]; Hirota et al. [15]; Rokni [16]; Hongxing [17]; Yang and Hwang [18]; Park [19]; Zhang et al. [20]; Zheng et al. [21]; Su and Da Silva Neto [22]; Saidi and Sundén [23]; Rokni [24]; Valencia [25]; Sharatchandra and Rhode [26]; Campo et al. [27]; Rokni and Sundén [28]; Yang and Ebadian [29] and others). The experimental work of Melling and Whitelaw [5] shows detailed characteristics of turbulent flow in a rectangular duct where they used a laser-Doppler anemometer to report the axial development of the mean velocity, secondary mean velocity, etc. Nakayama et al. [6] show the analysis of the fully developed flow field in ducts of rectangular and trapezoidal cross-sections using a finitedifference method with the model of Launder and Ying [13]. On the other hand, Hirota et al. [15] present an experimental work on the turbulent heat transfer in a square duct; they show detailed characteristics of turbulent flow and temperature field. Likewise, Rokni [16] carried out a comparison of four different turbulence models for predicting the turbulent Reynolds stresses and three turbulent heat flux models for ducts square. In the literature on turbulence modeling it is well known that Linear Eddy Viscosity Models (LEVM) can give rise to inaccurate
0735-1933/$ – see front matter. Crown Copyright © 2011 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2011.04.005
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Nomenclature specific heat at constant pressure duct height hydraulic diameter, Dh = 4.A/Pe = 2.L.D/(L + D) pressure gradient at the z direction (longitudinal axis) factor of friction Moody and friction coefficient of Fanning, respectively kf fluid thermal conductivity L duct width Nu Nusselt number Pe perimeter Pk turbulence production term Re Reynolds number T temperature Tb internal flow bulk temperature TWm wall mean temperature T1, T2, T3 and T4 temperature distributions at duct wall (bottom, right side, top and left side) Ub internal flow bulk velocity U, V and W averaged velocity in the Cartesian plane in the direction x, y and z; respectively y+ dimensionless wall distance cp D Dh dP/dz f, Cf
Greek symbols α thermal diffusivity, α = kf/(ρ.cp) η distance normal to the wall ϕ dimensionless temperature distribution μ dynamic viscosity μt turbulent viscosity ρ density σt turbulent Prandtl number
predictions for the Reynolds normal stresses and so that it does not have the ability to predict secondary flows in non-circular ducts due to its isotropic nature (LEVM). In spite of that, they are one of the most popular models in the engineering due to its simplicity, good numerical stability and it can be applied to a wide variety of flows. Thus, the Nonlinear Eddy Viscosity Model (NLEVM) represents a progress of the classical LEVM which permits inequality of the Reynolds normal stresses, a necessary condition for calculating turbulence-driven secondary flow in non-circular ducts within the relative cost of a two-equation formulation. The Reynolds Stress Model (RSM), also called the second order or second moment closure model is very accurate in the calculation of mean flow properties and all Reynolds stresses for many simple and more complex flows including wall jets, asymmetric channel and non-circular duct flows and curved flows, however it presents disadvantages, such as, very large computing costs. For calculating the turbulent heat flux, the Simple Eddy Diffusivity (SED) and Generalized Gradient Diffusion Hypothesis (GGDH) models have been adopted and investigated. Most of the works presented in the literature show results assuming constant temperatures on the wall. However, in many engineering applications the heat flux and surface temperature are non-uniform around the duct, becoming important the knowledge of the variation of the conductance around the duct, Kays and Crawford [30]. According to developments performed by Garcia [31], it is possible to carry out analysis with non-uniform wall temperature boundary conditions. In this case, it is necessary to define a value that represents the mean wall temperatures in a given cross section, TWm. In the present work important results in square cross-section ducts considering the fluids air and water and the influence of the nonuniform wall temperature distribution on the turbulent convective thermal field, are presented.
845
2. Mathematical formulation 2.1. Governing equations The equation systems to be solved are the Reynolds Averaged Navier Stokes (RANS), continuity equation, Eq. (1), momentum equation, Eq. (2), and energy equation, Eq. (3). ∂ U =0 ∂xj j
ð1Þ
" !# ∂Uj ∂ 1 ∂P 1 ∂ ∂Ui 1 ∂ Ui Uj = − + μ + −ρui′uj′ ð2Þ + ρ ∂xi ρ ∂xj ρ ∂xj ∂xj ∂xi ∂xj " # ∂ 1 ∂ μ ∂Tf ′ ′ UT = + −ρuj t ρ ∂xj Pr ∂xj ∂xj j f
ð3Þ
It is considered fully developed turbulent flow and heat transfer and the following hypotheses have been utilized: steady state, condition of nonslip on the wall, and fluid with constant properties. The turbulent Reynolds stress −ρui′uj′ and the turbulent heat flux −ρuj′t ′ are modeled and solved by solution of algebraic or differential expressions. 2.2. Turbulence models for Reynolds stresses 2.2.1. Nonlinear Eddy Viscosity Model (NLEVM) The NLEVM model to reproduce the tensions of Reynolds presents the inclusion of nonlinear terms to the basic constitutive equation trying to capture the sensitivity of the curvature of the stream lines. This model is based on the initial proposal of Speziale [32]. The Reynolds average equations, Eqs. ((1)–(3)) are applied for the device of Fig. 1(a) and (b). The components of velocities U, V represent the secondary flow and the axial velocity component W, the velocity of the main flow. The equations of transport in tensorial form for the turbulent kinetic energy, κ, and the rate of dissipation ε, respectively, they are given by: Ui
∂k ∂ μ t ∂k = + Pk −ε ∂xi ∂xi ρσk ∂xi
ð4Þ
Ui
∂ε ∂ μ t ∂ε ε ε2 = + c1 Pk −c2 f2 : k k ∂xi ∂xi ρσε ∂xi
ð5Þ
The symbols Pk and μ t, represent the rate of the turbulent kinetic energy production and the turbulent viscosity, respectively, and thus, we have: Pk = τij
∂Ui ; ∂xj
2
μ t = cμ fμ ρ
k : ε
ð6Þ
In the present work for NLEVM, the formulations of Low Reynolds Number will be assumed for wall treatment. The damping functions f2 and fμ are observed in Eqs. (5) and (6) and shown in Table 1. These functions and the constant C1 and C2 have been used together with the equations κ–ε, the subscribed letter P refers to the nodal point near to the wall. Thus UP and kP are the values of the velocity and kinetic energy in this point, respectively. The constant cμ , c1, c2, σk and σε adopt the values of 0.09; 1.5; 1.9; 1.4 and 1.3; respectively. The new constitutive relation for the tensions of Reynolds in the model NLEVM, assumed in the thesis of Assato (2001), is given by: NL L k 1 τij = μ t Sij + c1NL μ t Sik Skj − Skl Skl δij : ð7Þ ε 3 This expression shows that the second term of the right side of Eq. (7), represents the nonlinear term added the original constitutive
846
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Fig. 1. (a) Fully developed turbulent flow in rectangular ducts, (b) Rectangular duct — system of reference and transversal section.
relation. This quadratic term represents the degree of anisotropy between the normal tensions of Reynolds, which makes it possible to predict the presence of the secondary flow in non-circular ducts. The values of c1NL proposed by Speziale [32] are equal to 1.68. Here, c1NL will be analyzed and will adopt values for the formulation of Low Reynolds Numbers. The tensions of Reynolds, normal and shear, are presented in Eqs. (8), (9) and (10), are expressed as: " # k 1 ∂W 2 2 ∂W 2 τxx = c1NL μ t − ; ε 3 ∂x 3 ∂y " # k 1 ∂W 2 2 ∂W 2 − τyy = c1NL μ t ε 3 ∂y 3 ∂x τxy = c1NL μ t
k ∂W ∂W ; ε ∂x ∂y
τxz = μ t
∂W ∂W ; τyz = μ t : ∂x ∂y
ð8Þ
∂ ∂ ρuk ui′uj′ ðbÞ ρui′uj′ ðaÞ + ∂t ∂xk ∂ ∂ ∂ ρui′uj′uk′ + p δkj uj′ + δik uj′ ðcÞ + μ u′u′ ðdÞ =− ∂xk ∂xk ∂xk i j ! ! ð12Þ ∂uj′ ∂uj ∂u′i ∂u ð gÞ −ρ ui′uk′ + uj′uk′ i ðeÞ−ρβ gi uj′θ ð f Þ + p + ∂xk ∂xk ∂xj ∂xi −2μ
ð9Þ
The following differences for the normal tensions of Reynolds are presented and have been observed for this type of flow. " # k ∂W 2 ∂W 2 − τyy −τxx = c1NL μ t ε ∂y ∂x
equations of transport. The individual Reynolds tensions are utilized to close the average Reynolds equations of the momentum. This model has shown superiority regarding the models of two equations in complex flows that involve swirl, rotation, etc. The exact equations of transport for the Reynolds tensions, ρui′uj′, can be written as follows:
Table 1 Low Reynolds Number formulation for wall treatment using the NLEVM model.
ð10Þ fu
In such a way in Eq. (6), including the derivatives above the tensions of Reynolds, the turbulence production term, is expressed as: Pk = τxz
∂W ∂W + τyz : ∂x ∂y
Low Reynolds Numbers formulation proposed by Abe, Nagano and Kondoh [33] 8 2 2 !2 3 9 = n h i o2 < 5 Þ0:25 4− k =vε 5 1− exp − ðvε14v 1+ exp 0:75 : ; 2 200 ðk =vεÞ 2 " #)2 8 2 !2 39 < = k =vε ðvεÞ0:25 5 1− exp − 1−0:3 exp4− : ; 3:1v 6:5
( f2
ð11Þ
2.2.2. Reynolds Stress Model (RSM) The most complex turbulence is the Reynolds Stress Model (RSM), also called second order, it involves calculations of the tensions of Reynolds in an individual form, ρui′uj′, used for this differential
∂ui′ ∂u′j ′ εikm + ui′um ′ εjkm ðiÞ + Sð jÞ: ðhÞ−2ρΩk uj′um ∂xk ∂xk
σk σε C1 C2 τw
1.4 1.3 1.5 1.9 ∂UP μ ∂yP
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Where the respective letters represent: (a) local derivative of the time; (b) Cij≡ convection; (c) DT,ij≡ turbulent diffusion; (d) DL,ij≡ molecular diffusion; (e) Pij≡ production term of tensions; (f) Gij≡ buoyancy production term; (g) ϕij≡ term of pressure–tension (redistribution); (h) εij≡ term of dissipation; (i) Fij≡ term production for the rotation system; and (j) Sj≡ source term. The terms of the exact equations, presented previously,Cij, DL,ij, Pij and Fij do not require modeling. However, the terms DT,ij, Gij, ϕij and εij need to be modeled to close the equations. For the present analysis, the model LRR [34] is chosen, which assumes that the correlation velocity pressure is a linear function of the anisotropy tensor LRR in the phenomenology of the redistribution, ϕij. For the treatment of the wall, it is also assumed the Low Reynolds Numbers and the periodic conditions, Rokni [16]. This model had been simulated in the commercial code Fluent 6.3. 2.3. Turbulence models for turbulent heat flux
and ∫ W:Tf :dA Tb =
μ t ∂Tf σT ∂xj
∂Tf k = −ρCt u′j u′k : ε ∂xk
2.3.3. Dimensionless Energy Equation for SED and GGDH Models For a given cross section of area A, it is possible to define a mean velocity Ub and a bulk temperature Tb, express as:
1 Ub = ∬W:dx:dy A
ð15Þ
Table 2 Fluid properties (at 300 K) for internal flow (Moran et al. [36]).
Specific heat, cp [kJ/kg.K] Specific mass ρ [kg/m3] Thermal conductivity kf [W/m.K] Thermal diffusivity α [m²/s] Dynamic viscosity μ [N.s/m²]
ð16Þ
∂Tf ∂x
+V
∂Tf ∂y
+W
∂Tf
"
∂Tf ∂ −u′ t ′ α ∂z ∂x ∂x !# −
∂Tf ∂ −w′ t ′ α ∂z ∂z
! +
∂Tf ∂ −v′ t ′ α ∂y ∂y
= 0:
!
ð18Þ
The properties of the fluid are shown in Table 2, and the following considerations are applied to obtain the variables in dimensionless form: X=
x ; Dh
ð19Þ
Y=
y Dh
ð20Þ
ð14Þ
The constant Ct, assuming the value of 0.3, is adopted. The main advantage of this model is that it considers the anisotropic behavior of the fluid heat transport in ducts (Table 2).
Fluid
1 ∬W:Tf :dx:dy: A:Ub
It is possible to develop a formula similar to Kays and Crawford [30], the new expression for the turbulent energy equation, is presented as:
+
!
=
1 L 1 D 1 L 1 D :∫0 T1 ð0; yÞ:dy+ :∫0 T2 ðx; 0Þ:dx+ :∫0 T3 ðD; yÞ:dy + :∫0 T4 ðx; LÞ:dx L D L D : = 2ðL + DÞ
ð13Þ
2.3.2. Generalized Gradient Diffusion Hypothesis (GGDH) Daly and Harlow [35] present the following formulation to the turbulent heat flux:
ρu′j t ′
Ub :A
ð17Þ
U ρu′j t ′ = −
A
Kays and Crawford [30] developed an applicable formulation to rectangular cross section ducts. They considered the boundary conditions with prescribed uniform wall temperatures at the cross section, and along the duct length. According to developments performed by Garcia [31], it is possible to carry out analysis with non-uniform wall temperature boundary conditions. In this case, it is necessary to define a value that represents the mean wall temperatures in a given cross section, TWm, given as:
TWm
2.3.1. Simple Eddy Diffusivity (SED) This method is based on the Boussinesq viscosity model. The turbulent diffusivity for the energy equation can be expressed as: αt = μt=ρσt , where the turbulent Prandtl number σt needs to be given. The SED model assumes that the turbulent Prandtl number is constant in the entire region, for the air σt it assumes values of 0.89, independent on the wall proximity effect.
847
Air
Water
1005.7
4179.0
1.1614
997.606
2.63 × 10− 2
61.3 × 10− 2
225 × 10− 7
6.15 × 10− 7
1.846 × 10− 5
88.5 × 10− 5
and α: TWm −Tf : ϕ= dTb Ub :D2h : dz
ð21Þ
Replacing Eqs. (13), (14), (19) up to (21) in Eq. (18), dimensionless energy equation for SED and GGDH becomes, respectively:
1 ∂ ∂ϕ 1 ∂ ∂ϕ D ∂ϕ ∂ϕ ðα + αt Þ ðα + αt Þ U + − h +V α ∂X α ∂Y α ∂X ∂Y ∂X ∂Y W ϕ =− ð22Þ UB ϕB ∂ϕ 1 ∂ ∂ϕ 0 0 k 0 0 k α + Ct u u + α + Ct v v ε ∂X α ∂Y ε ∂Y D ∂ϕ ∂ϕ W ϕ C ∂ 0 0 k ∂ϕ uv − t +V =− − h U UB ϕB ε ∂Y α α ∂X ∂X ∂Y C ∂ 0 0 k ∂ϕ vu − t ε ∂X α ∂Y
1 ∂ α ∂X
ð23Þ
848
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
The fluid temperature field Tf can be replaced by Tb and Eq. (21) can be expressed as: ϕb =
α:ðTWm −Tb Þ ; dTb Ub :D2h : dz
ð24Þ
and D2h :Ub dTb : :ϕ : ð25Þ α dz b From Eq. (21), Eq. (26) is obtained, and applying this in Eq. (16), one gets Eq. (27): Tb = TWm −
Tf = TWm −
2 ϕ:Ub :Dh :
α
dTb dz
ϕb =
D2h :∬W:ϕ:dX:dY: A:Ub
ð28Þ
It is possible to compute the heat transfer rate per unit length on the wall surface, q′, Eq. (29). It depends on values for TWm, Tb, and on the average heat convection coefficient, h. From fluid enthalpy derivative gradient [dhb = cp.dTb], the heat transfer rate per unit length in the fluid, q′f, can be expressed by Eq. (30). ;
q = Pe :h:ðTWm −Tb Þ;
ð29Þ
and
qf′ = ρ:Ub :A:cp : ;
ð26Þ
dTb dz
ð30Þ
When Eq. (29) is made equal to Eq. (30), integrating two cross sections (inlet, namely z1, and outlet, namely z2), it is possible to develop the resulting expression, Eq. (31).
and
Tb = TWm −
Replacing Eq. (27) in Eq. (24), and using also Eqs. (19) and (20), the dimensionless bulk temperature is obtained:
D2h dTb : :∬W:ϕ:dx:dy: α:A dz
ð27Þ
− Tbz2 = TWm − TWm −Tbz1 :e
Pe 2 α : :Nu:ðz1 −z2 Þ 2A Ub
ð31Þ
Fig. 2. (a) Grid 120 × 120 for numerical simulation. (b) Secondary flow contours and comparisons of the axial mean velocity with Melling and Whitelaw [5] for water utilizing NLEVM model.
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
From Eq. (31), the bulk temperature longitudinal (z-axis) variation ΔTb is obtained. It is done by “cutting” the duct into a lot of segments and applying the numerical method to find ΔTb at each finite cross section. For a given bulk temperature at the duct inlet section (Tb1), after solving the equation system, duct outlet bulk temperature (Tb2) is computed from that expression, Eq. (31). The dimensionless boundary conditions are given by the following equations:
ϕð0; Y Þ =
α:ðTWm −T1 Þ α:ðTWm −T2 Þ ; and ϕð X; 0Þ = dTb dTb 2 Ub :Dh : Ub :D2h : dz dz
ð32Þ
When considering uniform wall temperature, Eqs. (32) and (33) are equal to zero, and for these particular conditions, it is possible to notice that these boundary conditions are not functions of dTb/dz. That simplification becomes equal to the one studied by Patankar [37]. Eqs. (24), (25) and (26), as well as the boundary conditions from Eqs. (32) and (33), form a set of differential equations, in which ϕ and dTb/dz parameters are unknown. When that equation system is solved, it is possible to obtain Tf. 2.3.4. Additional equations Additional equations were utilized for the calculation of the factor of friction Moody, f; coefficient of friction of Fanning, Cf; Prandtl law; local Nusselt number for the Low Reynolds formulation (Rokni [16]), Nuxp and correlation of Gnielinsky, Nu, respectively; these equations are given by:
(a) 10 9
Numerical Simulation
8
f≡
Model Non Linear k-e Model RSM
7 6
Cf (x10-3)
α:ðTWm −T3 Þ D ; ϕ =Dh ; Y = dTb Ub :D2h : dz
849
L α:ðTWm −T4 Þ : ð33Þ and ϕ X; = dTb Dh Ub :D2h : dz
−ðdP = dzÞ:Dh ρUB2 =2 ;
Cf =
ð34Þ
f ; 4
ð35Þ
pffiffiffi 1 pffiffiffi = 2 log Re f −0:8; f ðT −TP Þ Nuxp = Dh w ηðTw −Tb Þ
5
4
ð36Þ ð37Þ
Methods Experimental and Correlations 3
Prandtl Friction Law Harnett et al. Leutheusser Launder e Ying Lund
(a)
1
1 0.65
0.8
0.65
0.70
0.70
2 2
3
4
5
6
7 8 9 10
20
30
0.80
0.6
Re (x104)
103
0.90 0.95
(b)
0.6
0.85
0.85 0.90
0.4
0.4
0.95
0.2
0.2
Numerical Simulation
0 -1
SED Model GGDH Model
0.8
0.75 0.80
0.75
-0.8
-0.6
-0.4
-0.2
Numerical Prediction
0
0.2
0.4
0.6
0.8
0 1
Experimental
Nu/Pr0.4
(b) 102
Methods experimental - Correlations Lowdermilk et al. Gnielinsky Brundrett & Burroughs
101
104
105
106
Re Fig. 3. (a) Friction coefficient for fully developed flow, (b) Nusselt number dependence on Reynolds number for fully developed flow.
Fig. 4. (a) Results (RSM–SED) numerical and (Hirota et al. [15]) experimental for mean temperature (uniform wall temperature). (b) Fluid temperature with non-uniform wall temperature (Case I).
850
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
obtained (dTb/dz b tolerance). This is the end of the second iterative loop;
and 2
3 ðf = 8ÞðRe−1000ÞPr 4 1 5: Nu = 1 1 + 12:7ðf =8Þ =2 Pr =2 −1
ð38Þ
3. Numerical implementation After applying the method of finite differences to the algebraic equations, to obtain the temperature fields, the following five steps indicate the developed methodology in the numerical solution (Garcia [31]): Step 1: To define the function value of the non-uniform temperatures in the walls of the duct TWm = f(T1(0,y), T2(x,0), T3(D,y), T4(x,L)), what that can be expressed by a Fourier expansion; Step 2: To obtain velocity field and estimated values for Ub, TWm and dTb/dz; Step 3: Equations for the boundary conditions are evaluated (Eqs. (32) and (33)); Step 4: Dimensionless energy equation (temperature field, Ø) at Eq. (24) is solved and Øb is computed according to Eq. (29), until convergence is obtained (Øb b tolerance). This is the end of the first iterative loop; Step 5: A value for dTb/dz is computed in accordance with Eq. (25). Boundary conditions are updated (step 3) to obtain a solution for the new temperature field (step 4), until convergence is
(a)
For all steps, “tolerance of 10− 7” is the value to be accomplished by the convergence criteria, which is applicable to Øb (dimensionless bulk temperature), dTb/dz and Ø (dimensionless temperature field). 4. Results and discussion 4.1. Fluid flow and heat transfer field Fig. 2(a) shows the utilized grid (120 × 120) in the numerical simulation for the formulations of Low Reynolds, Fig. 2(b) represents the secondary flow contours and comparisons of the velocity profile (NLEVM, Assato [8]) with the experimental work of the Melling and Whitelaw [5] for fluid water and Re = 42,000. The predicted distributions of the friction coefficient (NLEVM and RSM) and Nusselt number (SED and GGDH) dependence on Reynolds number for fully developed flow and heat transfer in a square duct are shown in Fig. 3(a) and (b), respectively. In Fig. 4(a), the comparisons of the results (RSM–SED) numerical with the experimental for temperature profile (wall constant temperature) (TWm − Tf)/(TWm − TC) with fluid air and Re = 65,000 (Hirota [15]) are shown, while Fig. 4(b) shows the variation of the temperature profile with non-uniform wall temperature: south = 400 K, north = 373 K, east = 393 K, and west = 353 K; it called of Case I. Already Fig. 5(a) shows the variation of the temperature profile with non-uniform wall temperature, this it is represented by means of functions sine (Case II), south = (350–20Sin (ζ))K, north = (400–50Sin(ζ))K, east = (330 + 20Sin(ζ))K, and west = (350 + 50Sin(ζ))K. Where ζ is function of the radians (0–π/ 2) and i,j (points number of the grid in the direction x and y, respectively). Fig. 5(b) represents the behavior of the Tb and DTb/dz for different square cross-sectional sections along of the direction of the main flow, according to Eq. (31).
T Bul k (K)
(b)
380
18
370
16
360
14
350
12 10
340 330
Case II
8
Twm Tb(Z) (DTb/dz)
320
6
Case Constant Tw
310
The results presented for the friction factor and Nusselt number in function of a large range of the Reynolds number for uniform wall temperature present good agreement with the experimental works and correlation of the literature, (Fig. 3(a), (b)) using the turbulent convective heat transfer proposed. Figs. 4(b) and 5(a) show the new results investigated in the present work, note a distortion of the temperatures field and as consequence the variation of the Nusselt number caused mainly by the distribution of the non-uniform wall temperature (Cases I and II, with fluid air and Re = 65,000, respectively). Most applications can be approximated by the functions sine and cosine in the wall, but we are able to resolve by means of the methodology presented, any peripheral heat flux variation that can be expressed by a Fourier expansion (Kays and Crawford [30]). Fig. 5(b) shows the comparisons of the behavior of the curves Tb and DTb/dz to the long of the main direction of the flow for Case II and Case uniform wall temperature. These results can be helpful in the project of thermal devices as in heat transfer and secondary flows in cavities, seals, canal of gas turbines and others.
4 Twm Tb(Z) (DTb/dz)
300 290 0
DTb/dz
5. Conclusions
0.5
1
1.5
2
2.5
3
3.5
2 0 4
Z (m) Fig. 5. (a) Fluid temperature with non-uniform wall temperature (Case II). (b) Behavior Tb for different square cross-sectional sections and cases throughout of the direction of the main flow.
References [1] Z.H. Qin, R.H. Plecther, Large eddy simulation of turbulent heat transfer in a rotating square duct, International Journal Heat Fluid Flow 27 (2006) 371–390. [2] J. Nikuradse, 1926, Untersuchung uber die Geschwindigkeitsverteilung in turbulenten Stromungen, Diss. Göttingen, VDI - forschungsheft 281. [3] F.B. Gessner, A.F. Emery, A Reynolds stress model for turbulent corner flows — part I: development of the model, Journal Fluids Eng. 98 (1976) 261–268. [4] F.B. Gessner, J.K. Po, A Reynolds stress model for turbulent corner flows — part II: comparison between theory and experiment, Journal Fluids Eng. 98 (1976) 269–277.
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851 [5] A. Melling, J.H. Whitelaw, Turbulent flow in a rectangular duct, Journal Fluid Mechanical 78 (1976) 289–315. [6] A. Nakayama, W.L. Chow, D. Sharma, Calculation of fully development turbulent flows in ducts of arbitrary cross-section, Journal Fluid Mechanical 128 (1983) 199–217. [7] H.K. Myon, T. Kobayashi, Numerical Simulation of Three Dimensional Developing Turbulent Flow in a Square Duct with the Anisotropic κ–ε Model, Advances in Numerical Simulation of Turbulent Flows ASME, Fluids Engineering Conference, 1991. Vol.117, Portland, United States of America, 1991, pp. 17–23. [8] M, Assato, Análise numérica do escoamento turbulento em geometrias complexas usando uma formulação implícita, Doctoral Thesis, Departamento de Engenharia Mecânica, Instituto Tecnológico de Aeronáutica — ITA, São José dos campos — SP, Brazil, 2001. [9] M. Assato, M.J.S. De Lemos, Turbulent flow in wavy channels simulated with nonlinear models and a new implicit formulation, Numerical Heat Transfer Part A Applications. 56 (4) (2009) 301–324. [10] D. Home, M.F. Lightstone, M.S. Hamed, Validation of DES-SST based turbulence model for a fully developed turbulent channel flow problem, Numerical Heat Transfer Part A Applications 55 (4) (2009) 337–361. [11] D.D. Luo, C.W. Leung, T.L. Chan, W.O. Wong, Simulation of turbulent flow and forced convection in a triangular duct with internal ribbed surfaces, Numerical Heat Transfer Part A Applications 48 (5) (2005) 447–459. [12] S. Ergin, M. Ota, H. Yamaguchi, Numerical study of periodic turbulent flow through a corrugated duct, Numerical Heat Transfer Part A Applications 40 (2) (2001) 139–156. [13] B.E. Launder, W.M. Ying, Prediction of flow and heat transfer in ducts of square cross section, Proc. Inst. Mech. Eng. 187 (1973) 455–461. [14] A.F. Emery, P.K. Neighbors, F.B. Gessner, The numerical prediction of developing turbulent flow and heat transfer in a square duct, Journal Heat Transfer 102 (1980) 51–57. [15] M. Hirota, H. Fujita, H. Yokosawa, H. Nakai, H. Itoh, Turbulent heat transfer in a square duct, International Journal Heat and fluid flow 18 (1997) 170–180. [16] M. Rokni, Numerical investigation of turbulent fluid flow and heat transfer in complex duct, Doctoral Thesis, Department of Heat and Power Engineering. Lund Institute of Technology, Sweden, 1998. [17] Y. Hongxing, Numerical study of forced turbulent heat convection in a straight square duct, International Journal of Heat and Mass Transfer 52 (2009) 3128–3136. [18] Y.T. Yang, M.L. Hwang, Numerical simulation of turbulent fluid flow and heat transfer characteristics in a rectangular porous channel with periodically spaced heated blocks, Numerical Heat Transfer Part A Applications 54 (8) (2008) 819–836. [19] T.S. Park, Numerical study of turbulent flow and heat transfer in a convex channel of a calorimetric rocket chamber, Numerical Heat Transfer Part A Applications 45 (10) (2004) 1029–1047. [20] J. Zhang, L. Dong, L. Zhou, S. Nieh, Simulation of swirling turbulent flows and heat transfer in a annular duct, Numerical Heat Transfer Part A Applications 44 (6) (2003) 591–609.
851
[21] B. Zheng, C.X. Lin, M.A. Ebadian, Combined turbulent forced convection and thermal radiation in a curved pipe with uniform wall temperature, Numerical Heat Transfer Part A Applications 44 (2) (2003) 149–167. [22] J. Su, A.J. Da Silva Neto, Simultaneous estimation of inlet temperature and wall heat flux in turbulent circular pipe flow, Numerical Heat Transfer Part A Applications 40 (7) (2001) 751–766. [23] A. Saidi, B. Sundén, Numerical simulation of turbulent convective heat transfer in square ribbed ducts, Numerical Heat Transfer Part A Applications 38 (1) (2001) 67–88. [24] M. Rokni, A new low-Reynolds version of an explicit algebraic stress model for turbulent convective heat transfer in ducts, Numerical Heat Transfer Part B Fundamentals 37 (3) (2000) 331–363. [25] A. Valencia, Turbulent flow and heat transfer in a channel with a square bar detached from the wall, Numerical Heat Transfer Part A Applications 37 (3) (2000) 289–306. [26] M.C. Sharatchandra, D.L. Rhode, Turbulent flow and heat transfer in staggered tube banks with displaced tube rows, Numerical Heat Transfer Part A Applications 31 (6) (1997) 611–627. [27] A. Campo, K. Tebeest, U. Lacoa, J.C. Morales, Application of a finite volume based method of lines to turbulent forced convection in circular tubes, Numerical Heat Transfer Part A Applications 30 (5) (1996) 503–517. [28] M. Rokni, B. Sundén, Numerical investigation of turbulent forced convection in ducts with rectangular and trapezoidal cross section area by using different turbulence models, Numerical Heat Transfer Part A Applications 30 (4) (1996) 321–346. [29] G. Yang, M.A. Ebadian, Effect of Reynolds and Prandtl numbers on turbulent convective heat transfer in a three-dimensional square duct, Numerical Heat Transfer Part A Applications 20 (1) (1991) 111–122. [30] W.M. Kays, M. Crawford, Convective Heat and Mass Transfer, McGraw-Hill, New York, USA, 1980, pp. 250–252. [31] E.C. Garcia, Condução, convecção e radiação acopladas em coletores e radiadores solares, Doctor degree thesis, ITA - Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, 1996. [32] C.G. Speziale, On non-linear k-l and k-e models of turbulence, Journal of Fluid Mechanics 178 (1987) 459–475. [33] K. Abe, T. Kondoh, Y. Nagano, An improved k-ɛ model for prediction of turbulent flows with separation and reattachment, Trans. JSME 58 (Series B) (1992) 3003–3010. [34] B.E. Launder, G.J. Reece, W. Rodi, Progress in the development of a Reynolds-Stress turbulence closure, Journal of Fluid Mechanics 68 (3) (1975) 537–566. [35] B.J. Daly, F.H. Harlow, Transport equations in turbulence, Physical Fluids 13 (1970) 2634–2649. [36] M.J. Moran, D.P. Dewitt, B.R. Munson, Introdução à Engenharia de Sistemas Térmicos: Termodinâmica, Mecânica dos Fluidos e Transferência de Calor, LTC Ed. Rio de Janeiro-RJ, Brasil, 2005. [37] S.V. Patankar, Computation of Conduction and Duct Flow Heat Transfer, Innovative Research, Maple Grove, USA, 1991.
Contents lists available at ScienceDirect
International Communications in Heat and Mass Transfer j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / i c h m t
Forced turbulent heat convection in a square duct with non-uniform wall temperature☆ G.A. Rivas ⁎, E.C. Garcia, M. Assato Divisão de Engenharia Mecânica, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil
a r t i c l e
i n f o
Available online 20 April 2011 Keywords: Heat transfer Numerical methods Forced turbulent heat convection Non-uniform wall temperature Square ducts
a b s t r a c t In the present work, the numerical simulation to calculate the problem of the turbulent convection with nonuniform wall temperature in a square cross-section duct was adopted. To solve this problem some assumptions for the flow, such as: the condition of fully developed turbulence and incompressible flow have been assumed. The methodology of the dimensionless energy equation was used to calculate the fluid temperature field in the square cross-section in function of the non-uniform wall temperatures prescribed. Numerical simulations were done using two different turbulent models to resolve the momentum equations and two more models to resolve the energy equation. The models of turbulence k–ε Nonlinear Eddy Viscosity Model (NLEVM) and the Reynolds Stress Model (RSM) were used to determine the turbulent intensities as well as the profiles of axial and secondary mean velocities. The turbulence model RSM was simulated using a commercial software. The thermal field was determined from other two models: Simple Eddy Diffusivity (SED), based in the hypothesis of the constant turbulent Prandtl number; and Generalized Gradient Diffusion Hypothesis (GGDH). In this last model, as the turbulent heat transfer depends on the shear tensions, the anisotropy is considered. These two last equation models of the energy equation of the fluid have been implemented in FORTRAN, a code of programming. The performances of the models were evaluated by validating them based in the experimental and numerical results published in the literature. Two important parameters of great interest in engineering are presented: the friction factor and the Nusselt number. The results of this investigation allow the evaluation of the behavior of the turbulent flow and convective heat fluxes for different square cross-sectional sections throughout the direction of the main flow, which is mainly influenced by the temperature distribution in the wall. Crown Copyright © 2011 Published by Elsevier Ltd. All rights reserved.
1. Introduction Square ducts are widely used in heat transfer devices. For instance in compact heat exchangers, gas turbine cooling systems, cooling channels in combustion chambers and nuclear reactors. The forced turbulent heat convection in a square duct is one of the fundamental problems in the thermal science and fluid mechanics. Recently, Qin and Plecther [1] showed that the Prandtl's secondary flow of the second kind has a significant effect on the transport of heat and momentum as revealed by the recent Large Eddy Simulation (LES) technique. Several experimental and numerical studies have been conducted on turbulent flow though a non-circular duct, specifically, (Nikuradse [2]; Gessner and Emery [3]; Gessner and Po [4]; Melling and Whitelaw [5]; Nakayama et al. [6]; Myon and Kobayashi [7]; Assato [8]; Assato and De Lemos [9]; Home et al. [10]; Luo et al. [11]; Ergin et al. [12] and others). Similarly important works in the turbulent heat convection ☆ Communicated by W.J. Minkowycz. ⁎ Corresponding author. E-mail addresses: [email protected] (G.A. Rivas), [email protected] (E.C. Garcia).
were developed (Launder and Ying [13]; Emery et al. [14]; Hirota et al. [15]; Rokni [16]; Hongxing [17]; Yang and Hwang [18]; Park [19]; Zhang et al. [20]; Zheng et al. [21]; Su and Da Silva Neto [22]; Saidi and Sundén [23]; Rokni [24]; Valencia [25]; Sharatchandra and Rhode [26]; Campo et al. [27]; Rokni and Sundén [28]; Yang and Ebadian [29] and others). The experimental work of Melling and Whitelaw [5] shows detailed characteristics of turbulent flow in a rectangular duct where they used a laser-Doppler anemometer to report the axial development of the mean velocity, secondary mean velocity, etc. Nakayama et al. [6] show the analysis of the fully developed flow field in ducts of rectangular and trapezoidal cross-sections using a finitedifference method with the model of Launder and Ying [13]. On the other hand, Hirota et al. [15] present an experimental work on the turbulent heat transfer in a square duct; they show detailed characteristics of turbulent flow and temperature field. Likewise, Rokni [16] carried out a comparison of four different turbulence models for predicting the turbulent Reynolds stresses and three turbulent heat flux models for ducts square. In the literature on turbulence modeling it is well known that Linear Eddy Viscosity Models (LEVM) can give rise to inaccurate
0735-1933/$ – see front matter. Crown Copyright © 2011 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2011.04.005
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Nomenclature specific heat at constant pressure duct height hydraulic diameter, Dh = 4.A/Pe = 2.L.D/(L + D) pressure gradient at the z direction (longitudinal axis) factor of friction Moody and friction coefficient of Fanning, respectively kf fluid thermal conductivity L duct width Nu Nusselt number Pe perimeter Pk turbulence production term Re Reynolds number T temperature Tb internal flow bulk temperature TWm wall mean temperature T1, T2, T3 and T4 temperature distributions at duct wall (bottom, right side, top and left side) Ub internal flow bulk velocity U, V and W averaged velocity in the Cartesian plane in the direction x, y and z; respectively y+ dimensionless wall distance cp D Dh dP/dz f, Cf
Greek symbols α thermal diffusivity, α = kf/(ρ.cp) η distance normal to the wall ϕ dimensionless temperature distribution μ dynamic viscosity μt turbulent viscosity ρ density σt turbulent Prandtl number
predictions for the Reynolds normal stresses and so that it does not have the ability to predict secondary flows in non-circular ducts due to its isotropic nature (LEVM). In spite of that, they are one of the most popular models in the engineering due to its simplicity, good numerical stability and it can be applied to a wide variety of flows. Thus, the Nonlinear Eddy Viscosity Model (NLEVM) represents a progress of the classical LEVM which permits inequality of the Reynolds normal stresses, a necessary condition for calculating turbulence-driven secondary flow in non-circular ducts within the relative cost of a two-equation formulation. The Reynolds Stress Model (RSM), also called the second order or second moment closure model is very accurate in the calculation of mean flow properties and all Reynolds stresses for many simple and more complex flows including wall jets, asymmetric channel and non-circular duct flows and curved flows, however it presents disadvantages, such as, very large computing costs. For calculating the turbulent heat flux, the Simple Eddy Diffusivity (SED) and Generalized Gradient Diffusion Hypothesis (GGDH) models have been adopted and investigated. Most of the works presented in the literature show results assuming constant temperatures on the wall. However, in many engineering applications the heat flux and surface temperature are non-uniform around the duct, becoming important the knowledge of the variation of the conductance around the duct, Kays and Crawford [30]. According to developments performed by Garcia [31], it is possible to carry out analysis with non-uniform wall temperature boundary conditions. In this case, it is necessary to define a value that represents the mean wall temperatures in a given cross section, TWm. In the present work important results in square cross-section ducts considering the fluids air and water and the influence of the nonuniform wall temperature distribution on the turbulent convective thermal field, are presented.
845
2. Mathematical formulation 2.1. Governing equations The equation systems to be solved are the Reynolds Averaged Navier Stokes (RANS), continuity equation, Eq. (1), momentum equation, Eq. (2), and energy equation, Eq. (3). ∂ U =0 ∂xj j
ð1Þ
" !# ∂Uj ∂ 1 ∂P 1 ∂ ∂Ui 1 ∂ Ui Uj = − + μ + −ρui′uj′ ð2Þ + ρ ∂xi ρ ∂xj ρ ∂xj ∂xj ∂xi ∂xj " # ∂ 1 ∂ μ ∂Tf ′ ′ UT = + −ρuj t ρ ∂xj Pr ∂xj ∂xj j f
ð3Þ
It is considered fully developed turbulent flow and heat transfer and the following hypotheses have been utilized: steady state, condition of nonslip on the wall, and fluid with constant properties. The turbulent Reynolds stress −ρui′uj′ and the turbulent heat flux −ρuj′t ′ are modeled and solved by solution of algebraic or differential expressions. 2.2. Turbulence models for Reynolds stresses 2.2.1. Nonlinear Eddy Viscosity Model (NLEVM) The NLEVM model to reproduce the tensions of Reynolds presents the inclusion of nonlinear terms to the basic constitutive equation trying to capture the sensitivity of the curvature of the stream lines. This model is based on the initial proposal of Speziale [32]. The Reynolds average equations, Eqs. ((1)–(3)) are applied for the device of Fig. 1(a) and (b). The components of velocities U, V represent the secondary flow and the axial velocity component W, the velocity of the main flow. The equations of transport in tensorial form for the turbulent kinetic energy, κ, and the rate of dissipation ε, respectively, they are given by: Ui
∂k ∂ μ t ∂k = + Pk −ε ∂xi ∂xi ρσk ∂xi
ð4Þ
Ui
∂ε ∂ μ t ∂ε ε ε2 = + c1 Pk −c2 f2 : k k ∂xi ∂xi ρσε ∂xi
ð5Þ
The symbols Pk and μ t, represent the rate of the turbulent kinetic energy production and the turbulent viscosity, respectively, and thus, we have: Pk = τij
∂Ui ; ∂xj
2
μ t = cμ fμ ρ
k : ε
ð6Þ
In the present work for NLEVM, the formulations of Low Reynolds Number will be assumed for wall treatment. The damping functions f2 and fμ are observed in Eqs. (5) and (6) and shown in Table 1. These functions and the constant C1 and C2 have been used together with the equations κ–ε, the subscribed letter P refers to the nodal point near to the wall. Thus UP and kP are the values of the velocity and kinetic energy in this point, respectively. The constant cμ , c1, c2, σk and σε adopt the values of 0.09; 1.5; 1.9; 1.4 and 1.3; respectively. The new constitutive relation for the tensions of Reynolds in the model NLEVM, assumed in the thesis of Assato (2001), is given by: NL L k 1 τij = μ t Sij + c1NL μ t Sik Skj − Skl Skl δij : ð7Þ ε 3 This expression shows that the second term of the right side of Eq. (7), represents the nonlinear term added the original constitutive
846
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Fig. 1. (a) Fully developed turbulent flow in rectangular ducts, (b) Rectangular duct — system of reference and transversal section.
relation. This quadratic term represents the degree of anisotropy between the normal tensions of Reynolds, which makes it possible to predict the presence of the secondary flow in non-circular ducts. The values of c1NL proposed by Speziale [32] are equal to 1.68. Here, c1NL will be analyzed and will adopt values for the formulation of Low Reynolds Numbers. The tensions of Reynolds, normal and shear, are presented in Eqs. (8), (9) and (10), are expressed as: " # k 1 ∂W 2 2 ∂W 2 τxx = c1NL μ t − ; ε 3 ∂x 3 ∂y " # k 1 ∂W 2 2 ∂W 2 − τyy = c1NL μ t ε 3 ∂y 3 ∂x τxy = c1NL μ t
k ∂W ∂W ; ε ∂x ∂y
τxz = μ t
∂W ∂W ; τyz = μ t : ∂x ∂y
ð8Þ
∂ ∂ ρuk ui′uj′ ðbÞ ρui′uj′ ðaÞ + ∂t ∂xk ∂ ∂ ∂ ρui′uj′uk′ + p δkj uj′ + δik uj′ ðcÞ + μ u′u′ ðdÞ =− ∂xk ∂xk ∂xk i j ! ! ð12Þ ∂uj′ ∂uj ∂u′i ∂u ð gÞ −ρ ui′uk′ + uj′uk′ i ðeÞ−ρβ gi uj′θ ð f Þ + p + ∂xk ∂xk ∂xj ∂xi −2μ
ð9Þ
The following differences for the normal tensions of Reynolds are presented and have been observed for this type of flow. " # k ∂W 2 ∂W 2 − τyy −τxx = c1NL μ t ε ∂y ∂x
equations of transport. The individual Reynolds tensions are utilized to close the average Reynolds equations of the momentum. This model has shown superiority regarding the models of two equations in complex flows that involve swirl, rotation, etc. The exact equations of transport for the Reynolds tensions, ρui′uj′, can be written as follows:
Table 1 Low Reynolds Number formulation for wall treatment using the NLEVM model.
ð10Þ fu
In such a way in Eq. (6), including the derivatives above the tensions of Reynolds, the turbulence production term, is expressed as: Pk = τxz
∂W ∂W + τyz : ∂x ∂y
Low Reynolds Numbers formulation proposed by Abe, Nagano and Kondoh [33] 8 2 2 !2 3 9 = n h i o2 < 5 Þ0:25 4− k =vε 5 1− exp − ðvε14v 1+ exp 0:75 : ; 2 200 ðk =vεÞ 2 " #)2 8 2 !2 39 < = k =vε ðvεÞ0:25 5 1− exp − 1−0:3 exp4− : ; 3:1v 6:5
( f2
ð11Þ
2.2.2. Reynolds Stress Model (RSM) The most complex turbulence is the Reynolds Stress Model (RSM), also called second order, it involves calculations of the tensions of Reynolds in an individual form, ρui′uj′, used for this differential
∂ui′ ∂u′j ′ εikm + ui′um ′ εjkm ðiÞ + Sð jÞ: ðhÞ−2ρΩk uj′um ∂xk ∂xk
σk σε C1 C2 τw
1.4 1.3 1.5 1.9 ∂UP μ ∂yP
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
Where the respective letters represent: (a) local derivative of the time; (b) Cij≡ convection; (c) DT,ij≡ turbulent diffusion; (d) DL,ij≡ molecular diffusion; (e) Pij≡ production term of tensions; (f) Gij≡ buoyancy production term; (g) ϕij≡ term of pressure–tension (redistribution); (h) εij≡ term of dissipation; (i) Fij≡ term production for the rotation system; and (j) Sj≡ source term. The terms of the exact equations, presented previously,Cij, DL,ij, Pij and Fij do not require modeling. However, the terms DT,ij, Gij, ϕij and εij need to be modeled to close the equations. For the present analysis, the model LRR [34] is chosen, which assumes that the correlation velocity pressure is a linear function of the anisotropy tensor LRR in the phenomenology of the redistribution, ϕij. For the treatment of the wall, it is also assumed the Low Reynolds Numbers and the periodic conditions, Rokni [16]. This model had been simulated in the commercial code Fluent 6.3. 2.3. Turbulence models for turbulent heat flux
and ∫ W:Tf :dA Tb =
μ t ∂Tf σT ∂xj
∂Tf k = −ρCt u′j u′k : ε ∂xk
2.3.3. Dimensionless Energy Equation for SED and GGDH Models For a given cross section of area A, it is possible to define a mean velocity Ub and a bulk temperature Tb, express as:
1 Ub = ∬W:dx:dy A
ð15Þ
Table 2 Fluid properties (at 300 K) for internal flow (Moran et al. [36]).
Specific heat, cp [kJ/kg.K] Specific mass ρ [kg/m3] Thermal conductivity kf [W/m.K] Thermal diffusivity α [m²/s] Dynamic viscosity μ [N.s/m²]
ð16Þ
∂Tf ∂x
+V
∂Tf ∂y
+W
∂Tf
"
∂Tf ∂ −u′ t ′ α ∂z ∂x ∂x !# −
∂Tf ∂ −w′ t ′ α ∂z ∂z
! +
∂Tf ∂ −v′ t ′ α ∂y ∂y
= 0:
!
ð18Þ
The properties of the fluid are shown in Table 2, and the following considerations are applied to obtain the variables in dimensionless form: X=
x ; Dh
ð19Þ
Y=
y Dh
ð20Þ
ð14Þ
The constant Ct, assuming the value of 0.3, is adopted. The main advantage of this model is that it considers the anisotropic behavior of the fluid heat transport in ducts (Table 2).
Fluid
1 ∬W:Tf :dx:dy: A:Ub
It is possible to develop a formula similar to Kays and Crawford [30], the new expression for the turbulent energy equation, is presented as:
+
!
=
1 L 1 D 1 L 1 D :∫0 T1 ð0; yÞ:dy+ :∫0 T2 ðx; 0Þ:dx+ :∫0 T3 ðD; yÞ:dy + :∫0 T4 ðx; LÞ:dx L D L D : = 2ðL + DÞ
ð13Þ
2.3.2. Generalized Gradient Diffusion Hypothesis (GGDH) Daly and Harlow [35] present the following formulation to the turbulent heat flux:
ρu′j t ′
Ub :A
ð17Þ
U ρu′j t ′ = −
A
Kays and Crawford [30] developed an applicable formulation to rectangular cross section ducts. They considered the boundary conditions with prescribed uniform wall temperatures at the cross section, and along the duct length. According to developments performed by Garcia [31], it is possible to carry out analysis with non-uniform wall temperature boundary conditions. In this case, it is necessary to define a value that represents the mean wall temperatures in a given cross section, TWm, given as:
TWm
2.3.1. Simple Eddy Diffusivity (SED) This method is based on the Boussinesq viscosity model. The turbulent diffusivity for the energy equation can be expressed as: αt = μt=ρσt , where the turbulent Prandtl number σt needs to be given. The SED model assumes that the turbulent Prandtl number is constant in the entire region, for the air σt it assumes values of 0.89, independent on the wall proximity effect.
847
Air
Water
1005.7
4179.0
1.1614
997.606
2.63 × 10− 2
61.3 × 10− 2
225 × 10− 7
6.15 × 10− 7
1.846 × 10− 5
88.5 × 10− 5
and α: TWm −Tf : ϕ= dTb Ub :D2h : dz
ð21Þ
Replacing Eqs. (13), (14), (19) up to (21) in Eq. (18), dimensionless energy equation for SED and GGDH becomes, respectively:
1 ∂ ∂ϕ 1 ∂ ∂ϕ D ∂ϕ ∂ϕ ðα + αt Þ ðα + αt Þ U + − h +V α ∂X α ∂Y α ∂X ∂Y ∂X ∂Y W ϕ =− ð22Þ UB ϕB ∂ϕ 1 ∂ ∂ϕ 0 0 k 0 0 k α + Ct u u + α + Ct v v ε ∂X α ∂Y ε ∂Y D ∂ϕ ∂ϕ W ϕ C ∂ 0 0 k ∂ϕ uv − t +V =− − h U UB ϕB ε ∂Y α α ∂X ∂X ∂Y C ∂ 0 0 k ∂ϕ vu − t ε ∂X α ∂Y
1 ∂ α ∂X
ð23Þ
848
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
The fluid temperature field Tf can be replaced by Tb and Eq. (21) can be expressed as: ϕb =
α:ðTWm −Tb Þ ; dTb Ub :D2h : dz
ð24Þ
and D2h :Ub dTb : :ϕ : ð25Þ α dz b From Eq. (21), Eq. (26) is obtained, and applying this in Eq. (16), one gets Eq. (27): Tb = TWm −
Tf = TWm −
2 ϕ:Ub :Dh :
α
dTb dz
ϕb =
D2h :∬W:ϕ:dX:dY: A:Ub
ð28Þ
It is possible to compute the heat transfer rate per unit length on the wall surface, q′, Eq. (29). It depends on values for TWm, Tb, and on the average heat convection coefficient, h. From fluid enthalpy derivative gradient [dhb = cp.dTb], the heat transfer rate per unit length in the fluid, q′f, can be expressed by Eq. (30). ;
q = Pe :h:ðTWm −Tb Þ;
ð29Þ
and
qf′ = ρ:Ub :A:cp : ;
ð26Þ
dTb dz
ð30Þ
When Eq. (29) is made equal to Eq. (30), integrating two cross sections (inlet, namely z1, and outlet, namely z2), it is possible to develop the resulting expression, Eq. (31).
and
Tb = TWm −
Replacing Eq. (27) in Eq. (24), and using also Eqs. (19) and (20), the dimensionless bulk temperature is obtained:
D2h dTb : :∬W:ϕ:dx:dy: α:A dz
ð27Þ
− Tbz2 = TWm − TWm −Tbz1 :e
Pe 2 α : :Nu:ðz1 −z2 Þ 2A Ub
ð31Þ
Fig. 2. (a) Grid 120 × 120 for numerical simulation. (b) Secondary flow contours and comparisons of the axial mean velocity with Melling and Whitelaw [5] for water utilizing NLEVM model.
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
From Eq. (31), the bulk temperature longitudinal (z-axis) variation ΔTb is obtained. It is done by “cutting” the duct into a lot of segments and applying the numerical method to find ΔTb at each finite cross section. For a given bulk temperature at the duct inlet section (Tb1), after solving the equation system, duct outlet bulk temperature (Tb2) is computed from that expression, Eq. (31). The dimensionless boundary conditions are given by the following equations:
ϕð0; Y Þ =
α:ðTWm −T1 Þ α:ðTWm −T2 Þ ; and ϕð X; 0Þ = dTb dTb 2 Ub :Dh : Ub :D2h : dz dz
ð32Þ
When considering uniform wall temperature, Eqs. (32) and (33) are equal to zero, and for these particular conditions, it is possible to notice that these boundary conditions are not functions of dTb/dz. That simplification becomes equal to the one studied by Patankar [37]. Eqs. (24), (25) and (26), as well as the boundary conditions from Eqs. (32) and (33), form a set of differential equations, in which ϕ and dTb/dz parameters are unknown. When that equation system is solved, it is possible to obtain Tf. 2.3.4. Additional equations Additional equations were utilized for the calculation of the factor of friction Moody, f; coefficient of friction of Fanning, Cf; Prandtl law; local Nusselt number for the Low Reynolds formulation (Rokni [16]), Nuxp and correlation of Gnielinsky, Nu, respectively; these equations are given by:
(a) 10 9
Numerical Simulation
8
f≡
Model Non Linear k-e Model RSM
7 6
Cf (x10-3)
α:ðTWm −T3 Þ D ; ϕ =Dh ; Y = dTb Ub :D2h : dz
849
L α:ðTWm −T4 Þ : ð33Þ and ϕ X; = dTb Dh Ub :D2h : dz
−ðdP = dzÞ:Dh ρUB2 =2 ;
Cf =
ð34Þ
f ; 4
ð35Þ
pffiffiffi 1 pffiffiffi = 2 log Re f −0:8; f ðT −TP Þ Nuxp = Dh w ηðTw −Tb Þ
5
4
ð36Þ ð37Þ
Methods Experimental and Correlations 3
Prandtl Friction Law Harnett et al. Leutheusser Launder e Ying Lund
(a)
1
1 0.65
0.8
0.65
0.70
0.70
2 2
3
4
5
6
7 8 9 10
20
30
0.80
0.6
Re (x104)
103
0.90 0.95
(b)
0.6
0.85
0.85 0.90
0.4
0.4
0.95
0.2
0.2
Numerical Simulation
0 -1
SED Model GGDH Model
0.8
0.75 0.80
0.75
-0.8
-0.6
-0.4
-0.2
Numerical Prediction
0
0.2
0.4
0.6
0.8
0 1
Experimental
Nu/Pr0.4
(b) 102
Methods experimental - Correlations Lowdermilk et al. Gnielinsky Brundrett & Burroughs
101
104
105
106
Re Fig. 3. (a) Friction coefficient for fully developed flow, (b) Nusselt number dependence on Reynolds number for fully developed flow.
Fig. 4. (a) Results (RSM–SED) numerical and (Hirota et al. [15]) experimental for mean temperature (uniform wall temperature). (b) Fluid temperature with non-uniform wall temperature (Case I).
850
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851
obtained (dTb/dz b tolerance). This is the end of the second iterative loop;
and 2
3 ðf = 8ÞðRe−1000ÞPr 4 1 5: Nu = 1 1 + 12:7ðf =8Þ =2 Pr =2 −1
ð38Þ
3. Numerical implementation After applying the method of finite differences to the algebraic equations, to obtain the temperature fields, the following five steps indicate the developed methodology in the numerical solution (Garcia [31]): Step 1: To define the function value of the non-uniform temperatures in the walls of the duct TWm = f(T1(0,y), T2(x,0), T3(D,y), T4(x,L)), what that can be expressed by a Fourier expansion; Step 2: To obtain velocity field and estimated values for Ub, TWm and dTb/dz; Step 3: Equations for the boundary conditions are evaluated (Eqs. (32) and (33)); Step 4: Dimensionless energy equation (temperature field, Ø) at Eq. (24) is solved and Øb is computed according to Eq. (29), until convergence is obtained (Øb b tolerance). This is the end of the first iterative loop; Step 5: A value for dTb/dz is computed in accordance with Eq. (25). Boundary conditions are updated (step 3) to obtain a solution for the new temperature field (step 4), until convergence is
(a)
For all steps, “tolerance of 10− 7” is the value to be accomplished by the convergence criteria, which is applicable to Øb (dimensionless bulk temperature), dTb/dz and Ø (dimensionless temperature field). 4. Results and discussion 4.1. Fluid flow and heat transfer field Fig. 2(a) shows the utilized grid (120 × 120) in the numerical simulation for the formulations of Low Reynolds, Fig. 2(b) represents the secondary flow contours and comparisons of the velocity profile (NLEVM, Assato [8]) with the experimental work of the Melling and Whitelaw [5] for fluid water and Re = 42,000. The predicted distributions of the friction coefficient (NLEVM and RSM) and Nusselt number (SED and GGDH) dependence on Reynolds number for fully developed flow and heat transfer in a square duct are shown in Fig. 3(a) and (b), respectively. In Fig. 4(a), the comparisons of the results (RSM–SED) numerical with the experimental for temperature profile (wall constant temperature) (TWm − Tf)/(TWm − TC) with fluid air and Re = 65,000 (Hirota [15]) are shown, while Fig. 4(b) shows the variation of the temperature profile with non-uniform wall temperature: south = 400 K, north = 373 K, east = 393 K, and west = 353 K; it called of Case I. Already Fig. 5(a) shows the variation of the temperature profile with non-uniform wall temperature, this it is represented by means of functions sine (Case II), south = (350–20Sin (ζ))K, north = (400–50Sin(ζ))K, east = (330 + 20Sin(ζ))K, and west = (350 + 50Sin(ζ))K. Where ζ is function of the radians (0–π/ 2) and i,j (points number of the grid in the direction x and y, respectively). Fig. 5(b) represents the behavior of the Tb and DTb/dz for different square cross-sectional sections along of the direction of the main flow, according to Eq. (31).
T Bul k (K)
(b)
380
18
370
16
360
14
350
12 10
340 330
Case II
8
Twm Tb(Z) (DTb/dz)
320
6
Case Constant Tw
310
The results presented for the friction factor and Nusselt number in function of a large range of the Reynolds number for uniform wall temperature present good agreement with the experimental works and correlation of the literature, (Fig. 3(a), (b)) using the turbulent convective heat transfer proposed. Figs. 4(b) and 5(a) show the new results investigated in the present work, note a distortion of the temperatures field and as consequence the variation of the Nusselt number caused mainly by the distribution of the non-uniform wall temperature (Cases I and II, with fluid air and Re = 65,000, respectively). Most applications can be approximated by the functions sine and cosine in the wall, but we are able to resolve by means of the methodology presented, any peripheral heat flux variation that can be expressed by a Fourier expansion (Kays and Crawford [30]). Fig. 5(b) shows the comparisons of the behavior of the curves Tb and DTb/dz to the long of the main direction of the flow for Case II and Case uniform wall temperature. These results can be helpful in the project of thermal devices as in heat transfer and secondary flows in cavities, seals, canal of gas turbines and others.
4 Twm Tb(Z) (DTb/dz)
300 290 0
DTb/dz
5. Conclusions
0.5
1
1.5
2
2.5
3
3.5
2 0 4
Z (m) Fig. 5. (a) Fluid temperature with non-uniform wall temperature (Case II). (b) Behavior Tb for different square cross-sectional sections and cases throughout of the direction of the main flow.
References [1] Z.H. Qin, R.H. Plecther, Large eddy simulation of turbulent heat transfer in a rotating square duct, International Journal Heat Fluid Flow 27 (2006) 371–390. [2] J. Nikuradse, 1926, Untersuchung uber die Geschwindigkeitsverteilung in turbulenten Stromungen, Diss. Göttingen, VDI - forschungsheft 281. [3] F.B. Gessner, A.F. Emery, A Reynolds stress model for turbulent corner flows — part I: development of the model, Journal Fluids Eng. 98 (1976) 261–268. [4] F.B. Gessner, J.K. Po, A Reynolds stress model for turbulent corner flows — part II: comparison between theory and experiment, Journal Fluids Eng. 98 (1976) 269–277.
G.A. Rivas et al. / International Communications in Heat and Mass Transfer 38 (2011) 844–851 [5] A. Melling, J.H. Whitelaw, Turbulent flow in a rectangular duct, Journal Fluid Mechanical 78 (1976) 289–315. [6] A. Nakayama, W.L. Chow, D. Sharma, Calculation of fully development turbulent flows in ducts of arbitrary cross-section, Journal Fluid Mechanical 128 (1983) 199–217. [7] H.K. Myon, T. Kobayashi, Numerical Simulation of Three Dimensional Developing Turbulent Flow in a Square Duct with the Anisotropic κ–ε Model, Advances in Numerical Simulation of Turbulent Flows ASME, Fluids Engineering Conference, 1991. Vol.117, Portland, United States of America, 1991, pp. 17–23. [8] M, Assato, Análise numérica do escoamento turbulento em geometrias complexas usando uma formulação implícita, Doctoral Thesis, Departamento de Engenharia Mecânica, Instituto Tecnológico de Aeronáutica — ITA, São José dos campos — SP, Brazil, 2001. [9] M. Assato, M.J.S. De Lemos, Turbulent flow in wavy channels simulated with nonlinear models and a new implicit formulation, Numerical Heat Transfer Part A Applications. 56 (4) (2009) 301–324. [10] D. Home, M.F. Lightstone, M.S. Hamed, Validation of DES-SST based turbulence model for a fully developed turbulent channel flow problem, Numerical Heat Transfer Part A Applications 55 (4) (2009) 337–361. [11] D.D. Luo, C.W. Leung, T.L. Chan, W.O. Wong, Simulation of turbulent flow and forced convection in a triangular duct with internal ribbed surfaces, Numerical Heat Transfer Part A Applications 48 (5) (2005) 447–459. [12] S. Ergin, M. Ota, H. Yamaguchi, Numerical study of periodic turbulent flow through a corrugated duct, Numerical Heat Transfer Part A Applications 40 (2) (2001) 139–156. [13] B.E. Launder, W.M. Ying, Prediction of flow and heat transfer in ducts of square cross section, Proc. Inst. Mech. Eng. 187 (1973) 455–461. [14] A.F. Emery, P.K. Neighbors, F.B. Gessner, The numerical prediction of developing turbulent flow and heat transfer in a square duct, Journal Heat Transfer 102 (1980) 51–57. [15] M. Hirota, H. Fujita, H. Yokosawa, H. Nakai, H. Itoh, Turbulent heat transfer in a square duct, International Journal Heat and fluid flow 18 (1997) 170–180. [16] M. Rokni, Numerical investigation of turbulent fluid flow and heat transfer in complex duct, Doctoral Thesis, Department of Heat and Power Engineering. Lund Institute of Technology, Sweden, 1998. [17] Y. Hongxing, Numerical study of forced turbulent heat convection in a straight square duct, International Journal of Heat and Mass Transfer 52 (2009) 3128–3136. [18] Y.T. Yang, M.L. Hwang, Numerical simulation of turbulent fluid flow and heat transfer characteristics in a rectangular porous channel with periodically spaced heated blocks, Numerical Heat Transfer Part A Applications 54 (8) (2008) 819–836. [19] T.S. Park, Numerical study of turbulent flow and heat transfer in a convex channel of a calorimetric rocket chamber, Numerical Heat Transfer Part A Applications 45 (10) (2004) 1029–1047. [20] J. Zhang, L. Dong, L. Zhou, S. Nieh, Simulation of swirling turbulent flows and heat transfer in a annular duct, Numerical Heat Transfer Part A Applications 44 (6) (2003) 591–609.
851
[21] B. Zheng, C.X. Lin, M.A. Ebadian, Combined turbulent forced convection and thermal radiation in a curved pipe with uniform wall temperature, Numerical Heat Transfer Part A Applications 44 (2) (2003) 149–167. [22] J. Su, A.J. Da Silva Neto, Simultaneous estimation of inlet temperature and wall heat flux in turbulent circular pipe flow, Numerical Heat Transfer Part A Applications 40 (7) (2001) 751–766. [23] A. Saidi, B. Sundén, Numerical simulation of turbulent convective heat transfer in square ribbed ducts, Numerical Heat Transfer Part A Applications 38 (1) (2001) 67–88. [24] M. Rokni, A new low-Reynolds version of an explicit algebraic stress model for turbulent convective heat transfer in ducts, Numerical Heat Transfer Part B Fundamentals 37 (3) (2000) 331–363. [25] A. Valencia, Turbulent flow and heat transfer in a channel with a square bar detached from the wall, Numerical Heat Transfer Part A Applications 37 (3) (2000) 289–306. [26] M.C. Sharatchandra, D.L. Rhode, Turbulent flow and heat transfer in staggered tube banks with displaced tube rows, Numerical Heat Transfer Part A Applications 31 (6) (1997) 611–627. [27] A. Campo, K. Tebeest, U. Lacoa, J.C. Morales, Application of a finite volume based method of lines to turbulent forced convection in circular tubes, Numerical Heat Transfer Part A Applications 30 (5) (1996) 503–517. [28] M. Rokni, B. Sundén, Numerical investigation of turbulent forced convection in ducts with rectangular and trapezoidal cross section area by using different turbulence models, Numerical Heat Transfer Part A Applications 30 (4) (1996) 321–346. [29] G. Yang, M.A. Ebadian, Effect of Reynolds and Prandtl numbers on turbulent convective heat transfer in a three-dimensional square duct, Numerical Heat Transfer Part A Applications 20 (1) (1991) 111–122. [30] W.M. Kays, M. Crawford, Convective Heat and Mass Transfer, McGraw-Hill, New York, USA, 1980, pp. 250–252. [31] E.C. Garcia, Condução, convecção e radiação acopladas em coletores e radiadores solares, Doctor degree thesis, ITA - Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil, 1996. [32] C.G. Speziale, On non-linear k-l and k-e models of turbulence, Journal of Fluid Mechanics 178 (1987) 459–475. [33] K. Abe, T. Kondoh, Y. Nagano, An improved k-ɛ model for prediction of turbulent flows with separation and reattachment, Trans. JSME 58 (Series B) (1992) 3003–3010. [34] B.E. Launder, G.J. Reece, W. Rodi, Progress in the development of a Reynolds-Stress turbulence closure, Journal of Fluid Mechanics 68 (3) (1975) 537–566. [35] B.J. Daly, F.H. Harlow, Transport equations in turbulence, Physical Fluids 13 (1970) 2634–2649. [36] M.J. Moran, D.P. Dewitt, B.R. Munson, Introdução à Engenharia de Sistemas Térmicos: Termodinâmica, Mecânica dos Fluidos e Transferência de Calor, LTC Ed. Rio de Janeiro-RJ, Brasil, 2005. [37] S.V. Patankar, Computation of Conduction and Duct Flow Heat Transfer, Innovative Research, Maple Grove, USA, 1991.