Exact Graetz problem solution by using hypergeometric function

Exact Graetz problem solution by using hypergeometric function

Ali Belhocine*Wan Z.W. Omar 

Faculty of Mechanical Engineering, University of Sciences and the Technology of Oran, L.P 1505 El-MNAOUER, USTO 31000 ORAN Algeria

Faculty of Mechanical Engineering, University Teknologi Malaysia, 81310 UTM Skudai, Alaysia

Corresponding Author Email: 
30 June 2017
| Citation



This paper proposes an exact solution of the classical Graetz problem in terms of an infinite series represented by a nonlinear partial differential equation considering two space variables, two boundary conditions and one initial condition. The mathematical derivation is based on the method of separation of variables whose several stages were illustrated to reach the solution of the Graetz problem.A MATLAB code was used to compute the eigenvalues of the differential equation as well as the coefficient series. In addition, the analytical solution was compared to the numerical values obtained previously by Shah and London. It is important to note that the analytical solution is in good agreement with published numerical data.


graetz problem, sturm-liouville problem, hypergeometric function, heat transfer

1. Introduction

The solutions of one or more partial differential equations (PDEs), which are subjected to relatively simple limits, can be tackled either by analytical or numerical approach. There are two common techniques available to solve PDEs analytically, namely the variable separation and combination of variables.

A heat exchanger is a device used for transferring heat from one fluid to another. The fluid may not be allowed to mix by separating them by a solid wall or they may be in direct contact. They are operated in numerous industries such as power generation, petroleum refineries, chemical and processing plants and HVACs [1]. The heat transfer process associated with natural convection is extensively involved in numerous engineering applications due to its diverse applications in geophysics, nuclear reactor system, energy efficient buildings, cooling of electronics system, solar system etc [2]. Convection is one of the heat transfer modes in addition to conduction and radiation, this transfer type can be arises between solid and flowing fluid. Heat transfer convection is divided into: forced convection caused by external forces like pumps and fans, and natural (free) convection when the motion is due only to the temperature difference between the wall and the fluid [3].

The Graetz problem describes the temperature (or concentration) field in fully developed laminar flow in a circular tube where the wall temperature (or concentration) profile is a step-function [4]. The simple version of the Graetz problem was initially neglecting axial diffusion, considering simple wall heating conditions (isothermal and isoflux), using simple geometric cross-section (either parallel plates or circular channels), and also neglecting fluid flow heating effects, which can be generally denoted as Classical Graetz Problem [5]. Min et al. [6] presented an exact solution for a Graetz problem with axial diffusion and flow heating effects in a semi-infinite domain with a given inlet condition. Later, the Graetz series solution was further improved by Brown [7]. Ebadian and Zhang [8] analyzed the convective heat transfer properties of a hydrodynamically, fully developed viscous flow in a circular tube. Lahjomri and Oubarra [9] investigated a new method of analysis and an improved solution for the extended Graetz problem of heat transfer in a conduit. An extensive list of contributions related to this problem may be found in the papers of Papoutsakis et al.[10] and Liou and Wang [11]. In addition, the analytical solution proposed efficiently resolves the singularity and this methodology allows extension to other problems such as the Hartmann flow [12], conjugated problems [13] and other boundary conditions. Recently, Belhocine [14] developed a mathematical model to solve the classic problem of Graetz using two numerical approaches, the orthogonal collocation method and the method of Crank-Nicholson.

In this paper, the Graetz problem that consists of two differential partial equations will be solved using separation of variables method. The Kummer equation is employed to identify the confluent hypergeometric functions and its properties in order to determine the eigenvalues of the infinite series which appears in the proposed analytical solution. In addition, the exact analytical solution presented in this work was validated by the numerical value data previously published by Shah and London.

2.The Heat Equation in Cylindrical Coordinates

The general equation for heat transfer in cylindrical coordinates developed by Bird, Stewart and Lightfoot [15] is as follows;

$u_{Z} \frac{\partial T}{\partial z}=\frac{k}{\rho C_{p}} \nabla^{2} T$      (1)

$\rho C_{p}\left(\frac{\partial T}{\partial t}+u_{r} \frac{\partial T}{\partial r}+\frac{u_{\theta}}{r} \frac{\partial T}{\partial \theta}+u_{z} \frac{\partial T}{\partial z}\right)=k\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} T}{\partial \theta^{2}}+\frac{\partial^{2} T}{\partial z^{2}}\right]$

