Heavy Oil Laminar Flow in Corrugated Ducts: A Numerical Study Using the Galerkin-Based Integral Method

: Fluid ﬂow in pipes plays an important role in di ﬀ erent areas of academia and industry. Due to the importance of this kind of ﬂow, several studies have involved circular cylindrical pipes. This paper aims to study fully developed internal laminar ﬂow through a corrugated cylindrical duct, using the Galerkin-based integral method. As an application, we present a study using heavy oil with a relative density of 0.9648 (14.6 ◦ API) and temperature-dependent viscosities ranging from 1715 to 13000 cP. Results for di ﬀ erent ﬂuid dynamics parameters, such as the Fanning friction factor, Reynolds number, shear stress, and pressure gradient, are presented and analyzed based on the corrugation number established for each section and aspect ratio of the pipe.


Introduction
The oil industry is a major player in the energy sector worldwide. The increased demand for light oil reserves in recent decades has led to their depletion. As a result, heavy oil is increasingly a topic of interest, receiving considerable attention regarding its efficient transportation. Although there is a large amount of productive crude oil reserves, heavy oil transportation is a complex and expensive task due to its high viscosity, which requires a large amount of energy to pump [1].
One of the major concerns in the production and transportation of heavy oil is related to load loss or pressure drops, which are associated with high flow costs. This occurs due to the fluid friction effect on the inner walls of the duct, making its transport difficult, and increasing the refining cost [2].

Mathematical Formulation
In this work, we choose the cylindrical coordinate system (r, θ, z) for the mathematical formulation of the physical problem, and the following assumptions are considered: (a) The laminar flow is fully developed, isothermal, single-phase, and steady-state; (b) The cross-sectional area of the tube, along the z-axis, is constant; (c) The properties of the fluid, thermal and physical, are considered constant; (d) There is a condition of not-slipping on the duct wall.

The Geometry
The geometry to be analyzed is a corrugated cylindrical duct. For a duct with corrugated geometry we consider, in polar coordinates, the parameterization of the cross-section, r(θ) = a + b sin(N θ), where a and b are shown in Figure 1 and N is the number of corrugations, sine-type. The relationship a b is the so-called aspect ratio, Energies 2020, 13, x FOR PEER REVIEW 3 of 15 and/or velocity in ducts, we can cite works with circular, rectangular, isosceles triangular, right triangular, and annuli cross-sections [2224].

Mathematical Formulation
In this work, we choose the cylindrical coordinate system ( , , ) for the mathematical formulation of the physical problem, and the following assumptions are considered: a) The laminar flow is fully developed, isothermal, single-phase, and steady-state; b) The cross-sectional area of the tube, along the z-axis, is constant; c) The properties of the fluid, thermal and physical, are considered constant; d) There is a condition of not-slipping on the duct wall.

The Geometry
The geometry to be analyzed is a corrugated cylindrical duct. For a duct with corrugated geometry we consider, in polar coordinates, the parameterization of the cross-section, ( ) = + ( ), where and are shown in Figure 1 and is the number of corrugations, sine-type. The relationship is the so-called aspect ratio, Thus, considering the cross-section of the duct with the parameterization as reported before, the cross-sectional area of the duct is given as follows: For the perimeter of the cross-section of the duct, we use the definition of arc length as follows:

Momentum Equation
The momentum equation for the domain and assumptions considered in a cylindrical coordinate system are given by: In the equation, ( , ) is the velocity component along the z -axis and annulled in the boundary, μ is the dynamic velocity of the fluid, and is the pressure gradient.
By considering the dimensionless variables Thus, considering the cross-section of the duct with the parameterization as reported before, the cross-sectional area of the duct is given as follows: For the perimeter of the cross-section of the duct, we use the definition of arc length as follows:

Momentum Equation
The momentum equation for the domain and assumptions considered in a cylindrical coordinate system are given by: Energies 2020, 13, 1363 4 of 13 In the equation, u(r, θ) is the velocity component along the z-axis and annulled in the boundary, µ is the dynamic velocity of the fluid, and dp dz is the pressure gradient. By considering the dimensionless variables and where L = a is the characteristic length, we can write Equation (3) in the dimensionless form and apply it to the hydrodynamically fully developed flow as follows: The normalized dimensionless mean velocity of the fluid can be obtained as follows: where u m and W m are the mean and dimensionless mean velocity, respectively.

