Heat and Mass Transfers in a Two-Phase Stratified Turbulent Fluid Flow in a Geothermal Pipe with Chemical Reaction

Eric M. Nyariki¹, Mathew N. Kinyanjui², Jeconia O. Abonyo²

¹Department of Mathematics, Pan African Institute of Basic Science, Technology and Innovation, Nairobi, Kenya.
²Department of Pure and Applied Mathematics, Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya.

DOI: 10.4236/jamp.2023.112030                

Abstract

This research focused on the study of heat and mass transfers in a two-phase stratified turbulent fluid flow in a geothermal pipe with chemical reaction. The derived non-linear partial differential equations governing the flow were solved using the Finite Difference Method. The effects of various physical parameters on the concentration, skin friction, heat, and mass transfers have been determined. Analysis of the results obtained indicated that the coefficient of skin friction decreased with an increase in Reynolds number and solutal Grasholf number, the rate of heat transfer increased with an increase in Eckert number, Prandtl number, and angle of inclination, and the rate of mass transfer increased with increase in Reynolds number, Chemical reaction parameter and angle of inclination. The findings would be useful to engineers in designing and maintaining geothermal pipelines more effectively.

Keywords

Two-Phase, Turbulence, Non-Newtonian, Inclination, Heat Transfer, Mass Transfer, Stratified

1. Introduction

In two-phase systems like power plants, nuclear reactors, oil and gas pipelines, and refrigeration equipment, a stratified flow regime is regarded as one of the basic and simple flows. When brine and steam are present in the fluid flow, two-phase flows take place in geothermal pipelines. When brine is positioned at the bottom and steam is positioned at the top of a geothermal pipe, the fluid flow is said to be stratified.

In vertical and horizontal pipelines, many flow regimes, such as plug, stratified, and annular, can occur. [1] studied two-phase flow in geothermal wells using several models to establish flow regimes. Results of models of pressure drop were compared with experimental values. They discovered that stratified wavy flow depended on the flow velocity and the void fraction, which was the most prevalent flow regime in a horizontal pipe. Furthermore, [2] investigated the flow regime transitions and flow structure of downward two-phase flow in large-diameter pipes. The following three flow regimes—cap-bubbly, churn turbulent, and annular flow—were observed in a downward flow, whereas pseudo-slug, plug, and stratified flow regimes were seen in a horizontal section.

Mass transfer is the movement of molecules within a mixture from a high-concentration area to a low-concentration area. The study of thermal energy transfer between physical systems when there is a temperature difference is known as heat transfer. Several researchers have reported on heat and mass fluxes in various available literature. When there is a temperature difference between bodies, heat transfer occurs between the different parts of the same body. Mass transfer results due to species concentration gradient in the geothermal fluid. In order to predict temperature and pressure profiles, two-phase flow and heat transfer models were used by [3] to study two-phase fluid flow and heat transfer in vertical geothermal, oil, and gas wells. The governing equations of continuity, momentum, and energy were solved while taking into account a bubbly flow regime. According to the results, an increase in mass flow rate caused a rise in temperature. [4] employed a drift flux approach method to examine two-phase fluid and heat flow in geothermal wells. Modeling heat transfer from the well-bore fluid to the surroundings involved treating the well-bore as a heat sink. One-dimensional flow and an annular flow regime with no slip approach were taken into account. Analysis was done on heat transfer resistances caused by different well-bore components. The findings showed that greater heat transfer caused a decrease in enthalpy. Despite the analysis demonstrating that various models functioned similarly, the newly presented model was superior.

Modeling of a single and two-phase fluid flow with silica scaling deposition in geothermal wells was done by [5] . Temperature was a key factor in the formulation of the model that they suggested to predict fluid movement, heat transfer, and silica deposition. The solubility of silica in the solution decreased, according to the results, as the temperature decreased. [6] examined an angular two-phase flow with heat and mass transfers in a geothermal pipe. An annular flow regime, non-linear viscosity, and a two-dimensional incompressible laminar flow were taken into account. The nonlinear governing equations were solved using the BVP4c collocation method. The findings showed that a rise in Reynolds number had an impact on the coefficient of skin friction, mass, and heat transfer.

[7] examined a two-phase flow in order to comprehend heat transmission and phase transition processes in a scaled geothermal well-bore. For the vapor and liquid phases, the mass, momentum, and energy equations were solved using a six-equation model that is accessible in ANSYS Fluent. It was observed that the scaling zone had a high mass transfer rate. The impact of Reynolds numbers on the hydrodynamic and thermal parameters of a turbulent flow with heat transfer in an inclined circular channel was examined by [8] . The two-phase mixture model, the finite volume method, and the second-order difference scheme were used to solve the governing equations. According to the suggestion of [9] , the κ-ε turbulence model was employed to simulate the turbulence. The results demonstrated that as the Reynolds number rose, the coefficient of skin friction dropped but the rate of heat transfer coefficient increased. In a turbulent flow via a tube, heat transmission was examined by [10] . The turbulent fluid flow model developed by Spalart Allmaras was used to examine the movement of fluid inside the tube. The results showed that the rate of heat transfer increased when both the Prandtl number and the Reynolds number increased. [11] examined the effects of chemical reaction and a uniform heat source on the Soret-Dufour and thermophoresis on MHD mixed convective mass and heat transfers of a semi-infinite plate. With an increase in the chemical reaction parameter, the rate of heat transfer and the coefficient of skin friction decreased. [12] studied a continuous three-dimensional flow with heat transfer of nanofluidic due to the stretched sheet in the presence of a heat source and a magnetic field. The Runge Kutta-Fehlberg scheme was used to numerically solve the governing equations. The findings showed that as the thermal conductivity value increased, temperature and velocity rose.

Different models can be used to simulate the turbulent fluid flow that occurs in geothermal pipelines. A two-phase stratified smooth turbulent flow through a circular pipe was studied by [13] . A low Reynolds turbulence model was used to solve a stable 2-dimensional momentum equation. Governing equations for non-linear systems were solved using the Newton-Raphson method. The findings of the research revealed that the profile of both κ and ε near the interface was higher in the gas phase than it was in the liquid phase. It was also reported that the interface acted as a moving wall and reduced fluid flow resistance between it and the pipe wall.

[14] investigated the void fraction and two-phase unstable turbulent flow patterns in both vertical and horizontal pipes. The Volume of Fluid (VOF) model and the ANSYS FLUENT program were used to examine unstable turbulent fluid flow. The κ-ε (realizable) model was utilized to solve the slug and bubble flow regimes while the RNG κ-ε model was used to solve turbulent fluid flow with annular and churn flow regimes. They deduced from the findings that the air superficial velocity affected how the flow patterns transitioned. [15] used computational fluid dynamics to study steady 3-dimensional turbulent fluid flow in a conduit. Two models were compared, κ-ω and κ-ε. The results showed that κ-ε provided better center-line velocities approximations and was adequate for estimating turbulent flow in a pipe.

The aforementioned investigations on geothermal pipes have usually ignored inclined geothermal pipes, variable viscosity, variable thermal conductivity, and chemical reaction effects. In geothermal plants, heat and mass transport result from variations in temperature and concentration.

The rationale behind the current work is covered in the novel aspects that follow. Firstly, formulation of a two-phase, stratified, viscous, non-Newtonian, temperature-dependent, incompressible fluid flow with heat and mass transfers in a geothermal pipe. Due to the topography, the geothermal pipe is inclined at a slight inclination α. Variable thermal conductivity is also taken into account. The second step is to use MATLAB to implement the results of solving non-linear partial differential equations using the finite difference method. Lastly, ascertaining how skin friction, heat, and mass transfers are impacted by flow parameters.

2. Mathematical Model

Let us consider a two-phase stratified turbulent fluid flow with heat and mass transfers in a geothermal pipe inclined at a slight angle. Buoyancy forces are considered because of the pipe’s slope. The gaseous phase is at the top of the pipe, and the liquid phase is at the bottom as shown in Figure 1. The flow is assumed to be unsteady, incompressible, and turbulent. The momentum equation takes into account a non-linear viscosity that relies on tangential direction and temperature in both the liquid and gaseous phases. Thermal conductivity, which is a function of temperature, is taken into account in the energy equation. The model’s cylindrical coordinates are ( $r\mathrm{,}\theta \mathrm{,}z$ ), where z is along the pipe axis.

Using the aforementioned assumptions as a guide, the precise governing equations in cylindrical coordinates as presented by [16] are as follows:

The equation of continuity:

$\frac{\partial \left(w\right)}{\partial z}+\frac{1}{r}\frac{\partial \left(v\right)}{\partial \theta }+\frac{1}{r}\frac{\partial \left(ru\right)}{\partial r}=0,$ (1)

Figure 1. Geometry of the flow problem.

where w, v and u are velocities in the z, θ and r directions respectively.

Momentum equation in the radial direction:

$\begin{array}{l}\rho \left(\frac{\partial u}{\partial t}+w\frac{\partial u}{\partial z}+\frac{v}{r}\frac{\partial u}{\partial \theta }+u\frac{\partial u}{\partial r}-\frac{v^2}{r}\right)\\ =\mu \left[-\frac{u}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)\right]-\frac{\partial p}{\partial r}+\mu \left[\frac{{\partial }^{2}u}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}-\frac{2}{r^2}\frac{\partial v}{\partial \theta }\right]\\ +\frac{\partial \mu }{\partial \theta }\left(\frac{1}{r^2}\frac{\partial u}{\partial \theta }+\frac{1}{r}\frac{\partial v}{\partial r}-\frac{v}{r^2}\right)+\rho g\cos \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(T-T_{\infty }\right)+{\beta }_c^{\mathrm{<mo>\; *\; </mo>}}\left(C-C_{\infty }\right)\right]\mathrm{,}\end{array}$ (2)

where μ denotes the dynamic viscosity, ρ denotes the fluid’s temperature, α denotes the angle of inclination, and $g$ is the acceleration brought on by gravity. The coefficients of thermal and mass expansion are denoted by ${\beta }_T^{\mathrm{<mo>\; *\; </mo>}}$ and ${\beta }_C^{\mathrm{<mo>\; *\; </mo>}}$ , respectively.

Equation of momentum in the θ-direction:

$\begin{array}{l}\rho \left[\frac{\partial v}{\partial t}+w\frac{\partial v}{\partial z}+\frac{v}{r}\frac{\partial v}{\partial \theta }+\frac{uv}{r}+u\frac{\partial v}{\partial r}\right]\\ =\frac{\partial \mu }{\partial \theta }\left[2\frac{u}{r^2}+\frac{2}{r^2}\frac{\partial v}{\partial \theta }\right]+\mu \left[\frac{1}{r^2}\frac{{\partial }^{2}v}{\partial {\theta }^2}-\frac{v}{r^2}+\frac{{\partial }^{2}v}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial v}{\partial r}\right)+\frac{2}{r^2}\frac{\partial u}{\partial \theta }\right]\mathrm{,}\end{array}$ (3)

Equation of momentum equation in the z-direction:

$\begin{array}{l}\rho \left[\frac{\partial w}{\partial t}+w\frac{\partial w}{\partial z}+\frac{v}{r}\frac{\partial w}{\partial \theta }+u\frac{\partial w}{\partial r}\right]\\ =\frac{\partial \mu }{\partial \theta }\left[\frac{1}{r}\frac{\partial v}{\partial z}+\frac{1}{r^2}\frac{\partial w}{\partial \theta }\right]+{\mu }^{\star }\left[\frac{{\partial }^{2}w}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial w}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}w}{\partial {\theta }^2}\right]\\ +\rho g\sin \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(T-T_{\infty }\right)+{\beta }_c^{\mathrm{<mo>\; *\; </mo>}}\left(C-C_{\infty }\right)\right]\mathrm{,}\end{array}$ (4)

Equation for energy:

$\begin{array}{l}\rho C_p\left[\frac{v}{r}\frac{\partial T}{\partial \theta }+w\frac{\partial T}{\partial z}+u\frac{\partial T}{\partial r}+\frac{\partial T}{\partial t}\right]\\ =\frac{\partial {\kappa }^{\star }}{\partial T}{\left(\frac{\partial T}{\partial r}\right)}^2+\frac{1}{r^2}\frac{\partial {\kappa }^{\star }}{\partial T}{\left(\frac{\partial T}{\partial \theta }\right)}^2+\frac{\partial {\kappa }^{\star }}{\partial T}{\left(\frac{\partial T}{\partial z}\right)}^2\\ +{\kappa }^{\star }\left(\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)+\Phi \mathrm{,}\end{array}$ (5)

where ${\kappa }^{\star }$ is the fluid’s thermal conductivity and $C_p$ is the specific heat at constant pressure.

The equation for concentration:

$\begin{array}{l}u\frac{\partial C}{\partial r}+w\frac{\partial C}{\partial z}+\frac{v}{r}\frac{\partial C}{\partial \theta }+\frac{\partial C}{\partial t}\\ =D_c\left[\frac{1}{r^2}\frac{{\partial }^{2}C}{\partial {\theta }^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial C}{\partial r}\right)+\frac{{\partial }^{2}C}{\partial z^2}\right]-k_r\left(C-C_{\infty }\right)\mathrm{,}\end{array}$ (6)

where $D_c$ denotes the diffusion coefficient while $k_r$ denotes the chemical reaction parameter.

Since the fluid is non-Newtonian, the power law model, as used by [17] , has been utilized to characterize the viscosity;

${\mu }^{\star }={\mu }_0{f}^{a-1},f=f\left(\theta \right),$ (7)

where ${\mu }_0$ represents the flow consistency index and a represents the power law index.

Due to the two-phase nature of the fluid flow under consideration in this work, a modified piece-wise function has been employed according to [18] to define the flow;

$f\left(\theta \mathrm{,}T\right)=\{\begin{array}{l}{\theta }^{\eta }\left[1+\lambda \left(T-T_{\infty }\right)\right],\text{for}\lambda <0,\eta \ne 0;\\ \frac{{\theta }^{\eta }}{1+\lambda \left(T-T_{\infty }\right)},\text{for}\lambda >0,\eta \ne 0.\end{array}$ (8)

η and λ are constants, f represents the velocity gradient. $\lambda <0$ and $\lambda >0$ represent the viscosity of gas and liquid phases respectively.

According to [19] , thermal conductivity ${\kappa }^{\star }$ is temperature-dependent and is described as

${\kappa }^{\star }={\kappa }_{\infty }\left[1+\frac{\omega \left(T-T_{\infty }\right)}{\Delta T}\right]\mathrm{,}$ (9)

where ω is a tiny parameter that can take positive or negative values depending on whether it’s a liquid or a gas. The fluid’s thermal conductivity is ${\kappa }_{\infty }$ .

2.1. Reynolds Decomposition and Turbulence Modeling

Equations (7)-(9) are substituted into (1), (2), (3), (4), (5) and (6) in order to simplify them. For turbulent flows, the simplified governing equations obtained are decomposed into time average and fluctuating parts (velocity, temperature, pressure, and concentration parameters which are important for the analysis in engineering). To derive the Reynolds Averaged Navier Stokes equations (RANS), the momentum equations are replaced by the mean and fluctuating velocities, then each equation is multiplied by the fluctuating velocities in their respective directions. The decomposed equations are time-averaged again and then Reynolds decomposition rules are applied.

According to [20] , the Reynolds stresses are stated in the mean flow quantities in order to solve the RANS equations. The modeling statements utilized are as follows:

$-\overline{u^{\prime }{u}^{\prime }}={\nu }_T\left(\left(\nabla u\right)+{\left(\nabla u\right)}^{\text{T}}\right)-\frac{2}{3}k{\delta }_{ij}\mathrm{,}$ (10)

${\delta }_{ij}$ is the Kronecker delta

$\begin{array}{l}-\overline{u^{\prime }{T}^{\prime }}={\phi }_T\frac{\partial \bar{T}}{\partial r},-\overline{v^{\prime }{T}^{\prime }}=\frac{{\phi }_T}{r}\frac{\partial \bar{T}}{\partial \theta },-\overline{w^{\prime }{T}^{\prime }}={\phi }_T\frac{\partial \bar{T}}{\partial z},\\ -\overline{u^{\prime }{C}^{\prime }}={\phi }_c\frac{\partial \bar{C}}{\partial r},-\overline{v^{\prime }{C}^{\prime }}=\frac{{\phi }_c}{r}\frac{\partial \bar{C}}{\partial \theta },-\overline{w^{\prime }{C}^{\prime }}={\phi }_c\frac{\partial \bar{C}}{\partial z}.\end{array}$ (11)

The Reynolds Navier Stokes equations are subjected to modeling statements (10) and (11), and only the final form is provided;

Averaged momentum equation for the gas phase in the radial direction:

$\begin{array}{l}\bar{u}\frac{\partial \bar{u}}{\partial r}+\bar{w}\frac{\partial \bar{u}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{u}}{\partial \theta }+\frac{\partial \bar{u}}{\partial t}-\frac{\bar{v}\bar{v}}{r}\\ =\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\left[-\frac{2}{r^2}\frac{\partial \bar{v}}{\partial \theta }+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \bar{u}}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\bar{u}}{\partial {\theta }^2}-\frac{\bar{u}}{r^2}+\frac{{\partial }^2\bar{u}}{\partial z^2}\right]\\ -\frac{1}{\rho }\frac{\partial p}{\partial r}+\frac{{\mu }_0}{\rho }\eta \left(a-1\right)f^{a-2}{\theta }^{a-1}\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]\left[\frac{1}{r^2}\frac{\partial \bar{u}}{\partial \theta }-\frac{\bar{v}}{r^2}+\frac{1}{r}\frac{\partial \bar{v}}{\partial r}\right]\\ +g\cos \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{T}-T_{\infty }\right)+{\beta }_c^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{C}-C_{\infty }\right)\right]\mathrm{,}\end{array}$ (12)

where $\bar{u}\mathrm{,}\bar{v}\mathrm{,}\bar{w}$ are the time-averaged velocities in the radial, tangential and axial directions respectively, $\bar{T}$ is the time average temperature, $\bar{C}$ is the time average concentration, $\bar{p}$ is the time average pressure and ${\nu }_T$ is the eddy viscosity.

Averaged momentum equation for the liquid phase in the radial direction:

$\begin{array}{l}\bar{w}\frac{\partial \bar{u}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{u}}{\partial \theta }+\bar{u}\frac{\partial \bar{u}}{\partial r}-\frac{\bar{v}\bar{v}}{r}+\frac{\partial \bar{u}}{\partial t}\\ =\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\left[-\frac{\bar{u}}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \bar{u}}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\bar{u}}{\partial {\theta }^2}-\frac{2}{r^2}\frac{\partial \bar{v}}{\partial \theta }+\frac{{\partial }^2\bar{u}}{\partial z^2}\right]\\ -\frac{1}{\rho }\frac{\partial p}{\partial r}+\frac{{\mu }_0}{\rho }\frac{\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}}{1+\lambda \left(\bar{T}-T_{\infty }\right)}\left[\frac{1}{r^2}\frac{\partial \bar{u}}{\partial \theta }+\frac{1}{r}\frac{\partial \bar{v}}{\partial r}-\frac{\bar{v}}{r^2}\right]\\ +g\cos \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{T}-T_{\infty }\right)+{\beta }_C^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{C}-C_{\infty }\right)\right]\mathrm{,}\end{array}$ (13)

The averaged momentum equation for the liquid phase in the tangential direction:

$\begin{array}{l}\frac{\bar{u}\bar{v}}{r}+\bar{u}\frac{\partial \bar{v}}{\partial r}+\bar{w}\frac{\partial \bar{v}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{v}}{\partial \theta }+\frac{\partial \bar{v}}{\partial t}\\ =\frac{{\mu }_0}{\rho }\frac{\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}}{1+\lambda \left(T-T_{\infty }\right)}\left[2\frac{\bar{u}}{r^2}+\frac{2}{r^2}\frac{\partial \bar{v}}{\partial \theta }\right]\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\\ \times \left[-\frac{\bar{v}}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \bar{v}}{\partial r}\right)+\frac{2}{r^2}\frac{\partial \bar{u}}{\partial \theta }+\frac{{\partial }^2\bar{v}}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^2\bar{v}}{\partial {\theta }^2}\right]-\frac{1}{\rho r}\frac{\partial p}{\partial \theta }\mathrm{,}\end{array}$ (14)

The averaged momentum equation for the gas phase in tangential direction:

$\begin{array}{l}\frac{\bar{u}\bar{v}}{r}+\bar{u}\frac{\partial \bar{v}}{\partial r}+\bar{w}\frac{\partial \bar{v}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{v}}{\partial \theta }+\frac{\partial \bar{v}}{\partial t}\\ =\frac{{\mu }_0}{\rho }\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}\left[1+\lambda \left(T-T_{\infty }\right)\right]\left[\frac{2}{r^2}\frac{\partial \bar{v}}{\partial \theta }+2\frac{\bar{u}}{r^2}\right]-\frac{1}{\rho r}\frac{\partial p}{\partial \theta }\\ +\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\left[\frac{2}{r^2}\frac{\partial \bar{u}}{\partial \theta }-\frac{\bar{v}}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \bar{v}}{\partial r}\right)+\frac{{\partial }^2\bar{v}}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^2\bar{v}}{\partial {\theta }^2}\right]\mathrm{,}\end{array}$ (15)

The averaged momentum equation for the liquid phase in the axial direction:

$\begin{array}{l}\bar{w}\frac{\partial \bar{w}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{w}}{\partial \theta }+\bar{u}\frac{\partial \bar{w}}{\partial r}+\frac{\partial \bar{w}}{\partial t}\\ =\frac{\eta \left(a-1\right){\mu }_0{f}^{a-2}{\theta }^{\eta -1}}{\rho r^2\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]}\left(r\frac{\partial \bar{v}}{\partial z}+\frac{\partial \bar{w}}{\partial \theta }\right)-\frac{1}{\rho }\frac{\partial p}{\partial z}\end{array}$

$\begin{array}{l}+\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\left[\frac{1}{r^2}\frac{\partial \bar{w}}{\partial {\theta }^2}+\frac{1}{r}\frac{\partial }{\partial }\left(r\frac{\partial \bar{w}}{\partial r}\right)+\frac{{\partial }^2\bar{w}}{\partial z^2}\right]\\ +g\sin \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{T}-T_{\infty }\right)+{\beta }_c^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{C}-C_{\infty }\right)\right]\mathrm{,}\end{array}$ (16)

The averaged momentum equation for the gas phase in the axial direction:

$\begin{array}{l}\frac{\partial \bar{w}}{\partial t}+\bar{w}\frac{\partial \bar{w}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{w}}{\partial \theta }+\bar{u}\frac{\partial \bar{w}}{\partial r}\\ =-\frac{1}{\rho }\frac{\partial p}{\partial z}+\frac{\eta \left(a-1\right){\mu }_0{f}^{a-2}{\theta }^{\eta -1}}{\rho r^2}\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]\left(\frac{\partial \bar{w}}{\partial \theta }+r\frac{\partial \bar{v}}{\partial z}\right)\\ +\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+{\nu }_T\right)\left[\frac{1}{r}\frac{\partial }{\partial }\left(r\frac{\partial \bar{w}}{\partial r}\right)+\frac{{\partial }^2\bar{w}}{\partial z^2}+\frac{1}{r^2}\frac{\partial \bar{w}}{\partial {\theta }^2}\right]\\ +g\sin \alpha \left[{\beta }_T^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{T}-T_{\infty }\right)+{\beta }_c^{\mathrm{<mo>\; *\; </mo>}}\left(\bar{C}-C_{\infty }\right)\right]\end{array}$ (17)

Averaged energy equation:

$\begin{array}{l}\frac{\partial \bar{T}}{\partial t}+\bar{w}\frac{\partial \bar{T}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{T}}{\partial \theta }+\bar{u}\frac{\partial \bar{T}}{\partial r}\\ =-\frac{{\kappa }^{\star }}{\rho C_p}\frac{\omega }{\Delta T}\left[{\left(\frac{1}{r^2}{\left(\frac{\partial \bar{T}}{\partial \theta }\right)}^2+\frac{\partial \bar{T}}{\partial r}\right)}^2+{\left(\frac{\partial \bar{T}}{\partial z}\right)}^2\right]\\ +\left(\frac{{\kappa }^{\star }}{\rho C_p}\left[1+\omega \frac{\bar{T}-T_{\infty }}{\Delta T}\right]+{\phi }_T\right)\left(\frac{1}{r}\frac{\partial \bar{T}}{\partial r}+\frac{{\partial }^2\bar{T}}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2\bar{T}}{\partial {\theta }^2}+\frac{{\partial }^2\bar{T}}{\partial z^2}\right)\\ +\frac{\Phi }{\rho C_p}+\frac{\epsilon }{C_p}\mathrm{,}\end{array}$ (18)

$\begin{array}{c}\Phi =2{\mu }^{\star }\left[{\left(\frac{\partial \bar{u}}{\partial r}\right)}^2+{\left(\frac{\partial \bar{w}}{\partial z}\right)}^2+{\left(\frac{\bar{u}}{r}+\frac{1}{r}\frac{\partial \bar{v}}{\partial \theta }\right)}^2\right]\\ +{\mu }^{\star }\left[{\left(\frac{1}{r}\frac{\partial \bar{u}}{\partial \theta }-\frac{\bar{v}}{r}+\frac{\partial \bar{v}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{w}}{\partial \theta }+\frac{\partial \bar{v}}{\partial z}\right)}^2+{\left(\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial r}\right)}^2\right]\mathrm{,}\end{array}$ (19)

where ${\phi }_T$ is turbulent thermal diffusivity and $\Phi$ is the viscous dissipation.

Averaged concentration equation:

$\begin{array}{l}\frac{\partial \bar{C}}{\partial t}+\bar{w}\frac{\partial \bar{C}}{\partial z}+\frac{\bar{v}}{r}\frac{\partial \bar{C}}{\partial \theta }+\bar{u}\frac{\partial \bar{C}}{\partial r}\\ =\left(D_c+{\phi }_c\right)\left[\frac{1}{r^2}\frac{{\partial }^2\bar{C}}{\partial {\theta }^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \bar{C}}{\partial r}\right)+\frac{{\partial }^2\bar{C}}{\partial z^2}\right]-k_r\left(\bar{C}-C_{\infty }\right)\mathrm{,}\end{array}$ (20)

where $k_r$ is a chemical reaction parameter, $D_c$ is a diffusion coefficient, and ${\phi }_c$ is the turbulent mass diffusivity.

The $\varkappa \text{-}\epsilon$ turbulence model is used to calculate the turbulent viscosity. Given in (21) is the relationship between turbulent viscosity ( ${\nu }_T$ ), turbulent kinetic energy ( $\varkappa$ ), and dissipation (ε).

${\nu }_T=C_{\mu }\frac{f_{\mu }{\varkappa }^2}{\epsilon }\mathrm{,}$ (21)

where $C_{\mu }=0.09$ and $f_{\mu }$ are empirical constant and turbulent model functions respectively.

Derivation of $\varkappa$ and ε equations was explained in [21] . The momentum equations are first replaced with the mean and fluctuating velocities, and then each equation is multiplied by the fluctuating velocities in the corresponding directions to produce the turbulent kinetic energy equation. Time-averaged Reynolds decomposition rules are used together with the decomposed momentum equations. To form a single equation, the three equations obtained from the processes are added. The $\varkappa$ equation is formed by this summation. The $\varkappa$ -equation according to [21] and [22] is given as a vector in (22):

$\rho \frac{\partial \varkappa }{\partial t}+\rho \left[\bar{u}\cdot \nabla \varkappa \right]=\nabla \cdot \left[\left(\mu +\frac{{\mu }_T}{{\sigma }_{\varkappa }}\right)\nabla \varkappa \right]+{\mu }_T\varphi -\rho \epsilon \mathrm{,}$ (22)

where the local change in turbulent kinetic energy is represented by the first item on the left and the convection of turbulent kinetic energy is represented by the second term. The diffusion of turbulent energy is represented by the first term on the right-hand side, the production of turbulent kinetic energy is represented by the second term, and the dissipation of turbulent kinetic energy is represented by the last term. ${\sigma }_{\varkappa }$ is the turbulent Prandtl number for $\varkappa$ -equation.

In vector form, the ε-equation is presented as follows, according to [13] :

$\begin{array}{l}\rho \left[\bar{u}\cdot \nabla \epsilon \right]+\rho \frac{\partial \epsilon }{\partial t}\\ =\nabla \cdot \left[\left(\mu +\frac{{\mu }_T}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]+f_1{C}_{\epsilon \mathrm{,1}}\frac{\epsilon }{\varkappa }\varphi -f_2\rho C_{\epsilon \mathrm{,2}}\frac{{\epsilon }^2}{\varkappa }\mathrm{,}\end{array}$ (23)

where the wall model functions $f_2$ and $f_1$ are utilized to adjust the $\varkappa \text{-}\epsilon$ model’s eddy viscosity behavior.

Turbulent kinetic energy and dissipation equations can be written in cylinder coordinates using Equations (23) and (22). Substituting Equations (19) and (8) into (23) and (22), only the final form is given.

The $\varkappa$ -equation for gas phase:

$\begin{array}{l}\frac{\partial \varkappa }{\partial t}+\bar{u}\frac{\partial \varkappa }{\partial r}+\frac{\bar{v}}{r}\frac{\partial \varkappa }{\partial \theta }+\bar{w}\frac{\partial \varkappa }{\partial z}\\ =\frac{1}{\rho r^2}{\mu }_0\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]\frac{\partial \varkappa }{\partial \theta }\\ +\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+\frac{{\nu }_T}{{\sigma }_{\varkappa }}\right)\left[\frac{1}{r}\frac{\partial }{\partial r}\left(\frac{\partial \varkappa }{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\varkappa }{\partial {\theta }^2}+\frac{{\partial }^2\varkappa }{\partial z^2}\right]\\ +2{\nu }_T\left[{\left(\frac{\bar{u}}{r}+\frac{1}{r}\frac{\partial \bar{v}}{\partial \theta }\right)}^2+{\left(\frac{\partial \bar{u}}{\partial r}\right)}^2+{\left(\frac{\partial \bar{w}}{\partial z}\right)}^2\right]\\ +{\nu }_T\left[{\left(\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{u}}{\partial \theta }-\frac{\bar{v}}{r}+\frac{\partial \bar{v}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{w}}{\partial \theta }+\frac{\partial \bar{v}}{\partial z}\right)}^2\right]-\epsilon \mathrm{,}\end{array}$ (24)

The $\varkappa$ -equation for liquid phase:

$\begin{array}{l}\frac{\partial \varkappa }{\partial t}+\bar{u}\frac{\partial \varkappa }{\partial r}+\frac{\bar{v}}{r}\frac{\partial \varkappa }{\partial \theta }+\bar{w}\frac{\partial \varkappa }{\partial z}\\ =\frac{1}{r^2}\frac{{\mu }_0}{\rho }\frac{\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}}{1+\lambda \left(\bar{T}-T_{\infty }\right)}\frac{\partial \varkappa }{\partial \theta }+\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+\frac{{\nu }_T}{{\sigma }_{\varkappa }}\right)[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \varkappa }{\partial r}\right)\\ +\frac{{\partial }^2\varkappa }{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^2\varkappa }{\partial {\theta }^2}]+2{\nu }_T\left[{\left(\frac{1}{r}\frac{\partial \bar{v}}{\partial \theta }+\frac{\bar{u}}{r}\right)}^2+{\left(\frac{\partial \bar{u}}{\partial r}\right)}^2+{\left(\frac{\partial \bar{w}}{\partial z}\right)}^2\right]\\ +{\nu }_T\left[{\left(\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{u}}{\partial \theta }+\frac{\partial \bar{v}}{\partial r}-\frac{\bar{v}}{r}\right)}^2\right]+{\nu }_T\left[{\left(\frac{1}{r}\frac{\partial \bar{w}}{\partial \theta }+\frac{\partial \bar{v}}{\partial z}\right)}^2\right]-\epsilon \mathrm{,}\end{array}$ (25)

The ε-equation for gas phase:

$\begin{array}{l}\frac{\partial \epsilon }{\partial t}+\bar{u}\frac{\partial \epsilon }{\partial r}+\bar{w}\frac{\partial \epsilon }{\partial z}+\frac{\bar{v}}{r}\frac{\partial \epsilon }{\partial \theta }\\ =\frac{{\mu }_0\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}}{r^2\rho }\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]\frac{\partial \epsilon }{\partial \theta }+\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+\frac{{\nu }_T}{{\sigma }_{\epsilon }}\right)[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \epsilon }{\partial r}\right)\\ +\frac{{\partial }^2\epsilon }{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^2\epsilon }{\partial {\theta }^2}]+2f_1{c}_{\epsilon \mathrm{,1}}\frac{\epsilon }{\varkappa }{\nu }_T\left[{\left(\frac{\partial \bar{u}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{v}}{\partial \theta }+\frac{\bar{u}}{r}\right)}^2+{\left(\frac{\partial \bar{u}}{\partial z}\right)}^2\right]\\ +f_1{c}_{\epsilon \mathrm{,1}}\frac{\epsilon }{\varkappa }{\nu }_T\left[{\left(\frac{1}{r}\frac{\partial \bar{u}}{\partial \theta }+\frac{\partial \bar{v}}{\partial r}-\frac{\bar{v}}{r}\right)}^2+{\left(\frac{\partial \bar{v}}{\partial z}+\frac{1}{r}\frac{\partial \bar{w}}{\partial \theta }\right)}^2+{\left(\frac{\partial \bar{w}}{\partial r}+\frac{\partial \bar{u}}{\partial z}\right)}^2\right]\\ -f_2{c}_{\epsilon \mathrm{,2}}\frac{{\epsilon }^2}{\varkappa }\mathrm{,}\end{array}$ (26)

The liquid phase ε-equation:

$\begin{array}{l}\frac{\partial \epsilon }{\partial t}+\bar{w}\frac{\partial \epsilon }{\partial z}+\frac{\bar{v}}{r}\frac{\partial \epsilon }{\partial \theta }+\bar{u}\frac{\partial \epsilon }{\partial r}\\ =\frac{{\mu }_0\eta \left(a-1\right)f^{a-2}{\theta }^{\eta -1}}{r^2\rho \left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]}\frac{\partial \epsilon }{\partial \theta }+\left(\frac{{\mu }_0{f}^{a-1}}{\rho }+\frac{{\nu }_T}{{\sigma }_{\epsilon }}\right)[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \epsilon }{\partial r}\right)+\frac{{\partial }^2\epsilon }{\partial z^2}\\ +\frac{1}{r^2}\frac{{\partial }^2\epsilon }{\partial {\theta }^2}]+2f_1{c}_{\epsilon \mathrm{,1}}\frac{\epsilon }{\varkappa }{\nu }_T\left[{\left(\frac{\partial \bar{u}}{\partial r}\right)}^2+{\left(\frac{\partial \bar{w}}{\partial z}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{v}}{\partial \theta }+\frac{\bar{u}}{r}\right)}^2\right]\\ +f_1{c}_{\epsilon \mathrm{,1}}\frac{\epsilon }{\varkappa }{\nu }_T\left[{\left(-\frac{\bar{v}}{r}+\frac{1}{r}\frac{\partial \bar{u}}{\partial \theta }+\frac{\partial \bar{v}}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial \bar{w}}{\partial \theta }+\frac{\partial \bar{v}}{\partial z}\right)}^2+{\left(\frac{\partial \bar{u}}{\partial z}+\frac{\partial \bar{w}}{\partial r}\right)}^2\right]\\ -f_2{c}_{\epsilon \mathrm{,2}}\frac{{\epsilon }^2}{\varkappa }\mathrm{,}\end{array}$ (27)