$+2 \mu\left\{\left(\frac{\partial u_{r}}{\partial r}\right)^{2}+\left[\frac{1}{r}\left(\frac{\partial u_{\theta}}{\partial \theta}+u_{r}\right)\right]^{2}+\left(\frac{\partial u_{z}}{\partial z}\right)^{2}\right\}+\mu\left\{\left(\frac{\partial u_{\theta}}{\partial z}+\frac{1}{r} \frac{\partial u_{z}}{\partial \theta}\right)^{2}+\left(\frac{\partial u_{z}}{\partial r}+\frac{\partial u_{r}}{\partial z}\right)^{2}\right\}+\mu\left[\frac{1}{r} \frac{\partial u_{r}}{\partial \theta}+r \frac{\partial}{\partial r}\left(\frac{u_{\theta}}{r}\right)\right]^{2}$       (2)

Considering that the flow is steady, laminar and fully developed flow (Re < 2400), and if the thermal equilibrium had already been established in the flow, then $\frac{\partial T}{\partial t}=0$ . The dissipation of energy would also be negligible. Other physical properties would also be constant and would not vary with temperature such as ρ, µ, Cp, k

This assumption also implies incompressible Newtonian flow.

Axisymmetric temperature field $\left(\frac{\partial T}{\partial \theta}=0\right)$ , where we are using the symbol θ for the polar angle.

By applying the above assumptions, Equation (2) can be written as follows:

$u_{Z} \frac{\partial T}{\partial z}=\frac{k}{\rho C_{p}}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r}\right)\right]$        (3)

Given that the flow is fully developed laminar flow (Poiseuille flow), then the velocity profile would have followed the parabolic distribution across the pipe section, represented by

$u_{Z}=2 \overline{u}\left[1-\left(\frac{r}{R}\right)^{2}\right]$          (4)

where 2$\overline{u}$ is the Maximum velocity existing at the centerline

By replacing the speed term in Equation (3), we get:

$2 \overline{u}\left[1-\left(\frac{r}{R}\right)^{2}\right] \frac{\partial T}{\partial z}=\frac{k}{\rho C_{p}}\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r}\right)\right]$        (5)

Solving the equations requires the boundary conditions as set in Figure 1, thus

$\left.\begin{array}{ll}@{z=0} & {T=T_{0}} \\ @{r=0} & {\frac{\partial T}{\partial r}=0} \\ @{r=R} & {T=T_{\omega}}\end{array}\right\}$      (6)

It is more practical to study the problem with standardized variables from 0 to 1. To do this, new variables without dimension (known as adimensional) are introduced, defined as $\theta=\frac{T_{\omega}-T}{T_{\omega}-T_{0}}$ , $x=\frac{r}{R}$ and $y=\frac{z}{L}$ . The substitution of the adimensional variables in Equation (5) gives

$2 \overline{u}\left[1-\frac{x^{2} R^{2}}{R^{2}}\right] \frac{\left(T_{0}-T_{w}\right)}{L} \frac{\partial \theta}{\partial y}=\frac{k}{\rho c_{p}}\left[\frac{1}{x R} \frac{\left(T_{0}-T_{w}\right)}{R} \frac{\partial \theta}{\partial x}+\frac{\left(T_{0}-T_{w}\right)}{R^{2}} \frac{\partial^{2} \theta}{\partial x^{2}}\right]$    (7)

After making the necessary arrangements and simplifications, the following simplified equation is obtained.

$\left(1-x^{2}\right) \frac{\partial \theta}{\partial y}=\frac{k L}{\rho C_{P} 2 \overline{u} R^{2}}\left[\begin{array}{cc}{\frac{1}{x}} & {\frac{\partial \theta}{\partial x}+\frac{\partial^{2} \theta}{\partial x^{2}}} \\ {} & {}\end{array}\right]$        (8)

where the term $\frac{2 \overline{u} R \rho C_{p}}{k}$ is the dimensional number known as the Peclet number (Pe), which in fact is the Reynolds number divided by the Prandtl number. In steady state condition, the partial differential equation resulting from this, in the adimensional form can be written as follows:

$\left(1-x^{2}\right) \frac{\partial \theta}{\partial y}=\frac{L}{R P e}\left[\frac{1}{x} \frac{\partial}{\partial x}\left(x \frac{\partial \theta}{\partial x}\right)\right]$      (9)

This equation, if subjected to the new boundary conditions, would be transformed to the followings:

$\mathrm@{z}=0, T=T_{0} \Rightarrow(a \mathrm{y}=0, \theta=1$

$\mathrm@{r}=0, \frac{\partial T}{\partial r}=0 \Rightarrow \widehat{\mathrm{Q}} \mathrm{x}=0, \frac{\partial \theta}{\partial x}=0 \Rightarrow \theta(0, y) \neq 0$

$@\mathrm{r}=\mathrm{R}, T=T_{\omega} \Rightarrow \widehat{\alpha} \mathrm{x}=1, \theta=0 \Rightarrow \theta(1, y)=0$

It is hereby proposed that the separation of variables method could be applied, to solve Equation (9).

3. Analytical Solution Using Separation of Variables Method

As a good model problem, we consider the steady state heat transfer of fluid in a fully developed laminar flow through a circular pipe. The fluid enters at z=0 at a temperature of T0 and the pipe walls are maintained at a constant temperature of Tω.We will write the differential equation for the temperature distribution as a function of r and z , and then express this in a dimensionless form and identify the important dimensionless parameters. Heat generation in the pipe due to the viscous dissipation is neglected, and a Newtonian fluid is assumed. Also, we neglect the changes in viscosity in the temperature variation. A sketch of the system is shown below.


Figure 1. Schematics of the classical Graetz problem and the coordinate system

In both qualitative and numerical methods, the dependence of solutions on the parameters plays an important role, and there are always more difficulties when there are more parameters. We describe a technique that changes variables so that the new variables are “dimensionless”. This technique will lead to a simple form of the equation with fewer parameters. Let the Graetz problem is given by the following governing equation

$\left.1-x^{2}\right) \frac{\partial \theta}{\partial y}=\frac{L}{P e R}\left[\frac{1}{x} \frac{\partial \theta}{\partial x}+\frac{\partial^{2} \theta}{\partial x^{2}}\right]$        (10)

where the initial conditions are as follows;

$\begin{array}{ll}{\mathrm{IC} :} & {y=0, \quad \theta=1} \\ {\mathrm{BCl} :} & {x=0 \quad, \quad \frac{\partial \theta}{\partial x}=0} \\ {\mathrm{BC} 2 :} & {x=1, \quad \theta=0}\end{array}$

Introducing dimensionless variables [16], as follows:

$x=\xi=\frac{r}{R}$      (11)

$\zeta=\frac{k z}{\rho c_{p} v_{\max } R^{2}}$    (12)

$y=\frac{z}{L}$      (13)

By substituting Equation (13) into Equation (12) then it becomes:

$\zeta=\frac{k y L}{\rho c_{p} v_{\max } R^{2}}$       (14)

Knowing that $v_{\max }=2 \overline{u}$


$\zeta=\frac{k y L}{2 \overline{u} \rho c_{p} R^{2}}=\frac{y L}{\frac{2 \overline{u} \rho c_{p} R}{k} R}$       (15)

Notice that the term $\frac{2 \overline{u} \rho c_{p} R}{k}$  in Equation (15) is similar to the Peclet number, P.

Thus, Equation (15) can be written as

$\zeta=\frac{y L}{P e R}$     (16)

Based on Equations (11)-(16), ones can write the following expressions;

$\frac{\partial \theta}{\partial x}=\frac{\partial \theta}{\partial \xi}$    (17)

$\frac{\partial^{2} \theta}{\partial x^{2}}=\frac{\partial^{2} \theta}{\partial \xi^{2}}$      (18)

$\frac{\partial \theta}{\partial y}=\frac{\partial \zeta}{\partial y} \cdot \frac{\partial \theta}{\partial \zeta}=\frac{L}{P e R} \frac{\partial \theta}{\partial \zeta}$       (19)

Now, by replacing Equations (17)-(19) into Equation (10) the governing equation becomes:

$\frac{L}{P e R}\left(1-\xi^{2}\right) \frac{\partial \theta}{\partial \zeta}=\frac{L}{P e R}\left[\frac{1}{\xi} \frac{\partial \theta}{\partial \xi}+\frac{\partial^{2} \theta}{\partial \xi^{2}}\right]$       (20)

