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

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.

Share and Cite:

Nyariki, E. , Kinyanjui, M. and Abonyo, J. (2023) Heat and Mass Transfers in a Two-Phase Stratified Turbulent Fluid Flow in a Geothermal Pipe with Chemical Reaction. Journal of Applied Mathematics and Physics, 11, 484-513. doi: 10.4236/jamp.2023.112030.

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 , θ , 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:

( w ) z + 1 r ( v ) θ + 1 r ( r u ) 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:

ρ ( u t + w u z + v r u θ + u u r v 2 r ) = μ [ u r 2 + 1 r r ( r u r ) ] p r + μ [ 2 u z 2 + 1 r 2 2 u θ 2 2 r 2 v θ ] + μ θ ( 1 r 2 u θ + 1 r v r v r 2 ) + ρ g cos α [ β T * ( T T ) + β c * ( C C ) ] , (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 β T * and β C * , respectively.

Equation of momentum in the θ-direction:

ρ [ v t + w v z + v r v θ + u v r + u v r ] = μ θ [ 2 u r 2 + 2 r 2 v θ ] + μ [ 1 r 2 2 v θ 2 v r 2 + 2 v z 2 + 1 r r ( r v r ) + 2 r 2 u θ ] , (3)

Equation of momentum equation in the z-direction:

ρ [ w t + w w z + v r w θ + u w r ] = μ θ [ 1 r v z + 1 r 2 w θ ] + μ [ 2 w z 2 + 1 r r ( r w r ) + 1 r 2 2 w θ 2 ] + ρ g sin α [ β T * ( T T ) + β c * ( C C ) ] , (4)

Equation for energy:

ρ C p [ v r T θ + w T z + u T r + T t ] = κ T ( T r ) 2 + 1 r 2 κ T ( T θ ) 2 + κ T ( T z ) 2 + κ ( 1 r T r + 2 T r 2 + 1 r 2 2 T θ 2 + 2 T z 2 ) + Φ , (5)

where κ is the fluid’s thermal conductivity and C p is the specific heat at constant pressure.

The equation for concentration:

u C r + w C z + v r C θ + C t = D c [ 1 r 2 2 C θ 2 + 1 r r ( r C r ) + 2 C z 2 ] k r ( C C ) , (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;

μ = μ 0 f a 1 , f = f ( θ ) , (7)

where μ 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 ( θ , T ) = { θ η [ 1 + λ ( T T ) ] , for λ < 0 , η 0 ; θ η 1 + λ ( T T ) , for λ > 0 , η 0. (8)

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

According to [19] , thermal conductivity κ is temperature-dependent and is described as

κ = κ [ 1 + ω ( T T ) Δ T ] , (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 κ .

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:

u u ¯ = ν T ( ( u ) + ( u ) T ) 2 3 k δ i j , (10)

δ i j is the Kronecker delta

u T ¯ = φ T T ¯ r , v T ¯ = φ T r T ¯ θ , w T ¯ = φ T T ¯ z , u C ¯ = φ c C ¯ r , v C ¯ = φ c r C ¯ θ , w C ¯ = φ c C ¯ z . (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:

u ¯ u ¯ r + w ¯ u ¯ z + v ¯ r u ¯ θ + u ¯ t v ¯ v ¯ r = ( μ 0 f a 1 ρ + ν T ) [ 2 r 2 v ¯ θ + 1 r r ( r u ¯ r ) + 1 r 2 2 u ¯ θ 2 u ¯ r 2 + 2 u ¯ z 2 ] 1 ρ p r + μ 0 ρ η ( a 1 ) f a 2 θ a 1 [ 1 + λ ( T ¯ T ) ] [ 1 r 2 u ¯ θ v ¯ r 2 + 1 r v ¯ r ] + g cos α [ β T * ( T ¯ T ) + β c * ( C ¯ C ) ] , (12)

where u ¯ , v ¯ , w ¯ are the time-averaged velocities in the radial, tangential and axial directions respectively, T ¯ is the time average temperature, C ¯ is the time average concentration, p ¯ is the time average pressure and ν T is the eddy viscosity.

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

w ¯ u ¯ z + v ¯ r u ¯ θ + u ¯ u ¯ r v ¯ v ¯ r + u ¯ t = ( μ 0 f a 1 ρ + ν T ) [ u ¯ r 2 + 1 r r ( r u ¯ r ) + 1 r 2 2 u ¯ θ 2 2 r 2 v ¯ θ + 2 u ¯ z 2 ] 1 ρ p r + μ 0 ρ η ( a 1 ) f a 2 θ η 1 1 + λ ( T ¯ T ) [ 1 r 2 u ¯ θ + 1 r v ¯ r v ¯ r 2 ] + g cos α [ β T * ( T ¯ T ) + β C * ( C ¯ C ) ] , (13)

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

u ¯ v ¯ r + u ¯ v ¯ r + w ¯ v ¯ z + v ¯ r v ¯ θ + v ¯ t = μ 0 ρ η ( a 1 ) f a 2 θ η 1 1 + λ ( T T ) [ 2 u ¯ r 2 + 2 r 2 v ¯ θ ] ( μ 0 f a 1 ρ + ν T ) × [ v ¯ r 2 + 1 r r ( r v ¯ r ) + 2 r 2 u ¯ θ + 2 v ¯ z 2 + 1 r 2 2 v ¯ θ 2 ] 1 ρ r p θ , (14)

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

u ¯ v ¯ r + u ¯ v ¯ r + w ¯ v ¯ z + v ¯ r v ¯ θ + v ¯ t = μ 0 ρ η ( a 1 ) f a 2 θ η 1 [ 1 + λ ( T T ) ] [ 2 r 2 v ¯ θ + 2 u ¯ r 2 ] 1 ρ r p θ + ( μ 0 f a 1 ρ + ν T ) [ 2 r 2 u ¯ θ v ¯ r 2 + 1 r r ( r v ¯ r ) + 2 v ¯ z 2 + 1 r 2 2 v ¯ θ 2 ] , (15)

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

w ¯ w ¯ z + v ¯ r w ¯ θ + u ¯ w ¯ r + w ¯ t = η ( a 1 ) μ 0 f a 2 θ η 1 ρ r 2 [ 1 + λ ( T ¯ T ) ] ( r v ¯ z + w ¯ θ ) 1 ρ p z

+ ( μ 0 f a 1 ρ + ν T ) [ 1 r 2 w ¯ θ 2 + 1 r ( r w ¯ r ) + 2 w ¯ z 2 ] + g sin α [ β T * ( T ¯ T ) + β c * ( C ¯ C ) ] , (16)

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

w ¯ t + w ¯ w ¯ z + v ¯ r w ¯ θ + u ¯ w ¯ r = 1 ρ p z + η ( a 1 ) μ 0 f a 2 θ η 1 ρ r 2 [ 1 + λ ( T ¯ T ) ] ( w ¯ θ + r v ¯ z ) + ( μ 0 f a 1 ρ + ν T ) [ 1 r ( r w ¯ r ) + 2 w ¯ z 2 + 1 r 2 w ¯ θ 2 ] + g sin α [ β T * ( T ¯ T ) + β c * ( C ¯ C ) ] (17)

Averaged energy equation:

T ¯ t + w ¯ T ¯ z + v ¯ r T ¯ θ + u ¯ T ¯ r = κ ρ C p ω Δ T [ ( 1 r 2 ( T ¯ θ ) 2 + T ¯ r ) 2 + ( T ¯ z ) 2 ] + ( κ ρ C p [ 1 + ω T ¯ T Δ T ] + φ T ) ( 1 r T ¯ r + 2 T ¯ r 2 + 1 r 2 2 T ¯ θ 2 + 2 T ¯ z 2 ) + Φ ρ C p + ε C p , (18)

Φ = 2 μ [ ( u ¯ r ) 2 + ( w ¯ z ) 2 + ( u ¯ r + 1 r v ¯ θ ) 2 ] + μ [ ( 1 r u ¯ θ v ¯ r + v ¯ r ) 2 + ( 1 r w ¯ θ + v ¯ z ) 2 + ( u ¯ z + w ¯ r ) 2 ] , (19)

where φ T is turbulent thermal diffusivity and Φ is the viscous dissipation.

Averaged concentration equation:

C ¯ t + w ¯ C ¯ z + v ¯ r C ¯ θ + u ¯ C ¯ r = ( D c + φ c ) [ 1 r 2 2 C ¯ θ 2 + 1 r r ( r C ¯ r ) + 2 C ¯ z 2 ] k r ( C ¯ C ) , (20)

where k r is a chemical reaction parameter, D c is a diffusion coefficient, and φ c is the turbulent mass diffusivity.

The ϰ - ε turbulence model is used to calculate the turbulent viscosity. Given in (21) is the relationship between turbulent viscosity ( ν T ), turbulent kinetic energy ( ϰ ), and dissipation (ε).

ν T = C μ f μ ϰ 2 ε , (21)

where C μ = 0.09 and f μ are empirical constant and turbulent model functions respectively.

Derivation of ϰ 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 ϰ equation is formed by this summation. The ϰ -equation according to [21] and [22] is given as a vector in (22):

ρ ϰ t + ρ [ u ¯ ϰ ] = [ ( μ + μ T σ ϰ ) ϰ ] + μ T ϕ ρ ε , (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. σ ϰ is the turbulent Prandtl number for ϰ -equation.

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

ρ [ u ¯ ε ] + ρ ε t = [ ( μ + μ T σ ε ) ε ] + f 1 C ε ,1 ε ϰ ϕ f 2 ρ C ε ,2 ε 2 ϰ , (23)

where the wall model functions f 2 and f 1 are utilized to adjust the ϰ - ε 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 ϰ -equation for gas phase:

ϰ t + u ¯ ϰ r + v ¯ r ϰ θ + w ¯ ϰ z = 1 ρ r 2 μ 0 η ( a 1 ) f a 2 θ η 1 [ 1 + λ ( T ¯ T ) ] ϰ θ + ( μ 0 f a 1 ρ + ν T σ ϰ ) [ 1 r r ( ϰ r ) + 1 r 2 2 ϰ θ 2 + 2 ϰ z 2 ] + 2 ν T [ ( u ¯ r + 1 r v ¯ θ ) 2 + ( u ¯ r ) 2 + ( w ¯ z ) 2 ] + ν T [ ( u ¯ z + w ¯ r ) 2 + ( 1 r u ¯ θ v ¯ r + v ¯ r ) 2 + ( 1 r w ¯ θ + v ¯ z ) 2 ] ε , (24)

The ϰ -equation for liquid phase:

ϰ t + u ¯ ϰ r + v ¯ r ϰ θ + w ¯ ϰ z = 1 r 2 μ 0 ρ η ( a 1 ) f a 2 θ η 1 1 + λ ( T ¯ T ) ϰ θ + ( μ 0 f a 1 ρ + ν T σ ϰ ) [ 1 r r ( r ϰ r ) + 2 ϰ z 2 + 1 r 2 2 ϰ θ 2 ] + 2 ν T [ ( 1 r v ¯ θ + u ¯ r ) 2 + ( u ¯ r ) 2 + ( w ¯ z ) 2 ] + ν T [ ( u ¯ z + w ¯ r ) 2 + ( 1 r u ¯ θ + v ¯ r v ¯ r ) 2 ] + ν T [ ( 1 r w ¯ θ + v ¯ z ) 2 ] ε , (25)

The ε-equation for gas phase:

ε t + u ¯ ε r + w ¯ ε z + v ¯ r ε θ = μ 0 η ( a 1 ) f a 2 θ η 1 r 2 ρ [ 1 + λ ( T ¯ T ) ] ε θ + ( μ 0 f a 1 ρ + ν T σ ε ) [ 1 r r ( r ε r ) + 2 ε z 2 + 1 r 2 2 ε θ 2 ] + 2 f 1 c ε ,1 ε ϰ ν T [ ( u ¯ r ) 2 + ( 1 r v ¯ θ + u ¯ r ) 2 + ( u ¯ z ) 2 ] + f 1 c ε ,1 ε ϰ ν T [ ( 1 r u ¯ θ + v ¯ r v ¯ r ) 2 + ( v ¯ z + 1 r w ¯ θ ) 2 + ( w ¯ r + u ¯ z ) 2 ] f 2 c ε ,2 ε 2 ϰ , (26)

The liquid phase ε-equation:

ε t + w ¯ ε z + v ¯ r ε θ + u ¯ ε r = μ 0 η ( a 1 ) f a 2 θ η 1 r 2 ρ [ 1 + λ ( T ¯ T ) ] ε θ + ( μ 0 f a 1 ρ + ν T σ ε ) [ 1 r r ( r ε r ) + 2 ε z 2 + 1 r 2 2 ε θ 2 ] + 2 f 1 c ε ,1 ε ϰ ν T [ ( u ¯ r ) 2 + ( w ¯ z ) 2 + ( 1 r v ¯ θ + u ¯ r ) 2 ] + f 1 c ε ,1 ε ϰ ν T [ ( v ¯ r + 1 r u ¯ θ + v ¯ r ) 2 + ( 1 r w ¯ θ + v ¯ z ) 2 + ( u ¯ z + w ¯ r ) 2 ] f 2 c ε ,2 ε 2 ϰ , (27)

where σ ϰ , c ε ,1 , σ ε , c ε ,2 and c μ are [23] model constants, which are listed as c μ = 0.09 , σ ϰ = 1 , c ε ,1 = 1.44 , σ ε = 1.3 c ε ,2 = 1.92 .

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

f μ = [ 1 exp ( 0.0165 R ϰ ) ] 2 ( 1 + 20.5 R t ) , (28)

f 1 = ( 1 + 0.05 f μ ) 3 , (29)

f 2 = 1 exp ( R t 2 ) , (30)

where R t and R ϰ 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 = D + 2 :

u ¯ = w ¯ = v ¯ = 0 , T ¯ = T w , ϰ = 0 , C ¯ = C w , ε = 0

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

C ¯ = C , u ¯ = U , T ¯ = T , v ¯ = U , ϰ r = 0 , w ¯ = U , ε r = 0.

r = D 2 at the wall of the pipe

u ¯ = w ¯ = 0 , ϰ = 0 , C ¯ = C w , ε = 0 , v ¯ = 0 , T ¯ = T w

z = 0 at the entrance

w ¯ = U , C ¯ = C , v ¯ = 0 , u ¯ = 0 , T ¯ = T ,

ϰ = 3 2 ( U I ) 2 , I = 0.16 R e 1 8 , ε = C μ 3 4 ϰ 3 2 0.07 D ,

where I is the turbulent intensity.

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

u ¯ z = 0 , v ¯ z = 0 , w ¯ z = 0 , T ¯ z = 0 , C ¯ z = 0 , ϰ z = 0 , ε 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 ( r ) = 2 τ w ρ U 2 , (31)

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

τ w = μ ( u ¯ r ) r = D 2 . (32)

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

C f ( r ) = 2 μ ρ U 2 ( u ¯ r ) r = D 2 . (33)

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

C f ( r ) = 2 μ 0 ρ U 2 θ η ( a 1 ) [ 1 + λ ( T ¯ T ) ] a 1 ( u ¯ r ) r = D 2 , (34)

C f ( r ) = 2 μ 0 ρ U 2 θ η ( a 1 ) [ 1 + λ ( T ¯ T ) ] a 1 ( u ¯ r ) r = 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 = D q w κ ( T w T ) , (36)

where q w is heat flux given as follows;

q w = κ ( T ¯ r ) r = D 2 . (37)

N u = D T w T ( T ¯ r ) r = D 2 . (38)

Sherwood number is defined as:

S h = D C w C ( C ¯ r ) r = 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:

u = u ¯ U , v = v ¯ U , w = w ¯ U , t = t U D , p = p ¯ ρ U 2 , θ = θ , z = z D , ν T = ν T D U , r = r D , T = T ¯ T T w T , φ T = φ T D U , ε = ε D U 3 , ϰ = ϰ U 2 , φ c = φ c D U , C = C ¯ C C w C . (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:

u t + w u z + v r u θ + u u r v 2 r = ( θ η ( a 1 ) [ 1 + λ T ] a 1 R e + ν T ) [ u r 2 + 1 r r ( r u r ) + 2 u z 2 + 1 r 2 2 u θ 2 2 r 2 v θ ] p r + η ( a 1 ) θ η ( a 2 ) [ 1 + λ T ] a 1 R e × [ v r 2 + 1 r 2 u θ + 1 r v r ] + cos α ( T G r ( T ) R e 2 + C G r ( c ) R e 2 ) , (41)

where R e = U D ν is Reynolds Number, G r ( T ) = g β T ( T w T ) ν 2 is Thermal Grasholf Number and G r ( c ) = g β c ( C w C ) ν 2 is solutal Grasholf Number.

Dimensionless momentum equation for the liquid phase in radial direction:

u t + w u z + v r u θ + u u r v 2 r = ( θ η ( a 1 ) R e [ 1 + λ T ] a 1 + ν T ) [ 2 r 2 v θ + 1 r r ( r u r ) + 2 u z 2 + 1 r 2 2 u θ 2 u r 2 ] p r + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T ] a 1 × [ v r 2 + 1 r 2 u θ + 1 r v r ] + cos α ( T G r ( T ) R e 2 + C G r ( C ) R e 2 ) , (42)

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

v t + w v z + v r v θ + u v r + v u r = 2 η ( n 1 ) θ η ( a 2 ) [ 1 + λ T ] a 1 R e [ 1 r 2 v θ + u r 2 ] 1 r p θ + ( θ η ( a 1 ) [ 1 + λ T ] a 1 R e + ν T ) [ 2 v z 2 + 2 r 2 u θ + 1 r r ( r v r ) + 1 r 2 2 v θ 2 v r 2 ] , (43)

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

v t + w v z + v r v θ + u v r + v u r = η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T ] a 1 [ 2 u r 2 + 2 r 2 v θ ] 1 r p θ + ( θ η ( a 1 ) R e [ 1 + λ T ] a 1 + ν T ) × [ 1 r r ( r v r ) + 2 v z 2 + 1 r 2 2 v θ 2 v r 2 + 2 r 2 u θ ] , (44)

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

w t + w w z + v r w θ + u w r = η ( a 1 ) θ η ( a 2 ) r 2 R e [ 1 + λ T ] a 1 ( w θ + r v z ) p z + ( θ η ( a 1 ) R e [ 1 + λ T ] a 1 + ν T ) × [ 1 r r ( r w r ) + 2 w z 2 + 1 r 2 w θ 2 ] + sin α [ T G r ( T ) R e 2 + C G r ( C ) R e 2 ] , (45)

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

w t + w w z + v r w θ + u w r = η ( a 1 ) θ η ( a 2 ) [ 1 + λ T ] a 1 R e 1 r 2 ( r v z + w θ ) p z + ( θ η ( a 1 ) [ 1 + λ T ] a 1 R e + ν T * ) [ 1 r w r + 2 w r 2 + 2 w z 2 + 1 r 2 w θ 2 ] + sin α [ T G r ( T ) R e 2 + C G r ( C ) R e 2 ] , (46)

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

T t + w T z + v r T θ + u T r = ω P r R e [ ( T r ) 2 + 1 r 2 ( T θ ) 2 + ( T z ) 2 ] + ( 1 + ω T P r R e + φ T ) × ( 2 T z 2 + 1 r T r + 1 r 2 2 T θ 2 + 2 T r 2 ) + E c R e θ η ( a 1 ) [ 1 + λ T ] a 1 × { 2 [ ( u r ) 2 + ( w z ) 2 + ( u r + 1 r v θ ) 2 ] + ( v r + 1 r u θ v r ) 2 } + E c R e θ η ( a 1 ) [ 1 + λ T ] a 1 [ ( v z + 1 r w θ ) 2 + ( w r + u z ) 2 ] + E c ε , (47)

where E c = U 2 C p ( T w T ) is Eckert Number and P r = μ C p κ * is Prandtl Number.

Dimensionless form of the gas phase energy equation:

T t + w T z + v r T θ + u T r = ω P r R e [ ( T r ) 2 + ( T z ) 2 + 1 r 2 ( T θ ) 2 ] + ( 1 + ω T P r R e + φ T )

× ( 1 r T r + 1 r 2 2 T θ 2 + 2 T r 2 + 2 T z 2 ) + E c θ η ( a 1 ) [ 1 + λ T ] a 1 R e × { 2 [ ( w z ) 2 + ( u r ) 2 + ( u r + 1 r v θ ) 2 ] + ( v r + 1 r u θ * + v r ) 2 } + E c θ η ( a 1 ) [ 1 + λ T ] a 1 R e [ ( 1 r w θ + v z ) 2 + ( u z + w r ) 2 ] + E c ε , (48)

Dimensionless concentration equation:

C t + w C z + v r C θ + u C r = ( 1 S c R e + φ c ) [ 1 r 2 2 C θ 2 + 1 r C r + 2 C r 2 + 2 C z 2 ] k c R e C (49)

where S c = ν D c is the Schmidt Number and k c = k r ν U 2 is the Chemical Reaction Parameter.

Dimensionless ϰ -equation for liquid phase:

ϰ t + u ϰ r + v r ϰ θ + w ϰ z = η ( a 1 ) θ η ( a 2 ) r 2 R e [ 1 + λ T ] a 1 ϰ θ + ( θ η ( a 1 ) R e [ 1 + λ T ] a 1 + ν T σ ϰ ) [ 1 r ϰ r + 2 ϰ r 2 + 2 ϰ z 2 + 1 r 2 2 ϰ θ 2 ] + 2 ν T [ ( u r ) 2 + ( w z ) 2 + ( u r + 1 r v θ ) 2 ] + ν T [ ( 1 r u θ v r + v r ) 2 + ( 1 r w θ + v z ) 2 + ( u z + w r ) 2 ] ε , (50)

Dimensionless ϰ -equation for gas phase:

ϰ t + u ϰ r + v r ϰ θ + w ϰ z = η ( a 1 ) θ η ( a 2 ) [ 1 + λ T ] a 1 r 2 R e ϰ θ + ( θ η ( a 1 ) [ 1 + λ T ] a 1 R e + ν T σ ϰ ) × ( 1 r ϰ r + 2 ϰ * z * 2 + 2 ϰ r 2 + 1 r 2 2 ϰ θ 2 ) + 2 ν T [ ( u r + 1 r v θ ) 2 + ( u r ) 2 + ( w z ) 2 ] + ν T [ ( 1 r u θ v r + v r ) 2 + ( u z + w r ) 2 + ( v z + 1 r w θ ) 2 ] ε (51)

Dimensionless ε-equation for liquid phase:

ε t + w ε z + v r ε θ + u ε r = η ( a 1 ) θ η ( a 2 ) r 2 R e [ 1 + λ T ] a 1 ε θ + ( θ η ( a 1 ) R e [ 1 + λ T ] a 1 + ν T σ ε ) [ 1 r * ε r * + 2 ε r 2 + 2 ε z 2 + 1 r 2 2 ε θ 2 ] + 2 f 1 c ε ,1 ε ϰ ν T [ ( u r ) 2 + ( 1 r v θ + u r ) 2 + ( w z ) 2 ] + f 1 c ε ,1 ε ϰ ν T [ ( 1 r u θ v r + v r ) 2 + ( u z + w r ) 2 + ( 1 r w θ + v z ) 2 ] f 2 c ε ,2 ε 2 ϰ . (52)

ν T = c μ f μ ϰ 2 ε (53)

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

C f ( r ) = 2 θ η ( a 1 ) [ 1 + λ T ] a 1 R e ( u r ) r = 0.5 , (54)

C f ( r ) = 2 θ η ( a 1 ) R e [ 1 + λ T ] a 1 ( u r ) r = 0.5 , (55)

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

N u = ( T r ) r = 0.5 . (56)

S h = ( T r ) r = 0.5 (57)

Dimensionless Boundary Conditions

Pipe wall, r = 0.5 +

v = 0 , w = 0 , u = 0 , ϰ = 0 , C = 1 , ε = 0 , T = 1

The fluid interface, r = 0

u = 1 , T = 0 , v = 1 , C = 0 , w = 1 , ϰ r = 0 , ε r = 0

r = 0.5 , wall of the pipe

u = 0 , w = 0 , v = 0 , ϰ = 0 , ε = 0 , T = 1 , C = 1

Inlet, z = 0

w = 1, v = 0, u = 0, T = 0, C = 0, ϰ = 3 2 I 2 , I = 0.016 R e 1 8 , ε = C μ 3 4 ϰ 3 2 0.07 .

Outlet z =

u z = 0 , C z = 0 , w z = 0 , ϰ z = 0 , v z = 0 , T z = 0 , ε z = 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 ( r , θ , z , t ) is described at the point ( i , n , k , s ) and u ( r , θ , z , t ) as u i , n , 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 Δ r , z k = z 0 + k Δ z , θ n = θ 0 + n Δ θ and t s = t 0 + s Δ t where s = 1 , 2 , 3 , 4 , , S , i = 1 , 2 , 3 , 4 , , I and k = 1,2,3,4 , , Z , n = 1 , 2 , 3 , 4 , , N .

The forward finite difference scheme with respect to time is

u t = u i , n , k s + 1 u i , n , k s Δ t (58)

The central finite difference approximation scheme with respect to space

u r = u i + 1 , n , k s u i 1 , n , k s 2 Δ r , (59)

2 u r 2 = u i + 1 , n , k s 2 u i , n , k m + U i 1 , n , k s ( Δ r ) 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:

u i , n , k s + 1 u i , n , k s Δ t + u i , n , k s u i + 1, n , k s u i 1, j , k s 2 Δ r + v i , n , k s r i u i , n + 1, k s u i , n 1, k s 2 Δ θ + w i , n , k s u i , n , k + 1 s u i , n , k 1 s 2 Δ z ( v i , n , k s ) 2 r i = p r + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T ) ( 1 r i u i + 1, n , k s u i 1, n , k s 2 Δ r + u i + 1, n , k s 2 u i , n , k s + u i 1, n , k s ( Δ r ) 2 u i , n , k s r i 2 + u i , n + 1, k s 2 u i , n , k s + u i , n 1, k s ( Δ θ ) 2 ) + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T ) ( u i , n , k + 1 s 2 u i , n , k s + u i , n , k 1 s ( Δ z ) 2 2 r i 2 v i , n + 1, k s v i , j 1, k m 2 Δ θ )

+ 1 R e η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , n , k s ] a 1 ( 1 r i 2 u i , n + 1, k s u i , n 1, k s 2 Δ θ v i , n , k s r i 2 + 1 r i v i + 1, n , k s v i 1, n , k s 2 Δ r ) + cos α ( T i , n , k s G r ( T ) R e 2 + C i , n , k s G r ( c ) R e 2 ) , (61)

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

u i , n , k s + 1 = u i , n , k s Δ t 2 Δ r ( u i + 1, n , k s u i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( u i , n , k + 1 s u i , n , k 1 s ) w i , n , k s + Δ t r i ( v i , n , k s ) 2 Δ t p r + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T ) [ Δ t 2 r i Δ r ( u i + 1, n , k s u i 1, n , k s )

+ Δ t ( Δ r ) 2 ( u i + 1, n , k s 2 u i , n , k s + u i 1, n , k s ) Δ t r i 2 u i , n , k s + Δ t ( Δ θ ) 2 ( u i , n + 1, k s 2 u i , n , k s + u i , n 1, k s ) + Δ t ( Δ z ) 2 ( u i , n , k + 1 s 2 u i , n , k s + u i , n , k 1 s ) Δ t r i 2 Δ θ ( v i , n + 1, k s V i , n 1, k s ) ] + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T i , n , k s ] a 1 ( Δ t 2 r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) Δ t r i 2 V i , n , k s + Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) ) + Δ t cos α ( T i , n , k s G r ( T ) R e 2 + C i , n , k s G r ( c ) R e 2 ) , (62)