Solution Methodology
In this work, we have employed the Galerkin-based integral (GBI) method on the MAPLE 17 platform for finding the numerical solution of Equation (6). The solution is approximated as a linear combination of a set of base functions. Thus, for Equation (6), we can write: The bases functions f 1 , f 2 , . . . , f n are linearly independent functions and satisfy the same homogeneous boundary conditions of W.
Using Equation (8) and applying the GBI method in Equation (6), we obtain the following matrix system: with and The coefficients d 1 , d 2 , . . . , d n can be determined as follows: In this way, the mean velocity is given by: Energies 2020, 13, 1363 5 of 13 One of the flow parameters commonly used in practice is the Fanning friction factor, f , defined by: where τ w = d h 4 dp dz represents the shear stress and d h = 4A c P the hydraulic diameter. The Reynolds number is defined as: With these results, we can determine values for the Poiseuille number, which is defined as the product of the Reynolds number and the friction factor as follows: with D h = d h L as the dimensionless hydraulic diameter. Considering the following parameterization in polar coordinates and the set of bases functions of the form, we can write the base functions as follows: where the subscript i represents the i-th term of the set of base functions in Equation (18).

Heavy Oil Flow Application
The oil is classified as heavy with • API = 14.6, a relative density of 0.9648 (20/4 • C), and the oil viscosity is reported in Table 1. Thus, we have the value for the fluid density ρ = 952.17 kg/m 3 . The cross-section of the duct does not change along the longitudinal direction (z-axis), and we consider an oil mean velocity of u m = 1 m/s.

Results and Discussion
In this section, first, we present the numerical results for the Poiseuille number f Re obtained for the different aspect ratio, β, and the number of corrugations, N, of the corrugated cross-section duct. The effect of the number of corrugations in the shape of the pipe can be observed in Figure 2.
Energies 2020, 13, x FOR PEER REVIEW 6 of 15 To obtain computational gains, a series approximation is performed for the function that represents the parameterization of the duct cross-section. For example, with =0.08 and =12, the following approximation can be used: 1 + 0.08 sin (12 ) where is in radian. However, no changes are made to the area values and perimeter of the duct. This approximation can be given at least in the first quadrant, because there is a periodicity in this type of geometry. The approximation reported in the series (Equation (20)) can be seen in Figure 3. To obtain computational gains, a series approximation is performed for the function that represents the parameterization of the duct cross-section. For example, with b a = 0.08 and N = 12, the following approximation can be used: where θ is in radian. However, no changes are made to the area values and perimeter of the duct. This approximation can be given at least in the first quadrant, because there is a periodicity in this type of geometry. The approximation reported in the series (Equation (20)) can be seen in Figure 3.  The values for different cases are presented in Table 2 and compared with the values reported in the literature [13]. Hu and Chang [11] employed the conformal mapping method and Green's functions. It is also important to note that unlike the literature, which uses numerical approximations for the perimeter value, in our work we use the arc length function. The good choice of the parameterization and the number of base functions created a good agreement between the predicted values for this research and the results reported in the literature. In some cases, the choice of the number of base functions from 3 to 52, as can be seen in Table 3, produces a significant increase in the computational time. The minimum computational time for some values was on average 1 minute, while the average maximum time was 2 hours and 15 minutes. Further, it is also observed that the numerical solution with the GBI method is in good agreement with the values reported in the literature. When the aspect ratio value tends to zero, this approximates the predicted Poiseuille number value for the values =16 obtained for a circular cross-section pipe, validating the mathematical procedure used in this work.
The solution of Equation (7) makes it possible to find the dimensionless velocity distribution, as illustrated in Figure 4:  Table 2 and in Figure 5, for better visualization of this parameter.
In the second part of the results, we present values for the Fanning friction factor, Reynolds number, shear stress, and pressure gradient in corrugated cross-sectional ducts with different aspect The f Re values for different cases are presented in Table 2 and compared with the values reported in the literature [13]. Hu and Chang [11] employed the conformal mapping method and Green's functions. It is also important to note that unlike the literature, which uses numerical approximations for the perimeter value, in our work we use the arc length function. The good choice of the parameterization and the number of base functions created a good agreement between the predicted values for this research and the results reported in the literature. In some cases, the choice of the number of base functions from 3 to 52, as can be seen in Table 3, produces a significant increase in the computational time. The minimum computational time for some values was on average 1 minute, while the average maximum time was 2 hours and 15 minutes. Further, it is also observed that the numerical solution with the GBI method is in good agreement with the values reported in the literature. When the aspect ratio value tends to zero, this approximates the predicted Poiseuille number value for the values f Re = 16 obtained for a circular cross-section pipe, validating the mathematical procedure used in this work.
The solution of Equation (7) makes it possible to find the dimensionless velocity distribution, as illustrated in Figure 4:  The values for different cases are presented in Table 2 and compared with the values reported in the literature [13]. Hu and Chang [11] employed the conformal mapping method and Green's functions. It is also important to note that unlike the literature, which uses numerical approximations for the perimeter value, in our work we use the arc length function. The good choice of the parameterization and the number of base functions created a good agreement between the predicted values for this research and the results reported in the literature. In some cases, the choice of the number of base functions from 3 to 52, as can be seen in Table 3, produces a significant increase in the computational time. The minimum computational time for some values was on average 1 minute, while the average maximum time was 2 hours and 15 minutes. Further, it is also observed that the numerical solution with the GBI method is in good agreement with the values reported in the literature. When the aspect ratio value tends to zero, this approximates the predicted Poiseuille number value for the values =16 obtained for a circular cross-section pipe, validating the mathematical procedure used in this work.
The solution of Equation (7) makes it possible to find the dimensionless velocity distribution, as illustrated in Figure 4:  Table 2 and in Figure 5, for better visualization of this parameter.
In the second part of the results, we present values for the Fanning friction factor, Reynolds number, shear stress, and pressure gradient in corrugated cross-sectional ducts with different aspect  Table 2 and in Figure 5, for better visualization of this parameter.
In the second part of the results, we present values for the Fanning friction factor, Reynolds number, shear stress, and pressure gradient in corrugated cross-sectional ducts with different aspect ratios changing from β = 0.02 to 0.06, and the number of corrugations N = 8, 12, 16 and 24. In Tables 4-11, these results are shown.   When analyzing the results presented in Tables 4-11, it is observed that the flow of the considered oil has the following behavior: (a) the higher the temperature, the lower the oil viscosity, and (b) the higher the Reynolds number, the lower the friction factor, shear stress, and pressure gradient. The flow remains laminar in all variations of the pipe geometry. Also, when comparing    Table 5. Values of the fluid dynamic parameters obtained in the simulations with β = 0.02 and N = 12.  Table 7. Values of the fluid dynamic parameters obtained in the simulations with β = 0.02 and N = 24.   Table 9. Values of the fluid dynamic parameters obtained in the simulations with β = 0.06 and N = 12. When analyzing the results presented in Tables 4-11, it is observed that the flow of the considered oil has the following behavior: (a) the higher the temperature, the lower the oil viscosity, and (b) the higher the Reynolds number, the lower the friction factor, shear stress, and pressure gradient. The flow remains laminar in all variations of the pipe geometry. Also, when comparing values of the aspect ratio, for example, β = 0.02 and 0.06, and considering the same number of corrugations N, it was observed that the higher the aspect ratio, the greater the pressure gradient. Similarly, if we set a value for the aspect ratio, β, the pressure gradient increases as the number of corrugations N increases. Such statements can be seen in Figure 6.
(c) (d) When analyzing the results presented in Tables 4-11, it is observed that the flow of the considered oil has the following behavior: (a) the higher the temperature, the lower the oil viscosity, and (b) the higher the Reynolds number, the lower the friction factor, shear stress, and pressure gradient. The flow remains laminar in all variations of the pipe geometry. Also, when comparing values of the aspect ratio, for example, = 0.02 and 0.06, and considering the same number of corrugations , it was observed that the higher the aspect ratio, the greater the pressure gradient. Similarly, if we set a value for the aspect ratio, , the pressure gradient increases as the number of corrugations increases. Such statements can be seen in Figure 6.

