Abstract
In this study, two-phase non-Newtonian turbulent fluid flow in an inclined geothermal pipe with chemical reaction was considered. The governing nonlinear partial differential equations derived were solved numerically using Finite Difference Method. Influence of flow parameters on the temperature, concentration and velocity profiles were analyzed graphically. From the mathematical analysis of the model, it was established that the chemical reaction parameter significantly influenced the concentration distribution in both gaseous and liquid phases. Findings further revealed that decreasing the chemical reaction parameter resulted in decreased concentration of the geothermal fluid, which causes corrosion of geothermal pipes. These findings provide important information to engineers and researchers in making better decisions in terms of design, sizing and maintenance of flow systems in geothermal pipes.
1. Introduction
Corrosion of geothermal pipes is one of the causes of inefficiency in production of electricity. It is important for the geothermal plants to minimize the losses arising from these corrosion so as to ensure maximum production of electricity to serve adequately the growing population. The geothermal fluids pick up metals and minerals which lead to corrosion and scaling of geothermal pipes and tanks, which these fluids come into contact with.
Two-phase flows occur in geothermal pipes due to the variation of density of the fluid. Two-phase turbulent stratified flows find practical applications in geothermal power plants, nuclear reactors, oil and gas pipelines and refrigeration equipment. A turbulent fluid flow in a pipe represents real-life flow frequently encountered in engineering applications like plants for generating power. For power generation to be optimized through a stable operation in a geothermal power plant, it is important to secure the flow channel in pipe conduits [1]. Ojiambo et al. [2] investigated two-phase Jeffrey Hammel flow in a geothermal pipe. They studied a two dimensional, incompressible lamina flow by considering viscosity as a non-linear function of temperature, annular flow regime and silica deposition on the geothermal pipes. The resulting governing equations of the flow were solved using BVP4c collocation method. Their findings revealed that Reynolds’ number had a significant effect in the gaseous phase as compared to the liquid phase and there was an increase on the temperature gradient due to the increase in Eckert number.
Different flow regimes occur in horizontal and vertical pipes e.g annular, plug and stratified flow regimes. Palsson et al. [3] examined the behavior of a two phase flow from geothermal wells where different models were used to determine two phase flow regimes. Pressure drop models were compared with the actual measurements. They concluded that the most common flow regime in a horizontal pipe is stratified wavy flow which depends on the flow velocity and the void fraction. Also, Li et al. [4] examined a flow structure and flow regime transitions of downward two-phase flow in large diameter pipes. For the downward flow, the following three flow regimes were observed: cap-bubbly, churn turbulent and annular flow. While for the horizontal section, pseudo-slug , stratified and plug flow regimes were observed. They found out that transition between annular flow and churn turbulent flow occurred at a certain superficial liquid velocity.
Eghbali et al. [5] investigated a two-phase fluid flow and heat transfer in a vertical geothermal, oil and gas wells where heat transfer and two-phase flow models were used to predict temperature and pressure profiles. A bubbly flow regime was also considered and the equations of continuity, momentum and energy were solved. For a two phase system of air-water, the production pressure was higher compared to a single phase system. Hou et al. [6] used a modified heat transfer model to predict the temperature of a two-phase geothermal fluid flow in the production pipelines. To predict the behavior of the flow, a Computational Fluid Dynamics (CFD) model was used. The equations of continuity, momentum and energy were considered and the flow was a compressible, one dimensional turbulent flow. The turbulent kinetic energy and the velocity profiles increased across the flow meter but the temperature profile decreased across the flow meter.
Zhao et al. [7] investigated a geothermal two-phase flow. The Seventh power law was used for the new proposed void fraction correlation model for the two-phase velocity distribution. In geothermal pipes, the new void fraction correlation gave a good agreement between the measured and predicted two phase pressure drops. A mathematical modeling of silica scaling in geothermal wells was carried out by [8]. They came up with a mathematical model where the two-phase pressure drop well-bore fluid temperature and solubility temperature correlation were integrated. Based on their results, precipitation of silica scaling increased with decrease in temperature.
The fluid flow in geothermal pipes is turbulent and different models can be used to model the turbulence. Duan et al. [9] carried out a numerical analysis of a two phase stratified smooth turbulent flow through a circular pipe. They studied a steady, two-dimensional momentum equation and employed a low Reynolds turbulence model. The non-linear equations were solved using the Newtons-Raphson Method. The findings of the study showed that the profiles of and near the interface were higher in the gas phase compared to the liquid phase. It was also observed that between the interface and the near wall of the pipe, the interface caused less fluid flow resistance and it behaved like a moving wall. Abdulwahid et al. [10] analyzed unsteady two phase turbulent flow patterns and void fraction in a vertical and horizontal pipes where ANSYS FLUENT program with Volume of fluid(VOF) model was used to investigate unsteady turbulent flow. RNG model was also used to solve the turbulent fluid flow with annular and churn flow regimes, whereas (realizable) model was used to solve slug and bubble flow regimes. They concluded from the results that transition among the flow patterns depended on the air superficial velocity. Menge [11] investigated a steady 3D turbulent flow in a pipe using Computational Fluid Dynamics. They compared two models ( and models) and found out that model gave better approximations for the center-line velocities and is suitable to predict turbulent flow in a pipe.
Corrosion of geothermal pipes plays a significant part in the long term operation and stability of the Geothermal power plants. The chemical problems can be managed through chemical modeling. Effects of chemical reactions have been reported in literature by various studies. Akter et al. [12] investigated a boundary layer mass transfer through an inclined plate with the effects of thermal diffusion and chemical reaction. They found out that there was a decrease in concentration distribution and velocity distribution as the reaction parameter was increased, but the temperature profile remained unchanged for different values of the reaction parameter. Mondal et al. [13] investigated Soret-Dufour and thermophoresis on magneto-hydrodynamic mixed convective mass and heat transfers of a semi-infinite plate in the presence of chemical reaction and uniform heat source. They found out that there was retardation in concentration distribution at the boundary layer due to increase in the thermophoretic number. The concentration distribution decreased with increase in the Schmidt number. The concentration profiles decreased with increase in chemical parameter. Rasool et al. [14] considered an incompressible mixed convection of a second grade nanofluidic viscous flow past a heated vertical Riga plate. The findings showed that concentration of the nanoparticles was enhanced by the chemical reaction parameter.
The above studies on geothermal pipes have generally neglected chemical reaction effects, variable thermal conductivity and inclined geothermal pipe. In many geothermal plants, chemical reaction arise due to salts/metals present in the geothermal fluid.
The motivation for the present work covers the following novel aspects: Firstly, formulation of a viscous, Non-Newtonian, incompressible, two-phase stratified turbulent fluid flow in a geothermal pipe with chemical reaction. Geothermal pipe is inclined at an angle due to the terrain. Thermal conductivity is temperature dependent and a first order chemical reaction is considered. Secondly, the solution of the non-linear PDEs using Finite Difference Method and its implementation in MATLAB. Finally, analysis of the effect of flow parameters on the flow variables were examined graphically. Table 1 represents nomenclature used in the study.
2. Mathematical Model
A stratified two layer flow of an incompressible fluid in a Geothermal pipe inclined at an angle as shown in Figure 1 is considered. Due to the inclination of the pipe, buoyancy forces are taken into account. The liquid phase is placed at the bottom of the pipe while the gaseous phase is at the top. The fluid flow is turbulent and unsteady. A non-linear viscosity which is a function of temperature and tangential direction in both the liquid phase and the gaseous phase is taken into consideration in the momentum equation. In the energy equation, a temperature dependent thermal conductivity is considered together with the effects of viscous dissipation. The flow system is modelled using cylindrical polar coordinate and the - axis is taken along the pipe axis.
From the governing equations presented in [15] and also taking into consideration above assumptions, the specific equations governing the flow are expressed as follows:
Continuity equation
Momentum equations
Equation (2) represents the momentum equation in the radial direction.
Equation (3) represents the momentum equation in the -direction.
The momentum equation in the direction is given by (4).
Energy equation
Concentration equation
From the equations above, is the coefficient of viscosity, is density of the fluid, is the pressure, is gravitational force, is the angle of inclination, is coefficient of thermal expansion, is coefficient of mass expansion, is diffusion coefficient, is the chemical reaction parameter, is the coefficient of thermal conductivity and is the specific heat at constant pressure.
Using the power law model as presented in [16], the viscosity can be expressed as follows:
From equation (7), is the flow consistency index and represents velocity gradient.
Using the piece-wise function for a two-phase flow as presented in [2], a modified version of the function is given as:
Where represents the viscosity for gaseous phase and the liquid phase. and are constants. Thermal conductivity as presented in [17] is a variable which depends on temperature and it is defined as: where is the thermal conductivity of the fluid and is a small parameter which takes negative values for gases and positive values for liquids.
2.1. Reynolds Decomposition and Turbulence Modeling
The fluid flow under consideration is turbulent. Before decomposing equations (1), (2), (3), (4), (5), (6), we substitute equations (9), (8) and (7) into equations (1), (2), (3), (4), (5), (6) to simplify them and then apply Reynolds decomposition rules to obtain the Reynolds Averaged Navier Stokes equations(RANS). For the RANS equations to be solved, the Reynolds stresses are expressed in the mean flow quantities and the following modeling equations are used [18]:
is the Kronecker delta;
Substituting the modeling equations (10) and (11) into RANS equations, we have the following final equations: where is the time average velocity in the radial direction, is the time average velocity in the -direction, is time-average velocity in the axial direction, is the time-average pressure and is the time-average temperature and is the time-average concentration. where is the viscous dissipation.
2.2. Turbulence Kinetic Energy and Dissipation equations
According to [19], to derive the Turbulent Kinetic Energy (TKE) the following steps are taken: Firstly, the momentum equations are replaced with the mean and fluctuating velocities; secondly, each equation is multiplied by the fluctuating velocities in their respective directions; thirdly, take time average of the decomposed equations; and lastly, apply the Reynolds decomposition rules. The equations are put together by adding them to form a single equation and this summation forms the turbulent kinetic energy equation. According to [19, 20], the TKE equation in vector form is given as:
The first term on the LHS of (21) represents the local change of the turbulent kinetic energy, second term is the convection of turbulent kinetic energy. The first term on the RHS represents diffusion of turbulent energy, the second term represents the production of turbulent kinetic energy and the last term represents the dissipation of turbulent kinetic energy. is the turbulent Prandtl number for kinetic energy.
The dissipation equation in vector form as given in [9] is: where and are turbulent model wall functions used to correct the behavior of eddy viscosity.
After substituting (19) and (8) into (22) and (21) and simplifying, the final forms are given as:
Equations (23) and (24) represent TKE for the gaseous and liquid phases respectively.
Equations (25) and (26) represent viscous dissipation equations for the gaseous and liquid phases respectively. , , , and are model constants taken from [21]. Their values are given as: , , , , .
The relationship between dissipation(), turbulent kinetic energy() and turbulent viscosity() is given as:
As a result of turbulent eddies, there is damping near the walls. The following functions have been introduced as presented in the work of [21]: where , and is the shortest distance to the interface or wall. and are turbulent Reynolds numbers.
2.3. Boundary Conditions
With the assumption that there was no slip condition at the walls of the pipe, the following boundary conditions were applied successfully by [9, 19]:
At ;
At (the interface is treated as a moving wall);
At ;
At ;
At , the gradients of all variables in the flow direction are zero as suggested by [22];
2.4. Non-dimensionalization
In order to solve the problem, equations (12), (13), (14), (15), (16), (17), (18), (20), (23), (24), (25) and (26) are converted to non-dimensionless forms. The following non-dimensional variables have been employed for non-dimensionalization as presented in [23]:
Substituting these non-dimensional variables into equations (12), (13), (14), (15), (16), (17), (18), (20), (23), (24), (25) and (26), we obtain the following dimensionless governing equations: equations (35) and (36) represent the dimensionless form of the momentum equation in the radial direction for the gaseous and liquid phase respectively.
equations (37) and (38) represent the dimensionless form of the momentum equation in the -direction for the gaseous and liquid phase respectively. equations (39) and (40) represent the dimensionless form of the momentum equation in the axial direction for the liquid and gaseous phase respectively. equation (41) represents the dimensionless form of the Energy equation. equation (42) represents the dimensionless form of the Concentration equation. equation (43) represents the dimensionless form of the turbulent kinetic energy equation. equation (44) represents the dimensionless form of the dissipation equation. equation (45) represents the dimensionless form of the relationship between dissipation, turbulent kinetic energy and turbulent viscosity.
The following are non-dimensional numbers obtained from the equations above: Reynolds Number, ; Thermal Grasholf Number, ; Solutal Grasholf Number, ; Schmidt Number, ; Chemical Reaction Parameter, ; Eckert Number, and Prandtl Number, .
2.4.1. Dimensionless Boundary Conditions
At
At
At
At
At
3. Methodology
Equations (35)–(44) are non-linear PDEs, hence cannot be solved analytically. An explicit Finite Difference Method (FDM) has been used to solve the non-linear partial differential equations governing the fluid flow. Finite Difference methods are accurate, stable, of rapid convergence and simple in solving Partial differential equations. In implementing the method, the partial derivatives in PDEs are replaced by approximations based on Taylor series expansions near the point of interest. When FDM is applied, the continuous domain is discretized and the differential terms of the equation are converted to linear algebraic equations.
The domain is discretized using a cylindrical mesh. Let be described at the point and for simplicity the point is written as and as . Assume we have points along , points along , points along and points along directions to form the mesh. Let the step size be along , be along , along and along .
, , and where , , and
For our set of equations, a second order Finite difference approximation(forward in time and central in space) scheme would be used. These approximations are given as:
Substituting equations (51)–(53) into equations (35)–(44), we have the following set of finite difference equations: equations (54) and (55) are radial momentum equations in discrete form for the liquid and gaseous phases respectively. equations (56) and (57) represent tangential momentum equation in discrete form for the liquid and gaseous phases respectively. equations (58) and (59) represent axial momentum equations in discrete form for liquid and gas phases respectively. equation (60) represents energy equation in discrete form. equation (61) represents Concentration equation in discrete form. equations (63) and (62) are Turbulent Kinetic Energy equations in discrete form for liquid and gaseous phase respectively.
Equations (64) and (65) are dissipation equations in discrete form for gas and liquid phases respectively. Equations (54)–(65) are rearranged in a way that the quantities at time are on the left hand side (since they depend explicitly on the time at ), before being implemented in MATLAB.
4. Results and Discussions
This section deals with the physical explanation of the effect of the various flow parameters on the velocity and temperature profiles. Computations were carried out for: (Air-gas phase) and (water-liquid phase).
Figure 2 displays the effect of Reynolds number on velocity profile. It is important to note that as the Reynolds number increases the velocity distribution for both the gaseous and liquid phase increase. Moreover, the momentum boundary layer thickness increases. Physically, increase in the Reynolds number results in the reduction of the viscous forces which oppose the fluid flow. Therefore, there is an increase in the velocity profile.
(a) Gaseous phase
(b) Liquid Phase
Figure 3 depicts the influence of Solutal Grasholf number on the velocity profile. It is worth noting that the velocity profiles decrease as solutal Grashof number increases for both phases. Physically, increase in the solutal Grasholf Number leads to an increase in the viscous forces which resist the fluid flow. Consequently, the velocity profile decreases.
(a) Gaseous phase
(b) Liquid Phase
Effect of Thermal Grasholf number on velocity profile is analyzed in Figure 4 for both liquid and gas phases. There is an increase in velocity as Thermal Grashof number increases. An increase in Thermal Grasholf number leads to enhancement of buoyancy forces and there is also an addition of thermal energy into the molecules of the fluid which loosen up the inter-molecular forces of the fluid particles. Thus, there is an increase in the velocity profiles.
(a) Gaseous phase
(b) Liquid Phase
Figure 5 shows the effect of the angle of inclination, , on the velocity profiles for both phases. There is a decrease in the velocity distribution as the angle of inclination increases. As increases, the magnitude of decreases. Thus, the effect of buoyancy forces decrease along the radial direction. Consequently, there are decreased velocity profiles along the radial direction.
(a) Gaseous phase
(b) Liquid Phase
Effect of Reynolds number on the temperature profiles for both gas and liquid phases are displayed in Figure 6. It is important to note that temperature decreases with increase in Reynolds number for the gaseous phase but it increases with increase in Reynolds number for the liquid phase. As increases the viscous forces become less important. The viscous forces are directly proportional to temperature for gases but inversely proportional to temperature in liquids. Therefore, an increase in Reynolds number results in an increase in temperature profile for the liquid phase but a decrease in temperature profile for the gaseous phase.
(a) Gaseous phase
(b) Liquid Phase
Figure 7 displays the effect of Eckert number on the temperature profiles. Here, temperature profiles increase for both phases as Eckert number increases. As increases, the kinetic energy of the molecules increases, hence increased vibration of molecules which leads to conversion of kinetic energy to heat energy. Consequently, an increase in Eckert number leads to an increase in the temperature of the fluid. There is also an additional heating caused by viscous dissipation.
(a) Gaseous phase
(b) Liquid Phase
From Figure 8, an increase in leads to increase in the temperature profiles for both phases. An increase in results in an increase in the thermal boundary layer thickness. Thus, there is increased temperature distribution.
(a) Gaseous phase
(b) Liquid Phase
Figure 9 portrays the effect of Thermal Grasholf Number on the temperature field. Temperature distribution decreases with increase in Thermal Grasholf number. Increase in Thermal Grasholf number results in an enhancement of buoyancy forces. Buoyancy forces are inversely proportional to temperature of the fluid. Consequently, increase in Thermal Grasholf number leads to reduction in temperature distribution.
(a) Gaseous phase
(b) Liquid Phase
Figure 10 displays the effect of Schmidt number on concentration profiles for the gas and liquid phases respectively. It is worth noting that increase in Schmidt number leads to a decrease in concentration profile. A high number implies a decrease in molecular diffusivity which causes a reduction in the concentration boundary layer thickness. Consequently, increase in Schmidt number results to a decrease in the concentration.
(a) Gaseous phase
(b) Liquid Phase
Effect of Chemical Reaction Parameter on concentration is analyzed in Figure 11 for gas and liquid phases. An increase in Chemical reaction parameter increases the concentration boundary layer thickness. Thus, there is an increase in concentration distribution as chemical reaction parameter increases.
(a) Gaseous phase
(b) Liquid Phase
4.1. Validation Of Results
In validating the presented results, the computed results are compared with the results obtained by [9] as shown in Figure 12 which proved to be in good agreement. According to [9] the horizontal velocity profile of the gaseous phase was symmetrical to the center-line which were in good agreement with the experimental data by [24] which confirmed the validity of the method.
(a) Velocity profile for Gaseous phase (Present work)
![(b) Duan [9] velocity profile](https://static.hindawi.com/articles/jam/volume-2022/7617017/figures/7617017.fig.0012b.jpg)
(b) Duan [9] velocity profile
5. Conclusion
In this paper, two-phase turbulent fluid flow in a geothermal pipe with chemical reaction was investigated and effects of various flow parameters discussed. The results for various variations of flow parameters were presented in graphical form. The following important deductions were got from the results: (i)Fluid flow was enhanced with increase in Reynolds number and Thermal Grasholf number;(ii)Higher values of Reynolds number resulted in reduction of temperature profile for gas phase but enhancement of temperature profile for the liquid phase;(iii)Higher values of Eckert number and angle of inclination led to increase of temperature of the fluid;(iv)The velocity profile reduced for higher values of Mass Grasholf Number;(v)Higher values of Thermal Grasholf Number resulted in reduction of temperature of the fluid;(vi)Higher values of Schmidt Number resulted in reduction of concentration of the fluid;(vii)Concentration was enhanced with increase in Chemical reaction Parameter.
For future research, a two-phase stratified wavy flow in presence of an inclined variable magnetic field can be studied. Also a different turbulence modeling approach can be considered.
Data Availability
No data were used to support this study
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
The author(s) extend their appreciation to the Pan African University for Basic Sciences and Technology for funding this research.