- Email: [email protected]
Modeling flow and heat transfer in tubes using a fast CFD formulation
Computers and Chemical Engineering 25 (2001) 713– 722 www.elsevier.com/locate/compchemeng
Modeling flow and heat transfer in tubes using a fast CFD formulation E.R.L. Mercado, V.C. Souza, R. Guirardello, J.R. Nunhez * Faculdade de Engenharia Quı´mica, UNICAMP, CP 6066, 13083 -970 Campinas, SP, Brazil Received 10 May 2000; accepted 5 January 2001
Abstract An approach to study turbulent flow and conjugate heat transfer in tubes is proposed in this work. Instead of using the usual finite element or finite volume methods, this formulation applies a different technique that calculates both for the flow and heat transfer. It discretizes the flow in the radial direction using a fourth order finite differences method, which is more accurate than the traditional second order schemes. Using this technique, a system composed of several ordinary differential equations for the temperature and a set of linear equations for the velocities and pressure gradient is obtained. The equations are then integrated in the axial direction using a fourth order Runge–Kutta method. The values of viscosity, density and thermal conductivity are dependent on temperature, which makes the model suitable for the calculation of high temperature gradients, as in the case of refinery fired heaters. The turbulence is taken into account using a zero order turbulence model. A very interesting feature of this formulation is that for the laminar case the method is non-iterative, which makes it more desirable and faster than conventional computational fluid dynamics (CFD) packages. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Turbulent flow; Conjugate heat transfer and fluid flow; Computational fluid dynamics
Nomenclature C. p dP/dz H K kt K lm L P P( R R T T( Tw 6q Vr 6¯ r
Fluid thermal capacity Pressure drop Heat transfer coefficient Thermal conductivity Turbulent conductivity Constant= −(dP( /dz) Mixing length Tube length Pressure Average pressure (in time) Radial position Tube radius Temperature of fluid Average fluid temperature (in time) Wall temperature Angular velocity Radial velocity Average radial velocity
* Corresponding author. E-mail address: [email protected] (J.R. Nunhez). 0098-1354/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 6 7 2 - X
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
714
Vz 6¯ z Z Dr Dz v vt z
Axial velocity Average axial velocity (in time) Axial position Radial distance Axial distance between points Fluid viscosity Turbulent viscosity Fluid density
1. Introduction The heating and cooling inside tubes have been among the most important processes in the engineering field such as petrochemical fired heaters and petrol cracking. The applications are uncountable. The modeling of these processes, based on the conservation of mass, momentum and energy associated with its boundary conditions, normally lead to a set of partial differential equations with no analytical solution. Many details of the flow cannot be captured by experimental analysis, therefore, numerical procedures are needed to have a deeper understanding of these processes. Travello and Dias (1984) have developed a model to solve the energy equation in tubes under laminar flow. The fluid is incompressible and it is assumed a parabolic profile for the axial velocity. The fluid properties are constant and the work is primarily concerned to analyze axial conduction at the tube walls. Faghri and Sparrow (1980) also studied laminar flow in tubes. They were also interested in determining axial conduction at the wall of the tubes. They assume a parabolic profile for the axial velocity and the properties are considered independent on temperature. Barozzi and Pagliarini (1985) developed a method combining the finite element method with the superposition principle to solve simultaneously the momentum and energy equations in order to analyze axial conduction at the walls. The physical properties are also assumed to be constant. Bilir (1994) analyzed the conjugated heat transfer in tubes for laminar flow. He assumed that the wall temperature is constant and he also analyzed axial conduction at the walls. The physical properties again are assumed to be constant. Martinuzzi and Pollard (1989a,b) compared six turbulence models in tubes and interesting detail about the numerical methods were explained. They arrived at the conclusion that the low Reynolds number k– m model better predicts turbulence in tubes. However, physical properties are also constant. In order to get more insight of how these systems behave, a computational fluid dynamics model has been
developed to calculate the velocities, pressure and temperature. The flow can be either laminar or turbulent. Also, as the fluid properties are dependent on temperature, the model is suitable to investigate the flow inside petrochemical fired heaters. The hypotheses of the model are. The tube wall temperature is assumed to be constant, which is acceptable for the design of fired heaters (Wimpress, 1963). The flow is steady-state. It is considered that the inlet flow is already fully developed, because the flow comes from a previous tube section joined to the tube portion under consideration. This assumption is very important. Otherwise the radial velocity would not be small and the model approximations would not apply to the problem. Radial velocities are small in comparison to the axial velocities. The flow is symmetric about the axial axis (there is no variation in the angular direction and angular velocity is negligible). Physical properties are dependant on temperature, but not on pressure. There is no phase change and the fluid is a Newtonian liquid.
2. Modeling The exact model equations for tubes are well known and can be found in textbooks such as (Bird, Stewart & Lightfoot, 1982). Under the hypothesis described above, a study of the order of magnitude was applied to these equations and several terms are negligible in comparison to others. The resulting equations are described as follows, for the laminar and turbulent axisymmetric model. It is important to notice that the radial velocities can be calculated from the continuity equation, since in this equation both terms are of the same order of magnitude.
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
2.1. Laminar flow
Tube wall (r= R).
2.1.1. Momentum conser6ation Radial direction.
Vr = 0
(P =0 (r
−k
(P 1 ( (6 0= − + r ·v z (r (z r (r
(T 1 ( (T = r ·k (z r (r (r
(T = − h(Tw − T) (r
(12)
(2)
2.1.2. Energy conser6ation z · C. p · 6z
(11)
(1)
Axial direction.
715
(3)
Tube entrance (z= 0). At the tube entrance, it is considered that the flow is already fully developed. It is as if the fluid came from a previous tube section joined to the tube portion under consideration. As already mentioned, this is important to be considered since the model does not apply for situations where radial velocities are high. A constant temperature profile is assumed for the temperature since the wall temperature starts acting on the fluid only at the section being considered:
2.1.3. Continuity equation (13)
T= To ( 1 ( (z · r · 6r)+ (z · 6z ) = 0 (z r (r
(4)
=
2.2. Turbulent flow
(P( =0 (r
(5)
Axial direction.
(P( 1 ( (6¯ 0= − + r(v +vt) z (z r (r (r
(6)
z · 6z · 2 · y · r dr
(14)
(T( 1 ( (T( = r(k + kt) (z r (r (r
1 ( ( (z · r · 6¯ r)+ (z · 6¯ z ) = 0 r (r (z
In order to evaluate the turbulent viscosity and thermal conductivity the mixing length hypothesis is assumed (Prandtl theory). The mixing length theory assumes that: Turbulent viscosity:
) ) (6¯ z (r
(15)
Turbulent thermal conductivity:
kt = z · C. p · l 2m (7)
2.2.3. Continuity equation
) ) (6¯ z (r
(16)
In order to evaluate vt and kt, it is necessary to know the mixing length (lm). There are several expressions for tubes. The one chosen is by (Rodi, 1984):
(8)
lm r =0.14−0.08 R R
2
− 0.06
r R
4
(17)
2.4. Numerical calculation
2.3. Boundary conditions
Symmetry line (r= 0).
R
vt = z · l 2m
2.2.2. Energy conser6ation
&
0
2.2.1. Momentum conser6ation Radial direction.
z · C. p · 6z
Mass conservation (05 z5L).
(6z =0 (r
(9)
(T =0 (r
(10)
The fourth order finite differences method is applied to the radial direction. This leads to a linear set of equations for the velocities and a set of first order differential equations for the temperature. The discretization in the radial direction is applied to (m+1) points and there are (n+ 1) points in the axial direction. Fig. 1 shows the mesh used for the calculations.
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
716
Fig. 1. Mesh used for the calculations.
2.4.1. Laminar flow Navier –Stokes (from Eq. (1), after integrating in the r direction):
2.4.2. Turbulent flow From the continuity equations and the assumptions made, it follows that P= P(z), so the following is true:
(6z r (P = i (r ri 2vi (z ri
−
)
)
(18)
A system of m +1 equations is generated. However, there are m+2 unknowns (m + 1 velocities and pressure drop). An additional equation is obtained through a mass balance in the cross section (Eq. (14)). Integration is performed using Simpson’s rule. The linear system is solved using Gaussian elimination. Energy conservation (from Eq. (2) and isolating (#T/#z)):
)
1 (T = (z ri zi · cˆpi · 6zi
) )
)n
(k k (T ( 2T + i +ki 2 (r ri ri (r ri (r ri
(19)
The set of equations is solved by considering the calculation of the axial velocities as a subrouting in the resolution of the ordinary differential equation given by (Eq. (19)). The differential equations are solved by the Runge –Kutta –Fehlberg method, where at each required calculation the axial velocity is obtained by (Eq. (18)) by solving a set of linear equations. Continuity equation: zi · 6ri +ri
)
)
( ( (z · 6r ) = −ri (z · 6z ) (r (z ri zi
(P( =K (z
(21)
From Eq. (6); after integrating: Kr (6¯ z =− (r 2
(v + vt)
(22)
From Eq. (15). − z · f(r)
d6¯ z dr
2
+v
d6¯ z K ·r + =0 dr 2
(23)
Therefore the velocities can be calculated by the following relation, which is one root of Eq. (23)
)
(6¯ z
v 2i + 2Krizi f(ri )− vi = (r ri − 2zi f(ri ) where:
f(r)= R 2 0.14− 0.08
r R
2
(24)
n
− 0.06
r R
4 2
(20)
Applying finite differences to Eq. (20), a system of m + 1 linear equations for the radial velocity is obtained and it is solved using Gaussian elimination, after calculating the values of the axial velocity. At the tube entrance, a uniform temperature profile is assumed. From this, all physical properties are calculated. The axial velocities at the different radial positions and the pressure drop are calculated according to Eq. (18) and a set of first order differential equations is obtained according to Eq. (19), which is calculated using the fourth order Runge – Kutta method. The procedure is repeated for the whole length of the tube. After evaluating all axial velocities, the radial velocities are calculated according to Eq. (20).
Fig. 2. Laminar temperature profile (Re= 1800).
(25)
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
717
A uniform temperature profile is given at the tube entrance, similar to the linear case. The physical properties are estimated and a set of linear equations is obtained for the velocities according to Eq. (24). The values of the axial velocities are determined as a subroutine in the resolution of the ordinary differential equation given by (Eq. (26)). However, since K appears non-linearly in (Eq. (24)), an iterative approach is required for each velocity profile. The total pressure drop (both for laminar and turbulent flows) is obtained by the numerical integration of the pressure drop along the tube length. Using Simpson’s numerical integration: −DP( =
&
L
K dz $
0
Dz [K + 4K2 + 2K3 + 4K4 + ··· 3 1 +4Kn − 1 + Kn ]
where
Fig. 3. Turbulent temperature profile (Re= 13 000).
K= −
)
)
) )
(ki + kti ) (k (T( (k + + t ri (r ri (r ri (r ri
)n
( 2T( (r 2 ri
+ (ki +kti )
)
The continuity equation:
)
dP( dz
'
(26)
( ( zi · 6¯ ri +ri (z · 6¯ r) = −ri (z · 6¯ z ) (r (z ri zi
An analytical solution, with an experimental result for the velocity near the wall is used for the axial velocities near the wall for the turbulent flow. This is a similar approach to the Prandtl near wall logarithmic law.
The energy equation is discretized as follows: (T( 1 = (z ri zicˆpi6¯ zi
(28)
6¯ zn = 12.85 (29) (27)
K ·R + 2z
'
K ·R (1/N · R) ln 0.32z 26(v/z) (2z/K · R)
Where n is the number of intervals in the radial direction.
Fig. 4. Average temperature at the tube outlet as a function of the Reynolds number.
718
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
3. Results and discussion All results are based on a case study of a tube with radius r= 0.1 m; wall temperature Tw =90°C and a temperature at the entrance Te =30°C. Since the temperature varies, it was set a reference temperature equivalent to the average between Te and Tw for the calculation of the Reynolds number. Fig. 2 shows the variation of the temperature in the axial direction for a Reynolds number of 1800, which is laminar flow. Fig. 3 shows the variation of the temperature in the axial direction for a Reynolds number of
13 000, which is turbulent flow. As expected, heat transfer is improved for the turbulent flow. It can be easily noticed that the tube centerline, for the laminar case, starts to get heated only after 20 m whereas, for the turbulent flow, it starts after about 7 m. It should be pointed out that both tubes have the same length, i.e. z =51 m, which means that the residence time for the turbulent flow is considerably smaller. Fig. 4 shows the variation of average radial temperature at the tube outlet for tube sections having the same residence time for different Reynolds numbers. It is clear that heat transfer improves as Reynolds number
Fig. 5. Temperature as a function of the radius (Re=2000).
Fig. 6. Temperature as a function of the radius (Re= 12 000).
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
increases. However, average temperature in the radial direction does not increase significantly under turbulent flow. Additionally, results show that only small elevations on the temperature is noticed for tubes longer than 20 m. The great temperature variations occur between 10 and 20 m. Figs. 5 and 6 show the variation of the temperature in some given sections as a function of the radius for laminar and turbulent flow, respectively. As expected, the turbulent profile is flatter when compared with the profile of the laminar flow. It can be noticed that it would be necessary a tube length of 2030 m to obtain a flat temperature profile (equal to the wall temperature) at the tube outlet. The turbulent flow, requires 281 m,
719
nearly a tenth of the value required for the laminar flow. Figs. 7 and 8 show the temperature profile along the axial axis for different radius for laminar and turbulent flow, respectively. In order to validate the method, it has been compared against the experimental data obtained by Haber (1966); Laufer (1952). The data obtained by Laufer is for air, which is a compressible fluid. However, the ratio between the air stream and sound velocity in the fluid (air) is low. 6 M= $ 0 c
Fig. 7. Axial temperature for several radial positions (Re= 2000).
Fig. 8. Axial temperature for several radial positions (Re=12 000).
720
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
where M is the Mach number, 6 the air velocity and c is the sound velocity at the flow temperature and pressure conditions. Since the Mach number is very small, it means that the air flow behavior is not much different from an incompressible flow. This analysis is presented in more detail by Sisson and Pitts (1988). Furthermore, the physical properties in the model presented in this work are allowed to vary with the temperature, but not with the pressure. Figs. 9 and 10 show the velocity profile as a function of the non-dimensional radius for 41 radial points and for Reynolds number of 50 000 and 500 000, respectively. The results are very close. Figs. 11 and 12 compare the axial velocity as a function of the non-dimensional radius for the experi-
mental data obtained by Haber. The fluid is water and the Reynolds numbers are 27 281 and 48 997, respectively. The comparisons are for 21 and 41 radial points. A good agreement between the model and experimental data is again obtained. Fig. 13 shows the pressure gradient as a function of the tube length. Fig. 14 compares the pressure drop as a function of the Reynolds number for the model proposed in this work and the data obtained by Friend and Metzner (1958) from empirical correlations for isothermal ideal tubes. The computational time spent for the laminar cases was about 1 s, whereas it took about 2 s for the turbulent case, using a Pentium II 400 MHz for the calculations.
Fig. 9. Comparison of the model and the experimental data by Laufer (Re= 50 000).
Fig. 10. Comparison of the model and the experimental data by Laufer (Re= 500 000).
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
721
Fig. 11. Comparison of the model and the experimental data by Haber (Re= 27 281).
Fig. 12. Comparison of the model and the experimental data by Haber (Re=48 997).
4. Conclusion A new and fast CFD method which is able to calculate conjugate heat transfer in tubes with physical properties depending on temperature is presented in this paper. The method is particularly suitable for the prediction of temperature and velocities in petrochemical fired heaters, where a considerable temperature gradient is present. Results show a very good agreement against experimental data and the pressure gradient, which can be critical in these calculations. The wall treatment is
very simple, nevertheless precise. An extension for the k–m model is under way since it is reported in the literature it is more suitable than the zero order turbulence models.
Acknowledgements The authors would like to thank CNPq and FAPESP for the grants received for this project.
722
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
Fig. 13. Pressure drop as a function of the tube length (Re= 2000 and L= 30 m).
Fig. 14. Pressure drop as a function of the Reynolds number.
References Barozzi, G. S., & Pagliarini, G. (1985). A method to solve conjugate heat transfer problems: the case of fully developed laminar flow in a pipe. Journal of Heat Transfer, 107, 77–83. Bilir, S. (1994). Laminar flow heat transfer in pipes including two-dimensional wall and fluid axial conduction. International Journal of Heat and Mass Transfer, 38 (11), 1619–1625. Bird, R. B., Stewart, W. E., Lightfoot, E. N. (1982). Transport phenomena. New York Faghri, M., & Sparrow, E. M. (1980). Simultaneous and fluid axial conduction in laminar-flow heat transfer. Journal of Heat Transfer, 102, 58 – 63. Haber, D.F., (1966). The mean velocity distribution in fully developed turbulent water flow in a straight circular pipe. University Microphilms Inc, Ann Arbor, Michigan. Laufer, J., (1952). The structure of turbulence in fully developed pipe flow. National Advisory Committee of Aeronautics, Report num-
ber 1174. National Aeronautics and Space Administration of the USA, Washington, DC, October 28. Martinuzzi, R. and Pollard, A. (1989a). Comparative study of turbulence in predicting turbulent pipe flow: algebraic stress and k-epsilon models of part I, AAIA Journal, 27(1). Martinuzzi, R. and Pollard, A. (1989b). Comparative study of turbulence models in predicting turbulent pipe flow part II: Algebraic stress and models, AIAA Journal, 27(12). Rodi, W. (1984). Turbulence modelling. In Turbulence models and their application in hydraulics. Germany. Sisson, L.E. and Pitts, D.R., (1988). Fenoˆmenos de Transporte, Editora Guanabara. Travello, J., & Dias, L. A. (1984). Temperature calculation in na incompressible permanent laminar fluid flow through a circular pipe with axial conduction and viscosity. International Journal of Heat Mass Transfer, Great Britain, 27 (6), 1183– 1187. Wimpress, R. N. (1963). Rating fired heaters. Hydrocarbon Processing and Petroleum Refiner, 42 (10), 115– 126.
Modeling flow and heat transfer in tubes using a fast CFD formulation E.R.L. Mercado, V.C. Souza, R. Guirardello, J.R. Nunhez * Faculdade de Engenharia Quı´mica, UNICAMP, CP 6066, 13083 -970 Campinas, SP, Brazil Received 10 May 2000; accepted 5 January 2001
Abstract An approach to study turbulent flow and conjugate heat transfer in tubes is proposed in this work. Instead of using the usual finite element or finite volume methods, this formulation applies a different technique that calculates both for the flow and heat transfer. It discretizes the flow in the radial direction using a fourth order finite differences method, which is more accurate than the traditional second order schemes. Using this technique, a system composed of several ordinary differential equations for the temperature and a set of linear equations for the velocities and pressure gradient is obtained. The equations are then integrated in the axial direction using a fourth order Runge–Kutta method. The values of viscosity, density and thermal conductivity are dependent on temperature, which makes the model suitable for the calculation of high temperature gradients, as in the case of refinery fired heaters. The turbulence is taken into account using a zero order turbulence model. A very interesting feature of this formulation is that for the laminar case the method is non-iterative, which makes it more desirable and faster than conventional computational fluid dynamics (CFD) packages. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Turbulent flow; Conjugate heat transfer and fluid flow; Computational fluid dynamics
Nomenclature C. p dP/dz H K kt K lm L P P( R R T T( Tw 6q Vr 6¯ r
Fluid thermal capacity Pressure drop Heat transfer coefficient Thermal conductivity Turbulent conductivity Constant= −(dP( /dz) Mixing length Tube length Pressure Average pressure (in time) Radial position Tube radius Temperature of fluid Average fluid temperature (in time) Wall temperature Angular velocity Radial velocity Average radial velocity
* Corresponding author. E-mail address: [email protected] (J.R. Nunhez). 0098-1354/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 6 7 2 - X
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
714
Vz 6¯ z Z Dr Dz v vt z
Axial velocity Average axial velocity (in time) Axial position Radial distance Axial distance between points Fluid viscosity Turbulent viscosity Fluid density
1. Introduction The heating and cooling inside tubes have been among the most important processes in the engineering field such as petrochemical fired heaters and petrol cracking. The applications are uncountable. The modeling of these processes, based on the conservation of mass, momentum and energy associated with its boundary conditions, normally lead to a set of partial differential equations with no analytical solution. Many details of the flow cannot be captured by experimental analysis, therefore, numerical procedures are needed to have a deeper understanding of these processes. Travello and Dias (1984) have developed a model to solve the energy equation in tubes under laminar flow. The fluid is incompressible and it is assumed a parabolic profile for the axial velocity. The fluid properties are constant and the work is primarily concerned to analyze axial conduction at the tube walls. Faghri and Sparrow (1980) also studied laminar flow in tubes. They were also interested in determining axial conduction at the wall of the tubes. They assume a parabolic profile for the axial velocity and the properties are considered independent on temperature. Barozzi and Pagliarini (1985) developed a method combining the finite element method with the superposition principle to solve simultaneously the momentum and energy equations in order to analyze axial conduction at the walls. The physical properties are also assumed to be constant. Bilir (1994) analyzed the conjugated heat transfer in tubes for laminar flow. He assumed that the wall temperature is constant and he also analyzed axial conduction at the walls. The physical properties again are assumed to be constant. Martinuzzi and Pollard (1989a,b) compared six turbulence models in tubes and interesting detail about the numerical methods were explained. They arrived at the conclusion that the low Reynolds number k– m model better predicts turbulence in tubes. However, physical properties are also constant. In order to get more insight of how these systems behave, a computational fluid dynamics model has been
developed to calculate the velocities, pressure and temperature. The flow can be either laminar or turbulent. Also, as the fluid properties are dependent on temperature, the model is suitable to investigate the flow inside petrochemical fired heaters. The hypotheses of the model are. The tube wall temperature is assumed to be constant, which is acceptable for the design of fired heaters (Wimpress, 1963). The flow is steady-state. It is considered that the inlet flow is already fully developed, because the flow comes from a previous tube section joined to the tube portion under consideration. This assumption is very important. Otherwise the radial velocity would not be small and the model approximations would not apply to the problem. Radial velocities are small in comparison to the axial velocities. The flow is symmetric about the axial axis (there is no variation in the angular direction and angular velocity is negligible). Physical properties are dependant on temperature, but not on pressure. There is no phase change and the fluid is a Newtonian liquid.
2. Modeling The exact model equations for tubes are well known and can be found in textbooks such as (Bird, Stewart & Lightfoot, 1982). Under the hypothesis described above, a study of the order of magnitude was applied to these equations and several terms are negligible in comparison to others. The resulting equations are described as follows, for the laminar and turbulent axisymmetric model. It is important to notice that the radial velocities can be calculated from the continuity equation, since in this equation both terms are of the same order of magnitude.
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
2.1. Laminar flow
Tube wall (r= R).
2.1.1. Momentum conser6ation Radial direction.
Vr = 0
(P =0 (r
−k
(P 1 ( (6 0= − + r ·v z (r (z r (r
(T 1 ( (T = r ·k (z r (r (r
(T = − h(Tw − T) (r
(12)
(2)
2.1.2. Energy conser6ation z · C. p · 6z
(11)
(1)
Axial direction.
715
(3)
Tube entrance (z= 0). At the tube entrance, it is considered that the flow is already fully developed. It is as if the fluid came from a previous tube section joined to the tube portion under consideration. As already mentioned, this is important to be considered since the model does not apply for situations where radial velocities are high. A constant temperature profile is assumed for the temperature since the wall temperature starts acting on the fluid only at the section being considered:
2.1.3. Continuity equation (13)
T= To ( 1 ( (z · r · 6r)+ (z · 6z ) = 0 (z r (r
(4)
=
2.2. Turbulent flow
(P( =0 (r
(5)
Axial direction.
(P( 1 ( (6¯ 0= − + r(v +vt) z (z r (r (r
(6)
z · 6z · 2 · y · r dr
(14)
(T( 1 ( (T( = r(k + kt) (z r (r (r
1 ( ( (z · r · 6¯ r)+ (z · 6¯ z ) = 0 r (r (z
In order to evaluate the turbulent viscosity and thermal conductivity the mixing length hypothesis is assumed (Prandtl theory). The mixing length theory assumes that: Turbulent viscosity:
) ) (6¯ z (r
(15)
Turbulent thermal conductivity:
kt = z · C. p · l 2m (7)
2.2.3. Continuity equation
) ) (6¯ z (r
(16)
In order to evaluate vt and kt, it is necessary to know the mixing length (lm). There are several expressions for tubes. The one chosen is by (Rodi, 1984):
(8)
lm r =0.14−0.08 R R
2
− 0.06
r R
4
(17)
2.4. Numerical calculation
2.3. Boundary conditions
Symmetry line (r= 0).
R
vt = z · l 2m
2.2.2. Energy conser6ation
&
0
2.2.1. Momentum conser6ation Radial direction.
z · C. p · 6z
Mass conservation (05 z5L).
(6z =0 (r
(9)
(T =0 (r
(10)
The fourth order finite differences method is applied to the radial direction. This leads to a linear set of equations for the velocities and a set of first order differential equations for the temperature. The discretization in the radial direction is applied to (m+1) points and there are (n+ 1) points in the axial direction. Fig. 1 shows the mesh used for the calculations.
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
716
Fig. 1. Mesh used for the calculations.
2.4.1. Laminar flow Navier –Stokes (from Eq. (1), after integrating in the r direction):
2.4.2. Turbulent flow From the continuity equations and the assumptions made, it follows that P= P(z), so the following is true:
(6z r (P = i (r ri 2vi (z ri
−
)
)
(18)
A system of m +1 equations is generated. However, there are m+2 unknowns (m + 1 velocities and pressure drop). An additional equation is obtained through a mass balance in the cross section (Eq. (14)). Integration is performed using Simpson’s rule. The linear system is solved using Gaussian elimination. Energy conservation (from Eq. (2) and isolating (#T/#z)):
)
1 (T = (z ri zi · cˆpi · 6zi
) )
)n
(k k (T ( 2T + i +ki 2 (r ri ri (r ri (r ri
(19)
The set of equations is solved by considering the calculation of the axial velocities as a subrouting in the resolution of the ordinary differential equation given by (Eq. (19)). The differential equations are solved by the Runge –Kutta –Fehlberg method, where at each required calculation the axial velocity is obtained by (Eq. (18)) by solving a set of linear equations. Continuity equation: zi · 6ri +ri
)
)
( ( (z · 6r ) = −ri (z · 6z ) (r (z ri zi
(P( =K (z
(21)
From Eq. (6); after integrating: Kr (6¯ z =− (r 2
(v + vt)
(22)
From Eq. (15). − z · f(r)
d6¯ z dr
2
+v
d6¯ z K ·r + =0 dr 2
(23)
Therefore the velocities can be calculated by the following relation, which is one root of Eq. (23)
)
(6¯ z
v 2i + 2Krizi f(ri )− vi = (r ri − 2zi f(ri ) where:
f(r)= R 2 0.14− 0.08
r R
2
(24)
n
− 0.06
r R
4 2
(20)
Applying finite differences to Eq. (20), a system of m + 1 linear equations for the radial velocity is obtained and it is solved using Gaussian elimination, after calculating the values of the axial velocity. At the tube entrance, a uniform temperature profile is assumed. From this, all physical properties are calculated. The axial velocities at the different radial positions and the pressure drop are calculated according to Eq. (18) and a set of first order differential equations is obtained according to Eq. (19), which is calculated using the fourth order Runge – Kutta method. The procedure is repeated for the whole length of the tube. After evaluating all axial velocities, the radial velocities are calculated according to Eq. (20).
Fig. 2. Laminar temperature profile (Re= 1800).
(25)
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
717
A uniform temperature profile is given at the tube entrance, similar to the linear case. The physical properties are estimated and a set of linear equations is obtained for the velocities according to Eq. (24). The values of the axial velocities are determined as a subroutine in the resolution of the ordinary differential equation given by (Eq. (26)). However, since K appears non-linearly in (Eq. (24)), an iterative approach is required for each velocity profile. The total pressure drop (both for laminar and turbulent flows) is obtained by the numerical integration of the pressure drop along the tube length. Using Simpson’s numerical integration: −DP( =
&
L
K dz $
0
Dz [K + 4K2 + 2K3 + 4K4 + ··· 3 1 +4Kn − 1 + Kn ]
where
Fig. 3. Turbulent temperature profile (Re= 13 000).
K= −
)
)
) )
(ki + kti ) (k (T( (k + + t ri (r ri (r ri (r ri
)n
( 2T( (r 2 ri
+ (ki +kti )
)
The continuity equation:
)
dP( dz
'
(26)
( ( zi · 6¯ ri +ri (z · 6¯ r) = −ri (z · 6¯ z ) (r (z ri zi
An analytical solution, with an experimental result for the velocity near the wall is used for the axial velocities near the wall for the turbulent flow. This is a similar approach to the Prandtl near wall logarithmic law.
The energy equation is discretized as follows: (T( 1 = (z ri zicˆpi6¯ zi
(28)
6¯ zn = 12.85 (29) (27)
K ·R + 2z
'
K ·R (1/N · R) ln 0.32z 26(v/z) (2z/K · R)
Where n is the number of intervals in the radial direction.
Fig. 4. Average temperature at the tube outlet as a function of the Reynolds number.
718
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
3. Results and discussion All results are based on a case study of a tube with radius r= 0.1 m; wall temperature Tw =90°C and a temperature at the entrance Te =30°C. Since the temperature varies, it was set a reference temperature equivalent to the average between Te and Tw for the calculation of the Reynolds number. Fig. 2 shows the variation of the temperature in the axial direction for a Reynolds number of 1800, which is laminar flow. Fig. 3 shows the variation of the temperature in the axial direction for a Reynolds number of
13 000, which is turbulent flow. As expected, heat transfer is improved for the turbulent flow. It can be easily noticed that the tube centerline, for the laminar case, starts to get heated only after 20 m whereas, for the turbulent flow, it starts after about 7 m. It should be pointed out that both tubes have the same length, i.e. z =51 m, which means that the residence time for the turbulent flow is considerably smaller. Fig. 4 shows the variation of average radial temperature at the tube outlet for tube sections having the same residence time for different Reynolds numbers. It is clear that heat transfer improves as Reynolds number
Fig. 5. Temperature as a function of the radius (Re=2000).
Fig. 6. Temperature as a function of the radius (Re= 12 000).
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
increases. However, average temperature in the radial direction does not increase significantly under turbulent flow. Additionally, results show that only small elevations on the temperature is noticed for tubes longer than 20 m. The great temperature variations occur between 10 and 20 m. Figs. 5 and 6 show the variation of the temperature in some given sections as a function of the radius for laminar and turbulent flow, respectively. As expected, the turbulent profile is flatter when compared with the profile of the laminar flow. It can be noticed that it would be necessary a tube length of 2030 m to obtain a flat temperature profile (equal to the wall temperature) at the tube outlet. The turbulent flow, requires 281 m,
719
nearly a tenth of the value required for the laminar flow. Figs. 7 and 8 show the temperature profile along the axial axis for different radius for laminar and turbulent flow, respectively. In order to validate the method, it has been compared against the experimental data obtained by Haber (1966); Laufer (1952). The data obtained by Laufer is for air, which is a compressible fluid. However, the ratio between the air stream and sound velocity in the fluid (air) is low. 6 M= $ 0 c
Fig. 7. Axial temperature for several radial positions (Re= 2000).
Fig. 8. Axial temperature for several radial positions (Re=12 000).
720
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
where M is the Mach number, 6 the air velocity and c is the sound velocity at the flow temperature and pressure conditions. Since the Mach number is very small, it means that the air flow behavior is not much different from an incompressible flow. This analysis is presented in more detail by Sisson and Pitts (1988). Furthermore, the physical properties in the model presented in this work are allowed to vary with the temperature, but not with the pressure. Figs. 9 and 10 show the velocity profile as a function of the non-dimensional radius for 41 radial points and for Reynolds number of 50 000 and 500 000, respectively. The results are very close. Figs. 11 and 12 compare the axial velocity as a function of the non-dimensional radius for the experi-
mental data obtained by Haber. The fluid is water and the Reynolds numbers are 27 281 and 48 997, respectively. The comparisons are for 21 and 41 radial points. A good agreement between the model and experimental data is again obtained. Fig. 13 shows the pressure gradient as a function of the tube length. Fig. 14 compares the pressure drop as a function of the Reynolds number for the model proposed in this work and the data obtained by Friend and Metzner (1958) from empirical correlations for isothermal ideal tubes. The computational time spent for the laminar cases was about 1 s, whereas it took about 2 s for the turbulent case, using a Pentium II 400 MHz for the calculations.
Fig. 9. Comparison of the model and the experimental data by Laufer (Re= 50 000).
Fig. 10. Comparison of the model and the experimental data by Laufer (Re= 500 000).
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
721
Fig. 11. Comparison of the model and the experimental data by Haber (Re= 27 281).
Fig. 12. Comparison of the model and the experimental data by Haber (Re=48 997).
4. Conclusion A new and fast CFD method which is able to calculate conjugate heat transfer in tubes with physical properties depending on temperature is presented in this paper. The method is particularly suitable for the prediction of temperature and velocities in petrochemical fired heaters, where a considerable temperature gradient is present. Results show a very good agreement against experimental data and the pressure gradient, which can be critical in these calculations. The wall treatment is
very simple, nevertheless precise. An extension for the k–m model is under way since it is reported in the literature it is more suitable than the zero order turbulence models.
Acknowledgements The authors would like to thank CNPq and FAPESP for the grants received for this project.
722
E.R.L. Mercado et al. / Computers and Chemical Engineering 25 (2001) 713–722
Fig. 13. Pressure drop as a function of the tube length (Re= 2000 and L= 30 m).
Fig. 14. Pressure drop as a function of the Reynolds number.
References Barozzi, G. S., & Pagliarini, G. (1985). A method to solve conjugate heat transfer problems: the case of fully developed laminar flow in a pipe. Journal of Heat Transfer, 107, 77–83. Bilir, S. (1994). Laminar flow heat transfer in pipes including two-dimensional wall and fluid axial conduction. International Journal of Heat and Mass Transfer, 38 (11), 1619–1625. Bird, R. B., Stewart, W. E., Lightfoot, E. N. (1982). Transport phenomena. New York Faghri, M., & Sparrow, E. M. (1980). Simultaneous and fluid axial conduction in laminar-flow heat transfer. Journal of Heat Transfer, 102, 58 – 63. Haber, D.F., (1966). The mean velocity distribution in fully developed turbulent water flow in a straight circular pipe. University Microphilms Inc, Ann Arbor, Michigan. Laufer, J., (1952). The structure of turbulence in fully developed pipe flow. National Advisory Committee of Aeronautics, Report num-
ber 1174. National Aeronautics and Space Administration of the USA, Washington, DC, October 28. Martinuzzi, R. and Pollard, A. (1989a). Comparative study of turbulence in predicting turbulent pipe flow: algebraic stress and k-epsilon models of part I, AAIA Journal, 27(1). Martinuzzi, R. and Pollard, A. (1989b). Comparative study of turbulence models in predicting turbulent pipe flow part II: Algebraic stress and models, AIAA Journal, 27(12). Rodi, W. (1984). Turbulence modelling. In Turbulence models and their application in hydraulics. Germany. Sisson, L.E. and Pitts, D.R., (1988). Fenoˆmenos de Transporte, Editora Guanabara. Travello, J., & Dias, L. A. (1984). Temperature calculation in na incompressible permanent laminar fluid flow through a circular pipe with axial conduction and viscosity. International Journal of Heat Mass Transfer, Great Britain, 27 (6), 1183– 1187. Wimpress, R. N. (1963). Rating fired heaters. Hydrocarbon Processing and Petroleum Refiner, 42 (10), 115– 126.