Abstract

In this paper, convective heat transfer in a circular pipe is investigated for the case of a thermally developing laminar flow, by considering an axially varying wall heat flux and heat losses to the surrounding. The Integral Transform Technique is used in combination with a boundary filter solution, aiming to improve the convergence behavior of the associated eigenfunctions expansion. The associated linear eigenvalue problem yields a transcendental equation and eigenfunctions in the form of confluent hypergeometric function. The developed analytical solution is computed with MATLAB symbolic computation offering a straightforward way for the evaluation of the integral terms. Solution verification is performed by comparing the computed analytical solution with the solution obtained with COMSOL Multiphysics, a commercial Finite Element software, for the case of a uniform and a variable wall heat flux. Finally, it is shown that the temperature distribution and the local Nusselt number quickly converge with a relatively small number of eigenfunctions for the prescribed exponential wall heat flux examined herein.

Keywords

Convective Heat Transfer, Integral Transform Technique, Nusselt Number, Confluent Hypergeometric Function

Share and Cite:

1. Introduction

Numerical methods for solving differential equations (partial or ordinary) occupy an important place in the numerical simulation of physical, chemical, and biological phenomena due to their flexibility to take into account non-linear problems and complex geometries. Various commercial numerical simulation software available make use of numerical methods such as the finite difference method, the finite volume method, or the finite element method. A common point between these methods is that they require a spatial and temporal mesh, and user experience is often an asset to obtain satisfactory solutions. However, the development of analytical and hybrid solution methods remains necessary to provide benchmark solutions and perform asymptotic analysis for several problems of practical interest [1]. In addition, these methods can have a relatively low computational cost and thus prove to be an interesting alternative, particularly for solving inverse problems, optimization problems, or simulation problems under uncertainties in which the solution of the differential equation requires several evaluations. Moreover, these methods allow better error control.

Internal convective heat transfer has several applications in engineering. For instance, in the analysis of heat exchangers, one seeks to determine the amount of heat gained or lost by a fluid stream. Achieving this objective involves modeling convective heat transfer followed by solving the differential equations obtained using analytical, numerical, or even hybrid methods. Hydrodynamically developed flow, with developing thermal boundary layer, often referred as the classical Graetz problem, describes the evolution of the temperature field during laminar flow in a circular pipe with a constant wall temperature [2]. On the other hand, several versions of the Graetz problem, often referred to as extensions of the Graetz problem, have been investigated by different authors [2]-[6]. These include consideration of axial diffusion, viscous dissipation forces, and other types of boundary conditions [2]-[6]. In practical applications of heat exchangers or solar collectors, non-uniform thermal boundary conditions with heat losses are likely to occur, thus requiring a special attention [7]-[13]. While in [6], an exact analytical solution for thermally developing tube flow with variable wall heat flux was considered under ideal condition, herein we provide a solution by accounting for heat losses to the surrounding. Moreover, as compared to [6], we make use of a different splitting procedure and the eigenvalue problem is treated analytically yielding confluent hypergeometric functions.

The objective of this work is thus, to develop an analytical solution for the steady state forced convection heat transfer in a circular pipe, for an hydrodynamically developed flow in laminar regime and developing thermal boundary layer, by considering an axially varying wall heat flux and heat losses to the surrounding. From such development, the computation of the Nusselt number is straightforward.

2. General Formalism

Mikhailov and Ozisik, [14] analyzed a wide range of heat and mass transfer problems encountered in different applications of practical interest. In order to treat these problems in a systematic and unified approach, they proceeded to the classification of the linear problems of mass and heat diffusion examined in seven (07) classes. This classification will subsequently be adopted by several authors [1] [13] [14]. Class 1 problems, included the diffusion of heat or mass in transient regime in a finite domain of arbitrary geometry. The mathematical formulation of this class of problems is given by Equations (1a)-(1c) below [14]:

$w\left(x\right)\frac{\partial T\left(x,t\right)}{\partial t}=\nabla k\left(x\right)\nabla T\left(x,t\right)-d\left(x\right)T\left(x,t\right)+P\left(x,t\right),x\in V,t>0$ (1a)

$T\left(x,0\right)=f\left(x\right),x\in V$ (1b)

$\gamma \left(x\right)T\left(x,t\right)+\beta \left(x\right)k\left(x\right)\frac{\partial T\left(x,t\right)}{\partial n}=\varphi \left(x,t\right),x\in S$ (1c)

