Similarity solution and Runge Kutta method to a thermal boundary layer model at the entrance region of a circular tube: The Lévêque Approximation

The heat equation in cylindrical coordinates

The general equation for heat transfer in cylindrical coordinates developed by Bird, Stewart, and Lightfoot (1960) is as follows;

The Graetz-Poiseuille flow problem

The Graetz problem is the problem of determining the steady-state temperature field in a fluid flowing in a circular tube when the wall of the tube is held at a uniform temperature and the fluid enters the tube at a different uniform temperature. The fluid enters the tube at a temperature T₀ and encounters a wall temperature at Tω, which can be larger or smaller than T₀. The Graetz problem considers the thermal entry of an incompressible fluid in a circular tube with a fixed velocity profile. The situation is illustrated in figure 1.

Figure 1: Illustration of Graetz problem.

After some hydrodynamic entry length, the velocity profile approaches a steady profile; that is, it ceases to change the downstream position. A fully developed parabolic velocity profile is established in a circular duct and remains unchanged over the length of the duct. There is a sudden jump in the wall temperature, and the fluid temperature is initially uniform at the upstream wall temperature. The thermal entry problem is to determine the behavior temperature profile as it changes to be uniform at the downstream wall temperature. Because the flow is incompressible, the velocity distribution does not depend on the varying temperatures.

1. The flow is steady, laminar and fully developed flow (Re < 2400)

2. Constant physical properties would also be and would not vary with temperature such as ρ, µ, Cp , k. This assumption also implies incompressible Newtonian flow.

3. Axisymmetric temperature field ,

4. Negligible viscous dissipation

5. For a fully developed hydrodynamic boundary layer, namely Hagen Poiseuille flow, the velocity profile is:

Where is the maximum velocity existing at the centerline, , and .

The energy equation is subject to the assumptions mentioned above, Eq. (1) can be written as follows:

where is the thermal diffusivity of the fluid and, our problem is subjected to the following boundary conditions: at the inlet of the tube ; at the wall of the tube and at the centerline is finite or

A suitable set of nondimentional variables is given as:

where the Péclet Number

After replacement of these variables, the nondimentional energy equation and the boundary conditions applied to the system will take the following form:

It is apparent that as the Péclet number (Pe) increases, the effect of the axial diffusion is significantly diminished. The boundary layer assumptions call for neglecting the axial diffusion altogether, yielding the following governing equation:

This is a linear equation whose solution can be determined by the method of separation of variables. This method yields an infinite hypergeometric series solution for the scaled temperature field:

Where , are the eigenvalues and , are the eigenfunctions of a proper Sturm-Liouville system. The coefficients can be obtained by using the orthogonality property of the eigenfunctions defined as follows:

The Lévêque Approximation

For all values of the axial position, the orthogonal function expansion solution obtained in the resolution of the classical Graetz problem is quite convergent, but the convergence is very slow as soon as one approaches the input tube. Indeed, for very long values of Z, the factor has become converged. Lévêque (1969) observed the thermal input region in a circular tube and developed a periodic solution, which is useful precisely where the expansion of the orthogonal function converges too slowly (figure 2).

Figure 2: Schematics of the Lévêque problem and the coordinate system.

We will presently develop the Lévêque solution based on the hypothesis that the thickness of the thermal boundary layer This assumption leads to the following reductions:

• In the radial conduction term, the curvature effects can be neglected. Thus, derivative is approximated by

• Since we have been captivated that by the velocity distribution in the thermal boundary layer, the velocity field is developed in a Taylor series in the distance measured from the tube wall and if we keep the first non-zero term

If we set x=R-r, the speed distribution will take the following form:

We know that the boundary conditions outside the thermal boundary layer are those of the fluid entering the tube, we will use the boundary condition instead of the centerline boundary condition used to arrive at the Graetz solution.

Governing Lévêque’s Equation

