# Two-Phase Turbulent Fluid Flow in a Geothermal Pipe with Chemical Reaction

Authors: Nyariki M. Eric; Mathew N. Kinyanjui; Jeconia O. Abonyo
Journal: Journal of Applied Mathematics
 (2022)

---

## 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.


---

## Body

## 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 k 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 k−ε (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.Table 1
Nomenclature.
SymbolMeaningSI UnitsgAcceleration due to gravityms−2pPressure forceNm−2r,θ,zCylindrical coordinate variablesu⟶Velocity vectorms−1FBody forces tensor[N]ρFluid densitykgm−3μCoefficient of viscositykgm−3κCoefficient of thermal conductivityWm−1k−1F⟶BBuoyancy forces[N]αAngle of inclinationCpSpecific heat at constant pressureJkg−1k−1DmDiffusion coefficientTTemperature of the fluid[K]T∞Ambient temperature of the fluid[K]tTime[s]μ0Flow consistency index[m2/s]μTEddy viscosityνTKinematic Eddy Viscosityk∗Coefficient of thermal conductivityσkTurbulent Prandtl number for kinetic EnergyαTTurbulent thermal diffusivityαcTurbulent mass diffusivityσεTurbulent Prandtl number for dissipation EnergyReReynolds numbergVelocity gradient


## 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 r,θ,z and the z- axis is taken along the pipe axis.Figure 1
Flow Configuration for turbulent two-phase flow inclined at an angle.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(1)1r∂ru∂r+1r∂v∂θ+∂w∂z=0.Momentum equations(2)ρ∂u∂t+u∂u∂r+vr∂u∂θ+w∂u∂z−v2r=−∂p∂r+μ1r∂∂rr∂u∂r−ur2+μ1r2∂2u∂θ2+∂2u∂z2−2r2∂v∂θ+∂μ∂θ1r2∂u∂θ−vr2+1r∂v∂r+ρgcosαβTT−T∞+βcC−C∞.Equation (2) represents the momentum equation in the radial direction.
(3)ρ∂v∂t+u∂v∂r+vr∂v∂θ+w∂v∂z+uvr=−1r∂p∂θ+∂μ∂θ2r2∂v∂θ+2ur2+μ−vr2+1r∂∂rr∂v∂r+1r2∂2v∂θ2+2r2∂u∂θ+∂2v∂z2.Equation (3) represents the momentum equation in the θ-direction.
(4)ρ∂w∂t+u∂w∂r+vr∂w∂θ+w∂w∂z=−∂p∂z+∂μ∂θ1r2∂w∂θ+1r∂v∂z+μ1r∂∂rr∂w∂r+1r2∂2w∂θ2+∂2w∂z2+ρgsinαβTT−T∞+βcC−C∞.The momentum equation in thez direction is given by (4).Energy equation(5)ρCp∂T∂t+u∂T∂r+vr∂T∂θ+w∂T∂z=∂κ∂T∂T∂r2+1r2∂κ∂T∂T∂θ2+∂κ∂T∂T∂z2+κT1r∂T∂r+∂2T∂r2+1r2∂2T∂θ2+∂2T∂z2+Φ.Concentration equation(6)∂C∂t+u∂C∂r+vr∂C∂θ+w∂C∂z=Dm1r∂∂rr∂C∂r+1r2∂2C∂θ2+∂2C∂z2−krC−C∞.From the equations above,μ is the coefficient of viscosity, ρ is density of the fluid, p is the pressure, g is gravitational force, α is the angle of inclination, βT is coefficient of thermal expansion, βC is coefficient of mass expansion, Dm is diffusion coefficient, Kr is the chemical reaction parameter, κ is the coefficient of thermal conductivity and Cp is the specific heat at constant pressure.Using the power law model as presented in [16], the viscosity can be expressed as follows:
(7)μ=μ0gn−1,g=gθ.From equation (7), μ0 is the flow consistency index and g 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:
(8)gθ,T=θs1+γT−T∞forγ<0,s≠0.θs1+γT−T∞forγ>0,s≠0.Whereγ<0 represents the viscosity for gaseous phase and γ>0 the liquid phase. γ and s are constants. Thermal conductivity κ as presented in [17] is a variable which depends on temperature and it is defined as:
(9)κ=κ∗1+bT−T∞ΔT,where κ∗ is the thermal conductivity of the fluid and b 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]:
(10)−u′u′¯=νT∇u+∇uT−23kδij,δij is the Kronecker delta;
(11)−u′T′¯=αT∂T¯∂r,−v′T′¯=αTr∂T¯∂θ,−w′T′¯=αT∂T¯∂z,−u′C′¯=αc∂C¯∂r,−v′C′¯=αcr∂C¯∂θ,−w′C′¯=αc∂C¯∂z.Substituting the modeling equations (10) and (11) into RANS equations, we have the following final equations:
(12)∂u¯∂t+u¯∂u¯∂r+v¯r∂u¯∂θ+w¯∂u¯∂z−v¯v¯r=−1ρ∂p¯∂r+μ0gn−1ρ+νT1r∂∂rr∂u¯∂r−u¯r2+1r2∂2u¯∂θ2+∂2u¯∂z2−2r2∂v¯∂θ+μ0ρsn−1gn−2θn−11+γT¯−T∞1r2∂u¯∂θ−v¯r2+1r∂v¯∂r+gcosαβTT¯−T∞+βcC¯−C∞,where u¯ is the time average velocity in the radial direction, v¯ is the time average velocity in the θ-direction, w¯ is time-average velocity in the axial direction, p¯ is the time-average pressure and T¯ is the time-average temperature and C¯ is the time-average concentration.
(13)∂u¯∂t+u¯∂u¯∂r+v¯r∂u¯∂θ+w¯∂u¯∂z−v¯v¯r=−1ρ∂p¯∂r+μ0gn−1ρ+νT1r∂∂rr∂u¯∂r−u¯r2+1r2∂2u¯∂θ2+∂2u¯∂z2−2r2∂v¯∂θ+μ0ρsn−1gn−2θn−11+γT¯−T∞1r2∂u¯∂θ−v¯r2+1r∂v¯∂r+gcosαβTT¯−T∞+βcC¯−C∞,(14)∂v¯∂t+u¯∂v¯∂r+v¯r∂v¯∂θ+w¯∂v¯∂z+uvr¯=−1ρr∂p¯∂θ+μ0ρsn−1gn−2θs−11+γT−T∞2r2∂v¯∂θ+2u¯r2μ0gn−1ρ+νT×1r∂∂rr∂v¯∂r−v¯r2+1r2∂2v¯∂θ2+∂2v¯∂z2+2r2∂u¯∂θ,(15)∂v¯∂t+u¯∂v¯∂r+v¯r∂v¯∂θ+w¯∂v¯∂z+uvr¯=−1ρr∂p¯∂θ+μ0ρsn−1gn−2θs−11+γT−T∞2r2∂v¯∂θ+2u¯r2+μ0gn−1ρ+νT1r∂∂rr∂v¯∂r−v¯r2+1r2∂2v¯∂θ2+∂2v¯∂z2+2r2∂u¯∂θ,(16)∂w¯∂t+u¯∂w¯∂r+v¯r∂w¯∂θ+w¯∂w¯∂z=1ρ∂p¯∂z+sn−1μ0gn−2θs−1ρr21+γT¯−T∞∂w¯∂θ+r∂v¯∂z+μ0gn−1ρ+νT1r∂∂r∂w¯∂r+1r2∂w¯∂θ2+∂2w¯∂z2+gsinαβTT¯−T∞+βcC¯−C∞,(17)∂w¯∂t+u¯∂w¯∂r+v¯r∂w¯∂θ+w¯∂w¯∂z=1ρ∂p¯∂z+sn−1μ0gn−2θs−1ρr21+γT¯−T∞×∂w¯∂θ+r∂v¯∂z+μ0gn−1ρ+νT1r∂∂r∂w¯∂r+1r2∂w¯∂θ2+∂2w¯∂z2+gsinαβTT¯−T∞+βcC¯−C∞,(18)∂T¯∂t+u¯∂T¯∂r+v¯r∂T¯∂θ+w¯∂T¯∂z=−κ∗ρCpbΔT∂T¯∂r2+1r2∂T¯∂θ2+∂T¯∂z2+κ∗ρCp1+bT¯−T∞ΔT+αT1r∂T¯∂r+∂2T¯∂r2+1r2∂2T¯∂θ2+∂2T¯∂z2+ϕρCp+εCp,(19)ϕ=2μ∂ur¯∂r2+1r∂uθ¯∂θ+ur¯r2+∂uz¯∂z2+μ1r∂ur¯∂θ+∂uθ¯∂r−uθ¯r2+∂uθ¯∂z+1r∂uz¯∂θ2+∂uz¯∂r+∂ur¯∂z2,where ϕ is the viscous dissipation.
(20)∂C¯∂t+u¯∂C¯∂r+v¯r∂C¯∂θ+w¯∂C¯∂z=Dm+αc1r∂∂rr∂C¯∂r+1r2∂2C¯∂θ2+∂2C¯∂z2−krC¯−C∞.