by eliminating $\frac{L}{P e R}$ the term, Equation (20) will be reduced to:

$\left(1-\xi^{2}\right) \frac{\partial \theta}{\partial \zeta}=\frac{1}{\xi} \frac{\partial \theta}{\partial \xi}+\frac{\partial^{2} \theta}{\partial \xi^{2}}$       (21)

The right term in Equation (21) can be simplified as follows:

$\frac{1}{\xi} \frac{\partial \theta}{\partial \xi}+\frac{\partial^{2} \theta}{\partial \xi^{2}}=\frac{1}{\xi} \frac{\partial}{\partial \xi}\left(\xi \frac{\partial \theta}{\partial \xi}\right)$      (22)

Finally, the equation that characterizes the Graetz problem can be written in the form of:  

$\left(1-\xi^{2}\right) \frac{\partial \theta}{\partial \zeta}=\frac{1}{\xi} \frac{\partial}{\partial \xi}\left(\xi \frac{\partial \theta}{\partial \xi}\right)$       (23)

Now, using an energy balance method in the cylindrical coordinates, Equation (23) can be decomposed into two ordinary differential equations. This is done by assuming constant physical properties of a fluid and neglecting axial conduction and in steady state. By imposing initial conditions as given below:

@at $\zeta=0, \quad \theta=1$

@at $\quad \xi=1, \quad \theta=0$

@at $\quad \xi=0, \quad \frac{\partial \theta}{\partial \xi}=0$

and dimensionless variables are defined by:

$\theta=\frac{T_{\omega}-T}{T_{\omega}-T_{0}}$  , $\xi=\frac{r}{r_{1}}$ and $\zeta=\frac{k z}{\rho c_{p} v_{\max } r_{1}^{2}}$

while the separation of variables method is given by

$\theta=Z(\zeta) R(\xi)$    (24)

Finally, Equation (23) can be expressed as follows:

$\frac{d Z}{Z}=-\beta^{2} \xi$       (25)


$\xi \frac{d^{2} R}{d \xi^{2}}+\frac{d R}{d \xi}+\beta^{2} \xi\left(1-\xi^{2}\right) R=0$      (26)

where $\beta^{2}$ is a positive real number and represents the intrinsic value of the system.

The solution of Equation (25) can be given as:

$Z=c_{1} e^{-\beta^{2} \zeta}$       (27)

where $\mathcal{C}_{1}$ is an arbitrary constant. In order to solve Equation (26), transformations of dependent and independent variables need to be made by taking:

(Ι) $v=\beta \xi^{2}$

(Π) $R(v)=e^{-v / 2} S(v)$ 

Thus, Equation (17) is now given by;

$v \frac{d^{2} S}{d v^{2}}+(1-v) \frac{d S}{d v}-\left(\frac{1}{2}-\frac{\beta}{4}\right) S=0$       (28)

Equation (19) is also called as confluent hypergeometric [17] and it is commonly known as the Kummer equation.

A homogeneous linear differential equation of the second order is given by

$y^{\prime \prime}+P(Z) y^{\prime}+Q(Z) y=0$       (29)

If P(Z) and Q(Z) admit a pole at point Z=Z0, it is possible to find a solution developed in the whole series provided that the limits on and exist.

The method of Frobenius seeks a solution in the form of

$y(Z)=Z^{\lambda} \sum_{n=0}^{\infty} a_{n} Z^{n}$    (30)

where, $\lambda$ is a coefficient to be determined whilst properties of the hypergeometric functions are defined by;

$F(\alpha, \beta, Z)=1+\frac{\alpha}{\beta} Z+\frac{\alpha(\alpha+1)}{\beta(\beta+1)} \frac{Z^{2}}{2 !}+\cdots+\frac{\alpha(\alpha+1)(\alpha+2) \cdots(\alpha+n-1)}{\beta(\beta+1)(\beta+2) \cdots(\theta+n-1)} \frac{Z^{n}}{n !}+\cdots+$          (31)

Using derivation against Z, the function is now become