$T\left(x,t\right)$ represents the temperature or concentration at a point $x$ in space and at a given time $t$ ; $P\left(x,t\right)$ and $\varphi (x,t)$ represent source terms; $w\left(x\right)$ , $k\left(x\right)$ , $d\left(x\right)$ are the thermophysical properties of the medium; depending on the values of $\gamma \left(x\right)$ and $\beta \left(x\right)$ , we can obtain boundary conditions of the Dirichlet, Neuman or Robin type.

One of the widely used analytical methods for solving heat diffusion problems is the separation of variables method. However, it is important to note that the source terms of the governing equation and the boundary conditions are non-separable, which makes it impossible to obtain an analytical solution by the separation of variables method. On the other hand, the classical integral transformation method used in this study is perfectly suited for this type of problem. It consists of representing the solution of the problem in the form of an expansion in eigenfunctions of the solution. The systematic and unified approach as described by [1] [14]-[16] allows the formal solution of class 1 problems to be written in a general form. The development of the formal solution is omitted; interested readers may consult references [1] [14]-[18].

Experiments have shown slow convergence of the eigenfunction expansion near non-homogeneous boundary conditions, and therefore requires more terms to ensure convergence. However, finding the solution of the problem in the form of a solution to a particular problem and a problem with homogeneous boundary conditions can circumvent this difficulty.

Assuming that,

$T\left(x,t\right)=T_s\left(x,t\right)+T^{*}\left(x,t\right)$ (2)

Then the formal solution of class 1 problems obtained by integral transformation is in the form:

$T\left(x,t\right)=T_s\left(x;t\right)+{\displaystyle \sum _{i=0}^{\infty }\left(1/N_i\right){\Psi }_i\left(x\right)\left({\overline{f}}_i{\text{e}}^{-{\mu }_i^{2}t}+{\displaystyle \underset{0}{\overset{t}{\int }}{\overline{g}}_i}\left(t^{\prime }\right){\text{e}}^{-{\mu }_i^2\left(t-t^{\prime }\right)}\text{d}{t}^{\prime }\right)}$ (3)

where ${\mu }_i$ and ${\Psi }_i$ denote the associated eigenvalues and eigenfunctions and are obtained from the solution of the eigenvalue problem and the associated eigenvalue equation. The associated eigenvalue problem is part of a class of problems known as the Sturm-Liouville problem, and is chosen to correspond to the homogeneous version of the problem after separation of variables. $T_s\left(x,t\right)$ and $T^{*}\left(x,t\right)$ denote the variables of the particular problem and the problem with homogeneous boundary conditions and are obtained respectively from the boundary value problems given by the following Equations (4a)-(4b) and Equations (5a)-(5c):

$\nabla k\left(x\right)\nabla T_s\left(x;t\right)-d\left(x\right)T_s\left(x;t\right)+P\left(x\right)=0,x\in V$ (4a)

$\gamma \left(x\right)T_s\left(x;t\right)+\beta \left(x\right)k\left(x\right)\frac{\partial T_s\left(x;t\right)}{\partial n}=\varphi \left(x,t\right),x\in S$ (4b)

It should be noted that in the particular problem, time is treated as a parameter and not as an independent variable, hence the notation «;»

$w\left(x\right)\frac{\partial T^{*}\left(x,t\right)}{\partial t}=\nabla k\left(x\right)\nabla T^{*}\left(x,t\right)-d\left(x\right)T^{*}\left(x,t\right)+P^{*}\left(x,t\right)=0;x\in V,t>0$ (5a)

$T^{*}\left(x,0\right)=f^{*}\left(x\right)=f\left(x\right)-T_p\left(x;0\right)$ (5b)

$\gamma \left(x\right)T^{*}\left(x,t\right)+\beta \left(x\right)k\frac{\partial T^{*}\left(x,t\right)}{\partial n}=0,\underset{˜}{x}\in S,t>0$ (5c)

The terms $N_i$ , ${\overline{f}}_i$ , ${\overline{g}}_i$ appearing in the formal solution given by Equation (3) are respectively the norm, the integral transform of the initial condition Equation (5b), the integral transform of the source term and are defined as,

$N_i={\displaystyle \underset{V}{\int }w\left(x\right)}{\left[{\Psi }_i\left(x\right)\right]}^2\text{d}V$ (6)