### 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:
(21)ρ∂κ∂t+u¯·ρ∇κ=∇·μ+μTσk∇κ+μTϕ−ρε.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. σk is the turbulent Prandtl number for kinetic energy.The dissipation equation in vector form as given in [9] is:
(22)ρ∂ε∂t+u¯·ρ∇ε=∇·μ+μTσε∇ε+f1Cε,1εkϕ−f2ρCε,2ε2k,where f1 and f2 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:
(23)∂κ∂t+u¯∂κ∂r+v¯r∂κ∂θ+w¯∂κ∂z=1ρr2μ0sn−1gn−2θs−11+γT¯−T∞∂κ∂θ+μ0gn−1ρ+νTσκ×1r∂∂r∂κ∂r+1r2∂2κ∂θ2+∂2κ∂z2+2νT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+νT1r∂u¯∂θ+∂v¯∂r−v¯r2+νT∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−ε,(24)∂κ∂t+u¯∂κ∂r+v¯r∂κ∂θ+w¯∂k∂z=1r2μ0ρsn−1gn−2θs−11+γT¯−T∞∂κ∂θ+μ0gn−1ρ+νTσκ1r∂∂rr∂κ∂r+1r2∂2κ∂θ2+∂2κ∂z2+2νT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+νT1r∂u¯∂θ+∂v¯∂r−v¯r2+νT∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−ε.Equations (23) and (24) represent TKE for the gaseous and liquid phases respectively.
(25)∂ε∂t+u¯∂ε∂r+v¯r∂ε∂θ+w¯∂ε∂z=1ρr2μ0sn−1gn−2θs−11+γT¯−T∞∂ε∂θ+μ0gn−1ρ+νTσε×1r∂∂rr∂ε∂r+1r2∂2ε∂θ2+∂2ε∂z2+2f1Cε,1εkνT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂u¯∂z2+f1Cε,1εkνT1r∂u¯∂θ+∂v¯∂r−v¯r2+f1νTCε,1εk∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−f2Cε,2ε2k,(26)∂ε∂t+u¯∂ε∂r+v¯r∂ε∂θ+w¯∂ε∂z=1r2μ0ρsn−1gn−2θs−11+γT¯−T∞∂ε∂θ+μ0gn−1ρ+νTσε×1r∂∂rr∂ε∂r+1r2∂2ε∂θ2+∂2ε∂z2+2f1Cε,1εkνT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+f1Cε,1εkνT1r∂u¯∂θ+∂v¯∂r−v¯r2+f1νTCε,1εk∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−f2Cε,2ε2k.Equations (25) and (26) represent viscous dissipation equations for the gaseous and liquid phases respectively. Cμ, Cε,1, Cε,2, σε and σk are model constants taken from [21]. Their values are given as: Cμ=0.09, Cε,1=1.44, Cε,2=1.92, σε=1.3, σk=1.The relationship between dissipation(ε), turbulent kinetic energy(k) and turbulent viscosity(νT) is given as:
(27)νT=Cμfμk2ε.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]:
(28)fμ=1−exp−0.0165Rek21+20.5Ret,f1=1+0.05fμ3,f2=1−exp−Ret2,where Rek=ρκ1/2dw/μt=κ1/2dw/νt, Ret=ρk2/μtε=κ2/νtε and dw is the shortest distance to the interface or wall. Ret and Rek 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]:Atr=D+/2;
(29)u=v=w=0,k=0,ε=0,T=Tw,C=Cw.Atr=0 (the interface is treated as a moving wall);
(30)C=C∞,T=T∞,u=u∞,v=u∞,w=u∞,κ=0,∂ε∂r=0.Atr=D−/2;
(31)u=v=w=0,k=0,ε=0,T=Tw,C=Cw.Atz=0;
(32)w=u∞,v=0,u=0,T=T∞,C=C∞.Atz=∞, the gradients of all variables in the flow direction are zero as suggested by [22];
(33)∂v∂z=0,∂u∂z=0,∂w∂z=0,∂κ∂z=0,∂ε∂z=0