where ${\sigma }_{\varkappa }$ , $c_{\epsilon \mathrm{,1}}$ , ${\sigma }_{\epsilon }$ , $c_{\epsilon \mathrm{,2}}$ and $c_{\mu }$ are [23] model constants, which are listed as $c_{\mu }=0.09$ , ${\sigma }_{\varkappa }=1$ , $c_{\epsilon \mathrm{,1}}=1.44$ , ${\sigma }_{\epsilon }=1.3$ $c_{\epsilon \mathrm{,2}}=1.92$ .

According to [23] , the following modeling functions have been incorporated because turbulent eddies cause damping close to the walls.

$f_{\mu }={\left[1-\exp \left(-0.0165R_{\varkappa }\right)\right]}^2\left(1+\frac{20.5}{R_t}\right)\mathrm{,}$ (28)

$f_1={\left(1+\frac{0.05}{f_{\mu }}\right)}^3\mathrm{,}$ (29)

$f_2=1-\exp \left(-R_t^2\right)\mathrm{,}$ (30)

where $R_t$ and $R_{\varkappa }$ are turbulent Reynolds numbers.

2.2. Boundary Conditions

The following are relevant boundary conditions to the problem under consideration as applied successfully by [13] and [21] . The geothermal pipe’s walls and the gas-liquid interface don’t have a slip condition.

At the pipe wall, $r=\frac{D^{+}}{2}$ :

$\bar{u}=\bar{w}=\bar{v}=0,\bar{T}=T_w,\varkappa =0,\bar{C}=C_w,\epsilon =0$

The fluid interface is regarded as a moving wall when $r=0$ .

$\bar{C}=C_{\infty },\bar{u}=U_{\infty },\bar{T}=T_{\infty },\bar{v}=U_{\infty },\frac{\partial \varkappa }{\partial r}=0,\bar{w}=U_{\infty },\frac{\partial \epsilon }{\partial r}=0.$

$r=\frac{D^{-}}{2}$ at the wall of the pipe

$\bar{u}=\bar{w}=0,\varkappa =0,\bar{C}=C_w,\epsilon =0,\bar{v}=0,\bar{T}=T_w$

$z=0$ at the entrance

$\bar{w}=U_{\infty },\bar{C}=C_{\infty },\bar{v}=0,\bar{u}=0,\bar{T}=T_{\infty },$

$\varkappa =\frac{3}{2}{\left(U_{\infty }I\right)}^2\mathrm{,}I=0.16R{e}^{-\frac{1}{8}}\mathrm{,}\epsilon =C_{\mu }^{\frac{3}{4}}\frac{{\varkappa }^{\frac{3}{2}}}{0.07D}\mathrm{,}$

where I is the turbulent intensity.

According to [24] , at the outlet $z=\infty$ , all variable gradients are zero in the axial direction i.e.,

$\frac{\partial \bar{u}}{\partial z}=0,\frac{\partial \bar{v}}{\partial z}=0,\frac{\partial \bar{w}}{\partial z}=0,\frac{\partial \bar{T}}{\partial z}=0,\frac{\partial \bar{C}}{\partial z}=0,\frac{\partial \varkappa }{\partial z}=0,\frac{\partial \epsilon }{\partial z}=0$

2.3. Local Nusselt Number, Coefficient of Skin Friction and Sherwood Number

The local skin friction coefficient, Nusselt number, and Sherwood number are quantities of interest for the current engineering problem in terms of practical application. The definition of the mean wall shear stress (skin friction coefficient) is:

$C_{f_{\left(r\right)}}=-\frac{2{\tau }_w}{\rho U_{\infty }^2},$ (31)

where the skin friction on the pipe wall is given by,

${\tau }_w=\mu {\left(\frac{\partial \bar{u}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (32)

Substituting Equations (32) into (31) we have;

$C_{f_{\left(r\right)}}=-\frac{2\mu }{\rho U_{\infty }^2}{\left(\frac{\partial \bar{u}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (33)

Since viscosity is a variable, by Equations (7) and (8) Equation (33) can be expressed as;

$C_{f_{\left(r\right)}}=-\frac{2{\mu }_0}{\rho U_{\infty }^2}{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]}^{a-1}{\left(\frac{\partial \bar{u}}{\partial r}\right)}_{r=\frac{D}{2}},$ (34)

$C_{f_{\left(r\right)}}=-\frac{2{\mu }_0}{\rho U_{\infty }^2}\frac{{\theta }^{\eta \left(a-1\right)}}{{\left[1+\lambda \left(\bar{T}-T_{\infty }\right)\right]}^{a-1}}{\left(\frac{\partial \bar{u}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (35)

Skin friction for the gaseous and liquid phases are represented by Equations (34) and (35), respectively.

Another crucial parameter of the current engineering problem is the local Nusselt number which is defined as:

$N_u=-\frac{D{q}_w^{\star }}{{\kappa }^{\star }\left(T_w-T_{\infty }\right)}\mathrm{,}$ (36)

where $q_w^{\star }$ is heat flux given as follows;

$q_w^{\star }={\kappa }^{\star }{\left(\frac{\partial \bar{T}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (37)

$Nu=-\frac{D}{T_w-T_{\infty }}{\left(\frac{\partial \bar{T}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (38)

Sherwood number is defined as:

$Sh=-\frac{D}{C_w-C_{\infty }}{\left(\frac{\partial \bar{C}}{\partial r}\right)}_{r=\frac{D}{2}}.$ (39)

3. Non-Dimensionalization

To solve the problem, the governing equations are converted to non-dimensionless forms. The following dimensionless quantities are introduced:

$\begin{array}{l}{u}^{\star }=\frac{\bar{u}}{U_{\infty }}\mathrm{,}{v}^{\star }=\frac{\bar{v}}{U_{\infty }}\mathrm{,}{w}^{\star }=\frac{\bar{w}}{U_{\infty }}\mathrm{,}{t}^{\star }=\frac{t{U}_{\infty }}{D}\mathrm{,}{p}^{\star }=\frac{\bar{p}}{\rho U_{\infty }^2},{\theta }^{\star }=\theta ,z^{\star }=\frac{z}{D},\\ {\nu }_T^{\star }=\frac{{\nu }_T}{D{U}_{\infty }},r^{\star }=\frac{r}{D},T^{\star }=\frac{\bar{T}-T_{\infty }}{T_w-T_{\infty }},{\phi }_T^{\star }=\frac{{\phi }_T}{D{U}_{\infty }},{\epsilon }^{\star }=\frac{\epsilon D}{U_{\infty }^3},{\varkappa }^{\star }=\frac{\varkappa }{U_{\infty }^2},\\ {\phi }_c^{\star }=\frac{{\phi }_c}{D{U}_{\infty }},C^{\star }=\frac{\bar{C}-C_{\infty }}{C_w-C_{\infty }}.\end{array}$ (40)

Using these non-dimensional values, the boundary conditions and the governing equations are then non-dimensionalized and represented as follows.

For the gas phase, the dimensionless momentum equation in the r direction is:

$\begin{array}{l}\frac{\partial u^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial u^{\star }}{\partial r^{\star }}-\frac{v^{\star 2}}{r^{\star }}\\ =\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}+{\nu }_T^{\star }\right)[-\frac{u^{\star }}{r^{\star 2}}+\frac{1}{r^{\star }}\frac{\partial }{\partial r^{\star }}\left(r^{\star }\frac{\partial u^{\star }}{\partial r^{\star }}\right)+\frac{{\partial }^2{u}^{\star }}{\partial z^{\star 2}}\\ +\frac{1}{r^{\star 2}}\frac{{\partial }^2{u}^{\star }}{\partial {\theta }^{\star 2}}-\frac{2}{r^{\star 2}}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}]-\frac{\partial p^{\star }}{\partial r^{\star }}+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}\\ \times \left[-\frac{v^{\star }}{r^{\star 2}}+\frac{1}{r^{\star 2}}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}+\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial r^{\star }}\right]+\cos \alpha \left(T^{\star }\frac{G_{r_{\left(T\right)}}}{R{e}^2}+C^{\star }\frac{G_{r_{\left(c\right)}}}{R_e^2}\right)\mathrm{,}\end{array}$ (41)

where $Re=\frac{U_{\infty }D}{{\nu }_{\infty }}$ is Reynolds Number, $G_{r_{\left(T\right)}}=\frac{g{\beta }_T^{\star }\left(T_w-T_{\infty }\right)}{{\nu }_{\infty }^2}$ is Thermal Grasholf Number and $G_{r_{\left(c\right)}}=\frac{g{\beta }_c^{\star }\left(C_w-C_{\infty }\right)}{{\nu }_{\infty }^2}$ is solutal Grasholf Number.

Dimensionless momentum equation for the liquid phase in radial direction:

$\begin{array}{l}\frac{\partial u^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial u^{\star }}{\partial r^{\star }}-\frac{v^{\star 2}}{r^{\star }}\\ =\left(\frac{{\theta }^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}+{\nu }_T^{\star }\right)[-\frac{2}{r^{\star 2}}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}+\frac{1}{r^{\star }}\frac{\partial }{\partial r^{\star }}\left(r^{\star }\frac{\partial u^{\star }}{\partial r^{\star }}\right)\\ +\frac{{\partial }^2{u}^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{u}^{\star }}{\partial {\theta }^{\star 2}}-\frac{u^{\star }}{r^{\star 2}}]-\frac{\partial p^{\star }}{\partial r^{\star }}+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}\\ \times \left[-\frac{v^{\star }}{r^{\star 2}}+\frac{1}{r^{\star 2}}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}+\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial r^{\star }}\right]+\cos \alpha \left(T^{\star }\frac{G_{r_{\left(T\right)}}}{R{e}^2}+C^{\star }\frac{G_{r_{\left(C\right)}}}{R{e}^2}\right)\mathrm{,}\end{array}$ (42)

Dimensionless momentum equation for the gas phase in the θ-direction:

$\begin{array}{l}\frac{\partial v^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial v^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial v^{\star }}{\partial r^{\star }}+\frac{v^{\star }{u}^{\star }}{r^{\star }}\\ =\frac{2\eta \left(n-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}\left[\frac{1}{r^{\star 2}}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}+\frac{u^{\star }}{r^{\star 2}}\right]-\frac{1}{r^{\star }}\frac{\partial p^{\star }}{\partial {\theta }^{\star }}\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}+{\nu }_T^{\star }\right)[\frac{{\partial }^2{v}^{\star }}{\partial z^{\star 2}}+\frac{2}{r^{\star 2}}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}+\frac{1}{r^{\star }}\frac{\partial }{\partial r^{\star }}\left(r^{\star }\frac{\partial v^{\star }}{\partial r^{\star }}\right)\\ +\frac{1}{r^{\star 2}}\frac{{\partial }^2{v}^{\star }}{\partial {\theta }^{\star 2}}-\frac{v^{\star }}{r^{\star 2}}]\mathrm{,}\end{array}$ (43)

Dimensionless momentum equation for the liquid phase in the θ-direction:

$\begin{array}{l}\frac{\partial v^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial v^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial v^{\star }}{\partial r^{\star }}+\frac{v^{\star }{u}^{\star }}{r^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}\left[2\frac{u^{\star }}{r^{\star 2}}+\frac{2}{r^{\star 2}}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}\right]-\frac{1}{r^{\star }}\frac{\partial p^{\star }}{\partial {\theta }^{\star }}+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}+{\nu }_T^{\star }\right)\\ \times \left[\frac{1}{r^{\star }}\frac{\partial }{\partial r^{\star }}\left(r^{\star }\frac{\partial v^{\star }}{\partial r^{\star }}\right)+\frac{{\partial }^2{v}^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{v}^{\star }}{\partial {\theta }^{\star 2}}-\frac{v^{\star }}{r^{\star 2}}+\frac{2}{r^{\star 2}}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}\right]\mathrm{,}\end{array}$ (44)