${\overline{f}}_i={\displaystyle \underset{V}{\int }w\left(x\right)f^{*}\left(x\right){\Psi }_i\left(x\right)\text{d}V}$ (7)

with,

$f^{*}\left(x\right)=f\left(x\right)-T_s\left(x;0\right)$ (7a)

and

${\overline{g}}_i\left(t^{\prime }\right)={\displaystyle \underset{V}{\int }{P}^{*}\left(x,t\right){\Psi }_i\left(x\right)\text{d}V}$ (8)

with,

$P^{*}\left(x,t\right)=-w\left(x\right)\frac{\partial T_s\left(x;t\right)}{\partial t}$ (8a)

3. Analytical Solution for Thermally Developing Pipe Flow with Axially Varying Wall Heat Flux with Heat Losses

3.1. Description of the Physical Problem and Mathematical Formulation

We consider steady-state heat transfer in a hydrodynamically developed flow and in thermal development within a circular pipe subjected to an axially variable wall heat flux, and with heat losses to the surrounding environment (see Figure 1). Furthermore, axial conduction and natural convection are negligible. The fluid enters the pipe at a uniform temperature T₀.

Figure 1. Physical problem.

Based on the assumptions made, the mathematical formulation of the problem is given by the following Equation:

$u_{fd}\left(r\right)\frac{\partial T\left(r,z\right)}{\partial z}=\frac{\alpha }{r}\frac{\partial }{\partial r}\left(r\frac{\partial T\left(r,z\right)}{\partial r}\right),0<r<r_w,z>0$ (9a)

Subject to the initial condition,

$T\left(r,0\right)=T_0,0<r<r_w,z=0$ (9b)

And the boundary conditions,

${\frac{\partial T}{\partial r}|}_{r=0}=0,z>0$ (9c)

${k\frac{\partial T}{\partial r}|}_{r=r_w}=q_w\left(z\right)-h\left(T\left(r_w,z\right)-T_{\infty }\right),z>0$ (9d)

where,

$q_w\left(z\right)=q_{w0}{\text{e}}^{-mz}$ (10a)

$u_{fd}\left(r\right)=2u_m\left(1-\frac{r^2}{r_w^2}\right)$ (10b)

T denotes the temperature, $q_w$ denotes the applied boundary heat flux, where m is a scale parameter of dimension m⁻¹, h the heat transfer coefficient between the wall and the surrounding medium at temperature $T_{\infty }$ .

3.2. Mathematical Formulation in Dimensionless Form

Equations (9) given above, expressed in dimensionless form, are given as follows:

$\left(1-R^2\right)\frac{\partial \Theta \left(R,Z\right)}{\partial Z}=\frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial \Theta \left(R,Z\right)}{\partial R}\right)$ (11a)

Subject to the initial condition

$\Theta \left(R,Z\right)=10\le R\le 1,Z=0$ (11b)

And the boundary conditions

${\frac{\partial \Theta \left(R,Z\right)}{\partial R}|}_{R=0}=0,Z>0$ (11c)

${\frac{\partial \Theta \left(R,Z\right)}{\partial R}|}_{R=1}=Q_w^{*}\left(Z\right)-Bi\Theta \left(R,Z\right),R=1,Z>0$ (11d)

where the following dimensionless parameters and variables were introduced:

$\Theta \left(R,Z\right)=\frac{T\left(r,z\right)-T_{\infty }}{T_0-T_{\infty }}$ (12a)

$R=\frac{r}{r_w}$ (12b)

$U_{fd}^{*}\left(R\right)=\frac{u_{fd}\left(r\right)}{2u_m}=1-R^2$ (12c)

$Z=\frac{\alpha z}{2u_m{r}_w^2}$ (12d)

$Q_w^{*}\left(Z\right)=\frac{q_w\left(z\right)r_w}{k\left(T_0-T_{\infty }\right)}=Q_{w0}^{*}{\text{e}}^{-MZ}$ (12e)

$Bi=\frac{h{r}_w}{k_{}}$ (12f)

Equations (12a)-(12f) respectively designate the following dimensionless variables and parameters: temperature, radial position, velocity profile, axial position, prescribed wall heat flux, Biot number. In Equation (12e), M denotes the dimensionless counterpart of the scale parameter m.

3.3. Solution of the Problem through the Integral Transformation Technique

We first proceed in establishing a correspondence between the problem under study given by Equations (11a)-(11d) with the general formalism given by Equations (1a)-(1c) in Section 2. By doing so, we can write the following for the variables and the geometry