### 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]:
(34)u∗=u¯u∞,v∗=v¯u∞,w∗=w¯u∞,t∗=tu∞D,p∗=p¯ρu∞2,θ∗=θ,z∗=zD,νT∗=νTDu∞,r∗=rD,T∗=T¯−T∞Tw−T∞,αT=αTDu∞,ε∗=εDu∞3,k∗=ku∞2,αc∗=αcDu∞,C∗=C¯−C∞Cw−C∞.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:
(35)∂u∗∂t∗+u∗∂u∗∂r∗+v∗r∗∂u∗∂θ∗+w∗∂u∗∂z∗−v∗2r∗=−∂p∗∂r+1Reθsn−11+γT∗ΔT+νT∗×1r∗∂∂r∗r∗∂u∗∂r∗−u∗r∗2+1r∗2∂2u∗∂θ∗2+∂2u∗∂z∗2−2r∗2∂v∗∂θ∗+1Resn−1θsn−1−11+γT∗ΔTn−11r∗2∂u∗∂θ∗−v∗r∗2+1r∗∂v∗∂r∗+cosαT∗GrTRe2+C∗GrCRe2,(36)∂u∗∂t∗+u∗∂u∗∂r∗+v∗r∗∂u∗∂θ∗+w∗∂u∗∂z∗−v∗2r∗=−∂p∗∂r+1Reθsn−11+γT∗ΔTn−1+νT∗×1r∗∂∂r∗r∗∂u∗∂r∗−u∗r∗2+1r∗2∂2u∗∂θ∗2+∂2u∗∂z∗2−2r∗2∂v∗∂θ∗+1Resn−1θsn−1−11+γT∗ΔTn−11r∗2∂u∗∂θ∗−v∗r∗2+1r∗∂v∗∂r∗+cosαT∗GrTRe2+C∗GrCRe2,equations (35) and (36) represent the dimensionless form of the momentum equation in the radial direction for the gaseous and liquid phase respectively.
(37)∂v∗∂t∗+u∗∂v∗∂r∗+v∗r∗∂v∗∂θ∗+w∗∂v∗∂z∗+u∗v∗r∗=−1r∗∂p∗∂θ∗+1Resn−1gn−2θs−11+γT∗Tw−T∞2r∗2∂v∗∂θ∗+2u∗r∗2+1Regn−1+νT∗1r∗∂∂r∗r∗∂v∗∂r∗−v∗r∗2+1r∗2∂2v∗∂θ∗2+∂2v∗∂z∗2+2r∗2∂u∗∂θ∗,(38)∂v∗∂t∗+u∗∂v∗∂r∗+v∗r∗∂v∗∂θ∗+w∗∂v∗∂z∗+u∗uθ∗r∗=−1r∗∂p∗∂θ∗+1Resn−1gn−2θs−11+γT∗Tw−T∞2r∗2∂v∗∂θ∗+2u∗r∗2+1Regn−1+νT∗×1r∗∂∂r∗r∗∂v∗∂r∗−v∗r∗2+1r∗2∂2v∗∂θ∗2+∂2v∗∂z∗2+2r∗2∂u∗∂θ∗,equations (37) and (38) represent the dimensionless form of the momentum equation in the θ-direction for the gaseous and liquid phase respectively.
(39)∂w∗∂t∗+u∗∂w∗∂r∗+v∗r∗∂w∗∂θ∗+w∗∂w∗∂z∗=−∂p∗∂z∗+1Resn−1gn−2θs−1r∗21+γT∗Tw−T∞∂w∗∂θ∗+r∗∂v∗∂z∗+1Regn−1+νT∗1r∗∂∂r∗r∗∂w∗∂r∗+1r∗2∂w∗∂θ∗2+∂2w∗∂z∗2+sinαT∗GrTRe2+C∗GrCRe2,(40)∂w∗∂t∗+u∗∂w∗∂r∗+v∗r∗∂w∗∂θ∗+w∗∂w∗∂z∗=−∂p∗∂z∗+1Re1r∗2sn−1gn−2θs−11+γT∗Tw−T∞∂w∗∂θ∗+r∗∂v∗∂z∗+1Regn−1+νT∗1r∗∂∂r∗r∗∂w∗∂r∗+1r∗2∂w∗∂θ∗2+∂2w∗∂z∗2+sinαT∗GrTRe2+C∗GrCRe2,equations (39) and (40) represent the dimensionless form of the momentum equation in the axial direction for the liquid and gaseous phase respectively.
(41)∂T∗∂t∗+u∗∂T∗∂r∗+v∗r∗∂T∗∂θ∗+w∗∂T∗∂z∗=−bPrRe∂T∗∂r∗2+1r∗2∂T∗∂θ∗2+∂T∗∂z∗2+1+bT∗PrRe+αT∗1r∗∂T∗∂r∗+∂2T∗∂r∗2+1r∗2∂2T∗∂θ∗2+∂2T∗∂z∗2+EcRegn−1×2∂u∗∂r∗2+1r∂v∗∂θ∗+u∗r∗2+∂w∗∂z∗2+1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+EcRegn−1∂v∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2+Ecε∗,equation (41) represents the dimensionless form of the Energy equation.
(42)∂C∗∂t∗+u∗∂C∗∂r∗+v∗r∗∂C∗∂θ∗+w∗∂C∗∂z∗=1ScRe+αc∗×1r∗∂∂r∗r∗∂C∗∂r∗+1r∗2∂2C∗∂θ∗2+∂2C∗∂z∗2−KcReC∗,equation (42) represents the dimensionless form of the Concentration equation.
(43)∂κ∗∂t∗+u∗∂κ∗∂r∗+v∗r∗∂κ∗∂θ∗+w∗∂κ∗∂z∗=1r∗21Resn−1θns−s−11+γT∗ΔTn−1∂k∗∂θ∗+1Reθsn−11+γT∗ΔTn−1+νT∗σk1r∗∂∂r∗r∗∂κ∗∂r∗+1r∗2∂2κ∗∂θ∗2+∂2κ∗∂z∗2+2νT∗∂u∗∂r∗2+1r∗∂v∗∂θ∗+u∗r2+∂w∗∂z∗2+νT∗1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+νT∗∂v∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2−ε∗,equation (43) represents the dimensionless form of the turbulent kinetic energy equation.
(44)∂ε∗∂t∗+u∗∂ε∗∂r∗+v∗r∗∂ε∗∂θ∗+w∗∂ε∗∂z∗=1r∗21Resn−1θns−s−11+γT∗ΔTn−1∂ε∗∂θ∗+1Reθsn−11+γT∗ΔTn−1+νT∗σε1r∗∂∂r∗r∗∂ε∗∂r∗+1r∗2∂2ε∗∂θ∗2+∂2ε∗∂z∗2+2f1Cε,1ε∗k∗νT∗∂u∗∂r∗2+1r∗∂v∗∂θ∗+u∗r∗2+∂w∗∂z∗2+f1Cε,1ε∗k∗νT∗1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+f1νT∗Cε,1ε∗k∗∂uθ∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2−f2Cε,2ε∗2k∗,equation (44) represents the dimensionless form of the dissipation equation.
(45)νT∗=Cμfμκ∗2ε∗,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,Re=u∞D/ν∞; Thermal Grasholf Number, GrT=gβTTw−T∞/ν∞2; Solutal Grasholf Number, Grc=gβcCw−C∞/ν∞2; Schmidt Number, Sc=ν∞/Dm; Chemical Reaction Parameter, Kc=Krν∞/u∞2; Eckert Number, Ec=U∞2/CpTw−T∞ and Prandtl Number, Pr=μ∞Cp/κ∗.


#### 2.4.1. Dimensionless Boundary Conditions

Atr∗=0.5(46)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atr∗=0(47)C∗=0,T∗=0,u∗=1,v∗=1,w∗=1,κ∗=0,∂ε∗∂r∗=0.Atr∗=−0.5(48)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atz=0(49)w∗=1,v∗=0,u∗=0,T=0,C=0.Atz=∞(50)∂u∗∂z∗=0,∂u∗∂z∗=0,∂w∗∂z∗=0,∂κ∗∂z∗=0,∂ε∗∂z∗=0.