u i + 1, n , k s , u i , n , k s and u i 1, n , k s hence u i , n , k s + 1 can be computed explicitly.

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

u i , n , k s + 1 = u i , n , k s Δ t 2 Δ r ( u i + 1, n , k s u i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( u i , j , k + 1 m u i , j , k 1 m ) w i , j , k m + Δ t r i ( v i , j , k m ) 2 Δ t p r + ( θ η ( a 1 ) [ 1 + λ T i , j , k m ] a 1 R e + ν T ) [ Δ t 2 r i Δ r ( u i + 1, n , k s u i 1, n , k s )

+ Δ t ( Δ r ) 2 ( u i + 1, n , k s 2 u i , n , k s + u i 1, n , k s ) Δ t r i 2 u i , n , k s + Δ t ( Δ θ ) 2 ( u i , n + 1, k s 2 u i , n , k s + u i , n 1, k s ) + Δ t ( Δ z ) 2 ( u i , n , k + 1 s 2 u i , n , k s + u i , n , k 1 s ) Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) ] + η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , n , k s ] a 1 R e ( Δ t 2 r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) Δ t r i 2 v i , n , k s + Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) ) + Δ t cos α ( T i , n , k s G r ( T ) R e 2 + C i , n , k s G r ( c ) R e 2 ) , (63)

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