$x=R;t=Z;T\left(x,t\right)=\Theta \left(R,Z\right)$ (13a-d)

$V=\left[0,1\right];S=\left[0,1\right];\text{d}V=R\text{d}R$ (13e-g)

while for the governing equation and the initial condition, we have

$w\left(x\right)=\left(1-R^2\right);k\left(x\right)=1;d\left(x\right)=0;P\left(x,t\right)=0;f\left(x\right)=1$ (14a-e)

with the associated boundary conditions:

for $R=0:\gamma \left(x\right)=0;\beta \left(x\right)=1;\varphi \left(x,t\right)=0$

for $R=1:\gamma \left(x\right)=Bi;\beta \left(x\right)=1;\varphi \left(x,t\right)=Q_w^{*}\left(Z\right)$ (15a-f)

The formal solution in this case is given as follow (see Equation (3)):

$\Theta \left(R,Z\right)={\displaystyle \sum _{i=1}^{\infty }\frac{1}{N_i}{\Psi }_i}\left(R\right)\left[{\overline{f}}_i{\text{e}}^{-{\mu }_i^{2}Z}+{\displaystyle {\int }_0^Z{\overline{g}}_i\left(Z^{\prime }\right){\text{e}}^{-{\mu }_i^2\left(Z-Z^{\prime }\right)}\text{d}{Z}^{\prime }}\right]+{\Theta }_s\left(R;Z\right)$ (16)

It is recalled here that we are seeking the solution of the problem given by Equations (11a)-(11d) in the form,

$\Theta \left(R,Z\right)={\Theta }_s\left(R;Z\right)+{\Theta }_h\left(R,Z\right)$ (17)

where the first and second term of the right-hand side represent respectively, the dimensionless temperature of the particular problem and of the homogeneous problem.

Here, the choice of a boundary condition filter is appropriate, and the objective is to make the boundary conditions homogeneous. We use a linear boundary condition filter in the form

${\Theta }_s\left(R;Z\right)=a\left(Z\right)R+b\left(Z\right)$ (18)

Substituting Equation (18) into the boundary conditions given by Equations (11c)-(11d), we obtain the following system of equations:

$\{\begin{array}{l}a\left(Z\right)=0\\ a\left(Z\right)+Bi\left(a\left(Z\right)+b\left(Z\right)\right)=Q_w^{*}\left(Z\right)\end{array}$ (19a)

This allows us to obtain,

$b\left(Z\right)=\frac{Q_w^{*}\left(Z\right)}{Bi}$ (19b)

The solution to the particular problem is then given by Equation (20):

${\Theta }_s\left(R;Z\right)=\frac{Q_w^{*}\left(Z\right)}{Bi}$ (20)

Or by Equation (21), after substitution of the dimensionless wall heat flux $Q_w\left(Z\right)$ by its expression

${\Theta }_s\left(R;Z\right)=\frac{Q_{w0}^{*}{\text{e}}^{-MZ}}{Bi}$ (21)

Substituting Equation (17) into Equations (11a)-(11d) and taking into account Equation (21), after rearrangement, we obtain the non-dimensional form of the homogeneous problem,

$\left(1-R^2\right)\frac{\partial {\Theta }_h\left(R,Z\right)}{\partial Z}=\frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial {\Theta }_h\left(R,Z\right)}{\partial R}\right)+\left(1-R^2\right)\frac{M{Q}_{w0}^{*}{\text{e}}^{-MZ}}{Bi}$ (22a)

${\Theta }_h\left(R,Z\right)=1-\frac{Q_{w0}^{*}}{Bi},0\le R\le 1,Z=0$ (22b)

${\frac{\partial {\Theta }_h\left(R,Z\right)}{\partial R}|}_{R=0}=0,Z>0$ (22c)

$\frac{\partial {\Theta }_h\left(R,Z\right)}{\partial R}+Bi{\Theta }_h\left(R,Z\right)=0,R=1,Z>0$ (22d)

where $Q_w^{*}\left(0\right)$ and $\partial Q_w^{*}\left(Z\right)/\partial Z$ where replaced by their expressions in Equations (22).

The auxiliary eigenvalue problem associated with Equations (22a)-(22d) is given by:

$\frac{{\text{d}}^2\Psi \left(R\right)}{\text{d}{R}^2}+\frac{1}{R}\frac{\text{d}\Psi \left(R\right)}{\text{d}R}+{\mu }^2\left(1-R^2\right)\Psi \left(R\right)=0$ (23a)

${\frac{\text{d}\Psi \left(R\right)}{\text{d}R}|}_{R=0}=0$ (23b)

$\frac{\text{d}\Psi \left(R\right)}{\text{d}R}+Bi\Psi \left(R\right)=0,R=1$ (23c)

Under the symmetry condition given by Equation (23b), the solution of Equation (23a) can be represented by [14] [19]:

$\Psi \left(R\right)=\Omega \left(R\right)\exp \left(-\frac{\mu R^2}{2}\right)$ (24)

Substituting Equation (24) into Equation (23a) and after rearrangement, we obtain:

$\frac{{\text{d}}^2\Omega \left(R\right)}{\text{d}{R}^2}+\left(\frac{1}{R}-2\mu R\right)\frac{\text{d}\Omega \left(R\right)}{\text{d}R}+\left({\mu }^2-2\mu \right)\Omega \left(R\right)=0$ (25)

By defining,

$z=x^2=\mu R^2$ (26)

Using Equation (26) in Equation (25), it comes:

$\frac{{\text{d}}^2\Omega \left(x^2\right)}{\text{d}{x}^2}+\left(\frac{1}{x}-2x\right)\frac{\text{d}\Omega \left(x^2\right)}{\text{d}x}-4\left(\frac{1}{2}-\frac{\mu }{4}\right)\Omega \left(x^2\right)=0$ (27)

The function that satisfies Equation (27) is the confluent Kummer hypergeometric function in the form defined by [14] [19]:

$\frac{{\text{d}}^2\Omega \left(x^2\right)}{\text{d}{x}^2}+\left(\frac{2c-1}{x}-2x\right)\frac{\text{d}\Omega \left(x^2\right)}{\text{d}x}-4a\Omega \left(x^2\right)=0$ (28a)

With,

$a=\frac{1}{2}-\frac{\mu }{4};c=1$ (28b)

$\Omega \left(R\right)=1F1\left(a,c;z\right)={\displaystyle \sum _{n=0}^{\infty }\frac{a_n}{c_n}}\frac{z^n}{n!}$ (29)

Then Equation (24) becomes:

$\Psi \left(R\right)=\exp \left(-\mu \frac{R^2}{2}\right)1F1\left(a,c;\mu R^2\right)$ (30a)

With,

$a=\frac{1}{2}-\frac{\mu }{4};c=1$ (30b)

Using Equation (30a) in the second boundary condition given by Equation (23c), we obtain after rearrangement, the following non-linear algebraic equation:

$1F1\left(a,c;{\mu }_i{R}^2\right)\left[Bi-{\mu }_i\right]+1F{1}^{\prime }\left(a,c;{\mu }_i{R}^2\right)=0,i=1,2,\cdots$ (31)

where $1F{1}^{\prime }$ is the derivative with respect to R of the confluent Kummer hypergeometric function evaluated at R = 1.

The ${\mu }_i$ are the eigenvalues and are obtained from the eigenvalue equation given by Equation (31), while the eigenfunctions are given by the following equation

${\Psi }_i\left(R\right)=\exp \left(-\frac{{\mu }_i{R}^2}{2}\right)1F1\left(a,c;{\mu }_i{R}^2\right)$ . (32)

At this stage, we know how to compute the eigenvalues and the eigenfunctions, but also the solution of the particular problem. In order, to compute the solution given by Equation (16), we need to develop the expressions of the terms N_i (norm of the eigenfunctions), ${\overline{f}}_i$ (integral transform of the initial condition) and ${\overline{g}}_i$ (integral transform of the source term), as done next.

Expression of the transformed initial condition

By using Equations (7) and the appropriate correspondences Equations (13)-(14), we can write,

${\overline{f}}_i={\displaystyle \underset{0}{\overset{1}{\int }}{U}_{fd}^{*}}\left(R\right)\left(\Theta \left(R,0\right)-{\Theta }_s\left(R,0\right)\right){\Psi }_i\left(R\right)R\text{d}R$ (33)

Substituting $U_{fd}^{*}$ , $\Theta \left(R,0\right)-{\Theta }_s\left(R;0\right)$ , and ${\Psi }_i$ by their expressions in Equation (33), and after rearrangement we obtain:

${\overline{f}}_i=\left(1-\frac{Q_{w0}^{*}}{Bi}\right){\displaystyle \underset{0}{\overset{1}{\int }}R\left(1-R^2\right)}\exp \left(-\frac{{\mu }_i{R}^2}{2}\right)1F1\left(a,c;{\mu }_i{R}^2\right)\text{d}R$ (34a)

With,

$a=\frac{1}{2}-\frac{{\mu }_i}{4};c=1$ (34b)

Expression of the transformed source term

Similarly, by using Equations (8) and the appropriate correspondences Equations (13)-(14), we can write,

${\overline{g}}_i\left(Z\right)={\displaystyle \underset{0}{\overset{1}{\int }}{P}^{*}}\left(R,Z\right){\Psi }_i\left(R\right)R\text{d}R$ (35a)

where,

$P^{*}\left(R,Z\right)=-U_{fd}^{*}\left(R\right)\frac{\partial {\Theta }_s}{\partial Z}$ (35b)

With ${\Theta }_s$ given by Equation (21). Then,

$P^{*}\left(R,Z\right)=\left(1-R^2\right)\frac{M{Q}_{w,0}^{*}{\text{e}}^{-MZ}}{Bi}$ (36)

By substituting Equations (30a) and (36) in Equation (35a), and after rearrangement we obtain:

${\overline{g}}_i\left(Z\right)=\frac{P{Q}_{w,0}^{*}{\text{e}}^{-pZ}}{Bi}{\displaystyle \underset{0}{\overset{1}{\int }}R\left(1-R^2\right)}\exp \left(-\frac{{\mu }_i{R}^2}{2}\right)1F1\left(a,c;{\mu }_i{R}^2\right)\text{d}R$ (37)

Expression of the norm

Finally, the norm is obtained by repeating the same procedure as follows:

$N_i={\displaystyle \underset{0}{\overset{1}{\int }}{U}_{fd}^{*}\left(R\right){\Psi }_i^2\left(R\right)R\text{d}R}$ (38)

Which can be rewritten as,

$N_i={\displaystyle \underset{0}{\overset{1}{\int }}R\left(1-R^2\right)}{\left(\exp \left(-\frac{{\mu }_i{R}^2}{2}\right)1F1\left(a,c;{\mu }_i{R}^2\right)\right)}^2\text{d}R$ (39)

To summarize, the solution of problem (11a)-(11d) is given by Equation (13) where the terms, ${\Psi }_i\left(R\right)$ , ${\overline{f}}_i$ and ${\overline{g}}_i\left(Z^{\prime }\right)$ are defined respectively by Equations (30), (34), (37) and (39). While the eigenvalues ${\mu }_i$ are roots of Equation (31).

The Nusselt number is easily computed with,

$Nu=\frac{\partial \Theta \left(R=1,Z\right)/\partial R}{\Theta \left(R=1,Z\right)-{\Theta }_b\left(Z\right)}$ (40a)

where,

${\Theta }_b\left(Z\right)=\frac{{\displaystyle \underset{0}{\overset{1}{\int }}R\left(1-R^2\right)\Theta \left(R,Z\right)\text{d}R}}{{\displaystyle \underset{0}{\overset{1}{\int }}R\left(1-R^2\right)\text{d}R}}$ (40b)

represents the dimensionless bulk temperature or mixing-cup temperature and the numerator is the dimensionless temperature gradient at the wall.

4. Results and Discussions

Equations (16), (32), (34), (37), and (39) represent the solution to the problem defined by (11a)-(11d). To illustrate the results, we consider the case of a flow of water in a circular pipe with a diameter D = 25 mm and a length of L = 12 m. The water enters the pipe at a temperature of 75˚C and the medium surrounding the pipe is at a temperature of 25˚C. The thermophysical properties of the water are taken at T = 95˚C as k = 0.68 W/m∙K, c_p = 4212 J/(kg K), ρ = 961.5 kg/m³ [20]. The water enters in the pipe at an average velocity of 0.027 m/s, which correspond to an average Reynolds number of 2185; thus, the flow is laminar. The prescribed wall heat flux is in the form of an exponentially decaying function with an amplitude q_w0 = 450 W/m² and a scale parameter m. By setting the value of the scale parameter m to 0 m⁻¹, 0.02 m⁻¹ or 0.2 m⁻¹, one obtains a uniform wall heat flux, a linear-like decaying wall heat flux, or a truly exponential decaying behavior. The Biot number is taken as 0.001. For these different three cases, the total heat flux is kept constant by normalizing the prescribed wall heat flux as follow,