## 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]:
(10)−u′u′¯=νT∇u+∇uT−23kδij,δij is the Kronecker delta;
(11)−u′T′¯=αT∂T¯∂r,−v′T′¯=αTr∂T¯∂θ,−w′T′¯=αT∂T¯∂z,−u′C′¯=αc∂C¯∂r,−v′C′¯=αcr∂C¯∂θ,−w′C′¯=αc∂C¯∂z.Substituting the modeling equations (10) and (11) into RANS equations, we have the following final equations:
(12)∂u¯∂t+u¯∂u¯∂r+v¯r∂u¯∂θ+w¯∂u¯∂z−v¯v¯r=−1ρ∂p¯∂r+μ0gn−1ρ+νT1r∂∂rr∂u¯∂r−u¯r2+1r2∂2u¯∂θ2+∂2u¯∂z2−2r2∂v¯∂θ+μ0ρsn−1gn−2θn−11+γT¯−T∞1r2∂u¯∂θ−v¯r2+1r∂v¯∂r+gcosαβTT¯−T∞+βcC¯−C∞,where u¯ is the time average velocity in the radial direction, v¯ is the time average velocity in the θ-direction, w¯ is time-average velocity in the axial direction, p¯ is the time-average pressure and T¯ is the time-average temperature and C¯ is the time-average concentration.
(13)∂u¯∂t+u¯∂u¯∂r+v¯r∂u¯∂θ+w¯∂u¯∂z−v¯v¯r=−1ρ∂p¯∂r+μ0gn−1ρ+νT1r∂∂rr∂u¯∂r−u¯r2+1r2∂2u¯∂θ2+∂2u¯∂z2−2r2∂v¯∂θ+μ0ρsn−1gn−2θn−11+γT¯−T∞1r2∂u¯∂θ−v¯r2+1r∂v¯∂r+gcosαβTT¯−T∞+βcC¯−C∞,(14)∂v¯∂t+u¯∂v¯∂r+v¯r∂v¯∂θ+w¯∂v¯∂z+uvr¯=−1ρr∂p¯∂θ+μ0ρsn−1gn−2θs−11+γT−T∞2r2∂v¯∂θ+2u¯r2μ0gn−1ρ+νT×1r∂∂rr∂v¯∂r−v¯r2+1r2∂2v¯∂θ2+∂2v¯∂z2+2r2∂u¯∂θ,(15)∂v¯∂t+u¯∂v¯∂r+v¯r∂v¯∂θ+w¯∂v¯∂z+uvr¯=−1ρr∂p¯∂θ+μ0ρsn−1gn−2θs−11+γT−T∞2r2∂v¯∂θ+2u¯r2+μ0gn−1ρ+νT1r∂∂rr∂v¯∂r−v¯r2+1r2∂2v¯∂θ2+∂2v¯∂z2+2r2∂u¯∂θ,(16)∂w¯∂t+u¯∂w¯∂r+v¯r∂w¯∂θ+w¯∂w¯∂z=1ρ∂p¯∂z+sn−1μ0gn−2θs−1ρr21+γT¯−T∞∂w¯∂θ+r∂v¯∂z+μ0gn−1ρ+νT1r∂∂r∂w¯∂r+1r2∂w¯∂θ2+∂2w¯∂z2+gsinαβTT¯−T∞+βcC¯−C∞,(17)∂w¯∂t+u¯∂w¯∂r+v¯r∂w¯∂θ+w¯∂w¯∂z=1ρ∂p¯∂z+sn−1μ0gn−2θs−1ρr21+γT¯−T∞×∂w¯∂θ+r∂v¯∂z+μ0gn−1ρ+νT1r∂∂r∂w¯∂r+1r2∂w¯∂θ2+∂2w¯∂z2+gsinαβTT¯−T∞+βcC¯−C∞,(18)∂T¯∂t+u¯∂T¯∂r+v¯r∂T¯∂θ+w¯∂T¯∂z=−κ∗ρCpbΔT∂T¯∂r2+1r2∂T¯∂θ2+∂T¯∂z2+κ∗ρCp1+bT¯−T∞ΔT+αT1r∂T¯∂r+∂2T¯∂r2+1r2∂2T¯∂θ2+∂2T¯∂z2+ϕρCp+εCp,(19)ϕ=2μ∂ur¯∂r2+1r∂uθ¯∂θ+ur¯r2+∂uz¯∂z2+μ1r∂ur¯∂θ+∂uθ¯∂r−uθ¯r2+∂uθ¯∂z+1r∂uz¯∂θ2+∂uz¯∂r+∂ur¯∂z2,where ϕ is the viscous dissipation.
(20)∂C¯∂t+u¯∂C¯∂r+v¯r∂C¯∂θ+w¯∂C¯∂z=Dm+αc1r∂∂rr∂C¯∂r+1r2∂2C¯∂θ2+∂2C¯∂z2−krC¯−C∞.


## 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:
(21)ρ∂κ∂t+u¯·ρ∇κ=∇·μ+μTσk∇κ+μTϕ−ρε.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. σk is the turbulent Prandtl number for kinetic energy.The dissipation equation in vector form as given in [9] is:
(22)ρ∂ε∂t+u¯·ρ∇ε=∇·μ+μTσε∇ε+f1Cε,1εkϕ−f2ρCε,2ε2k,where f1 and f2 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:
(23)∂κ∂t+u¯∂κ∂r+v¯r∂κ∂θ+w¯∂κ∂z=1ρr2μ0sn−1gn−2θs−11+γT¯−T∞∂κ∂θ+μ0gn−1ρ+νTσκ×1r∂∂r∂κ∂r+1r2∂2κ∂θ2+∂2κ∂z2+2νT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+νT1r∂u¯∂θ+∂v¯∂r−v¯r2+νT∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−ε,(24)∂κ∂t+u¯∂κ∂r+v¯r∂κ∂θ+w¯∂k∂z=1r2μ0ρsn−1gn−2θs−11+γT¯−T∞∂κ∂θ+μ0gn−1ρ+νTσκ1r∂∂rr∂κ∂r+1r2∂2κ∂θ2+∂2κ∂z2+2νT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+νT1r∂u¯∂θ+∂v¯∂r−v¯r2+νT∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−ε.Equations (23) and (24) represent TKE for the gaseous and liquid phases respectively.
(25)∂ε∂t+u¯∂ε∂r+v¯r∂ε∂θ+w¯∂ε∂z=1ρr2μ0sn−1gn−2θs−11+γT¯−T∞∂ε∂θ+μ0gn−1ρ+νTσε×1r∂∂rr∂ε∂r+1r2∂2ε∂θ2+∂2ε∂z2+2f1Cε,1εkνT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂u¯∂z2+f1Cε,1εkνT1r∂u¯∂θ+∂v¯∂r−v¯r2+f1νTCε,1εk∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−f2Cε,2ε2k,(26)∂ε∂t+u¯∂ε∂r+v¯r∂ε∂θ+w¯∂ε∂z=1r2μ0ρsn−1gn−2θs−11+γT¯−T∞∂ε∂θ+μ0gn−1ρ+νTσε×1r∂∂rr∂ε∂r+1r2∂2ε∂θ2+∂2ε∂z2+2f1Cε,1εkνT∂u¯∂r2+1r∂v¯∂θ+u¯r2+∂w¯∂z2+f1Cε,1εkνT1r∂u¯∂θ+∂v¯∂r−v¯r2+f1νTCε,1εk∂v¯∂z+1r∂w¯∂θ2+∂w¯∂r+∂u¯∂z2−f2Cε,2ε2k.Equations (25) and (26) represent viscous dissipation equations for the gaseous and liquid phases respectively. Cμ, Cε,1, Cε,2, σε and σk are model constants taken from [21]. Their values are given as: Cμ=0.09, Cε,1=1.44, Cε,2=1.92, σε=1.3, σk=1.The relationship between dissipation(ε), turbulent kinetic energy(k) and turbulent viscosity(νT) is given as:
(27)νT=Cμfμk2ε.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]:
(28)fμ=1−exp−0.0165Rek21+20.5Ret,f1=1+0.05fμ3,f2=1−exp−Ret2,where Rek=ρκ1/2dw/μt=κ1/2dw/νt, Ret=ρk2/μtε=κ2/νtε and dw is the shortest distance to the interface or wall. Ret and Rek 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]:Atr=D+/2;
(29)u=v=w=0,k=0,ε=0,T=Tw,C=Cw.Atr=0 (the interface is treated as a moving wall);
(30)C=C∞,T=T∞,u=u∞,v=u∞,w=u∞,κ=0,∂ε∂r=0.Atr=D−/2;
(31)u=v=w=0,k=0,ε=0,T=Tw,C=Cw.Atz=0;
(32)w=u∞,v=0,u=0,T=T∞,C=C∞.Atz=∞, the gradients of all variables in the flow direction are zero as suggested by [22];
(33)∂v∂z=0,∂u∂z=0,∂w∂z=0,∂κ∂z=0,∂ε∂z=0