$\frac{d}{d Z}[F(\alpha, \beta, Z)]=\frac{\alpha}{\beta}\left[1+\frac{(\alpha+1)}{(\beta+1)} Z+\frac{(\alpha+1)(\alpha+2)}{(\beta+1)(\beta+2)} \frac{Z^{2}}{2 !}\right. \left.+\cdots+\frac{\alpha(\alpha+1)(\alpha+2) \cdots(\alpha+n-1)}{\beta(\beta+1)(\beta+2) \cdots(\alpha+n-1)} \frac{Z^{n}}{n !}+\cdots+-\right]=\frac{\alpha}{\beta} F(\alpha+1, \beta+1, Z)$          (32)

From Equation (32), ones will get;

$\frac{d}{d \xi} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)=\left(-\beta_{n}^{2}\right) \xi\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n} \xi^{2}\right)$              (33)

Thus, the solution of Equation (26) can be obtained by:

$R=c_{2} e^{-\beta \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta}{4}, 1, \beta \xi^{2}\right)$    (34)

$F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n}\right)=0$        (35)

Where n = 1, 2, 3, ... and eigenvalues $\beta_{n}$ are the roots of Equation (35). Since the system is linear, the general solution can be determined using superposition approach:

$\theta=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \xi} e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)$     (36)

The constants in Equation (36) can be sought using orthogonality property of the Sturm-Liouville systems after the initial condition is being applied as stated below;

$C_{n}=\frac{\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)}{\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2}}\left[F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)\right]^{2}}$      (37)

The integral in the denominator of Equation (37) can be evaluated using numerical integration.

For the Graetz problem, it is noticed that;

$\left(\xi-\xi^{3}\right) \frac{\partial \theta}{\partial \zeta}=\frac{\partial}{\partial \xi}\left(\xi \frac{\partial \theta}{\partial \xi}\right)$        (38)

where is the function of the weight / $\beta_{n}$ eigenvalues

B.C  $\xi=1, \theta=0$

B.C  $\xi=0, \theta=1$

IC    $\zeta=0, \theta=1$

$\theta=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta} G_{n}(\xi)=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta} e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) G_{n}(\xi)=e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)$                            (39)

is the function of the weight

Sturm-Liouville problem.

$\frac{1}{\xi} \frac{d}{d \xi}\left(\xi \frac{d G_{n}}{d \xi}\right)-\left(1-\xi^{2}\right) \beta_{n}^{2} G_{n}=0$     (40)

$\frac{d G_{n}}{d \xi}=0$ for $\xi=0, G_{n}=0$ for $\xi=1$       (41)

IC $\quad \zeta=0, \theta=1$

$\theta(\zeta=0)=1=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)$     (42)

Relation of orthogonality

$\int_{0}^{1} W(x) Y_{i}(x) Y_{j}(x)=0,(i \neq j)$        (43)

$\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) e^{-\beta_{m} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{m}}{4}, 1, \beta_{m} \xi^{2}\right) d \xi=0$       (44)

$\frac{\partial \theta}{\partial \xi}=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta}\left(-\beta_{n} \xi\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)+\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta}\left(-\beta_{n}^{2}\right) \xi\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n} \xi^{2}\right)$   (45)

$\frac{\partial \theta}{\partial \zeta}=\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta}\left(-\beta_{n}^{2}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)$       (46)

By considering $F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n}\right)$ Equation (38) and onwards can be given as;

$\int_{0}^{1}\left(\xi-\xi^{3}\right) \frac{\partial \theta}{\partial \zeta} d \xi=\int_{0}^{1} \frac{\partial}{\partial \xi}\left(\xi \frac{\partial \theta}{\partial \xi}\right)=\left.\xi \frac{\partial \theta}{\partial \xi}\right|_{0} ^{1}$          (47)

$\int_{0}^{1}\left(\xi-\xi^{3}\right) \frac{\partial \theta}{\partial \zeta} d \xi=\left.\frac{\partial \theta}{\partial \xi}\right|_{\xi=1}$          (48)

$\int_{0}^{1}\left(\xi-\xi^{3}\right)\left[\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \xi}\left(-\beta_{n}^{2}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)\right] d \xi =\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \xi}\left(-\beta_{n}^{2}\right) e^{-\beta_{n} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n}\right)+\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \xi}\left(-\beta_{n}^{2}\right)\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)$ (49)

