1
The Graetz Problem
R. Shankar Subramanian
As a good model problem, we consider steady state heat transfer to fluid in steady flow through a
tube. The fluid enters the tube at a temperature
0
T and encounters a wall temperature at
w
T ,
which can be larger or smaller than
0
T . A simple version of this problem was first analyzed by
Graetz (1883). A sketch of the system is shown below.
Objective
To obtain the steady temperature distribution ( ) , T r z in the fluid, and to calculate the rate of heat
transfer from the wall to the fluid
Assumptions
1. Steady fully developed laminar flow; steady temperature field.
2. Constant physical properties , , ,
p
k C -- This assumption also implies incompressible
Newtonian flow.
3. Axisymmetric temperature field 0
T
.
4. Negligible viscous dissipation
Fluid at
0
T
r
z
R
( )
,
w
T R z T =
2
Velocity Field
Poiseuille Flow
0; 0
r
v v
= =
( )
2
0 2
1
z
r
v r v
R
| |
=
|
\ .
0
: v Maximum velocity existing at the centerline
Energy Equation
Subject to assumption (2), Equation (B.9.2) from Bird et al. (page 850) can be written as follows.
1 0 0
r
p
v v
T
C
t
= =
r
v +
v
T
r
+
2
2 2
3 4
1 1
z
T T
v
r z
T T
k r
r r r r
(
+ (
(
| |
= +
|
\ .
2
2 v
T
z
(
+ +
(
(
and therefore, simplified to
2 2
0 2 2
1
1
r T T T
v r
R z r r r z
| | ( | |
= +
| | (
\ .
\ .
where
( )
/
p
k C = is the thermal diffusivity of the fluid.
Boundary Conditions
Inlet: ( )
0
, 0 T r T =
Wall: ( ) ,
w
T R z T =
Centerline: ( ) 0, T z is finite
or ( ) 0, 0
T
z
r
3
Because of the appearance of the axial conduction term in the governing differential
equation, we should write another boundary condition in the z coordinate. But actually, the
inlet condition written above is incompatible with the inclusion of axial conduction in the
problem, because conduction will lead to some of the information about the step change in wall
temperature at the inlet to propagate backward. As we shall see shortly, well neglect axial
conduction, which will obviate the need for writing a second condition in the z coordinate.
Non-Dimensionalization
We shall use the following scheme for scaling (or non-dimensionalizing) the variables.
0
w
w
T T
T T
,
r
Y
R
= ,
z
Z
R Pe
= , where the Pclet Number
0
Rv
Pe
= .
This permits us to transform the governing differential equation and boundary conditions to the
following form.
( )
2
2
2 2
1 1
1 Y Y
Z Y Y Y Pe Z
| |
= +
|
\ .
( ) , 0 1 Y =
( ) 1, 0 Z =
( ) 0, Z is finite
or ( ) 0, 0 Z
Y
=
The Pclet Number
The Pclet number plays the same role in heat transport as the Reynolds number does in fluid
mechanics. First, we note that the Pclet number is the product of the Reynolds and Prandtl
numbers.
0 0
Re Pr
Rv Rv
Pe
= = =
The physical significance of the Pclet number can be inferred by recasting it slightly.
0
Rate of energy transport by convection
Rate of energy transport by conduction
p
v C T
Pe
T
k
R
= =
4
Note that the numerator represents the order of magnitude of the convective flux in the main
flow direction, whereas the denominator stands for the order of magnitude of the conduction flux
in the radial direction. If we wish to compare the rates of energy transport by these two
mechanisms in the same direction, we can multiply the Pclet number by / L R where L is a
characteristic length in the axial direction.
For large values of Pe , we can see that
2
1
1
Pe
. Therefore, in the scaled energy equation, the
term involving axial conduction can be safely neglected. Physically, there are two mechanisms
for transporting energy in the axial direction, namely, convection and conduction. Because the
Pclet number is large, we are able to neglect transport by conduction in comparison with
transport by convection. On the other hand, in the radial direction, there is only a single
mechanism for transport of energy, namely conduction. By performing calculations including
conduction in the axial direction, it has been established that it is safe to neglect axial conduction
for 100. Pe To learn about how to include axial conduction, you can consult the articles by
Davis (1973), Acrivos (1980), and Papoutsakis et al. (1980).
Let us make a sample calculation of the Pclet number for laminar flow heat transfer in a tube.
The thermal diffusivity of common liquids is typically in the range
2
7 7
10 2 10
m
s
, and well
use the larger limit. Choose
2
7
0
10 , 0.05 , 2 10
m m
R mm v
s s
= = =
This yields, 2, 500 Pe = , which is much larger than 100. We can check to see if the flow is
laminar by calculating the Reynolds number. If the fluid is water,
2
6
10
m
s
, which yields a
Prandtl number Pr 5
= = . Therefore, the Reynolds number is Re 500 = , which is comfortably
in the laminar flow regime.
The final version of the scaled energy equation is
( )
2
1
1 Y Y
Z Y Y Y
| |
=
|
\ .
We can solve this equation by separation of variables, because the boundary conditions in the
Y coordinate are homogeneous. The method of separation of variables yields an infinite series
solution for the scaled temperature field.
( ) ( )
2
1
,
n
Z
n n
n
Y Z A e Y
=
=
5
In the above solution, the functions ( )
n
Y are the characteristic functions or eigenfunctions of
a proper Sturm-Liouville system.
( )
2 2
1
1 0
d d
Y Y
Y dY dY
| |
+ =
|
\ .
( ) ( ) 0 0 or 0
d
dY
= is finite
( ) 1 0 =
The above ordinary differential equation for ( ) Y can be solved by applying the following
transformations to both the dependent and the independent variables (Lauwerier, 1951, Davis,
1973).
( ) ( )
2
2
X
X Y W X e Y = =
This leads to the following differential equation for ( ) W X .
( )
2
2
1
1 0
4 2
d W dW
X X W
dX dX
| |
+ + =
|
\ .
This is known as Kummers equation. It has two linearly independent solutions, but only one is
bounded at 0. X = Because ( ) 0 must be bounded, we must require that ( ) 0 W also remain
bounded. This rules out the singular solution, leaving us with the regular solution
( )
1
, 1,
2 4
W X c M X
| |
=
|
\ .
where c is an arbitrary multiplicative constant. The function ( ) , , M a b X is the confluent
hypergeometric function, or Kummer function, and is discussed in Chapter 13 of the Handbook
of Mathematical Functions by M. Abramowitz and I. A. Stegun, It is an extension of the
exponential function, and is written in the form of the following series.
( )
( )
( )
( ) ( )
( ) ( )
2
1
, , 1
1 2!
1 1
1 1 !
n
a a
a X
M a b X X
b b b
a a a n
X
b b b n n
+
= + + +
+
+ +
+ +
+ +
You can see that when a b = ,
( ) , ,
X
M a a X e =
6
Application of the boundary condition at the tube wall, ( ) 1 0 = , leads to the following
transcendental equation for the eigenvalues.
1
, 1, 0
2 4
M
| |
=
|
\ .
The above equation has infinitely many discrete solutions for , which we designate as
n
, with
n assuming positive integer values beginning from 1. Corresponding to each value
n
, there is
an eigenfunction ( )
n
Y given by
( ) ( )
2
2
2
n
Y
n n n
Y e W Y
=
The first few eigenvalues are reported in the table.
n
n
2
n
1 2.7044 7.3136
2 6.6790 44.609
3 10.673 113.92
4 14.671 215.24
5 18.670 348.57
Note that technically
{ }
2
n
is the set of eigenvalues, even though we use the term loosely to
designate { }
n
as that set for convenience.
The most important property of a proper Sturm-Liouville system is that the eigenfunctions are
orthogonal with respect to a weighting function that is specific to that system. In the present
case, the orthogonality property of the eigenfunctions can be stated as follows.
( ) ( ) ( )
1
2
0
1 0,
m n
Y Y Y Y dY m n =
Using this orthogonality property, it is possible to obtain a result for the coefficients in the
solution by separation variables.
( ) ( )
( ) ( )
1
2
0
1
2 2
0
1
1
n
n
n
Y Y Y dY
A
Y Y Y dY
7
The Heat Transfer Coefficient
The heat flux from the wall to the fluid, ( )
w
q z is a function of axial position. It can be
calculated directly by using the result
( ) ( ) ,
w
T
q z k R z
r
but as we noted earlier, it is customary to define a heat transfer coefficient ( ) h z via
( ) ( )( )
w w b
q z h z T T =
where the bulk or cup-mixing average temperature
b
T is introduced. The way to experimentally
determine the bulk average temperature is to collect the fluid coming out of the system at a given
axial location, mix it completely, and measure its temperature. The mathematical definition of
the bulk average temperature was given in an earlier section.
( )
0
0
2 ( ) ,
2 ( )
R
b R
rV r T r z dr
T
rV r dr
where the velocity field ( ) ( )
2 2
0
1 / V r v r R = . You can see from the definition of the heat
transfer coefficient that it is related to the temperature gradient at the tube wall in a simple
manner.
( )
( )
( )
,
w b
T
k R z
r
h z
T T
We can define a dimensionless heat transfer coefficient, which is known as the Nusselt number.
( )
( )
( )
1,
2
2
b
Z
hR
Y
Nu Z
k Z
= =
where
b
is the dimensionless bulk average temperature.
8
By substituting from the infinite series solution for both the numerator and the denominator, the
Nusselt number can be written as follows.
( )
( )
( ) ( )
2
2
1
1
2
1
0
1
2
4 1
n
n
Z
n
n
n
Z
n n
n
d
A e
dY
Nu Z
A e Y Y Y dY
=
=
The denominator can be simplified by using the governing differential equation for ( )
n
Y , along
with the boundary conditions, to finally yield the following result.
( )
( )
2
2
1
2
1
1
2 1
n
n
Z
n
n
n
Z
n
n
n n
d
A e
dY
Nu
d e
A
dY
=
=
=
We can see that for large Z , only the first term in the infinite series in the numerator, and
likewise the first term in the infinite series in the denominator, is important. Therefore, as
Z ,
2
1
3.656
2
Nu
= .
The sketch qualitatively illustrates the behavior of the Nusselt number as a function of
dimensionless axial position.
3.656
Nu
Dimensionless Axial Position Z
0
0
9
A similar analysis is possible in the case of a uniform wall flux boundary condition.
Extensions of the Graetz solution by separation of variables have been made in a variety of ways,
accommodating non-Newtonian flow, turbulent flow, and other geometries besides a circular
tube.
The Lvque Approximation
The orthogonal function expansion solution obtained above is convergent at all values of the
axial position, but convergence is very slow at the inlet is approached. The main reason for this
is the assistance provided by
2
n
Z
e
in accelerating convergence for sufficiently large values of
Z . Lvque (1928) considered the thermal entrance region in a tube and developed an
alternative solution, which is useful precisely where the orthogonal function expansion
converges too slowly.
We shall now construct the Lvque solution which is built on the assumption that the thickness
of the thermal boundary layer
t
R . This assumption leads to the following simplifications.
1. Curvature effects can be neglected in the radial conduction term. This means that the
derivative
1 T
r
r r r
| |
|
\ .
can be approximated by
2
2
1 T T
R
R r r r
| |
=
|
\ .
.
2. Because we are only interested in the velocity distribution within the thermal boundary layer,
we expand the velocity field in a Taylor series in distance measured from the tube wall and retain
the first non-zero term.
Defining x R r = , we can rewrite the velocity distribution as
( )
( )
2
2
0 0 0 2 2
1 2 2
z
R x
x x x
v r v v v
R R R R
| |
| |
= = |
|
|
\ .
\ .
Fluid at
0
T
r
z
R
( )
,
w
T R z T =
t
10
Recall that a power series obtained by any method is a Taylor series. The above approach is
simpler than working out the derivatives of ( )
z
v r in the x coordinate, evaluating them at the
wall, and constructing the Taylor series.
3. Because the conditions outside the thermal boundary layer are those in the fluid entering the
tube, we shall use the boundary condition ( )
0
T x T instead of the centerline boundary
condition employed in obtaining the Graetz solution.
Beginning with the simplified energy equation in which axial conduction has been neglected
already, and invoking the above assumptions, we have the following governing equation for the
temperature field.
2
0 2
2
x T T
v
R z x
=
where the chain rule has been used to transform the second derivative in r to the second
derivative in x .
The temperature field ( ) , T x z satisfies the following boundary conditions.
( )
0
, 0 T x T =
( ) 0,
w
T z T =
( )
0
, T z T =
We shall work with a dimensionless version of these equations. For consistency, we scale the
temperature and axial coordinate in the same manner as before.
0
w
w
T T
T T
z
Z
R Pe
=
We define a new scaled distance from the wall via / X x R = . The scaled governing equation
and boundary conditions are given below.
2
2
2X
Z X
=
( ) , 0 1 X =
( ) 0, 0 Z =
( ) , 1 Z =
11
The similarity of this governing equation and boundary conditions to those in the fluid
mechanical problem in which we solved for the velocity distribution between two plates when
one of them is held fixed and the other is moved suddenly is not a coincidence. For small values
of time in the fluid mechanical problem, we replaced the boundary condition at the top plate with
one at an infinite distance from the suddenly moved plate, and used the method of combination
of variables to solve the equations. It would be worthwhile for you to go back and review the
notes on combination of variables at this stage.
By invoking ideas very similar to those used in the fluid mechanical problem, we postulate that a
similarity solution exists for the temperature field in the present problem. That is, we assume
( ) ( ) , X Z F = where the similarity variable ( ) / X Z = . The variable ( ) Z represents the
scaled thermal boundary layer thickness, and is unknown at this stage. We make the necessary
transformations using the chain rule.
2
dF X d dF d dF
Z Z d dZ d dZ d
| |
= = =
|
\ .
1 dF dF
X X d d
= =
( )
2 2
2 2 2
1 1 1 1 dF dF d dF d F
X X Z d X d X d d d
(
( (
= = = =
(
( (
Using these results, the partial differential equation for ( ) , X Z is transformed to an ordinary
differential equation for ( ) F .
2
2 2
2
2 0
d F d dF
d dZ d
| |
+ =
|
\ .
It is evident that the similarity hypothesis will fail unless the quantity inside the parentheses is
required to be independent of Z , and therefore, a constant. For convenience, we set this
constant to 3/2. Therefore, we have an ordinary differential equation for ( ) F and another for
( ) Z .
2
2
2
3 0
d F dF
d d
+ =
2
3
2
d
dZ
=
12
To derive the boundary conditions on these functions, we must go to the boundary conditions on
( ) , X Z . In a straightforward way, we see that ( ) 0, 0 Z = yields ( ) 0 0 F = , and ( ) , 1 Z =
leads to ( ) 1 F = . The remaining (inlet) condition gives
( )
( )
, 0 1
0
X
X F
| |
= =
|
|
\ .
By choosing ( ) 0 0 = , this condition collapses into the condition ( ) 1 F = obtained already
from the boundary condition on the scaled temperature field as X . Summarizing the
boundary conditions on ( ) F and ( ) Z , we have
( ) ( ) 0 0, 1 F F = = , and
( ) 0 0 =
Integration yields the following solution for the scaled boundary layer thickness ( ) Z .
( )
1/ 3
9
2
Z Z
| |
=
|
\ .
The solution for ( ) F is
( )
( )
3
3
3
0
0
0
1
4/ 3
e d
F e d
e d
= =
Here, ( ) x represents the Gamma function, discussed in the Handbook of Mathematical
Functions by Abramowitz and Stegun. The numerical value of ( ) 4/ 3 1.120 .
Heat Transfer Coefficient
In the thermal entrance region, when the thermal boundary layer is thin, we can approximate the
bulk average temperature
b
T by the temperature of the fluid entering the tube
0
T . Therefore, we
define the heat transfer coefficient in this entrance region by
( ) ( )
0
,
w w
T
q k R z h T T
r
= =
13
Transforming to dimensionless variables, and defining a Nusselt number 2 / Nu hR k = , we
can write
( ) ( )
( )
( )
2
2 0, 0
dF
Nu Z Z
X Z d
= =
By substituting for ( ) Z and ( ) 0
dF
d
, we obtain the following approximate result for the
Nusselt number in the thermal entrance region.
( )
1/ 3
1/ 3
1.357
R
Nu Z Pe
z
| |
|
\ .
Comparison with the exact solution shows this result is a good approximation in the range
2500 50
Pe z Pe
R
| |
|
\ .
References
M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York
(1965).
A. Acrivos, The extended Graetz problem at low Pclet numbers, Appl. Sci. Res. 36, 35 (1980).
E.J. Davis, Can. J. Chem. Engg. 51, 562 (1973).
L. Graetz, Ueber die Wrmeleitungsfhigkeit von Flssigkeiten, Annalen der Physik und Chemie
18, 79 (1883).
H.A. Lauwerier, Appl. Sci. Res. A2, 184 (1951).
M.A. Lvque, Les lois de la transmission de chaleur par convection, Annales des Mines,
Memoires, Series 12, 13, 201-299, 305-362, 381-415 (1928) {as cited by J. Newman, Trans.
ASME J. Heat Transfer, 91, 177 (1969)}
E. Papoutsakis, D. Ramkrishna, and H.C. Lim, The extended Graetz problem with Dirichlet wall
boundary conditions, Appl. Sci. Res. 36, 13 (1980).