v i , n , k s + 1 = v i , n , k s Δ t 2 Δ r ( v i + 1, n , k s v i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( v i , n + 1, k s v i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( v i , n , k + 1 s v i , n , k + 1 s ) w i , n , k s Δ t r i u i , n , k s v i , n , k s Δ t r i p θ + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T i , n , k s ] a 1 × Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T ) [ Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) + Δ t ( Δ r ) 2 ( v i + 1, n , k s 2 v i , n , k s + v i 1, n , k s ) Δ t r i 2 v i , n , k s + Δ t r i 2 ( Δ θ ) 2 ( v i , n + 1, k s 2 v i , n , k s + v i , n 1, k s ) + Δ t ( Δ z ) 2 ( v i , n , k + 1 s 2 v i , n , k s + v i , n , k 1 s ) + Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) ] (64)

v i , n + 1, k s , v i , n , k s and v i , n 1, k s are known thus v i , n , k s + 1 can be computed explicitly.

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

v i , n , k s + 1 = v i , n , k s Δ t 2 Δ r ( v i + 1, n , k s v i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( v i , n + 1, k m v i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( v i , n , k + 1 s v i , n , k + 1 s ) w i , n , k s Δ t r i u i , n , k s v i , n , k s Δ t r i p θ + 1 R e [ η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , n , k s ] a 1 ] Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) + ( θ η ( a 1 ) [ 1 + λ T i , n , k s ] a 1 R e + ν T ) [ Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) + Δ t ( Δ r ) 2 ( v i + 1, n , k s 2 v i , n , k s + v i 1, n , k s ) Δ t r i 2 v i , n , k s + Δ t r i 2 ( Δ θ ) 2 ( v i , n + 1, k s 2 v i , n , k s + v i , n 1, k s ) + Δ t ( Δ z ) 2 ( v i , n , k + 1 s 2 v i , n , k s + v i , n , k 1 s ) + Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) ] , (65)

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