## 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]:
(34)u∗=u¯u∞,v∗=v¯u∞,w∗=w¯u∞,t∗=tu∞D,p∗=p¯ρu∞2,θ∗=θ,z∗=zD,νT∗=νTDu∞,r∗=rD,T∗=T¯−T∞Tw−T∞,αT=αTDu∞,ε∗=εDu∞3,k∗=ku∞2,αc∗=αcDu∞,C∗=C¯−C∞Cw−C∞.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:
(35)∂u∗∂t∗+u∗∂u∗∂r∗+v∗r∗∂u∗∂θ∗+w∗∂u∗∂z∗−v∗2r∗=−∂p∗∂r+1Reθsn−11+γT∗ΔT+νT∗×1r∗∂∂r∗r∗∂u∗∂r∗−u∗r∗2+1r∗2∂2u∗∂θ∗2+∂2u∗∂z∗2−2r∗2∂v∗∂θ∗+1Resn−1θsn−1−11+γT∗ΔTn−11r∗2∂u∗∂θ∗−v∗r∗2+1r∗∂v∗∂r∗+cosαT∗GrTRe2+C∗GrCRe2,(36)∂u∗∂t∗+u∗∂u∗∂r∗+v∗r∗∂u∗∂θ∗+w∗∂u∗∂z∗−v∗2r∗=−∂p∗∂r+1Reθsn−11+γT∗ΔTn−1+νT∗×1r∗∂∂r∗r∗∂u∗∂r∗−u∗r∗2+1r∗2∂2u∗∂θ∗2+∂2u∗∂z∗2−2r∗2∂v∗∂θ∗+1Resn−1θsn−1−11+γT∗ΔTn−11r∗2∂u∗∂θ∗−v∗r∗2+1r∗∂v∗∂r∗+cosαT∗GrTRe2+C∗GrCRe2,equations (35) and (36) represent the dimensionless form of the momentum equation in the radial direction for the gaseous and liquid phase respectively.
(37)∂v∗∂t∗+u∗∂v∗∂r∗+v∗r∗∂v∗∂θ∗+w∗∂v∗∂z∗+u∗v∗r∗=−1r∗∂p∗∂θ∗+1Resn−1gn−2θs−11+γT∗Tw−T∞2r∗2∂v∗∂θ∗+2u∗r∗2+1Regn−1+νT∗1r∗∂∂r∗r∗∂v∗∂r∗−v∗r∗2+1r∗2∂2v∗∂θ∗2+∂2v∗∂z∗2+2r∗2∂u∗∂θ∗,(38)∂v∗∂t∗+u∗∂v∗∂r∗+v∗r∗∂v∗∂θ∗+w∗∂v∗∂z∗+u∗uθ∗r∗=−1r∗∂p∗∂θ∗+1Resn−1gn−2θs−11+γT∗Tw−T∞2r∗2∂v∗∂θ∗+2u∗r∗2+1Regn−1+νT∗×1r∗∂∂r∗r∗∂v∗∂r∗−v∗r∗2+1r∗2∂2v∗∂θ∗2+∂2v∗∂z∗2+2r∗2∂u∗∂θ∗,equations (37) and (38) represent the dimensionless form of the momentum equation in the θ-direction for the gaseous and liquid phase respectively.
(39)∂w∗∂t∗+u∗∂w∗∂r∗+v∗r∗∂w∗∂θ∗+w∗∂w∗∂z∗=−∂p∗∂z∗+1Resn−1gn−2θs−1r∗21+γT∗Tw−T∞∂w∗∂θ∗+r∗∂v∗∂z∗+1Regn−1+νT∗1r∗∂∂r∗r∗∂w∗∂r∗+1r∗2∂w∗∂θ∗2+∂2w∗∂z∗2+sinαT∗GrTRe2+C∗GrCRe2,(40)∂w∗∂t∗+u∗∂w∗∂r∗+v∗r∗∂w∗∂θ∗+w∗∂w∗∂z∗=−∂p∗∂z∗+1Re1r∗2sn−1gn−2θs−11+γT∗Tw−T∞∂w∗∂θ∗+r∗∂v∗∂z∗+1Regn−1+νT∗1r∗∂∂r∗r∗∂w∗∂r∗+1r∗2∂w∗∂θ∗2+∂2w∗∂z∗2+sinαT∗GrTRe2+C∗GrCRe2,equations (39) and (40) represent the dimensionless form of the momentum equation in the axial direction for the liquid and gaseous phase respectively.
(41)∂T∗∂t∗+u∗∂T∗∂r∗+v∗r∗∂T∗∂θ∗+w∗∂T∗∂z∗=−bPrRe∂T∗∂r∗2+1r∗2∂T∗∂θ∗2+∂T∗∂z∗2+1+bT∗PrRe+αT∗1r∗∂T∗∂r∗+∂2T∗∂r∗2+1r∗2∂2T∗∂θ∗2+∂2T∗∂z∗2+EcRegn−1×2∂u∗∂r∗2+1r∂v∗∂θ∗+u∗r∗2+∂w∗∂z∗2+1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+EcRegn−1∂v∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2+Ecε∗,equation (41) represents the dimensionless form of the Energy equation.
(42)∂C∗∂t∗+u∗∂C∗∂r∗+v∗r∗∂C∗∂θ∗+w∗∂C∗∂z∗=1ScRe+αc∗×1r∗∂∂r∗r∗∂C∗∂r∗+1r∗2∂2C∗∂θ∗2+∂2C∗∂z∗2−KcReC∗,equation (42) represents the dimensionless form of the Concentration equation.
(43)∂κ∗∂t∗+u∗∂κ∗∂r∗+v∗r∗∂κ∗∂θ∗+w∗∂κ∗∂z∗=1r∗21Resn−1θns−s−11+γT∗ΔTn−1∂k∗∂θ∗+1Reθsn−11+γT∗ΔTn−1+νT∗σk1r∗∂∂r∗r∗∂κ∗∂r∗+1r∗2∂2κ∗∂θ∗2+∂2κ∗∂z∗2+2νT∗∂u∗∂r∗2+1r∗∂v∗∂θ∗+u∗r2+∂w∗∂z∗2+νT∗1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+νT∗∂v∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2−ε∗,equation (43) represents the dimensionless form of the turbulent kinetic energy equation.
(44)∂ε∗∂t∗+u∗∂ε∗∂r∗+v∗r∗∂ε∗∂θ∗+w∗∂ε∗∂z∗=1r∗21Resn−1θns−s−11+γT∗ΔTn−1∂ε∗∂θ∗+1Reθsn−11+γT∗ΔTn−1+νT∗σε1r∗∂∂r∗r∗∂ε∗∂r∗+1r∗2∂2ε∗∂θ∗2+∂2ε∗∂z∗2+2f1Cε,1ε∗k∗νT∗∂u∗∂r∗2+1r∗∂v∗∂θ∗+u∗r∗2+∂w∗∂z∗2+f1Cε,1ε∗k∗νT∗1r∗∂u∗∂θ∗+∂v∗∂r∗−v∗r∗2+f1νT∗Cε,1ε∗k∗∂uθ∗∂z∗+1r∗∂w∗∂θ∗2+∂w∗∂r∗+∂u∗∂z∗2−f2Cε,2ε∗2k∗,equation (44) represents the dimensionless form of the dissipation equation.
(45)νT∗=Cμfμκ∗2ε∗,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,Re=u∞D/ν∞; Thermal Grasholf Number, GrT=gβTTw−T∞/ν∞2; Solutal Grasholf Number, Grc=gβcCw−C∞/ν∞2; Schmidt Number, Sc=ν∞/Dm; Chemical Reaction Parameter, Kc=Krν∞/u∞2; Eckert Number, Ec=U∞2/CpTw−T∞ and Prandtl Number, Pr=μ∞Cp/κ∗.