Dimensionless momentum equation for the liquid phase in the z-direction:

$\begin{array}{l}\frac{\partial w^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial w^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial w^{\star }}{\partial r^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{r^{\star 2}Re{\left[1+\lambda T^{\star }\right]}^{a-1}}\left(\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+r^{\star }\frac{\partial v^{\star }}{\partial z^{\star }}\right)-\frac{\partial p^{\star }}{\partial z^{\star }}+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}+{\nu }_T^{\star }\right)\\ \times \left[\frac{1}{r^{\star }}\frac{\partial }{\partial r^{\star }}\left(r^{\star }\frac{\partial w^{\star }}{\partial r^{\star }}\right)+\frac{{\partial }^2{w}^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{\partial w^{\star }}{\partial {\theta }^{\star 2}}\right]+\sin \alpha \left[T^{\star }\frac{G_{r_{\left(T\right)}}}{R{e}^2}+C^{\star }\frac{G_{r_{\left(C\right)}}}{R{e}^2}\right]\mathrm{,}\end{array}$ (45)

Dimensionless momentum equation for the gas phase in the z direction:

$\begin{array}{l}\frac{\partial w^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial w^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial w^{\star }}{\partial r^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}\frac{1}{r^{\star 2}}\left(r^{\star }\frac{\partial v^{\star }}{\partial z^{\star }}+\frac{\partial w^{\star }}{\partial {\theta }^{\star }}\right)-\frac{\partial p^{\star }}{\partial z^{\star }}\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}+{\nu }_T^{\mathrm{<mo>\; *\; </mo>}}\right)\left[\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial r^{\star }}+\frac{{\partial }^2{w}^{\star }}{\partial r^{\star 2}}+\frac{{\partial }^2{w}^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{\partial w^{\star }}{\partial {\theta }^{\star 2}}\right]\\ +\sin \alpha \left[T^{\star }\frac{G_{r_{\left(T\right)}}}{R{e}^2}+C^{\star }\frac{G_{r_{\left(C\right)}}}{R{e}^2}\right]\mathrm{,}\end{array}$ (46)

The energy equation in its dimensionless form for the liquid phase:

$\begin{array}{l}\frac{\partial T^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial T^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial T^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial T^{\star }}{\partial r^{\star }}\\ =-\frac{\omega }{P_r{R}_e}\left[{\left(\frac{\partial T^{\star }}{\partial r^{\star }}\right)}^2+\frac{1}{r^{\star 2}}{\left(\frac{\partial T^{\star }}{\partial {\theta }^{\star }}\right)}^2+{\left(\frac{\partial T^{\star }}{\partial z^{\star }}\right)}^2\right]+\left(\frac{1+\omega T^{\star }}{P_{r}Re}+{\phi }_T^{\star }\right)\\ \times \left(\frac{{\partial }^2{T}^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star }}\frac{\partial T^{\star }}{\partial r^{\star }}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{T}^{\star }}{\partial {\theta }^{\star 2}}+\frac{{\partial }^2{T}^{\star }}{\partial r^{\star 2}}\right)+\frac{E_c}{Re}\frac{{\theta }^{\eta \left(a-1\right)}}{{\left[1+\lambda T^{\star }\right]}^{a-1}}\\ \times \left\{2\left[{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial w^{\star }}{\partial z^{\star }}\right)}^2+{\left(\frac{u^{\star }}{r^{\star }}+\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}\right)}^2\right]+{\left(\frac{\partial v^{\star }}{\partial r^{\star }}+\frac{1}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}-\frac{v^{\star }}{r^{\star }}\right)}^2\right\}\\ +\frac{E_c}{Re}\frac{{\theta }^{\eta \left(a-1\right)}}{{\left[1+\lambda T^{\star }\right]}^{a-1}}\left[{\left(\frac{\partial v^{\star }}{\partial z^{\star }}+\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}\right)}^2+{\left(\frac{\partial w^{\star }}{\partial r^{\star }}+\frac{\partial u^{\star }}{\partial z^{\star }}\right)}^2\right]+E_c{\epsilon }^{\star }\mathrm{,}\end{array}$ (47)

where $E_c=\frac{U_{\infty }^2}{C_p\left(T_w-T_{\infty }\right)}$ is Eckert Number and $P_r=\frac{{\mu }_{\infty }{C}_p}{{\kappa }^{\mathrm{<mo>\; *\; </mo>}}}$ is Prandtl Number.

Dimensionless form of the gas phase energy equation:

$\begin{array}{l}\frac{\partial T^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial T^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial T^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial T^{\star }}{\partial r^{\star }}\\ =-\frac{\omega }{P_{r}Re}\left[{\left(\frac{\partial T^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial T^{\star }}{\partial z^{\star }}\right)}^2+\frac{1}{r^{\star 2}}{\left(\frac{\partial T^{\star }}{\partial {\theta }^{\star }}\right)}^2\right]+\left(\frac{1+\omega T^{\star }}{P_r{R}_e}+{\phi }_T^{\star }\right)\end{array}$

$\begin{array}{l}\times \left(\frac{1}{r^{\star }}\frac{\partial T^{\star }}{\partial r^{\star }}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{T}^{\star }}{\partial {\theta }^{\star 2}}+\frac{{\partial }^2{T}^{\star }}{\partial r^{\star 2}}+\frac{{\partial }^2{T}^{\star }}{\partial z^{\star 2}}\right)+\frac{Ec{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}\\ \times \left\{2\left[{\left(\frac{\partial w^{\star }}{\partial z^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{u^{\star }}{r^{\star }}+\frac{1}{r}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}\right)}^2\right]+{\left(-\frac{v^{\star }}{r^{\star }}+\frac{1}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\mathrm{<mo>\; *\; </mo>}}}+\frac{\partial v^{\star }}{\partial r^{\star }}\right)}^2\right\}\\ +\frac{Ec{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}\left[{\left(\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+\frac{\partial v^{\star }}{\partial z^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{\partial w^{\star }}{\partial r^{\star }}\right)}^2\right]+E_c{\epsilon }^{\star }\mathrm{,}\end{array}$ (48)

Dimensionless concentration equation:

$\begin{array}{l}\frac{\partial C^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial C^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial C^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial C^{\star }}{\partial r^{\star }}\\ =\left(\frac{1}{S_{c}Re}+{\phi }_c^{\star }\right)\left[\frac{1}{r^{\star 2}}\frac{{\partial }^2{C}^{\star }}{\partial {\theta }^{\star 2}}+\frac{1}{r^{\star }}\frac{\partial C^{\star }}{\partial r^{\star }}+\frac{{\partial }^2{C}^{\star }}{\partial r^{\star 2}}+\frac{{\partial }^2{C}^{\star }}{\partial z^{\star 2}}\right]-k_{c}Re{C}^{\star }\end{array}$ (49)

where $Sc=\frac{{\nu }_{\infty }}{D_c}$ is the Schmidt Number and $k_c=\frac{k_r{\nu }_{\infty }}{U_{\infty }^2}$ is the Chemical Reaction Parameter.

Dimensionless $\varkappa$ -equation for liquid phase:

$\begin{array}{l}\frac{\partial {\varkappa }^{\star }}{\partial t^{\star }}+u^{\star }\frac{\partial {\varkappa }^{\star }}{\partial r^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial {\varkappa }^{\star }}{\partial {\theta }^{\star }}+w^{\star }\frac{\partial {\varkappa }^{\star }}{\partial z^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{r^{\star 2}Re{\left[1+\lambda T^{\star }\right]}^{a-1}}\frac{\partial {\varkappa }^{\star }}{\partial {\theta }^{\star }}+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T^{\star }\right]}^{a-1}}+\frac{{\nu }_T^{\star }}{{\sigma }_{\varkappa }}\right)[\frac{1}{r^{\star }}\frac{\partial {\varkappa }^{\star }}{\partial r^{\star }}+\frac{{\partial }^2{\varkappa }^{\star }}{\partial r^{2\star }}\\ +\frac{{\partial }^2{\varkappa }^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{\varkappa }^{\star }}{\partial {\theta }^{\star 2}}]+2{\nu }_T^{\star }\left[{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial w^{\star }}{\partial z^{\star }}\right)}^2+{\left(\frac{u^{\star }}{r}+\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}\right)}^2\right]\\ +{\nu }_T^{\star }\left[{\left(\frac{1}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}-\frac{v^{\star }}{r^{\star }}+\frac{\partial v^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+\frac{\partial v^{\star }}{\partial z^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{\partial w^{\star }}{\partial r^{\star }}\right)}^2\right]-{\epsilon }^{\star }\mathrm{,}\end{array}$ (50)

Dimensionless $\varkappa$ -equation for gas phase:

$\begin{array}{l}\frac{\partial {\varkappa }^{\star }}{\partial t^{\star }}+u^{\star }\frac{\partial {\varkappa }^{\star }}{\partial r^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial {\varkappa }^{\star }}{\partial {\theta }^{\star }}+w^{\star }\frac{\partial {\varkappa }^{\star }}{\partial z^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{r^{\star 2}{R}_e}\frac{\partial {\varkappa }^{\star }}{\partial {\theta }^{\star }}+\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}+\frac{{\nu }_T^{\star }}{{\sigma }_{\varkappa }}\right)\\ \times \left(\frac{1}{r^{\star }}\frac{\partial {\varkappa }^{\star }}{\partial r^{\star }}+\frac{{\partial }^2{\varkappa }^{\mathrm{<mo>\; *\; </mo>}}}{\partial z^{\mathrm{<mo>\; *\; </mo>2}}}+\frac{{\partial }^2{\varkappa }^{\star }}{\partial r^{\star 2}}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{\varkappa }^{\star }}{\partial {\theta }^{\star 2}}\right)\\ +2{\nu }_T^{\star }\left[{\left(\frac{u^{\star }}{r^{\star }}+\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial w^{\star }}{\partial z^{\star }}\right)}^2\right]\\ +{\nu }_T^{\star }\left[{\left(\frac{1}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}-\frac{v^{\star }}{r^{\star }}+\frac{\partial v^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{\partial w^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial v^{\star }}{\partial z^{\star }}+\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}\right)}^2\right]-{\epsilon }^{\star }\end{array}$ (51)

Dimensionless ε-equation for liquid phase:

$\begin{array}{l}\frac{\partial {\epsilon }^{\star }}{\partial t^{\star }}+w^{\star }\frac{\partial {\epsilon }^{\star }}{\partial z^{\star }}+\frac{v^{\star }}{r^{\star }}\frac{\partial {\epsilon }^{\star }}{\partial {\theta }^{\star }}+u^{\star }\frac{\partial {\epsilon }^{\star }}{\partial r^{\star }}\\ =\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{r^{\star 2}{R}_e{\left[1+\lambda T^{\star }\right]}^{a-1}}\frac{\partial {\epsilon }^{\star }}{\partial {\theta }^{\star }}+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}+\frac{{\nu }_T^{\star }}{{\sigma }_{\epsilon }}\right)[\frac{1}{r^{\mathrm{<mo>\; *\; </mo>}}}\frac{\partial {\epsilon }^{\star }}{\partial r^{\mathrm{<mo>\; *\; </mo>}}}+\frac{{\partial }^2{\epsilon }^{\star }}{\partial r^{\star 2}}\\ +\frac{{\partial }^2{\epsilon }^{\star }}{\partial z^{\star 2}}+\frac{1}{r^{\star 2}}\frac{{\partial }^2{\epsilon }^{\star }}{\partial {\theta }^{\star 2}}]+2f_1{c}_{\epsilon \mathrm{,1}}\frac{{\epsilon }^{\star }}{{\varkappa }^{\star }}{\nu }_T^{\star }\left[{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{1}{r^{\star }}\frac{\partial v^{\star }}{\partial {\theta }^{\star }}+\frac{u^{\star }}{r^{\star }}\right)}^2+{\left(\frac{\partial w^{\star }}{\partial z^{\star }}\right)}^2\right]\\ +f_1{c}_{\epsilon \mathrm{,1}}\frac{{\epsilon }^{\star }}{{\varkappa }^{\star }}{\nu }_T^{\star }\left[{\left(\frac{1}{r^{\star }}\frac{\partial u^{\star }}{\partial {\theta }^{\star }}-\frac{v^{\star }}{r^{\star }}+\frac{\partial v^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{\partial u^{\star }}{\partial z^{\star }}+\frac{\partial w^{\star }}{\partial r^{\star }}\right)}^2+{\left(\frac{1}{r^{\star }}\frac{\partial w^{\star }}{\partial {\theta }^{\star }}+\frac{\partial v^{\star }}{\partial z^{\star }}\right)}^2\right]\\ -f_2{c}_{\epsilon \mathrm{,2}}\frac{{\epsilon }^{\star 2}}{{\varkappa }^{\star }}\mathrm{.}\end{array}$ (52)

${\nu }_{T^{\star }}=c_{\mu }{f}_{\mu }\frac{{\varkappa }^{\star 2}}{{\epsilon }^{\star }}$ (53)

The equivalent dimensionless form of the Skin friction coefficient, the local Nusselt number, and the Sherwood number.

$C_{f_{\left(r\right)}}^{\star }=-\frac{2{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T^{\star }\right]}^{a-1}}{Re}{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}_{r=0.5}\mathrm{,}$ (54)