w i , n , k s + 1 = w i , n , k s Δ t 2 Δ r ( w i + 1, n , k s w i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( w i , n + 1, k s w i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( w i , n , k + 1 s w i , n , k 1 s ) w i , n , k s Δ t p z + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T i , n , k s ] a 1 × ( Δ t 2 r i 2 Δ θ ( w i , n + 1, k s w i , n 1, k s ) + Δ t 2 r i Δ z ( v i , n , k + 1 s v i , n , k + 1 s ) ) + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T ) [ Δ t 2 r i Δ r ( w i + 1, n , k s w i 1, n , k s ) + Δ t ( Δ r ) 2 ( w i + 1, n , k s 2 w i , n , k s + w i 1, n , k s ) + Δ t r i 2 ( Δ θ ) 2 ( w i , n + 1, k s 2 w i , n , k s + w i , n 1, k s ) + Δ t ( Δ z ) 2 ( w i , n , k + 1 s 2 w i , n , k s + w i , n , k 1 s ) ] + Δ t sin α [ T i , n , k s G r ( T ) R e 2 + C i , n , k s G r ( c ) R e 2 ] , (66)

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

w i , n , k s + 1 = w i , n , k s Δ t 2 Δ r ( w i + 1, n , k s w i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( w i , n + 1, k s w i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( w i , n , k + 1 s w i , n , k 1 s ) w i , n , k s Δ t p z + 1 R e [ η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , j , k m ] a 1 ] × ( Δ t 2 r i 2 Δ θ ( w i , n + 1, k s w i , n 1, k s ) + Δ t 2 r i Δ z ( v i , n , k + 1 s v i , n , k + 1 s ) ) + ( θ η ( a 1 ) [ 1 + λ T i , n , k s ] a 1 R e + ν T ) [ Δ t 2 r i Δ r ( w i + 1, n , k s w i 1, n , k s ) + Δ t ( Δ r ) 2 ( w i + 1, n , k s 2 w i , n , k s + w i 1, n , k s ) + Δ t r i 2 ( Δ θ ) 2 ( w i , n + 1, k s 2 w i , n , k s + w i , n 1, k s ) + Δ t ( Δ z ) 2 ( w i , n , k + 1 s 2 w i , n , k s + w i , n , k 1 s ) ] + Δ t sin α [ T i , n , k s G r ( T ) R e 2 + C i , n , k s G r ( c ) R e 2 ] , (67)

Energy equation in discrete form:

T i , n , k s + 1 = T i , n , k s Δ t 2 Δ r ( T i + 1, n , k s T i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( T i , n + 1, k s T i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( T i , n , k + 1 s T i , n , k 1 s ) w i , n , k s + ω P r R e [ Δ t 4 ( Δ r ) 2 ( T i + 1, n , k s T i 1, n , k s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( T i , n + 1, k s T i , n 1, k s 2 ) 2 + Δ t 4 ( Δ z ) 2 ( T i , n , k + 1 s T i , n , k 1 s ) 2 ] + [ 1 + ω T i , n , k s P r R e + φ T ] [ Δ t 2 r i Δ r ( T i + 1, n , k s T i 1, n , k s ) + Δ t ( Δ r ) 2 ( T i + 1, n , k s 2 T i , n , k s + T i 1, n , k s ) + Δ t r i 2 ( Δ ) 2 ( T i , n + 1, k s 2 T i , n , k s + T i , n 1, k s ) + Δ t ( Δ z ) 2 ( T i , n , k + 1 s 2 T i , n , k s + T i , n , k 1 s ) ]

+ E c R e θ η ( a 1 ) 1 + λ T i , n , k s { 2 [ Δ t 4 ( Δ r ) 2 ( u i + 1, n , k s u i 1, n , k s ) 2 + Δ t 4 r i ( Δ ) 2 ( v i , n + 1, k s v i , n 1, k s ) 2 + Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) u i , n , k s + Δ t r i 2 ( u i , n , k s ) 2 + Δ t 2 Δ z ( w i , n , k + 1 s w i , n , k 1 s ) 2 ] + Δ t 4 r i 2 ( Δ θ ) 2 ( u i , n + 1, k s u i , n 1, k s ) 2 + Δ t 2 r i Δ r Δ θ ( u i , n + 1, k s u i , n 1, k s ) ( v i + 1, n , k s v i 1, n , k s ) Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s + Δ t 4 ( Δ r ) 2 ( v i + 1, n , k s v i 1, n , k s ) 2

Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) v i , n , k s + Δ t r i 2 ( v i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( v i , n , k + 1 s v i , n , k 1 s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( w i , n + 1, k s w i , n 1, k s ) 2 + Δ t 2 r i Δ θ Δ z ( v i , n , k + 1 s v i , n , k 1 s ) ( w i , n + 1, k s w i , n 1, k s ) + Δ t 4 ( Δ r ) 2 ( w i + 1, n , k s w i 1, n , k s ) 2 + Δ t 4 Δ r Δ z ( w i + 1, n , k s w i 1, n , k s ) ( u i , n , k + 1 s u i , n , k 1 s ) + Δ t 4 ( Δ z ) 2 ( u i , n , k + 1 s u i , n , k 1 s ) 2 } + E c ε i , n , k s , (68)

Concentration equation in discrete form:

C i , n , k s + 1 = C i , n , k s Δ t 2 Δ r ( C i + 1, n , k s C i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( C i , n + 1, k s C i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( C i , n , k + 1 s C i , n , k 1 s ) w i , n , k s + ( 1 S c R e + φ c ) [ Δ t 2 r i Δ r ( C i + 1, n , k s C i 1, n , k s ) + Δ t ( Δ r ) 2 ( C i + 1, n , k s 2 C i , n , k s + C i 1, n , k s ) + Δ t r i 2 ( Δ θ ) 2 ( C i , n + 1, k s 2 C i , n , k s + C i , n 1, k s ) + Δ t ( Δ z ) 2 ( C i , n , k + 1 s 2 C i , n , k s + C i , n , k 1 s ) ] Δ t k c R e C i , n , k s , (69)

Discrete ϰ -equation for the liquid phase:

ϰ i , n , k s + 1 = ϰ i , n , k s Δ t 2 Δ r ( ϰ i + 1, n , k s ϰ i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( ϰ i , n + 1, k s ϰ i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( ϰ i , n , k + 1 s ϰ i , n , k 1 s ) w i , n , k s + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T i , n , k s ] a 1 Δ t 2 r i 2 Δ θ ( ϰ i , n + 1, k s ϰ i , n 1, k s ) + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T σ ϰ ) [ Δ t ( Δ r ) 2 ( ϰ i + 1, n , k s 2 ϰ i , n , k s + ϰ i 1, n , k s ) + Δ t 2 r i Δ r ( ϰ i + 1, n , k s ϰ i 1, n , k s ) + Δ t r i 2 ( Δ θ ) 2 ( ϰ i , n + 1, k s 2 ϰ i , n , k s + ϰ i , n 1, k s ) + Δ t ( Δ z ) 2 ( ϰ i , n , k + 1 s 2 ϰ i , n , k s + ϰ i , n , k 1 s ) ] + 2 ν T [ Δ t 4 ( Δ r ) 2 ( u i + 1, n , k s u i 1, n , k s ) 2

+ Δ t 4 r i 2 ( Δ θ ) 2 ( v i , n + 1, k s v i , n 1, k s ) 2 + Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) u i , n , k s + Δ t r i 2 ( u i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( w i , n , k + 1 s w i , n , k 1 s ) 2 ] + ν T { Δ t 4 r i 2 ( Δ θ ) 2 ( u i , n + 1, k s u i , n 1, k s ) 2 + Δ t 2 r i Δ r Δ θ ( u i , n + 1, k s u i , n 1, k s ) ( v i + 1, n , k s v i 1, n , k s ) Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s + Δ t 4 ( Δ r ) 2 ( v i + 1, n , k s v i 1, n , k s ) 2 Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) v i , n , k s + Δ t r i 2 ( v i , n , k s ) 2

+ Δ t 4 ( Δ z ) 2 ( v i , n , k + 1 s v i , n , k 1 s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( w i , n + 1, k s w i , n 1, k s ) 2 + Δ t 2 r i Δ θ Δ z ( v i , n , k + 1 s v i , n , k 1 s ) ( w i , n + 1, k s w i , n 1, k s ) + Δ t 4 ( Δ r ) 2 ( w i + 1, n , k s w i 1, n , k s ) 2 + Δ t 4 Δ r Δ z ( w i + 1, n , k s w i 1, n , k s ) ( u i , n , k + 1 s u i , n , k 1 s ) + Δ t 4 ( Δ z ) 2 ( u i , n , k + 1 s u i , n , k 1 s ) 2 } ε i , n , k s , (70)

Discrete ϰ -equation for the gas phase:

ϰ i , n , k s + 1 = ϰ i , n , k s Δ t 2 Δ r ( ϰ i + 1, n , k s ϰ i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( ϰ i , n + 1, k s ϰ i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( ϰ i , n , k + 1 s ϰ i , n , k 1 s ) w i , n , k s + η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , n , k s ] a 1 R e × Δ t 2 r i 2 Δ θ ( ϰ i , n + 1, k s ϰ i , n 1, k s ) + ( θ η ( a 1 ) [ 1 + λ T i , n , k s ] a 1 R e + ν T σ ϰ ) × [ Δ t ( Δ r ) 2 ( ϰ i + 1, n , k s 2 ϰ i , n , k s + ϰ i 1, n , k s ) + Δ t 2 r i Δ r ( ϰ i + 1, n , k s ϰ i 1, n , k s ) + Δ t r i 2 ( Δ θ ) 2 ( ϰ i , n + 1, k s 2 ϰ i , n , k s + ϰ i , n 1, k s ) + Δ t ( Δ z ) 2 ( ϰ i , n , k + 1 s 2 ϰ i , n , k s + ϰ i , n , k 1 s ) ]

+ 2 ν T [ Δ t 4 ( Δ r ) 2 ( u i + 1, n , k s u i 1, n , k s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( v i , n + 1, k s v i , n 1, k s ) 2 + Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) u i , n , k s + Δ t r i 2 ( u i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( w i , n , k + 1 s w i , n , k 1 s ) 2 ] + ν T { Δ t 4 r i 2 ( Δ θ ) 2 ( u i , n + 1, k s u i , n 1, k s ) 2 + Δ t 2 r i Δ r Δ θ ( u i , n + 1, k s u i , n 1, k s ) ( v i + 1, n , k s v i 1, n , k s ) Δ t r i 2 Δ θ ( U i , n + 1, k s U i , n 1, k s ) v i , n , k s + Δ t 4 ( Δ r ) 2 ( v i + 1, n , k s v i 1, n , k s ) 2

Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) v i , n , k s + Δ t r i 2 ( v i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( v i , n , k + 1 s v i , n , k 1 s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( w i , n + 1, k s w i , n 1, k s ) 2 + Δ t 2 r i Δ θ Δ z ( v i , n , k + 1 s v i , n , k 1 s ) ( w i , n + 1, k s w i , n 1, k s ) + Δ t 4 ( Δ r ) 2 ( w i + 1, n , k s w i 1, n , k s ) 2 + Δ t 4 Δ r Δ z ( w i + 1, n , k s w i 1, n , k s ) ( u i , n , k + 1 s u i , n , k 1 s ) + Δ t 4 ( Δ z ) 2 ( u i , n , k + 1 s u i , n , k 1 s ) 2 } ε i , n , k s , (71)

Discrete ε-equation for the liquid phase:

ε i , n , k s + 1 = ε i , n , k s Δ t 2 Δ r ( ε i + 1, n , k s ε i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( ε i , n + 1, k s ε i , n 1, k s ) v i , n , k m Δ t 2 Δ z ( ε i , n , k + 1 s ε i , n , k 1 s ) w i , n , k s + η ( a 1 ) θ η ( a 2 ) R e [ 1 + λ T i , j , k m ] a 1 Δ t 2 r i 2 Δ θ ( ε i , n + 1, k s ε i , n 1, k s ) + ( θ η ( a 1 ) R e [ 1 + λ T i , n , k s ] a 1 + ν T σ ε ) [ Δ t 2 r i Δ r ( ε i + 1, n , k s ε i 1, n , k s ) + Δ t ( Δ r ) 2 ( ε i + 1, n , k s 2 ε i , n , k s + ε i 1, n , k s ) + Δ t r i 2 ( θ ) 2 ( ε i , n + 1, k s 2 ε i , n , k s + ε i , n 1, k s ) + Δ t ( Δ z ) 2 ( ε i , n , k + 1 s 2 ε i , n , k s + ε i , n , k 1 s ) ] + 2 f 1 C ε ,1 ε i , n , k s ϰ i , n , k s ν T [ Δ t 4 ( Δ r ) 2 ( u i + 1, n , k s u i 1, n , k s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( v i , n + 1, k s v i , n 1, k s ) 2

+ Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) u i , n , k s + Δ t r i 2 ( u i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( w i , n , k + 1 s w i , n , k 1 s ) 2 ] + f 1 C ε ,1 ε i , j , k m ϰ i , n , k s ν T { Δ t 4 r i 2 ( Δ θ ) 2 ( u i , n + 1, k s u i , n 1, k s ) 2 + Δ t 2 r i Δ r Δ θ ( u i , n + 1, k s u i , n 1, k s ) ( v i + 1, n , k s v i 1, n , k s ) Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s + Δ t 4 ( Δ r ) 2 ( v i + 1, n , k s v i 1, n , k s ) 2 Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) v i , n , k s

+ Δ t r i 2 ( v i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( v i , n , k + 1 s v i , n , k 1 s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( w i , n + 1, k s w i , n 1, k s ) 2 + Δ t 2 r i Δ θ Δ z ( v i , n , k + 1 s v i , n , k 1 s ) ( w i , n + 1, k s w i , n 1, k s ) + Δ t 4 ( Δ r ) 2 ( w i + 1, n , k s w i 1, n , k s ) 2 + Δ t 4 Δ r Δ z ( w i + 1, n , k s w i 1, n , k s ) ( u i , n , k + 1 s u i , n , k 1 s ) + Δ t 4 ( Δ z ) 2 ( u i , n , k + 1 s u i , n , k 1 s ) 2 } Δ t f 2 C ε ,2 ( ε i , n , k s ) 2 ϰ i , n , k s (72)

Discrete ε-equation for the gas phase:

ε i , n , k s + 1 = ε i , n , k s Δ t 2 Δ r ( ε i + 1, n , k s ε i 1, n , k s ) u i , n , k s Δ t 2 r i Δ θ ( ε i , n + 1, k s ε i , n 1, k s ) v i , n , k s Δ t 2 Δ z ( ε i , n , k + 1 s ε i , n , k 1 s ) w i , n , k s + η ( a 1 ) θ η ( a 2 ) [ 1 + λ T i , n , k s ] a 1 R e × Δ t 2 r i 2 Δ θ ( ε i , n + 1, k s ε i , n 1, k s ) + ( θ η ( a 1 ) [ 1 + λ T i , n , k s ] a 1 R e + ν T σ ϰ ) × [ Δ t 2 r i Δ r ( ε i + 1, n , k s ε i 1, n , k s ) + Δ t ( Δ r ) 2 ( ε i + 1, n , k s 2 ε i , n , k s + ε i 1, n , k s ) + Δ t r i 2 ( θ ) 2 ( ε i , n + 1, k s 2 ε i , n , k s + ε i , n 1, k s ) + Δ t ( Δ z ) 2 ( ε i , n , k + 1 s 2 ε i , n , k s + ε i , n , k 1 s ) ]

+ 2 f 1 C ε ,1 ε i , n , k s ϰ i , n , k s ν T [ Δ t 4 ( Δ r ) 2 ( u i + 1, n , k s u i 1, n , k s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( v i , n + 1, k s v i , n 1, k s ) 2 + Δ t r i 2 Δ θ ( v i , n + 1, k s v i , n 1, k s ) u i , n , k s + Δ t r i 2 ( u i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( w i , n , k + 1 s w i , n , k 1 s ) 2 ] + f 1 C ε ,1 ε i , n , k s ϰ i , n , k s ν T { Δ t 4 r i 2 ( Δ θ ) 2 ( u i , n + 1, k s u i , n 1, k s ) 2 + Δ t 2 r i Δ r Δ θ ( u i , n + 1, k s u i , n 1, k s ) ( v i + 1, n , k s v i 1, n , k s ) Δ t r i 2 Δ θ ( u i , n + 1, k s u i , n 1, k s ) v i , n , k s + Δ t 4 ( Δ r ) 2 ( v i + 1, n , k s v i 1, n , k s ) 2 Δ t 2 r i Δ r ( v i + 1, n , k s v i 1, n , k s ) v i , n , k s

+ Δ t r i 2 ( v i , n , k s ) 2 + Δ t 4 ( Δ z ) 2 ( v i , n , k + 1 s v i , n , k 1 s ) 2 + Δ t 4 r i 2 ( Δ θ ) 2 ( w i , n + 1, k s w i , n 1, k s ) 2 + Δ t 2 r i Δ θ Δ z ( v i , n , k + 1 s v i , n , k