Starting from the reduced energy equation whose axial conduction has been neglected yet, and considering the said hypotheses, for the temperature field, we obtain the following governing equation

We have used the chain rule in order to transform the second derivative in r into a secondary derivative of x.

Boundary Conditions

The temperature field T(x, z),is subject to the following boundary conditions.

Non-Dimensionalization

Now, we will use dimensionless variables for the simplification of the equation. For this, we introduce the temperature and the axial coordinate of the following forms

The scaled governing equation from the wall via X=x/R and the boundary conditions are given as follows

Analytical Methodology for Problem Solving: Temperature Field and Thermal Boundary Layer

At the current problem, we are looking for a similarity solution for the temperature field, we assume, , where is the similarity variable and is an unknown variable that draws the thickness of the scaled thermal boundary layer. Using the chain rule, we will perform the following necessary transformations.

When we use these results, the partial differential equation for θ (X, Z), is converted to an ordinary differential equation for _^F(n)

We put the term in parentheses () a constant is equal to 3/2 because it is certain that the similarity hypothesis will fail unless this quantity is required to be independent of Z. Finally, we obtain an ordinary differential equation for and another for .

In order to derive the boundary conditions of these functions, it is enough to go to the limiting conditions on We notice that which implies F(0)=0 and which leads to The residual condition (input) gives us

By favoring this condition collapses into the condition obtained already from the boundary condition on the scaled temperature field as X → ∞. By synthesizing the boundary conditions on / we obtain:

From the equation, we can write

By integrating the two terms of the obtained equation

We arrive at the following expression

Where k and C are constants of the integral. By analogy, the following is drawn

Finally, the solution of the equation will take the following expression:

The function checks the initial condition for , and also considers the boundary condition for ; which implies

From where

So

The scaled boundary layer thickness is calculated by the integration, which gives the following solution

Finally, the solution of our differential equation takes the following form

Where is the Gamma function (Abramowitz and Stegun, 1965), a Matlab code was used to approximate the values of the integral and the function F(η) for each abscissa η.

Resolution of the problem using the Runge-Kutta 4th order method

The original ODE of our problem is defined as follows:

With η=

However, we wish to use the 4th order Runge-Kutta method, so I have the system:

Withand Now, we know that for two general 1st order ODE’s

The 4th order Runge-Kutta formula’s for a system of 2 ODE’s are:

Where

We typically need some inputs for the algorithm:

• A range that we want to do the calculations over: a≤t≤b, let’s use a=0, b=1

• The number of steps N, say N=10

• The step size h=(b−a)/N=1/10=0.1

The flowchart for the above process is shown in figure 3.

Figure 3: Flowchart of the RK-4 method for resolving the second ODE’s systems.

The main program is written in Fortran language to solve the Levèque approximation containing a system of two differential equations of order 1 using the fourth-order Runge Kutta method RK04. This program relies on a definition of two functions whose subroutine RK04 is called at each repetition of the loop that intervenes in the calculations. The code edited in the machine that was executed is illustrated in detail in figure 4.

Figure 4: Fortran code of Runge Kutta for set of first order differential equations.

Comparison between the analytical and the numerical results

Analytical and numerical results that describe solutions with both methods in this analysis are summarized in table 1.

Table 1: Results of the exact and the numerical solution.

Figure 5 shows a comparison between the resolution results of the equation predicted by the analytical method and the numerical data derived from the Fortran code, the two sets of results of which are plotted in the same figure. On the basis of figure 5, it can be seen that the two curves are fairly identical, while observing that the dimensionless temperature θ gradually and gradually increases to the abscissa Z = 0.7, then loops and arches a little, by varying its path until it reaches the position Z = 1.7 where it stabilizes at a constant value 0.79 along the tube until the outlet of the fluid stream. In the same figure, the derivative function which is physically interpreted as the variation of the thermal transfer coefficient (h) indicated in blue color which is a solution of our system of differential equation by the Runge Kutta method decreases exceptionally as it moves away from the inlet region and then reaches the value zero on the abscissa Z = 1.6 until the exit of the flow. Figure 5 indicates that the present analytical results are in good agreement with the numerical results and theoretical analysis calculated using a fourth order Runge Kutta method. However, with RK04 approach, the method provides an extremely accurate approximation.