$q_{w,norm}(z)=Q_{CHF}\frac{q_w\left(z\right)}{{\displaystyle \underset{A}{\int }{q}_w\left(z\right)dA}}$ (41)

where $Q_{CHF}$ represents the total heat flux for a uniform heating of the wall. The total heat flux for a uniform heating with heat flux density q_w0 = 450 W/m² is computed by multiplying the heat flux density with the lateral surface area of the pipe, and it amounts to 424.12 W. Similarly, the integral term of Equation (41) is obtained by integrating the prescribed wall heat flux density over the lateral surface area of the pipe.

In order to perform solution verification, the case of a prescribed uniform wall heat flux is first examined. It must be noticed that, such a case can be obtained by simply keeping all the parameters specified above, except the scale parameter m and the Biot number which are set to 0 and 10⁻⁵, respectively. By doing so, we obtain the case of a uniform wall heat flux of magnitude 450 W/m² with negligible heat losses to the surrounding. Table 1 presents the first five eigenvalues obtained by solving Equation (31) using the secant method for this case. The analytical solution was implemented on the MATLAB platform using symbolic computation. Figure 2(a) presents the temperature distribution obtained for this case with N = 20 eigenfunctions used in the eigenfunction expansion. As can be noticed, the temperature of the fluid increases as it flows along the pipe, and the fluid appears warmer close to the wall and get colder away of the wall.

Table 1. First five eigenvalues.

(a)

(b)

Figure 2. Temperature distribution (˚C) of the fluid stream for a prescribed uniform wall heat flux with (a) N = 20 eigenfunctions; (b) COMSOL.

In order to perform solution verification of the analytical solution, a numerical solution of the problem given by Equations (9)-(10) is obtained through the finite element software COMSOL Multiphysics with the same input parameters specified above. The heat transfer in fluids module is used for such purpose. The inflow and outflow settings of the later are used to specify the entrance temperature and a zero normal temperature gradient at the exit, while the prescribed wall heat flux is given in the form of a user defined function. A customized refined mesh built with roughly 20,434 elements composed of 19800 quads, 1260 edge elements and 4 vertex elements. The relative tolerance is 10⁻³ and correspond to the default solver settings.

Figure 2(b) presents the numerical solution obtained with COMSOL Multiphysics. As can be noticed in these figures, both numerical and analytical temperature fields are very similar. Table 2 illustrates the convergence behavior of the temperature eigenfunction expansion with an increasing number of terms used in the series solution given by Equations (16), (32), (34), (37), and (39), for selected positions in the fluid. As it can be noticed in this table, convergence is quickly achieved with a small number of terms. Also shown in this table, is the computed numerical solution. Again, one can note that both numerical and analytical solution are very close.

Table 2. Illustration of convergence behavior of temperature field eigenfunction expansion for uniform heating, with associated relative error (%).

Both temperature distributions obtained from the numerical and analytical solutions are used to compute the local Nusselt number and is given in Figure 3. It can be noticed in this figure, that both numerical and analytical solutions perfectly match at the graph scale. Moreover, the Nusselt number is shown to be very close to 4.36 in the fully developed region, which corresponds to the exact value for a prescribed uniform wall heat flux as given in classical heat transfer books [20] [21] and related papers [22] [23].

Figure 3. Variation of the local Nusselt number with the dimensionless axial position.

We now consider, the case of a prescribed wall heat flux in the form of an exponentially decaying function, by setting the scale parameter m = 0.2 and the Biot number to 10⁻³ in order to account for heat losses to the surrounding. The results obtained for the temperature field as a function of the number of eigenfunctions used in the eigenfunction expansion are presented in Figures 4(a)-(b). One can notice that, the temperature fields match very well. Moreover, by comparing these temperature fields with the temperature field obtained with COMSOL (see Figure 4(c)), one can infer the validity of the analytical solution developed herein.

(a)

(b)

(c)

Figure 4. Temperature distribution of the fluid stream obtained for a prescribed exponential wall heat flux for m = 0.2 with (a) N = 15; (b) N = 20 eigenfunctions; (c) COMSOL.