$C_{f_{\left(r\right)}}^{\star }=-\frac{2{\theta }^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T^{\star }\right]}^{a-1}}{\left(\frac{\partial u^{\star }}{\partial r^{\star }}\right)}_{r=0.5}\mathrm{,}$ (55)

Equations (54) and (55) represent, for the gas and liquid phases, the dimensionless form of skin friction.

$Nu=-{\left(\frac{\partial T^{\star }}{\partial r^{\star }}\right)}_{r^{\star }=0.5}.$ (56)

$Sh=-{\left(\frac{\partial T^{\star }}{\partial r^{\star }}\right)}_{r^{\star }=0.5}$ (57)

Dimensionless Boundary Conditions

Pipe wall, $r^{\star }={0.5}^{+}$

$v^{\star }=0,w^{\star }=0,u^{\star }=0,{\varkappa }^{\star }=0,C^{\star }=1,{\epsilon }^{\star }=0,T^{\star }=1$

The fluid interface, $r^{\star }\mathrm{<mo>\; =\; </mo>0}$

$u^{\star }=1,T^{\star }=0,v^{\star }=1,C^{\star }=0,w^{\star }=1,\frac{\partial {\varkappa }^{\star }}{\partial r^{\star }}=0,\frac{\partial {\epsilon }^{\star }}{\partial r^{\star }}=0$

$r^{\star }={0.5}^{-}$ , wall of the pipe

$u^{\star }=0,w^{\star }=0,v^{\star }=0,{\varkappa }^{\star }=0,{\epsilon }^{\star }=0,T^{\star }=1,C^{\star }=1$

Inlet, $z^{\star }=0$

$w^{\star }=\mathrm{1,}{v}^{\star }=\mathrm{0,}{u}^{\star }=\mathrm{0,}{T}^{\star }=\mathrm{0,}{C}^{\star }=\mathrm{0,}{\varkappa }^{\star }=\frac{3}{2}{I}^2\mathrm{,}I=0.016R{e}^{-\frac{1}{8}}\mathrm{,}{\epsilon }^{\star }=C_{\mu }^{\frac{3}{4}}\frac{{\varkappa }^{\frac{3}{2}}}{0.07}\mathrm{.}$

Outlet $z^{\star }=\infty$

$\frac{\partial u^{\star }}{\partial z^{\star }}=0,\frac{\partial C^{\star }}{\partial z^{\star }}=0,\frac{\partial w^{\star }}{\partial z^{\star }}=0,\frac{\partial {\varkappa }^{\star }}{\partial z^{\star }}=0,\frac{\partial v^{\star }}{\partial z^{\star }}=0,\frac{\partial T^{\star }}{\partial z^{\star }}=0,\frac{\partial {\epsilon }^{\star }}{\partial z^{\star }}=0$

4. Methodology

The governing Equations (41) to (52) cannot be solved analytically due to their complex nature. An explicit finite difference method has been employed to solve the coupled nonlinear PDEs that govern the flow. PDEs have been solved using the finite difference approach due to its accuracy, stability, and quick convergence. The derivatives of the dependent variables appearing in the PDEs are approximated by Finite differences. A cylindrical mesh has been used to discretize the domain. Let’s say that $u\left(r\mathrm{,}\theta \mathrm{,}z\mathrm{,}t\right)$ is described at the point $\left(i\mathrm{,}n\mathrm{,}k\mathrm{,}s\right)$ and $u\left(r\mathrm{,}\theta \mathrm{,}z\mathrm{,}t\right)$ as $u_{i\mathrm{,}n\mathrm{,}k}^s$ . In order to build the mesh, it’s assumed there exist S points along t, Z points along z, N points along θ, and I points along r.

Let the step size along r be Δr, along θ be Δθ, along z be Δz and along t be Δt. $r_i=r_0+i\Delta r$ , $z_k=z_0+k\Delta z$ , ${\theta }_n={\theta }_0+n\Delta \theta$ and $t_s=t_0+s\Delta t$ where $s=1,2,3,4,\cdots ,S$ , $i=1,2,3,4,\cdots ,I$ and $k=\mathrm{1,2,3,4},\cdots \mathrm{,}Z$ , $n=1,2,3,4,\cdots ,N$ .

The forward finite difference scheme with respect to time is

$\frac{\partial u}{\partial t}=\frac{u_{i,n,k}^{s+1}-u_{i,n,k}^s}{\Delta t}$ (58)

The central finite difference approximation scheme with respect to space

$\frac{\partial u}{\partial r}=\frac{u_{i+1,n,k}^s-u_{i-1,n,k}^s}{2\Delta r},$ (59)

$\frac{{\partial }^{2}u}{\partial r^2}=\frac{u_{i+1,n,k}^s-2u_{i,n,k}^m+U_{i-1,n,k}^s}{{\left(\Delta r\right)}^2}.$ (60)

Substituting Equations (58) to (60) into Equations (41) to (52) we have the following linear algebraic equations:

Discrete form of the momentum equation for the liquid phase in the radial direction:

$\begin{array}{l}\frac{u_{i\mathrm{,}n\mathrm{,}k}^{s+1}-u_{i\mathrm{,}n\mathrm{,}k}^s}{\Delta t}+u_{i\mathrm{,}n\mathrm{,}k}^s\frac{u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}j\mathrm{,}k}^s}{2\Delta r}+\frac{v_{i\mathrm{,}n\mathrm{,}k}^s}{r_i}\frac{u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s}{2\Delta \theta }\\ +w_{i\mathrm{,}n\mathrm{,}k}^s\frac{u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s}{2\Delta z}-\frac{{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2}{r_i}\\ =-\frac{\partial p}{\partial r}+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+{\nu }_T\right)(\frac{1}{r_i}\frac{u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s}{2\Delta r}\\ +\frac{u_{i+\mathrm{1,}n\mathrm{,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i-\mathrm{1,}n\mathrm{,}k}^s}{{\left(\Delta r\right)}^2}-\frac{u_{i\mathrm{,}n\mathrm{,}k}^s}{r_i^2}+\frac{u_{i\mathrm{,}n+\mathrm{1,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i\mathrm{,}n-\mathrm{1,}k}^s}{{\left(\Delta \theta \right)}^2})\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+{\nu }_T\right)\left(\frac{u_{i\mathrm{,}n\mathrm{,}k+1}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i\mathrm{,}n\mathrm{,}k-1}^s}{{\left(\Delta z\right)}^2}-\frac{2}{r_i^2}\frac{v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}j-\mathrm{1,}k}^m}{2\Delta \theta }\right)\end{array}$

$\begin{array}{l}+\frac{1}{R_e}\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}\left(\frac{1}{r_i^2}\frac{u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s}{2\Delta \theta }-\frac{v_{i\mathrm{,}n\mathrm{,}k}^s}{r_i^2}+\frac{1}{r_i}\frac{v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s}{2\Delta r}\right)\\ +\cos \alpha \left(T_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(T\right)}}{R{e}^2}+C_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(c\right)}}{R{e}^2}\right)\mathrm{,}\end{array}$ (61)

multiplying Equation (61) by Δt and rearranging so that the values at time s + 1 are on the left.

$\begin{array}{c}{u}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2-\Delta t\frac{\partial p}{\partial r}\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\end{array}$

$\begin{array}{l}+\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2}{u}_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{{\left(\Delta \theta \right)}^2}(u_{i\mathrm{,}n+\mathrm{1,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s\\ +u_{i\mathrm{,}n-\mathrm{1,}k}^s)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)-\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-V_{i\mathrm{,}n-\mathrm{1,}k}^s\right)]\\ +\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}(\frac{\Delta t}{2r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)-\frac{\Delta t}{r_i^2}{V}_{i\mathrm{,}n\mathrm{,}k}^s\\ +\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right))+\Delta t\cos \alpha \left(T_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(T\right)}}{R{e}^2}+C_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(c\right)}}{R{e}^2}\right)\mathrm{,}\end{array}$ (62)

$u_{i+\mathrm{1,}n\mathrm{,}k}^s\mathrm{,}{u}_{i\mathrm{,}n\mathrm{,}k}^s$ and $u_{i-\mathrm{1,}n\mathrm{,}k}^s$ hence $u_{i\mathrm{,}n\mathrm{,}k}^{s+1}$ can be computed explicitly.

Discrete form of the momentum equation for the gas phase in the radial direction:

$\begin{array}{c}{u}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(u_{i\mathrm{,}j\mathrm{,}k+1}^m-u_{i\mathrm{,}j\mathrm{,}k-1}^m\right)w_{i\mathrm{,}j\mathrm{,}k}^m+\frac{\Delta t}{r_i}{\left(v_{i\mathrm{,}j\mathrm{,}k}^m\right)}^2-\Delta t\frac{\partial p}{\partial r}\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}j\mathrm{,}k}^m\right]}^{a-1}}{R_e}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\end{array}$

$\begin{array}{l}+\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2}{u}_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{{\left(\Delta \theta \right)}^2}(u_{i\mathrm{,}n+\mathrm{1,}k}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s\\ +u_{i\mathrm{,}n-\mathrm{1,}k}^s)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-2u_{i\mathrm{,}n\mathrm{,}k}^s+u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)-\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)]\\ +\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}(\frac{\Delta t}{2r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)-\frac{\Delta t}{r_i^2}{v}_{i\mathrm{,}n\mathrm{,}k}^s\\ +\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right))+\Delta t\cos \alpha \left(T_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(T\right)}}{R{e}^2}+C_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(c\right)}}{R{e}^2}\right)\mathrm{,}\end{array}$ (63)

Momentum equation for the liquid phase in discrete form in the θ direction:

$\begin{array}{c}{v}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=v_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k+1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{r_i}{u}_{i\mathrm{,}n\mathrm{,}k}^s{v}_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{r_i}\frac{\partial p}{\partial \theta }+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}\\ \times \frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s+v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2}{v}_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}(v_{i\mathrm{,}n+\mathrm{1,}k}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s\\ +v_{i\mathrm{,}n-\mathrm{1,}k}^s)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s+v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)+\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)]\end{array}$ (64)

$v_{i\mathrm{,}n+\mathrm{1,}k}^s\mathrm{,}{v}_{i\mathrm{,}n\mathrm{,}k}^s$ and $v_{i\mathrm{,}n-\mathrm{1,}k}^s$ are known thus $v_{i\mathrm{,}n\mathrm{,}k}^{s+1}$ can be computed explicitly.

Momentum equation for the gas phase in discrete form in the θ direction:

$\begin{array}{c}{v}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=v_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^m-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k+1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{r_i}{u}_{i\mathrm{,}n\mathrm{,}k}^s{v}_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{r_i}\frac{\partial p}{\partial \theta }\\ +\frac{1}{R_e}\left[\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}\right]\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s+v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2}{v}_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}(v_{i\mathrm{,}n+\mathrm{1,}k}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s\\ +v_{i\mathrm{,}n-\mathrm{1,}k}^s)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-2v_{i\mathrm{,}n\mathrm{,}k}^s+v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)+\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)]\mathrm{,}\end{array}$ (65)

Momentum equation in z-direction for liquid phase in discrete form:

$\begin{array}{c}{w}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=w_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s-\Delta t\frac{\partial p}{\partial z}+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}\\ \times \left(\frac{\Delta t}{2r_i^2\Delta \theta }\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{2r_i\Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k+1}^s\right)\right)\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]+\Delta t\sin \alpha \left[T_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(T\right)}}{R_e^2}+C_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(c\right)}}{R_e^2}\right]\mathrm{,}\end{array}$ (66)

Momentum equation in z-direction for gaseous phase in discrete form:

$\begin{array}{c}{w}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=w_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s-\Delta t\frac{\partial p}{\partial z}+\frac{1}{R_e}\left[\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T_{i\mathrm{,}j\mathrm{,}k}^m\right]}^{a-1}\right]\\ \times \left(\frac{\Delta t}{2r_i^2\Delta \theta }\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{2r_i\Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k+1}^s\right)\right)\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}+{\nu }_T\right)[\frac{\Delta t}{2r_i\Delta r}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-2w_{i\mathrm{,}n\mathrm{,}k}^s+w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]+\Delta t\sin \alpha \left[T_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(T\right)}}{R_e^2}+C_{i\mathrm{,}n\mathrm{,}k}^s\frac{G_{r\left(c\right)}}{R_e^2}\right]\mathrm{,}\end{array}$ (67)

Energy equation in discrete form:

$\begin{array}{c}{T}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=T_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(T_{i+\mathrm{1,}n\mathrm{,}k}^s-T_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(T_{i\mathrm{,}n+\mathrm{1,}k}^s-T_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(T_{i\mathrm{,}n\mathrm{,}k+1}^s-T_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\omega }{P_r{R}_e}[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(T_{i+\mathrm{1,}n\mathrm{,}k}^s-T_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\\ +\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(T_{i\mathrm{,}n+\mathrm{1,}k}^s-T_{i\mathrm{,}n-\mathrm{1,}k}^{s}2\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(T_{i\mathrm{,}n\mathrm{,}k+1}^s-T_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]\\ +\left[\frac{1+\omega T_{i\mathrm{,}n\mathrm{,}k}^s}{P_r{R}_e}+{\phi }_T\right][\frac{\Delta t}{2r_i\Delta r}\left(T_{i+\mathrm{1,}n\mathrm{,}k}^s-T_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(T_{i+\mathrm{1,}n\mathrm{,}k}^s-2T_{i\mathrm{,}n\mathrm{,}k}^s+T_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{r_i^2{\left(\Delta \right)}^2}\left(T_{i\mathrm{,}n+\mathrm{1,}k}^s-2T_{i\mathrm{,}n\mathrm{,}k}^s+T_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(T_{i\mathrm{,}n\mathrm{,}k+1}^s-2T_{i\mathrm{,}n\mathrm{,}k}^s+T_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]\end{array}$

$\begin{array}{l}+\frac{E_c}{R_e}\frac{{\theta }^{\eta \left(a-1\right)}}{1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s}\{2[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4r_i{\left(\Delta \right)}^2}{\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(u_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{2\Delta z}{\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]\\ +\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2+\frac{\Delta t}{2r_i\Delta r\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ -\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\end{array}$

$\begin{array}{l}-\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\\ +\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2+\frac{\Delta t}{2r_i\Delta \theta \Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4\Delta r\Delta z}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\\ +\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\}+E_c{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\mathrm{,}\end{array}$ (68)

Concentration equation in discrete form:

$\begin{array}{c}{C}_{i\mathrm{,}n\mathrm{,}k}^{s+1}=C_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left(C_{i+\mathrm{1,}n\mathrm{,}k}^s-C_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left(C_{i\mathrm{,}n+\mathrm{1,}k}^s-C_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left(C_{i\mathrm{,}n\mathrm{,}k+1}^s-C_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\left(\frac{1}{S_c{R}_e}+{\phi }_c\right)[\frac{\Delta t}{2r_i\Delta r}\left(C_{i+\mathrm{1,}n\mathrm{,}k}^s-C_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta r\right)}^2}\left(C_{i+\mathrm{1,}n\mathrm{,}k}^s-2C_{i\mathrm{,}n\mathrm{,}k}^s+C_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}\left(C_{i\mathrm{,}n+\mathrm{1,}k}^s-2C_{i\mathrm{,}n\mathrm{,}k}^s+C_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta z\right)}^2}\left(C_{i\mathrm{,}n\mathrm{,}k+1}^s-2C_{i\mathrm{,}n\mathrm{,}k}^s+C_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]-\Delta t{k}_c{R}_e{C}_{i\mathrm{,}n\mathrm{,}k}^s\mathrm{,}\end{array}$ (69)

Discrete $\varkappa$ -equation for the liquid phase:

$\begin{array}{c}{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^{s+1}={\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left({\varkappa }_{i\mathrm{,}n\mathrm{,}k+1}^s-{\varkappa }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}\frac{\Delta t}{2r_i^2\Delta \theta }\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+\frac{{\nu }_T}{{\sigma }_{\varkappa }}\right)[\frac{\Delta t}{{\left(\Delta r\right)}^2}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{2r_i\Delta r}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{{\left(\Delta z\right)}^2}\left({\varkappa }_{i\mathrm{,}n\mathrm{,}k+1}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]+2{\nu }_T[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\end{array}$

$\begin{array}{c}+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2+\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(u_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2\\ +\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]+{\nu }_T\{\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta r\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ +\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2-\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2\end{array}$

$\begin{array}{c}+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta \theta \Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\\ +\frac{\Delta t}{4\Delta r\Delta z}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\}\\ -{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\mathrm{,}\end{array}$ (70)

Discrete $\varkappa$ -equation for the gas phase:

$\begin{array}{c}{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^{s+1}={\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left({\varkappa }_{i\mathrm{,}n\mathrm{,}k+1}^s-{\varkappa }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}\\ \times \frac{\Delta t}{2r_i^2\Delta \theta }\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}+\frac{{\nu }_T}{{\sigma }_{\varkappa }}\right)\\ \times [\frac{\Delta t}{{\left(\Delta r\right)}^2}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{2r_i\Delta r}\left({\varkappa }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\varkappa }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{r_i^2{\left(\Delta \theta \right)}^2}\left({\varkappa }_{i\mathrm{,}n+\mathrm{1,}k}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left({\varkappa }_{i\mathrm{,}n\mathrm{,}k+1}^s-2{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s+{\varkappa }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]\end{array}$

$\begin{array}{c}+2{\nu }_T[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(u_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]\\ +{\nu }_T\{\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2+\frac{\Delta t}{2r_i\Delta r\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ -\frac{\Delta t}{r_i^2\Delta \theta }\left(U_{i\mathrm{,}n+\mathrm{1,}k}^s-U_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\end{array}$

$\begin{array}{c}-\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\\ +\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2+\frac{\Delta t}{2r_i\Delta \theta \Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4\Delta r\Delta z}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\\ +\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\}-{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\mathrm{,}\end{array}$ (71)

Discrete ε-equation for the liquid phase:

$\begin{array}{c}{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^{s+1}={\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^m\\ -\frac{\Delta t}{2\Delta z}\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k+1}^s-{\epsilon }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}j\mathrm{,}k}^m\right]}^{a-1}}\frac{\Delta t}{2r_i^2\Delta \theta }\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\\ +\left(\frac{{\theta }^{\eta \left(a-1\right)}}{R_e{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}+\frac{{\nu }_T}{{\sigma }_{\epsilon }}\right)[\frac{\Delta t}{2r_i\Delta r}\left({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)+\frac{\Delta t}{{\left(\Delta r\right)}^2}({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\\ +{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s)+\frac{\Delta t}{r_i^2{\left(\theta \right)}^2}\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s+{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k+1}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s+{\epsilon }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]\\ +2f_1{C}_{\epsilon \mathrm{,1}}\frac{{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}{\nu }_T[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\end{array}$

$\begin{array}{c}+\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(u_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]\\ +f_1{C}_{\epsilon \mathrm{,1}}\frac{{\epsilon }_{i\mathrm{,}j\mathrm{,}k}^m}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}{\nu }_T\{\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta r\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ +\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2-\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\end{array}$

$\begin{array}{c}+\frac{\Delta t}{r_i^2}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta \theta \Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\\ +\frac{\Delta t}{4\Delta r\Delta z}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\}\\ -\Delta t{f}_2{C}_{\epsilon \mathrm{,2}}\frac{{\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}\end{array}$ (72)

Discrete ε-equation for the gas phase:

$\begin{array}{c}{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^{s+1}={\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2\Delta r}\left({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s-\frac{\Delta t}{2r_i\Delta \theta }\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ -\frac{\Delta t}{2\Delta z}\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k+1}^s-{\epsilon }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)w_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\eta \left(a-1\right){\theta }^{\eta \left(a-2\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}\\ \times \frac{\Delta t}{2r_i^2\Delta \theta }\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\left(\frac{{\theta }^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{R_e}+\frac{{\nu }_T}{{\sigma }_{\varkappa }}\right)\\ \times [\frac{\Delta t}{2r_i\Delta r}({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s)+\frac{\Delta t}{{\left(\Delta r\right)}^2}\left({\epsilon }_{i+\mathrm{1,}n\mathrm{,}k}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s+{\epsilon }_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\\ +\frac{\Delta t}{r_i^2{\left(\theta \right)}^2}\left({\epsilon }_{i\mathrm{,}n+\mathrm{1,}k}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s+{\epsilon }_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{{\left(\Delta z\right)}^2}\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k+1}^s-2{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s+{\epsilon }_{i\mathrm{,}n\mathrm{,}k-1}^s\right)]\end{array}$

$\begin{array}{l}+2f_1{C}_{\epsilon \mathrm{,1}}\frac{{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}{\nu }_T[\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(u_{i+\mathrm{1,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{r_i^2\Delta \theta }\left(v_{i\mathrm{,}n+\mathrm{1,}k}^s-v_{i\mathrm{,}n-\mathrm{1,}k}^s\right)u_{i\mathrm{,}n\mathrm{,}k}^s+\frac{\Delta t}{r_i^2}{\left(u_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(w_{i\mathrm{,}n\mathrm{,}k+1}^s-w_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2]\\ +f_1{C}_{\epsilon \mathrm{,1}}\frac{{\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}{\nu }_T\{\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta r\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)-\frac{\Delta t}{r_i^2\Delta \theta }\left(u_{i\mathrm{,}n+\mathrm{1,}k}^s-u_{i\mathrm{,}n-\mathrm{1,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\\ +\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2-\frac{\Delta t}{2r_i\Delta r}\left(v_{i+\mathrm{1,}n\mathrm{,}k}^s-v_{i-\mathrm{1,}n\mathrm{,}k}^s\right)v_{i\mathrm{,}n\mathrm{,}k}^s\end{array}$

$\begin{array}{l}+\frac{\Delta t}{r_i^2}{\left(v_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2+\frac{\Delta t}{4r_i^2{\left(\Delta \theta \right)}^2}{\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)}^2\\ +\frac{\Delta t}{2r_i\Delta \theta \Delta z}\left(v_{i\mathrm{,}n\mathrm{,}k+1}^s-v_{i\mathrm{,}n\mathrm{,}k-1}^s\right)\left(w_{i\mathrm{,}n+\mathrm{1,}k}^s-w_{i\mathrm{,}n-\mathrm{1,}k}^s\right)+\frac{\Delta t}{4{\left(\Delta r\right)}^2}{\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)}^2\\ +\frac{\Delta t}{4\Delta r\Delta z}\left(w_{i+\mathrm{1,}n\mathrm{,}k}^s-w_{i-\mathrm{1,}n\mathrm{,}k}^s\right)\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)+\frac{\Delta t}{4{\left(\Delta z\right)}^2}{\left(u_{i\mathrm{,}n\mathrm{,}k+1}^s-u_{i\mathrm{,}n\mathrm{,}k-1}^s\right)}^2\}\\ -\Delta t{f}_2{C}_{\epsilon \mathrm{,2}}\frac{{\left({\epsilon }_{i\mathrm{,}n\mathrm{,}k}^s\right)}^2}{{\varkappa }_{i\mathrm{,}n\mathrm{,}k}^s}\mathrm{.}\end{array}$ (73)

Space discretization for the local Nusselt number, Sherwood number, and coefficient of skin friction as provided in Equations (54)-(57) are performed based on the Finite Difference backward approximation and the results evaluated at $r=0.5$ as follows:

$C_{f_{\left(r\right)}}=-\frac{2\theta {\left(n\right)}^{\eta \left(a-1\right)}{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}{Re}\left[\frac{u_{i\mathrm{,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s}{\Delta r}\right]\mathrm{,}$ (74)

$C_{f_{\left(r\right)}}=-\frac{2\theta {\left(j\right)}^{\eta \left(a-1\right)}}{Re{\left[1+\lambda T_{i\mathrm{,}n\mathrm{,}k}^s\right]}^{a-1}}\left[\frac{u_{i\mathrm{,}n\mathrm{,}k}^s-u_{i-\mathrm{1,}n\mathrm{,}k}^s}{\Delta r}\right]\mathrm{,}$ (75)

$N_u=-\frac{T_{i\mathrm{,}n\mathrm{,}k}^s-T_{i-\mathrm{1,}n\mathrm{,}k}^s}{\Delta r}\mathrm{,}$ (76)

$Sh=-\frac{C_{i\mathrm{,}n\mathrm{,}k}^s-C_{i-\mathrm{1,}n\mathrm{,}k}^s}{\Delta r}\mathrm{.}$ (77)

5. Results and Discussions

This section explains how the different flow parameters affect concentration, skin friction, Nusselt number, and Sherwood number physically. The Prandtl number is chosen to be 0.71 for air and 6.9 for water.

5.1. Flow Parameter Effects on Concentration Distribution

For both phases, Figure 2 demonstrates that the concentration distribution increases as the Reynolds number increases. Since both the fluid velocity and temperature increase as Re increases, the molecules will diffuse at an extremely rapid rate. The concentration boundary layer increases as species diffusivity increases.

Figure 3 demonstrates that for both the gaseous phase and the liquid phase, concentration profiles increase as α, the angle of inclination, increases. The concentration boundary layer thickness grows as the angle of inclination increases, increasing concentration profiles.

The influence of the solutal Grasholf number on concentration is depicted in Figure 4 for the gas phase and the liquid phase, respectively. It is important to note that as the solutal Grasholf number increases, so does the concentration distribution. The concentration boundary layer thickness increases as the solutal Grasholf number rises.

Figure 2. Reynolds number’s effect on the concentration distribution.

Figure 3. Effect of inclination angle on concentration distribution.

Figure 4. Solutal Grasholf number’s effect on concentration distribution.

5.2. Effect of Flow Parameters on Heat, Mass, and Skin Friction Transfers

From Table 1, it is crucial to remember that a decrease in skin friction for both phases is caused by an increase in the Reynolds number, thermal Grasholf number,

Table 1. Local skin friction values for different physical parameters.

and singular Grasholf number. Increasing the Reynolds number causes the viscous forces to decrease. Because skin friction is inversely related to viscous forces, it becomes less significant as skin friction increases. An increase in inclination angle leads to an increase in skin friction coefficient.

The data in Table 2 were obtained by numerical simulations of various physical parameter values for the local Nusselt number. Increase in the Reynolds number (Re), Eckert number (Ec), Prandtl number (P_r), and angle of inclination (α) results in an increase in the rate of heat transfer. As the Reynolds number rises, the flow velocity rises as well, which in turn causes the temperature gradient to rise and the Nusselt number to rise. The thickness of the thermal boundary layer decreases as the angle of inclination, Prandtl number, Eckert number, and Reynolds number rise. The rate of heat transfer is enhanced by reducing the thermal boundary layer. As the power law index rises, the rate of heat transfer reduces. The thickness of the thermal boundary layer increases as the power law index rises. The rate of heat transfer is unaffected by the chemical reaction parameter.

For various values of the physical parameters, numerical values of the local Sherwood numbers for the gas and liquid phase are shown in Table 3. The rate of mass transfer increases with an increase in the chemical reaction parameter, Reynolds number, and angle of inclination but decreases with an increase in

Table 2. Nusselt number numerical values for physical parameters.

Table 3. Sherwood number numerical values for values of physical parameters.

Schmidt number. An increase in Reynolds number $R_e,k_c$ and α implies there is more interaction of the species concentration with the momentum boundary layer. Larger values of the Schmidt Number result in a slower rate of mass transfer. The rate of mass transfer is unaffected by the power law index, n.

6. Conclusion

This study examined the impact of various flow parameters on two-phase turbulent fluid flow in a geothermal pipe with a chemical reaction. The numerical findings for the Sherwood number, Nusselt number, and coefficient of skin friction have been tabulated. From the findings, the following significant points may be made:

· Concentration distribution could be enhanced with increase of Solutal Grasholf number, angle of inclination, and Reynolds number.

· An increase in inclination angle increases the coefficient of skin friction.

· Higher values of Reynolds number and solutal Grasholf number result in a reduction of skin friction. Hence, there is enhanced efficiency of the system.

· Heat transfer rate increases by increasing Eckert number, Prandtl number, Reynolds number, and angle of inclination.

· The rate of mass transfer is unaffected by the power law index.

· When the Reynolds number, Chemical reaction parameter, and angle of inclination rise, the amount of mass transfer also rises but falls as Schmidt number rises.

Acknowledgements

The author(s) extend their gratitude to the Pan African University Institute for Basic Sciences, Technology and Innovation for funding this research.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Palsson, H., Bergthòrsson, E.S. and Palsson, O.P. (2006) Estimation and Validation of Models Two Phase Flow from Geothermal Wells. 10th International Symposium on District Heating and Cooling, Reykjavik, 3-5 September 2016.
[2]	Li, Z., Wang, G., Yousaf, M., Yang, X. and Ishii, M. (2018) Flow Structure and Flow Regime Transitions of Downward Two-Phase Flow in Large Diameter Pipes. International Journal of Heat and Mass Transfer, 118, 812-822. https://doi.org/10.1016/j.ijheatmasstransfer.2017.11.037
[3]	Eghbali, S., Banks, J. and Nobes, D.S. (2021) A Numerical Study on Compositional Modeling of Two-Phase Fluid Flow and Heat Transfer in Vertical Wells. Journal of Petroleum Science and Engineering, 201, Article ID: 108400. https://doi.org/10.1016/j.petrol.2021.108400
[4]	Hasan, A.R. and Kabir, C.S. (2010) Modeling Two-Phase Fluid and Heat Flows in Geothermal Wells. Journal of Petroleum Science and Engineering, 71, 77-86. https://doi.org/10.1016/j.petrol.2010.01.008
[5]	Nizami, M., et al. (2016) Mathematical Modelling of Silica Scaling Deposition in Geothermal Wells. IOP Conference Series: Earth and Environmental Science, 42, Article ID: 012013. https://doi.org/10.1088/1755-1315/42/1/012013
[6]	Ojiambo, V., Kinyanjui, M. and Kimathi, M. (2018) A Mathematical Model of Angular Two-Phase Jeffery Hamel Flow in a Geothermal Pipe. International Journal of Advances in Applied Mathematics and Mechanics, 6, 1-13.
[7]	Chauhan, V., Saevarsdottir, G., Tesfahunegn, Y.A., Asbjornsson, E. and Gudjonsdottir, M. (2021) Computational Study of Two-Phase Flashing Flow in a Calcite Scaled Geothermal Wellbore. Geothermics, 97, Article ID: 102239. https://doi.org/10.1016/j.geothermics.2021.102239
[8]	Vahidinia, F. and Miri, M. (2015) Numerical Study of the Effect of the Reynolds Numbers on Thermal and Hydrodynamic Parameters of Turbulent Flow Mixed Convection Heat Transfer in an Inclined Tube. Strojniški vestnik-Journal of Mechanical Engineering, 61, 669-679. https://doi.org/10.5545/sv-jme.2015.2818
[9]	Launder, B.E., Spalding, D.B., et al. (1972) Lectures in Mathematical Models of Turbulence.
[10]	Asl, K.D. and Jalali, A. (2006) Analysis of Heat Transfer in Turbulent Flow through the Tube with Uniform Heat Flux.
[11]	Mondal, H., Pal, D., Chatterjee, S. and Sibanda, P. (2018) Thermophoresis and Soret-Dufour on MHD Mixed Convection Mass Transfer over an Inclined Plate with Non-Uniform Heat Source/Sink and Chemical Reaction. Ain Shams Engineering Journal, 9, 2111-2121. https://doi.org/10.1016/j.asej.2016.10.015
[12]	Tarakaramu, N., Satya Narayana, P.V. and Venkateswarlu, B. (2020) Numerical Simulation of Variable Thermal Conductivity on 3d Flow of Nanofluid over a Stretching Sheet. Nonlinear Engineering, 9, 233-243. https://doi.org/10.1515/nleng-2020-0011
[13]	Duan, J.M., Liu, H.S., Wang, N., et al. (2015) Hydro Dynamic Modeling of Stratified Smooth Two-Phase Turbulent Flow with Curved Interface through Circular Pipe. International Journal of Heat and Mass Transfer, 89, 1034-1043. https://doi.org/10.1016/j.ijheatmasstransfer.2015.05.093
[14]	Abdulwahid, M.A., Kareem, H.J. and Almudhaffar, M.A. (2017) Numerical Analysis of Two Phase Flow Patterns in Vertical and Horizontal Pipes. WSEAS Transactions on Fluid Mechanics, 12, 131-140.
[15]	Menge, B.K. (2015) Analysis of Turbulent Flow in a Pipe at Constant Reynolds Number Using Computational Fluid Dynamics. Ph.D. Thesis, Jomo Kenyatta University of Agriculture and Technology, Nairobi.
[16]	Vitturi, M. (2016) Navier-Stokes Equations in Cylindrical Coordinates.
[17]	Nagler, J. (2017) Jeffery-Hamel Flow of Non-Newtonian Fluid with Nonlinear Viscosity and Wall Friction. Applied Mathematics and Mechanics, 38, 815-830. https://doi.org/10.1007/s10483-017-2206-8
[18]	Ojiambo, V., Kinyanjui, M. and Kimathi, M. (2018) A Study of Two-Phase Jeffery Hamel Flow in a Geothermal Pipe.
[19]	Hayat, T., Asad, S. and Alsaedi, A. (2014) Flow of Variable Thermal Conductivity Fluid Due to Inclined Stretching Cylinder with Viscous Dissipation and Thermal Radiation. Applied Mathematics and Mechanics, 35, 717-728. https://doi.org/10.1007/s10483-014-1824-6
[20]	Celik, I.B. (1999) Introductory Turbulence Modeling Lecture Notes. West Virginia University, Morgantown.
[21]	Canli, E., Ate, A. and Bilir, S. (2020) Derivation of Dimensionless Governing Equations for Axisymmetric Incompressible Turbulent Flow Heat Transfer Based on Standard κ-ω Model. Afyon Kocatepe University Journal of Sciences and Engineering, 20, 1096-1111. https://doi.org/10.35414/akufemubid.821009
[22]	Kumar, A. (2016) Derivation of Turbulence Kinetic Energy.
[23]	Lam, C.K.G. and Bremhorst, K. (1981) A Modified Form of the κ-ω Model for Predicting Wall Turbulence. Journal of Fluids Engineering, 103, 456-460. https://doi.org/10.1115/1.3240815
[24]	Versteeg, H.K. and Malalasekera, W. (2007) An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education, London.