Figure 5: Comparison of exact and fourth-order Runge Kutta (RK4) numerical solutions.

Figure 6 depicts variation of thermal boundary layer thickness with axial position. The thickness of the boundary layer gradually elongates from zero moving in the direction of the flow of the fluid when the fluid enters the tube and eventually leads to the center of the tube and invades the whole of the pipe. The wall shear stress is the highest at the pipe inlet where the thickness of the boundary layer is smallest, and decreases gradually to the fully developed value. Therefore, the pressure drop is higher in the entrance regions of a pipe, and the effect of the entrance region is always to increase the average friction factor for the entire pipe. This increase may be significant for short pipes, but is negligible for long ones. It can be seen that directly at the wall is a thin layer where the velocity is considerably lower than it is at some distance from the wall. The thickness of this layer increases along the tube from front to back. The fluid velocity in a pipe changes from zero at the wall because of the no-slip condition to a maximum at the pipe center.

Figure 6: Thermal boundary layer thickness distribution by analytical method.

The heat transfer coefficient

Depending on the axial position, the heat flux from the wall to the fluid, we can calculate it directly using the following formula:

By usual notation, we define the heat transfer coefficient h (z) as follows:

Where is the bulk or cup-mixing average temperature. The mathematical definition of the bulk average temperature is given by the following expression:

Where is the velocity field. The heat transfer coefficient is related to the temperature gradient at the tube wall, we can estimate it as follows:

We can define the Nusselt number as a dimensionless heat transfer coefficient.

Where is the dimensionless bulk average temperature

When the thermal boundary layer is thin in the thermal entrance region, the average bulk temperature Tb can be approximated by the temperature of the fluid entering the tube T0.

As a result, the heat transfer coefficient in this input region is defined as:

Changing into dimensionless variables and defining a Nusselt number Nu=2hR/k we can write the following:

By substituting and , the approximate Nusselt number obtained in the thermal entrance region is given by the following form.

By comparing with the exact solution, we can draw that this result is a good approximation in the range:

Figure 7 exhibits the Nusselt number versus axial distance, Z obtained in the thermal entrance region, for different Reynolds numbers. As expected, the Nusselt number, Nu(Z), enhances by increasing the tube radius, and that this effect is magnified near the entrance. When Z is greater than a certaine distance all the plots become nearly flat, indicating a thermally fully-developed condition. Indeed, when fluid enters the tube with tube walls at a different temperature from the fluid temperature, thermal boundary layer starts growing. After some distance downstream (thermal entry length) thermallyfully developed condition is eventually reached.

Figure 7: Nusselt number in the thermal entrance region versus axial position for different tube radius.

Figure 8 shows the effect of the Péclet number on the Nusselt number at a various axial distance. It’s clearly seen, the Nusselt number increases with the increasing Péclet number. As can be seen, the Péclet number has a much more pronounced effect on the Nusselt values for positions near the tube entrance. However the curve exhibits the same overall behavior-larger Nu at small Z and more or less constant value of large Z. It is evident that the local Nusselt number decreased in the entry region as the thermal boundary layer developed. It reached a constant value independent of heat flux and Reynolds number in the thermally fully developed region. Hence, the heat transfer coefficient (h) is infinite in the beginning (boundary layers just building up), then decays exponentially to a constant value when flow is fully developed (thermally) and thereafter remains constant. The results show that the value of the Nusselt number starts high and decrease rapidly along the length of the tube.

Figure 8: Nusselt number in the thermal entrance region versus axial position with various Peclet numbers