The accuracy of the analytical solution is further illustrated in Table 3 through the convergence behavior of the temperature field eigenfunction expansion at selected positions within the fluid domain for an increasing number of terms. It is worth noting in this table that, these temperatures are consistent with the temperatures obtained with COMSOL Multiphysics.

Table 3. Illustration of convergence behavior of temperature field eigenfunction expansion for exponential heating (m = 0.2), with associated relative error (%).

In addition to the temperature distribution, one can compute the local Nusselt number. Figure 5 presents the local Nusselt number as a function of the axial position, together with the prescribed wall heat flux. As can be noticed in these figures, the local Nusselt number is higher near the entrance of the pipe and exhibits an asymptotic behavior for increased axial positions within the pipe. A similar behavior is observed for the prescribed uniform wall heat flux. Moreover, one can notice that, convergence of the computed local Nusselt number is achieved for the specified number of eigenfunctions. It is worth noting that, the values of the eigenvalues increase rapidly (see Table 1), reducing the contribution of higher order terms in the convergence of the summation.

(a)

(b)

Figure 5. (a) Prescribed exponential wall heat flux for m = 0.2; (b) Local Nusselt number as a function of the dimensionless axial position and the number of eigenfunctions.

Finally, we consider the case of a prescribed exponential wall heat flux with the scale parameter m set to 0.02. As can be seen in Figures 6(a)-(b), the temperature distributions are again physically meaningful and are very similar for both number of eigenfunctions used in the eigenfunction expansion. In fact, the fluid appears warmer near the wall, and get colder as the distance to the wall increases. It is worth noting that, the temperature distribution as well as, the maximum temperature is very similar to those obtained in the case of a uniform wall heat flux (see Figure 2). Moreover, as compared with the case of an exponential wall heat flux (m = 0.2), the present case exhibits a maximum temperature which is higher by 4˚C. This shows that, by manipulating the prescribed wall heat flux, one can reduce the maximum wall temperature while using the same total heat flux. Again, one can notice in Figure 6(c), that the numerical solution obtained with COMSOL is very similar to the analytical temperature fields given in Figures 6(a)-(b).

(a)

(b)

(c)

Figure 6. Temperature distribution of the fluid stream for a prescribed exponential wall heat flux for m = 0.02 obtained with (a) N = 15; (b) N = 20 eigenfunctions; (c) COMSOL.

Similarly, with the previous cases, the accuracy of the analytical solution is evidenced in Table 4. One can note, the convergence behavior of the temperature field eigenfunction expansion at selected positions within the fluid domain for an increasing number of terms. Moreover, these analytical temperatures are consistent with the temperatures obtained with the finite element software COMSOL Multiphysics.

Table 4. Illustration of convergence behavior of temperature field eigenfunction expansion for exponential heating (m = 0.02), with associated relative error (%).

Figure 7(a) presents the profile of the dimensionless exponential wall heat flux for m = 0.02, while the local Nusselt number is given in Figure 7(b). Again, one can notice that, convergence of the computed local Nusselt number is quickly achieved with a relatively small number of eigenfunctions.

(a)

(b)

Figure 7. (a) Prescribed exponential wall heat flux for m = 0.02; (b) Local Nusselt number as a function of the dimensionless axial position and the number of eigenfunctions.

5. Conclusion

In this work, an analytical solution was developed through the Integral Transform Technique, for the convective heat transfer in a circular pipe for a thermally developing and hydrodynamically developed laminar flow under an axially varying wall heat flux with heat losses to the surrounding medium. The solution of the original problem is sought as the sum of a homogenous problem and a non-homogenous one. A boundary condition filter is used to homogenize the boundary condition in order to accelerate the convergence. The solution of the homogeneous problem yields an eigenvalue problem for which confluent hypergeometric functions are the eigenfunctions. The solution of the problem is obtained for a prescribed number of eigenfunctions and convergence is quickly achieved for the different wall heat fluxes profiles, namely: uniform, exponential with two different values of the scale parameter.

Author Contributions

B. LAMIEN developed the analytical solution. B. LAMIEN, S. SANOGO and I. TOUGRI performed numerical simulations and interpretation. B. LAMIEN prepared the figures and drafted the manuscript. B. LAMIEN, S. SANOGO, I. TOUGRI and M. BEIDARI made critical revisions to the manuscript. All authors revised and approved the final manuscript.

Acknowledgements

The work was supported by the Ministry of Higher Education and Research of Burkina Faso.

Nomenclature

Greek Letters

Subscripts and Superscripts

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References