$\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta}\left(-\beta_{n}^{2}\right) \int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) d \xi =\sum_{n=1}^{\infty} C_{n} e^{-\beta_{n}^{2} \zeta}\left(-\beta_{n}^{2}\right)\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)$  (50)

By combining Equations (48), (49) and (50), the equation can be reduced to;

$\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) d \xi =\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)$        (51)

Let’s multiply Equation (10) by Equation (52) and then integrate Equation (53),

$\left(\xi-\xi^{3}\right) e^{-\beta_{m} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{m}}{4}, 1, \beta_{m} \xi^{2}\right)$       (52)

$\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{m} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{m}}{4}, 1, \beta_{m} \xi^{2}\right) d \xi =\sum_{n=1}^{\infty} C_{n} \int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{m} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{m}}{4}, 1, \beta_{m} \xi^{2}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) d \xi$        (53)

The outcomes of multiplication and integration process will produce the following:

(i) If (n ≠ m)the result is equal to zero (0)

(ii) If (n = m) the result is

$\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2} / 2} F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right) d \xi =C_{n} \int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2}}\left[F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)\right]^{2} d \xi$      (54)

Substituting Equation (51) into Equation (54), the equation becomes;

$\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)= C_{n} \int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi^{2}}\left[F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)\right]^{2} d \xi$        (55)

And the constants Cn can be obtained by;

$C_{n}=\frac{\left(\frac{1}{2}-\frac{1}{\beta_{n}}\right) e^{-\beta_{n} / 2} F\left(\frac{3}{2}-\frac{\beta_{n}}{4}, 2, \beta_{n}\right)}{\int_{0}^{1}\left(\xi-\xi^{3}\right) e^{-\beta_{n} \xi}\left[F\left(\frac{1}{2}-\frac{\beta_{n}}{4}, 1, \beta_{n} \xi^{2}\right)\right]^{2} d \xi}$       (56)

4. Results and Discussion

4.1. Evaluation of the first four eigenvalues and the constant Cn 

A few values of the series coefficients are given in Table 1 together with the corresponding eigenvalues. The results of the calculated values of the center temperature as a function of the axial coordinate ζ are also summarized in Table 2.

Table 1. Eigenvalues and constants for Graetz’s problem


Eigenvalues βn

Coefficient Cn






















Table 2. Results of the center temperature functions θ (ζ)


Temperature (θ)

$\theta(\zeta, 0)$

















































The leading term in the solution for the center temperature is there for:

$\theta(\zeta, 0) \approx 0.9774 e^{-2.704 \zeta} G(0)$       (57)

4.2. Graphical representation of the exact solution of the Gratez problem

The center temperature profile is shown in Figure 2 using five terms to sum the series. As seen in this figure, the value of dimensionless temperature (θ) decreases with increasing values of dimensionless axial position (ζ). Note that the five-term series solution is not accurate for ζ<0.05 More terms needed here for the series to converge.

Figure 2. Variation of dimensionless temperature profile (θ) with dimensionless axial distance (ζ)

4.3. Comparison between the analytical model and the previous model simulation results

In order to compare the previous numerical results carried out previously by Shah and London [4] with the analytical model of our heat transfer problem, we chose to present the results of the numerical distribution of temperature with the numerical solution approached by these authors which gives the best results. Figure 3 plots the comparison results. It is clear from Figure 3 that there is a good agreement between numerical results and center analytical solutions  of the Graetz problem. 


Figure 3. A comparison between the present analytical results with the solution given by Shah and London [4]

5. Conclusions

In this paper, an exact solution of the Graetz problem is successfully obtained using the method of separation of variables. The hypergeometric functions are employed in order to determine the eigenvalues and constants, Cn and later with a find solution for the Graetz problem. The mathematical method performed in this study can be applied to the prediction of the temperature distribution in steady state thermally laminar heat transfer based on the fully developed velocity for fluid flow through a circular tube. In future work extensions, we recommend performing the Graetz solution by separation of variables in a variety of ways of accommodating non-Newtonian flow, turbulent flow, and other geometries besides a circular tube. It will be also interesting to solve the equation of the Graetz problem using numerical methods such as finite difference or orthogonal collocation for comparison purposes with the proposed exact solution. These numerical works will be carried out and included in the upcoming papers.



parameter of confluent hypergeometric function


parameter of confluent hypergeometric function


