Fully developed pipe flow with heat transfer


Heat transfer is a process that is abundant in nature and extensively used for engineering applications, therefore a good understanding of the phenomenon allows to tackle various scientific and technological problems.  Numerical simulation of heating and cooling processes, if properly conducted, reduces development costs, improves safety and underlies optimization. The objective of today’s article is to present the fundamental principles of forced convection and to validate our QuickerSim CFD Toolbox against an analytical solution of a fully developed laminar pipe flow with constant flux value heat transfer.

Convective heat transfer

The transfer of heat can take place by one of three mechanisms – convection, conduction or radiation. Convection occurs when heat is transferred within a fluid by means of motion. Depending on the nature of this process, convection can be divided into two types – natural (free) convection and forced convection. The former takes place when the fluid’s motion is driven by buoyancy, that is an upthrust caused by a non-equilibrium state of gravitational forces, e.g. motion of a fluid in a cup heated at its bottom. The latter occurs when the fluid is forced to move by an external source such as a ventilator and this mechanism of heat transfer, due to its importance in engineering applications, is going to be the centerpiece of today’s article.

Laminar pipe flow 

The case under investigation is a fully developed, pressure-driven laminar pipe flow, also known as the Hagen-Poiseuille flow.  The laminarity of the flow allows to assume that there is no lateral disruption between layers of the fluid, therefore radial velocity components can be omitted. Furthermore, full development of the flow means that the streamwise velocity profile is independent of the axial coordinate. As previously mentioned, the flow is pressure-driven so the fluid’s axial velocity is dependent on the pipe pressure drop. The volumetric flow rate is given by the Hagen-Poiseuille equation:

Q=\frac{\Delta p \pi R^4}{ 8 \nu L}

Since in many engineering applications the pressure drop is seldom an input parameter but rather an unknown adapted to fulfill a system’s design requirements, we shall introduce a mean axial velocity based on the desired Reynolds number (Re) and the volumetric flow rate. Applying the above into the Hagen-Poiseuille equation yields the following formula:

u_z(r) = \frac{2 Q}{\pi R^2} (1 - \frac{r^2}{R^2})

where u_m=\frac{Q}{\pi R^2} is the mean velocity based on the volumetric flow rate. The velocity profile is a paraboloid with a maximum value of twice the mean velocity at the centerline and a value of zero at the wall. A section of it is presented below:

Heat transfer in a laminar pipe flow

The case is a solution of the thermal energy convection-diffusion equation:

\rho c(u_r \frac{\partial T}{\partial r}+u_z\frac{\partial T}{\partial z} = k[\frac{1}{r} \frac{\partial }{\partial r}(r\frac{\partial T}{\partial r}) + \frac{\partial ^2 T}{\partial z^2}]

The left-hand side of the equation represents convective heat transfer, that is heat transfered by the fluid’s motion. The radial velocity is equal to zero, so the first term of the left-hand side can be neglected. The right-hand side of the equation is responsible for thermal diffusion. Solving the equation requires a certain observation that allows to analyze the flow’s dynamic behaviour and thermal energy transport separately. Since the flow is laminar, we can assume that the dimensionless Eckert number, which represents the ratio between a flow’s kinetic energy and its heat transfer driving force, is small enough to disregard viscous dissipation. Therefore, the thermal energy equation can be supplemented with the velocity profile defined in the previous section.

A constant heat flux value condition implies that the temperature difference between the wall and the fluid is equal. However we already know that the temperature of the fluid is of non-constant value within the pipe. Therefore, we shall introduce a bulk mean temperature denoted by:

T_m =\frac{1}{u_m A} \int _A u_z T dA

Assuming that the local temperature gradient and the bulk mean temperature gradient in the streamwise direction are equal and of constant value, integration of the aforementioned thermal energy transport equation results in the following formula for radial temperature distribution:

T(r) = T_w - 2 \frac{u_m}{a} \frac{\partial T_m}{\partial z} R^2 (\frac{3}{16} + \frac{r^4}{16 R^4} -\frac{r^2}{4 R^2})

where a=k/\rho c is the thermal diffusivity coefficient (k stands for thermal conductivity and c is the fluid’s specific heat). The mean temperature gradient can be obtained by applying the desired volumetric flow rate Q and heat flux q to the heat conservation equation:

Q\rho c \frac{dT_m}{dz}=\pi Dq

To satisfy the constant wall flux condition, the value of the wall temperature has been coupled with the bulk mean temperature gradient.

T_w (z) = T_w(0) + \frac{\partial T_m}{\partial z} z


The objective of the case was to prove that compputational method is of the 2nd order. That means that the root mean square error of the simulated temperature profile drops quadratically with the number of grid elements. The analytical solution served both as a reference for calculating the error and as the inlet boundary condition for the simulation. Validation was performed for 6 succesively finer computational grids of a 0.5 metre-long pipe with a diameter of 0.1 metre within which water was flowing at a arbitralily chosen Reynolds number of 100 to ensure a fully laminar flow. The heat flux applied to the pipe walls was set as 150 \frac{W}{m^2}.

The simulated temperature profiles (dashed lines) have been compared with the analytical solution (solid line) at 0.8 of the pipe’s length.

The impact of mesh element size is already visible but to quantify the convergence rate, the root-mean-square error has to be plotted as a function of the linear number of grid elements. To verify whether the convergence is in fact of second order, a logarithmic scale was used for the graph axes. The rate of convergence can be compared to the plotted line which represents the appropriate function.

Final remarks

The validation case revealed that the convergence rate is almost of the second order. The most probable reason for this deviation is the grid structure. The mesh was built from tetrahedral elements and therefore grid resolution scaling doesn’t necessarily result in a uniform change of the number of elements in each direction. Furthermore, the core of the cylindrical domain, for geometrical reasons, had to be divided differently than the rest of the pipe interior, which contributed to the non-uniformity of the mesh. Another aspect that should be taken into consideration is the mesh quality. Due to hardware limitations the tetrahedral elements are quite skewed which also may have affected the overall simulation accuracy.  Nevertheless, our validation case can be confidently described as satisfactory.

Author: Robert Tykocki-Crow