### 2.4.1. Dimensionless Boundary Conditions

Atr∗=0.5(46)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atr∗=0(47)C∗=0,T∗=0,u∗=1,v∗=1,w∗=1,κ∗=0,∂ε∗∂r∗=0.Atr∗=−0.5(48)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atz=0(49)w∗=1,v∗=0,u∗=0,T=0,C=0.Atz=∞(50)∂u∗∂z∗=0,∂u∗∂z∗=0,∂w∗∂z∗=0,∂κ∗∂z∗=0,∂ε∗∂z∗=0.


## 2.4.1. Dimensionless Boundary Conditions

Atr∗=0.5(46)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atr∗=0(47)C∗=0,T∗=0,u∗=1,v∗=1,w∗=1,κ∗=0,∂ε∗∂r∗=0.Atr∗=−0.5(48)u∗=v∗=w∗=0,k∗=0,ε∗=0,T∗=1,C∗=1.Atz=0(49)w∗=1,v∗=0,u∗=0,T=0,C=0.Atz=∞(50)∂u∗∂z∗=0,∂u∗∂z∗=0,∂w∗∂z∗=0,∂κ∗∂z∗=0,∂ε∗∂z∗=0.


## 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. Letur,θ,z,t be described at the point ri,θj,zk,tm and for simplicity the point ri,θj,zk,tm is written as i,j,k,m and Ur,θ,z,t as ui,j,km. Assume we have I points along r, N points along θ, Q points along z and M points along t directions to form the mesh. Let the step size Δr be along r, Δθ be along θ, Δz along z and Δt along t.ri=R0+iΔr, θj=θ0+jΔθ, zk=z0+kΔz and tm=t0+mΔt where i=1,2,⋯⋯,I, j=1,2,⋯⋯,N, k=1,2,⋯⋯,Q and m=1,2,⋯⋯,MFor 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:(51)∂u∂t=ui,j,km+1−ui,j,kmΔt,(52)∂u∂r=ui+1,j,km−ui−1,j,km2Δr,(53)∂2u∂r2=ui+1,j,km−2ui,j,km+ui−1,j,kmΔr2.Substituting equations (51)–(53) into equations (35)–(44), we have the following set of finite difference equations:
(54)ui,j,km+1−ui,j,kmΔt+ui,j,kmui+1,j,km−ui−1,j,km2Δr+vi,j,kmriui,j+1,km−ui,j−1,km2Δθ+wi,j,kmui,j,k+1m−ui,j,k−1m2Δz−vi,j,km2ri=−∂p∂r+θsn−1Re1+γTi,j,kmΔTn−1+νT×1riui+1,j,km−ui−1,j,km2Δr+ui+1,j,km−2ui,j,km+ui−1,j,kmΔr2−ui,j,kmri2+ui,j+1,km−2ui,j,km+ui,j−1,kmΔθ2+θsn−1Re1+γTi,j,kmΔTn−1+νTui,j,k+1m−2ui,j,km+ui,j,k−1mΔz2−2ri2vi,j+1,km−vi,j−1,km2Δθ+1Resn−1θsn−1−11+γTi,j,kmΔTn−11ri2ui,j+1,km−ui,j−1,km2Δθ−vi,j,kmri2+1rivi+1,j,km−vi−1,j,km2Δr+cosαTi,j,kmGrTRe2+Ci,j,kmGrcRe2,(55)Ui,j,km+1−Ui,j,kmΔt+Ui,j,kmUi+1,j,km−Ui−1,j,km2Δr+Vi,j,kmriUi,j+1,km−Ui,j−1,km2Δθ+Wi,j,kmUi,j,k+1m−Wi,j,k−1m2Δz−Vi,j,km2ri=∂p∂r+θsn−11+γTi,j,kmΔTn−1Re+νT×1riUi+1,j,km−Ui−1,j,km2Δr+Ui+1,j,km−2Ui,j,km+Ui−1,j,kmΔr2−Ui,j,kmri2+Ui,j+1,km−2Ui,j,km+Ui,j−1,kmΔθ2+θsn−11+γTi,j,kmΔTn−1Re+νTUi,j,k+1m−2Ui,j,km+Ui,j,k−1mΔz2−2ri2Vi,j+1,km−Vi,j−1,km2Δθ+sn−1θsn−1−11+γTi,j,kmΔTn−1Re1ri2Ui,j+1,km−Ui,j−1,km2Δθ−Vi,j,kmri2+1riVi+1,j,km−Vi−1,j,km2Δr+cosαTi,j,kmGrTRe2+Ci,j,kmGrcRe2,equations (54) and (55) are radial momentum equations in discrete form for the liquid and gaseous phases respectively.
(56)Vi,j,km+1−Vi,j,kmΔt+Ui,j,kmVi+1,j,km−Vi−1,j,km2Δr+Vi,j,kmriVi,j+1,km−Vi,j−1,km2Δθ+Wi,j,kmVi,j+1,km−Vi,j−1,km2Δz+Ui,j,kmVi,j,kmri=−1ri∂p∂θ+1Resn−1θsn−1−11+γTi,j,kmΔTn−1×2ri2Vi,j+1,km−Vi,j−1,km2Δθ+θsn−1Re1+γTi,j,kmΔTn−1+νT1riVi+1,j,km−Vi−1,j,km2Δr+Vi+1,j,km−2Vi,j,km+Vi−1,j,kmΔr2−Vi,j,kmri2+1ri2Vi,j+1,km−2Vi,j,km+Vi,j−1,kmΔθ2+Vi,j,k+1m−2Vi,j,km+Vi,j,k−1mΔz2+2ri2Ui,j+1,km−Ui,j−1,km2Δθ,(57)Vi,j,km+1−Vi,j,kmΔt+Ui,j,kmVi+1,j,km−Vi−1,j,km2Δr+Vi,j,kmriVi,j+1,km−Vi,j−1,km2Δθ+Wi,j,kmVi,j+1,km−Vi,j−1,km2Δz+Ui,j,kmVi,j,kmri=−1ri∂p∂θ+sn−1θsn−1−11+γTi,j,kmΔTn−1Re×2ri2Vi,j+1,km−Vi,j−1,km2Δθ+θsn−11+γTi,j,kmΔTn−1Re+νT1riVi+1,j,km−Vi−1,j,km2Δr+Vi+1,j,km−2Vi,j,km+Vi−1,j,kmΔr2−Vi,j,kmri2+1ri2Vi,j+1,km−2Vi,j,km+Vi,j−1,kmΔθ2+Vi,j,k+1m−2Vi,j,km+Vi,j,k−1mΔz2+2ri2Ui,j+1,km−Ui,j−1,km2Δθ,equations (56) and (57) represent tangential momentum equation in discrete form for the liquid and gaseous phases respectively.
(58)Wi,j,km+1−Wi,j,kmΔt+Ui,j,kmWi+1,j,km−Wi−1,j,km2Δr+Vi,j,kmriWi,j+1,km−Wi,j−1,km2Δθ+Wi,j,kmWi,j,k+1m−Wi,j,k+1m2Δz=−∂p∂z+1Resn−1θsn−1−11+γTi,j,kmΔTn−1×1ri2Wi,j+1,km−Wi,j−1,km2Δθ+1riVi,j,k+1m−Vi,j,k+1m2Δz+θsn−1Re1+γTi,j,kmΔTn−1+νT1riWi+1,j,km−Wi−1,j,km2Δr+Wi+1,j,km−2Wi,j,km+Wi−1,j,kmΔr2+1ri2Wi,j+1,km−2Wi,j,km+Wi,j−1,kmΔθ2+Wi,j,k+1m−2Wi,j,km+Wi,j,k−1mΔz2+sinαTi,j,kmGrTRe2+Ci,j,kmGrcRe2,(59)Wi,j,km+1−Wi,j,kmΔt+Ui,j,kmWi+1,j,km−Wi−1,j,km2Δr+Vi,j,kmriWi,j+1,km−Wi,j−1,km2Δθ+Wi,j,kmWi,j,k+1m−Wi,j,k+1m2Δz=−∂p∂z+sn−1θsn−1−11+γTi,j,kmΔTn−1Re×1ri2Wi,j+1,km−Wi,j−1,km2Δθ+1riVi,j,k+1m−Vi,j,k+1m2Δz+θsn−11+γTi,j,kmΔTn−1Re+νT1riWi+1,j,km−Wi−1,j,km2Δr+Wi+1,j,km−2Wi,j,km+Wi−1,j,kmΔr2+1ri2Wi,j+1,km−2Wi,j,km+Wi,j−1,kmΔθ2+Wi,j,k+1m−2Wi,j,km+Wi,j,k−1mΔz2+sinαTi,j,kmGrTRe2+Ci,j,kmGrcRe2,equations (58) and (59) represent axial momentum equations in discrete form for liquid and gas phases respectively.
(60)Ti,j,km+1−Ti,j,kmΔt+Ui,j,kmTi+1,j,km−Ti−1,j,km2Δr+Vi,j,kmriTi,j+1,km−Ti,j−1,km2Δθ+Wi,j,kmTi,j+1,km−Ti,j−1,km2Δz=bRePrTi+1,j,km−Ti−1,j,km2Δr2+1ri2Ti,j+1,km−Ti,j−1,km2Δθ2+Ti,j,k+1m−Ti,j,k−1m2Δr2+1+bTi,j,kmRePr+αT×1riTi+1,j,km−Ti−1,j,km2Δr+Ti+1,j,km−2Ti,j,km+Ti−1,j,kmΔr2+1ri2Ti,j+1,km−2Ti,j,km+Ti,j−1,kmΔθ2+Ti,j,k+1m−2Ti,j,km+Ti,j,k−1mΔz2+EcReθsn−11+γTi,j,kmΔT2Ui+1,j,km−Ui−1,j,km2Δr2+1riVi,j+1,km−Vi,j−1,km2Δθ+Ui,j,kmri2+Wi,j,k+1m−Ui,j,k−1m2Δz2+1riUi,j+1,km−Ui,j−1,km2Δθ+Vi+1,j,km−Vi−1,j,km2Δr−Vi,j1,kmri2+Vi,j,k+1m−Vi,j,k−1m2Δz+1riWi,j+1,km−Wi,j−1,km2Δθ2+Wi+1,j,km−Wi−1,j,km2Δr+Ui,j,k+1m−Ui,j,k−1m2Δz2+Ecε,equation (60) represents energy equation in discrete form.
(61)Ci,j,km+1−Ci,j,kmΔt+Ui,j,kmCi+1,j,km−Ci−1,j,km2Δr+Vi,j,kmriCi,j+1,km−Ci,j−1,km2Δθ+Wi,j,kmCi,j,k+1m−Ci,j,k−1m2Δz=1ScRe+αc1riCi+1,j,km−Ci−1,j,km2Δr+Ci+1,j,km−2Ci,j,km+Ci−1,j,kmΔr2+1ri2Ci,j+1,km−2Ci,j,km+Ci,j−1,kmΔθ2+Ci,j,k+1m−2Ci,j,km+Ci,j,k−1mΔz2−ΔtKcReCi,j,km,equation (61) represents Concentration equation in discrete form.
(62)κi,j,km+1−κi,j,kmΔt+Ui,j,kmκi+1,j,km−κi−1,j,km2Δr+Vi,j,kmriκi,j+1,km−κi,j−1,km2Δθ+Wi,j,kmκi,j,k+1m−κi,j,k−1m2Δz=1ri2sn−1θsn−1−1Re1+γTi,j,kmΔn−1κi,j+1,km−κi,j−1,km2Δθ+θsn−1Re1+γTi,j,kmΔn−1+νTσkκi+1,j,km−2κi,j,km+κi−1,j,kmΔr2+1riκi+1,j,km−κi−1,j,km2Δr+1ri2κi,j+1,km−2κi,j,km+κi,j−1,kmΔθ2+κi,j,k+1m−2κi,j,km+κi,j,k−1mΔz2+2νTUi+1,j,km−Ui−1,j,km2Δr2+1riVi,j+1,km−Vi,j−1,km2Δθ+Ui,j,kmri2+Wi,j,k+1m−Ui,j,k−1m2Δz2+νT1riUi,j+1,km−Ui,j−1,km2Δθ+Vi+1,j,km−Vi−1,j,km2Δr−Vi,j,kmri2+Vi,j,k+1m−Vi,j,k−1m2Δz+1riWi,j+1,km−Wi,j−1,km2Δθ2+Wi+1,j,km−Wi−1,j,km2Δr+Ui,j,k+1m−Ui,j,k−1m2Δz2−ε,(63)κi,j,km+1−κi,j,kmΔt+Ui,j,kmκi+1,j,km−κi−1,j,km2Δr+Vi,j,kmriκi,j+1,km−κi,j−1,km2Δθ+Wi,j,kmκi,j,k+1m−κi,j,k−1m2Δz=sn−1θsn−1−11+γTi,j,kmΔn−1Reri2κi,j+1,km−κi,j−1,km2Δθ+θsn−11+γTi,j,kmΔn−1Re+νTσkκi+1,j,km−2κi,j,km+κi−1,j,kmΔr2+1riκi+1,j,km−κi−1,j,km2Δr+1ri2κi,j+1,km−2κi,j,km+κi,j−1,kmΔθ2+κi,j,k+1m−2κi,j,km+κi,j,k−1mΔz2+2νTUi+1,j,km−Ui−1,j,km2Δr2+1riVi,j+1,km−Vi,j−1,km2Δθ+Ui,j,kmri2+Wi,j,k+1m−Ui,j,k−1m2Δz2+νT1riUi,j+1,km−Ui,j−1,km2Δθ+Vi+1,j,km−Vi−1,j,km2Δr−Vi,j,kmri2+Vi,j,k+1m−Vi,j,k−1m2Δz+1riWi,j+1,km−Wi,j−1,km2Δθ2+Wi+1,j,km−Wi−1,j,km2Δr+Ui,j,k+1m−Ui,j,k−1m2Δz2−ε,equations (63) and (62) are Turbulent Kinetic Energy equations in discrete form for liquid and gaseous phase respectively.
(64)εi,j,km+1−εi,j,kmΔt+Ui,j,kmεi+1,j,km−εi−1,j,km2Δr+Vi,j,kmriεi,j+1,km−εi,j−1,km2Δθ+Wi,j,kmεi,j,k+1m−εi,j,k−1m2Δz=1ri2sn−1θsn−1−1Re1+γTi,j,kmΔn−1εi,j+1,km−εi,j−1,km2Δθ+θsn−1Re1+γTi,j,kmΔn−1+νTσk1riεi+1,j,km−εi−1,j,km2Δr+εi+1,j,km−2εi,j,km+εi−1,j,kmΔr2+1ri2εi,j+1,km−2εi,j,km+εi,j−1,kmΔθ2+εi,j,k+1m−2εi,j,km+εi,j,k−1mΔz2+2f1Cε,1εi,j,kmκi,j,kmνTUi+1,j,km−Ui−1,j,km2Δr2+1riVi,j+1,km−Vi,j−1,km2Δθ+Ui,j,kmri2+Wi,j,k+1m−Ui,j,k−1m2Δz2+f1Cε,1εi,j,kmκi,j,kmνT1riUi,j+1,km−Ui,j−1,km2Δθ+Vi+1,j,km−Vi−1,j,km2Δr−Vi,j,kmri2+Vi,j,k+1m−Vi,j,k−1m2Δz+1riWi,j+1,km−Wi,j−1,km2Δθ2+Wi+1,j,km−Wi−1,j,km2Δr+Ui,j,k+1m−Ui,j,k−1m2Δz2−f2Cε,2εi,j,km2κi,j,km,(65)εi,j,km+1−εi,j,kmΔt+Ui,j,kmεi+1,j,km−εi−1,j,km2Δr+Vi,j,kmriεi,j+1,km−εi,j−1,km2Δθ+Wi,j,kmεi,j,k+1m−εi,j,k−1m2Δz=1ri2sn−1θsn−1−11+γTi,j,kmΔn−1Reεi,j+1,km−εi,j−1,km2Δθ+θsn−11+γTi,j,kmΔn−1Re+νTσk1riεi+1,j,km−εi−1,j,km2Δr+εi+1,j,km−2εi,j,km+εi−1,j,kmΔr2+1ri2εi,j+1,km−2εi,j,km+εi,j−1,kmΔθ2+εi,j,k+1m−2εi,j,km+εi,j,k−1mΔz2+2f1Cε,1εi,j,kmκi,j,kmνTUi+1,j,km−Ui−1,j,km2Δr2+1riVi,j+1,km−Vi,j−1,km2Δθ+Ui,j,kmri2+Wi,j,k+1m−Ui,j,k−1m2Δz2+f1Cε,1εi,j,kmκi,j,kmνT1riUi,j+1,km−Ui,j−1,km2Δθ+Vi+1,j,km−Vi−1,j,km2Δr−Vi,j,kmri2+Vi,j,k+1m−Vi,j,k−1m2Δz+1riWi,j+1,km−Wi,j−1,km2Δθ2+Wi+1,j,km−Wi−1,j,km2Δr+Ui,j,k+1m−Ui,j,k−1m2Δz2−f2Cε,2εi,j,km2κi,j,km.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 m+1 are on the left hand side (since they depend explicitly on the time at m), 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:Pr=0.71(Air-gas phase) and Pr=7.0(water-liquid phase).Figure2 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.Figure 2
Effect of Reynolds Number on Velocity.
(a)
Gaseous phase(b)
Liquid PhaseFigure3 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.Figure 3
Effect of Grashof Solutal Number on Velocity.
(a)
Gaseous phase(b)
Liquid PhaseEffect of Thermal Grasholf number on velocity profile is analyzed in Figure4 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.Figure 4
Effects of Grashof Thermal Number on Velocity.
(a)
Gaseous phase(b)
Liquid PhaseFigure5 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 cosα decreases. Thus, the effect of buoyancy forces decrease along the radial direction. Consequently, there are decreased velocity profiles along the radial direction.Figure 5
Effects of Angle of Inclinationα on Velocity.
(a)
Gaseous phase(b)
Liquid PhaseEffect of Reynolds number on the temperature profiles for both gas and liquid phases are displayed in Figure6. 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 Re 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.Figure 6
EffectsRe on Temperature.
(a)
Gaseous phase(b)
Liquid PhaseFigure7 displays the effect of Eckert number on the temperature profiles. Here, temperature profiles increase for both phases as Eckert number increases. As Ec 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.Figure 7
Effect Eckert Number on Temperature.
(a)
Gaseous phase(b)
Liquid PhaseFrom Figure8, 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.Figure 8
Effect of Angle of Inclination on Temperature Profile.
(a)
Gaseous phase(b)
Liquid PhaseFigure9 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.Figure 9
Effect of Thermal Grasholf Number on Temperature Profile.
(a)
Gaseous phase(b)
Liquid PhaseFigure10 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 Sc 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.Figure 10
Effect of Schmidt Number on Concentration.
(a)
Gaseous phase(b)
Liquid PhaseEffect of Chemical Reaction Parameter on concentration is analyzed in Figure11 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.Figure 11
Effect of Chemical Reaction Parameter on Concentration.
(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.Figure 12
Comparison between Duan [9] velocity profile for gaseous phase and velocity profile for this study.
(a)
Velocity profile for Gaseous phase (Present work)(b)
Duan [9] velocity profile


## 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.Figure 12
Comparison between Duan [9] velocity profile for gaseous phase and velocity profile for this study.
(a)
Velocity profile for Gaseous phase (Present work)(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.



---

*Source: 7617017-2022-11-01.xml*