heat capacity


coefficient of solution defined in Equation

F (a;b;x)



standard confluent hypergeometric function thermal conductivity


length of the circular tube


radial direction of the cylindrical coordinates


radius of the circular tube


temperature of the fluid inside a circular tube


temperature of the fluid entering the tube


temperature of the fluid at the wall of the tube


maximum axial velocity of the fluid

Greek symbols




dimensionless axial direction


dimensionless temperature


dimensionless radial direction



density of the fluid

viscosity of the fluid


[1] Muhammad A., Aftab H., Syed H.A., Ali M.A., Muizz P.M. (2017). Simulation of corrugated plate heat exchanger for heat and flow analysis, International Journal of Heat and Technology, Vol. 35, No. 1, pp. 205-210. DOI: 10.18280/ijht.350127

[2] Manoj K.T., Rajsekhar P. (2017). Numerical analysis of natural convection in a triangular cavity with different configurations of hot wall, International Journal of Heat and Technology, Vol. 35, No. 1, pp. 11-18. DOI: 10.18280/ijht. 350102

[3] Lakhdar A., Yahia L., Khaled L. (2016). Numerical analysis of the parameters governing 3D laminar mixed convection flow in a rectangular channel with imposed wall flux density, International Journal of Heat and Technology, Vol. 34, No. 4, pp. 581-589. DOI: 10. 18280/ijht. 340405

[4] Shah R.K., London A.L.(1978). Laminar Flow Forced Convection in Ducts, Chap.V, Academic Press, New York.

[5] Braga N.R., de Barros L.S., Sphaier L.A. (2014). Generalized integral transform solution of extended Graetz problems with axial diffusion, ICCM2014 28-30th July, Cambridge, England, pp.1-14.

[6] Min T., Yoo J.Y., Choi H. (1997). Laminar convective heat transfer of a bingham plastic in a circular pipe—I. Analytical approach—thermally fully developed flow and thermally developing flow (the Graetz problem extended), International Journal of Heat and Mass Transfer, Vol. 40, No. 13, pp. 3025-3037. 

[7] Brown G.M. (1960). Heat or mass transfer in a fluid in laminar flow in a circular or flat conduit, AIChE Journal, Vol. 6, No. 2, pp. 179-183.

[8] Ebadian M.A., Zhang H.Y. (1989). An exact solution of extended Graetz problem with axial heat conduction, Interantional Journal of Heat and Mass Tranfer, Vol. 32, No. 9, pp. 1709-1717. 

[9] Lahjomri J., Oubarra A. (1999). Analytical solution of the Graetz problem with axial conduction, Journal of Heat Transfer,Vol. 121, No. 4, pp. 1078 -1083.

[10] Papoutsakis E., Damkrishna D., Lim H.C. (1980). The extended Graetz problem with Dirichlet wall boundary conditions, Applied Scientific Research, Vol. 36, pp. 13-34. 

[11] Liou C.T., Wang F.S. (1990). A computation for the boundary value problem of a double tube heat exchanger, Numerical Heat Transfer, Part A., Vol. 17, No. 1, pp. 109-125.

[12] Lahjomri J., Oubarra A., Alemany A. (2002). Heat transfer by laminar Hartmann flow in thermal entrance eregion with a step change in wall temperatures: the Graetz problem extended, International Journal of Heat and Mass Transfer ,Vol. 45, No. 5, pp. 1127-1148. 

[13] Fithen R.M., Anand N.K. (1988). Finite element analysis of conjugate heat transfer in axisymmetric pipe flows, Numerical Heat Transfer, Vol. 13, No. 2, pp. 189-203.

[14] Belhocine A. (2016). Numerical study of heat transfer in fully developed laminar flow inside a circular tube, International Journal of Advanced Manufacturing Technology, Vol. 85, No. 9, pp. 2681-2692.

[15] Bird R.B., Stewart W.E., Lightfoot E.N. (1960). Transport Phenomena, John Wiley and Sons, New York.

[16] Huang C.R., Matloz M., Dow-Pan W., Snyder W. (1984). Heat transfer to a laminar flow in a circular tube, AIChE Journal, Vol. 30 , No. 5, pp. 833-835.

[17] Slater L.J. (1960). Confluent Hypergeometric Functions, Cambridge University Press.