Conclusions
From the obtained results, the following conclusions can be cited: 1. The Galerkin-based integral method assisted by symbolic manipulation software is an effective tool for investigating the fully developed laminar flow fluid in a corrugated cross-section cylindrical duct. Fully developed velocity profiles for any duct shape can be predicted by this method.

Conclusions
From the obtained results, the following conclusions can be cited: 1.
The Galerkin-based integral method assisted by symbolic manipulation software is an effective tool for investigating the fully developed laminar flow fluid in a corrugated cross-section cylindrical duct. Fully developed velocity profiles for any duct shape can be predicted by this method.

2.
The Poiseuille number, f Re, obtained for different aspect ratios of the pipe was compared with the existing literature, and a good concordance was obtained in all studied cases. 3.
The appropriate number of base functions is different depending on the aspect ratio and the number of corrugations of the studied geometry. Such variation altered the computational time significantly. The smaller the number of base functions, the lower the computational time.

4.
The pressure gradient increases with the number of corrugations and the aspect ratio of the pipe and decreases with the increase in the fluid temperature.

5.
Values for the Reynolds number, Fanning friction factor, shear stress, and pressure gradient for heavy oil flow, with temperature-dependent viscosity, were presented. These results make it possible to better understand the behavior of heavy oil in corrugated cross-section ducts.