Computational Modeling and Dynamics of the Oil and Water Flow Using the Galerkin Approximation

Show more

1. Introduction

The two-phase flows are found widely in the nature and in a whole range of industrial applications. The two-phase material is often called the mixture [1] . In a single phase flow, the design parameter can be modeled easily but the presence of a secondary phase increases the complexity of the fluid flow. The two phase flow phenomena is still difficult to understand because of their complex behavior [2] . With the progress of science and technology, the study of mathematical models which may describe the two-phase flow phenomena has taken an increasing attention and becomes a challenge to the researcher [3] [4] [5] .

A variety of two-phase flows can occur depending on the combination of two-phases such as gas-liquid, liquid-solid, liquid-liquid and gas-solid. Generally, we see these two phase flows in transport of oil-gas mixtures in the pipeline system, air conditioning and refrigeration system, crystallization system, clay extraction, oil water mixture flow in pipelines and fuel combustion etc.

There are various studies which have been done on two phase flow. Aparallel streamline Upwind Petrov-Galerkin model using the Finite Element method for two-phase flow in a T-junction showed that the primary phase accumulates at the middle of the pipe and top of the wall while the secondary phase accumulates near the wall after the junction [6] . In an air-oil flow, it was found that oil was concentrating at the lower portion of pipe’s cross section due to gravity. A wave of oil and air mixture was present at the middle that covers the whole pipe cross section area and a new wave of the mixture was just entering the pipe segment [7] . A mathematical model of gas liquid incompressible fluid flow was introduced to explore the feature of fluid under certain centrifugal force in vertical centrifugal casting using the projection level set method and showed that the free surface near the wall was higher compared to other parts of the domain. The distribution of liquid phase calculated by the projection level set method coincided with the physical truth and the pressure distribution was found a parabola [8] . So many types of flow patterns were predicted at different phase velocities of the flow of two immiscible liquids through the horizontal pipelines using volume of fluid method. The slug flow appeared at low velocities while in annular flow, the total pressure of the mixture decreases with the increases of the oil velocity. It is found that the volume fraction of oil is maximum at the centre in a core annular flow and at the upper side of the pipe in case of stratified flow [2] . In a stratified oil-water two-phase flow the pressure gradient increases as the velocity increases [9] . The water phase accumulates at the bottom of the pipe for low velocities in a horizontal pipe [10] . Solid phase velocity profile is asymmetrical about the central axis and maximum velocity position moves toward the top of the pipe in a solid-liquid slurry flow [11] . The high velocity region is observed at the outer side of the pipe in a CFD simulation of two-phase flow in helical pipe using population balanced model [12] . The radial velocity of a dispersed phase is larger than that of a gaseous phase. Also lighter and smaller ash particles move faster than corundum particles of same mass [13] . For different inlet velocities, pressure loss was different in two-phase flow of vertical pipe and increasing gas velocity reduced the total pressure loss in the pipe [14] .

Nowadays CFD has become an essential tool in simulation studies. In our study, a two-phase flow model is solved by using the Finite Element approximation based on the Galerkin method. The two-phases are considered as oil and water respectively. The Finite Element method is a numerical technique to obtain an approximate solution to boundary value problems and the method has a distinct advantage that it allows naturally for handling complex arbitrary geometries so that it can be applied using irregular grid of various shapes. It also provides a set of functions which give the variation of differential equations between the grid points [15] .

The oil-water flow is very common in the petroleum industries. Generally the oil phase is transported in a multiphase flow condition as water and oil are normally produced together. The presence of water has a significant effect during the transportation of oil. The complex interfacial structure of oil water flow makes it difficult to predict the hydrodynamics of the fluid flow [10] [16] . A comparison was done for the different shape of photo bioreactor for Microalgae culture and we considered the tubular one is better [17] .

In literature, we found different models for two-phase flows including the Petrov-Galerkin method, the projection level set method, the volume of fluid method, and the population balanced model.

In our study, we considered a FEM model for two-phase flow using the Galerkin method and an oil-water mixture flow is taken to simulate according to our model. To investigate the velocity profile, the volume fraction and the pressure distribution of the oil-water two-phase flow through a circular tube we applied boundary and initial condition in the inlet, outlet and wall of a typical computational domain.

2. Mathematical Model Formulation

In this section, a mathematical model will be developed to discuss the two-phase flow. A computational domain will also be constructed to simulate the fluid flow. A suitable mesh design is considered as the computational domain. The initial and boundary conditions are for the oil-water two-phase flow in a phase field platform.

2.1. Governing Equations

In this study, we consider a two-phase incompressible Newtonian flow in which the phases are liquid and the flow mixture is laminar. Therefore, the governing equation of the problem becomes as follows along with Boundary Conditions (BC) and Initial conditions (IC) in the computational domain $\Omega $ ,

$\nabla \cdot u=0$ (1)

$\rho \left(\frac{\partial u}{\partial t}+u\cdot \nabla u\right)=-\nabla p+\nabla \cdot \left(\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)\right)+\rho g+{F}_{st}\text{\hspace{0.17em}}\text{\hspace{0.17em}},$ (2)

where $u$ denotes the velocity of the mixture, $\rho $ and $\mu $ are its density and viscosity respectively, $g$ is the gravity and ${F}_{st}$ is the surface tension force. The phase field variable [18] $\phi $ is governed by the following partial differential equation (PDE)

$\frac{\partial \phi}{\partial t}+u\cdot \nabla \phi =\text{}\text{\hspace{0.05em}}\nabla \cdot \gamma \nabla G$ , (3)

where G is the chemical potential.

The following equation is used to track the interface in the Phase field method obtained from the Cahn-Hilliard equation [19] which describes the two-phases.

$\frac{\partial \phi}{\partial t}+u\cdot \nabla \phi =\nabla \cdot \frac{\gamma \lambda}{{\epsilon}^{2}}\nabla \psi $ , (4)

$\psi =-\nabla \cdot {\epsilon}^{2}\nabla \phi +\left({\phi}^{2}-1\right)\phi +\left(\frac{{\epsilon}^{2}}{\lambda}\right)\frac{\partial f}{\partial \phi}$ , (5)

where $\phi $ is the dimensionless phase field variable, $\gamma $ is the mobility, $\epsilon $ is a

controlling interface parameter, $\lambda $ is the mixing energy density, the term $\frac{\partial f}{\partial \phi}$

denotes the $\phi $ derivative of external free energy and $\chi $ is the mobility tuning parameter. The density and viscosity of the mixture are the function of volume

fraction of water ${V}_{w}$ . The volume fraction of water and oil are ${V}_{w}=\frac{1+\phi}{2}$ and ${V}_{o}=\frac{1-\phi}{2}$ respectively. The density and dynamic viscosity of the two-phase

model are calculated over the interface according to

$\rho ={\rho}_{o}+\left({\rho}_{w}-{\rho}_{o}\right){V}_{w}$ (6)

$\mu ={\mu}_{o}+\left({\mu}_{w}-{\mu}_{o}\right){V}_{w}$ , (7)

where the subscript w and o are used for water and oil respectively. The surface tension force for the method is applied as a body force

${F}_{st}=G\nabla \phi $ . (8)

The chemical potential G is given by the following equation related to the phase field variable and the interface controlling parameter.

$G=\lambda \left[-{\nabla}^{2}\phi +\frac{\phi \left({\phi}^{2}-1\right)}{{\epsilon}^{2}}\right]$ . (9)

2.2. Initial and Boundary Conditions for the Simulation

As the mixture of oil and water is incompressible, the flow is assumed to be laminar two-phase flow where the volume fraction of water and oil are considered 0.91 and 0.09 respectively [20] . The phase field method gives better interface than other methods. In phase field method, the interfacial layer is governed by a phase field variable which is obtained from the Cahn-Hilliard equation. The COMSOL Multiphysics software is used for the CFD simulation. At the inlet we considered velocity $\left|u\right|={u}_{0}$ , a no slip boundary condition on the wall and the traction boundary condition were imposed at the outlet which are given by

$u=0$ (10)

$T=\sigma \cdot n=\stackrel{\xaf}{T}$ , where $\text{\hspace{0.17em}}\sigma =-pI+\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)$ (11)

2.3. Computational Domain and Mesh Design

Atypical circular tube domain as in Figure 1 is considered in our study with 8 m length and 0.05 m radius while the working volume and surface area are 0.0618 m^{3} and 2.518 m^{2} respectively. As the mesh design plays a significant role in accuracy of numerical results, we considered a suitable mesh design for our simulation with 81,143 elements and 619,748 degrees of freedom according to Figure 2.

3. Finite Element Approximation

The variational statement of a boundary value problem is the integral representation of the problem. We can obtain the variational statement by setting the total weighted residual error to zero. When the local residual error is multiplied by a weighting function and then integrating over the domain, the weighted residual error is obtained [15] . The aim is to find an approximate solution such that the total weighted residual error vanishes. If the weighting functions are same as the basis functions, the method is known as Galerkin method. The finite element method is used to calculate the basis functions for the Galerkin approximation of the boundary value problem.

3.1. Variational Statement

$r\left(x,t\right)=\rho \left(\frac{\partial u}{\partial t}+u\cdot \nabla u\right)+\nabla p-\nabla \cdot \left(\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)\right)-\rho g-{F}_{st}$ (12)

According to weighted residual method [15] there is $\left(u,p\right)\in V\times Q$ such that for every $t\in I$

Figure 1. The geometry of the computational flow domain.

(a) (b)

Figure 2. Mesh design of the computational domain (a) inlet and outlet (b) along the pipe.

$\left(r,v\right)=0\text{\hspace{1em}}\text{\hspace{0.17em}}\forall v\in V$ and $v=0$ on ${\Gamma}_{u}$ (13)

$\left(\nabla \cdot u,q\right)=0,\forall q\in Q$ (14)

$u\left(x,0\right)={u}_{0}$ in $\Omega $ (15)

$u=\stackrel{\xaf}{u}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\Gamma}_{u}$ (16)

where $V$ and $Q$ are velocity and pressure spaces and the inner product is defined by

$\left(a,b\right)={\displaystyle \underset{\Omega}{\int}a\cdot b}\text{\hspace{0.17em}}\text{d}\Omega \text{\hspace{0.17em}}$ (17)

The integral representation of (13) is presented as follows

$\underset{\Omega}{\int}\rho}\frac{Du}{Dt}\cdot v\text{d}\Omega +{\displaystyle \underset{\Omega}{\int}\left(-p\nabla \cdot v+\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)\nabla v\right)\text{\hspace{0.17em}}\text{d}\Omega}-{\displaystyle \underset{\Omega}{\int}\rho g\cdot v}\text{d}\Omega -{\displaystyle \underset{\Omega}{\int}{F}_{st}\cdot v\text{d}\Omega$ $={\displaystyle \underset{\partial \Omega}{\int}\stackrel{\xaf}{T}\cdot v\text{d}s}$ (18)

Since $\partial \Omega ={\Gamma}_{u}\cup {\Gamma}_{t}$ and $u$ is specified on ${\Gamma}_{u}$ and $v$ is chosen to be zero on ${\Gamma}_{u}$ .

The variational statement of the problem is,

$\begin{array}{l}\text{Find}\text{\hspace{0.17em}}\left(u,p\right)\in V\times Q\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{such}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{that}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{every}\text{\hspace{0.17em}}\text{\hspace{0.17em}}t\in I\\ \begin{array}{l}\left(\rho \frac{Du}{Dt},v\right)-\left(p,\nabla \cdot v\right)+\left(\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)\text{\hspace{0.17em}},\nabla v\right)-\left(\rho g,v\right)-\left({F}_{st},v\right)=b\left(\stackrel{\xaf}{T},v\right)\text{\hspace{0.17em}}\text{\hspace{1em}}\forall v\in {V}_{0}\\ \left(\nabla \cdot u,q\right)=0,\text{\hspace{1em}}\forall q\in Q\text{\hspace{1em}}u\left(x,0\right)={u}_{0}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\Omega \\ u=\stackrel{\xaf}{u}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\Gamma}_{u}\\ V=\left\{v\text{\hspace{0.17em}}|\text{\hspace{0.17em}}v\in {\left[{H}^{1}\left(\Omega \right)\right]}^{2}\right\},\text{}{V}_{0}=\{v|v\in V\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}v=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma}_{u}\}\\ Q=\left\{\beta \text{\hspace{0.17em}}|\text{\hspace{0.17em}}\beta \in {H}^{1}\left(\Omega \right)\right\}\end{array}\}\end{array}$
(19)

This variational statement is also known as the weak formulation of a two-phase Newtonian incompressible flow in a computational domain $\Omega $ .

3.2. The Galerkin Finite Element Approximation

To find the Finite Element approximation, let ${V}_{h}\subset V$ be a N-Dimensional subspace of $V$ with basis functions $\left\{{\varphi}_{1},{\varphi}_{2},\cdots ,{\varphi}_{N}\right\}$ . Approximating $v$ and q in (19) by

${v}_{h}={\displaystyle \underset{k=1}{\overset{N}{\sum}}{\varphi}_{k}{v}_{k}}$ and $q={\displaystyle \underset{p=1}{\overset{M}{\sum}}{\theta}_{p}}{q}_{p}$

$\begin{array}{l}{\displaystyle \underset{k=1}{\overset{N}{\sum}}\{\left(\rho \frac{Du}{Dt},{\varphi}_{k}\right)-\left(p,\text{\hspace{0.17em}}\nabla \cdot {\varphi}_{k}\right)+\left(\mu \left(\nabla u+\nabla {u}^{\text{T}}\right)\text{\hspace{0.17em}},\nabla {\varphi}_{k}\right)-\left(\rho g\text{\hspace{0.17em}},{\varphi}_{k}\right)-\left({F}_{st}\text{\hspace{0.17em}},{\varphi}_{k}\right)}\end{array}$ $-b\left(\stackrel{\xaf}{T},{\varphi}_{k}\right)\}{v}_{k}=0$ (20)

$\underset{p=1}{\overset{M}{\sum}}\left(\nabla \cdot u,\text{\hspace{0.17em}}{\theta}_{p}\right)}\text{\hspace{0.17em}}{q}_{p}=0\text{\hspace{1em}$ (21)

which can be written as

$\begin{array}{l}\left(\rho \frac{\partial u}{\partial t},{\varphi}_{k}\right)+\left(\rho u\cdot \nabla u,{\varphi}_{k}\right)-\left(p,\nabla \cdot {\varphi}_{k}\right)+\left(\mu \left(\nabla u+\nabla {u}^{\text{T}}\right),\nabla {\varphi}_{k}\right)\\ =\left(\rho g,{\varphi}_{k}\right)-\left({F}_{st},{\phi}_{k}\right)-b\left(\stackrel{\xaf}{T},{\varphi}_{k}\right)\end{array}$ (22)

$\left(\nabla \cdot u\text{\hspace{0.17em}},\text{\hspace{0.17em}}{\theta}_{p}\right)=0$ (23)

Now approximating $u$ and p respectively by

${u}_{h}={\displaystyle \underset{l=1}{\overset{N}{\sum}}{\varphi}_{l}{u}_{l}}$ and ${p}_{h}={\displaystyle \underset{p=1}{\overset{M}{\sum}}{\omega}_{p}{p}_{p}}$

Therefore, we get from (22) and (23) that

$\begin{array}{l}{\displaystyle \underset{l=1}{\overset{N}{\sum}}\left\{\left(\rho {\varphi}_{l},{\varphi}_{k}\right){\stackrel{\dot{}}{u}}_{l}+\left(\rho \text{\hspace{0.17em}}{\varphi}_{l}\cdot \nabla {\varphi}_{l},{\varphi}_{k}\right){u}_{l}+\left(\mu \nabla {\varphi}_{l},\nabla {\varphi}_{k}\right){u}_{l}+\left(\mu \nabla {\varphi}_{l}^{\text{T}},\nabla {\varphi}_{k}\right)\text{\hspace{0.17em}}{u}_{l}^{\text{T}}\right\}}\\ -{\displaystyle \underset{p=1}{\overset{M}{\sum}}\left({\omega}_{p},\nabla \cdot {\varphi}_{k}\right){p}_{p}=\left(\rho g,{\varphi}_{k}\right)+\left({F}_{st},{\varphi}_{k}\right)+b\left(\stackrel{\xaf}{T},{\varphi}_{k}\right)}\end{array}$ (24)

$\underset{k=1}{\overset{N}{\sum}}\left(\nabla \cdot {\varphi}_{k},{\omega}_{p}\right)}\text{\hspace{0.17em}}{u}_{k}=0$ (25)

which can be represented by

$\begin{array}{l}M\stackrel{\dot{}}{U}+AU-CP=F\\ -{C}_{1}^{\text{T}}{U}_{1}-{C}_{2}^{\text{T}}{U}_{2}-{C}_{3}^{\text{T}}{U}_{3}=0\end{array}$

using $M=\left({m}_{kl}\right)$ with ${m}_{kl}=\left(\rho {\varphi}_{k},{\varphi}_{l}\right)$

$D=\left({D}_{kl}\right)$ with ${D}_{kl}=\left(\rho {\varphi}_{l}\cdot \nabla {\varphi}_{l},{\varphi}_{k}\right)$

${K}_{ij}={K}_{ijkl}=\left(\mu \nabla {\varphi}_{k},\nabla {\varphi}_{l}\right)$

$A={K}_{ij}+D$

$C={C}_{ikp}$ with ${C}_{ikp}=\left({\omega}_{p},\nabla \cdot {\varphi}_{k}\right)$

$F={F}_{ik}$ with ${F}_{ik}=\left(\rho g,{\varphi}_{k}\right)+\left({F}_{st},{\varphi}_{k}\right)+b\left(\stackrel{\xaf}{T},{\varphi}_{k}\right)$

$k,l=\left(1,\text{\hspace{0.17em}}2,\text{\hspace{0.17em}}\cdots ,N\right)$ and $i,j=1,2,3$

4. Numerical Results and Discussion

In our study, COMSOL Multiphysics version 4.2 has been used to simulate the oil-water two-phase flows. Properties of the fluid phases and parameters values are given in Table 1 [9] and Table 2 [18] respectively. The zero normal stress is considered at the outlet in our simulation. We have analyzed the velocity magnitude, the volume fractions of fluids, the pressure and the shear rate distributions for this two-phase flow mixture. It is also notable that as we simulated a two-phase flow, the interface thickness is an important criterion to understand the general state of the mixture. In this part we discuss about the oil-water interfaces for different times.

Table 1. Physical properties of fluid phases.

The Figure 3(a) shows the slices of the velocity magnitude while the Figure 3(b) represents their corresponding line graph at three different cross sections of the horizontal pipe at $x=1\text{m}$ , $x=4\text{m}$ and $x=7\text{m}$ respectively. It is found that velocity magnitude is generally higher at the middle part of the tube and the lowest at the boundary of the tube. The velocity profiles form a parabolic shape showing the maximum velocity at the center.

Table 2. Parameters value used for simulation.

i) Velocity magnitude at x = 1 m ii) Velocity magnitude at x = 4 m iii) Velocity magnitude at x = 7 m(a)

(b)

Figure 3. (a): The velocity magnitudes at different cross sections of the computational domain. (b): The corresponding line graphs of the Figure 3(a) of the computational domain.

In case of fluid flow inside a tube domain the shear stress can be found on the wall of the domain. We investigated the shear rate distribution for the entire domain which is shown in Figure 4. It shows the shear rate decreases gradually with small fluctuations along the domain.

One of the important parameters used to characterize the two-phase flow is the volume fraction of different phases. Figure 5(a) and Figure 5(b) shows the volume fraction of the oil and the water phase respectively along the horizontal pipe. We see that in the mixture the volume of the oil phase sharply increased while the water phase decreased at the same ratio.

Figure 6 represents the pressure drop along the pipe. A linear pressure drop is observed from the inlet to the outlet though there is a sudden drop near the inlet. This indicates the velocity of the mixture is increased initially and it keeps remain until the end of the tube.

The Figure 7 represents the oil-water interfaces at different times of our simulation. At the beginning of the mixture, the two-phases are mixed completely

Figure 4. The shear rate distribution along the domain.

(a) The Volume Fraction of Oil (b) The Volume fraction of Water

Figure 5. Volume fraction of the oil and the water phases in the mixture.

Figure 6. The pressure distribution along the domain.

Figure 7. The oil-water interfaces at different times of the mixture.

except for a small change around zero of the value of the phase field variable. The phases have started to separate (Figure 7(b)) at t = 1 s. After that the oil and water phases have started to form different phases (Figure 7(c) and Figure 7(d)) according to the phase field variable and we observed the two-phases exist separately. Gradually, the phases formed large phase domains by creating an interface to separate the phases. This shows that the two fluids tend to separate into distinct phases.

5. Conclusion

A CFD model is proposed based on the Galerkin Finite Element approximation to simulate the two-phase flows. According to the Galerkin FEM, we obtained the variational statement of the fluid flow problem by setting the total weighted residual error to zero. A computational domain with 8 m length and 0.05 m radius respectively is constructed to simulate an oil-water two-phase flow. The velocity magnitude has been observed higher at the middle of the tube and the lowest at the boundary of the tube. The shear rate decreased with small fluctuation along the domain. The Oil phase volume fraction increased consequently the water phase decreased at the same ratio. From the inlet to the outlet, a linear pressure drop is observed. It is also notable that the fluids tend to separate into two distinct phases.

Acknowledgements

The authors gratefully acknowledge for the technical support to The Centre of excellence in Mathematics, Department of Mathematics, Mahidol University, Bangkok, Thailand.

References

[1] Drew, D.A. (1983) Mathematical Modeling of Two-Phase Flow. Annual Reviews on Fluid Mechanics, 15, 261-291.

https://doi.org/10.1146/annurev.fl.15.010183.001401

[2] Desamala, A.B., Dasamahapatra, A.K. and Mandal, T.K. (2014) Oil-Water Two-phase Flow Characteristics in Horizontal Pipeline A Comprehensive CFD Study. International journal of Chemical, Molecular, Nuclear, Materials and Metallurgical Engineering, World Academy of Science, Engineering and Technology, 8, 360-364.

[3] Liao, Q., Liu, D., Ye, D., Zhu, X. and Lee, D. (2011) Mathematical Modelling of Two-phase Flow and Transport in an Immobilized-Cell Photo Bioreactor. International Journal of Hydrogen Energy, 36, 13939-13948.

[4] Morin, A. (2012) Mathematical Modelling and Numerical Simulation of Two-phase Multi-component Flows of CO2 Mixtures in Pipes. PhD Thesis, Norwegian University of Science and Technology, Norway.

[5] Hudson, J. and Harris, D. (2006) Numerical Modeling of Eulerian Two-Phase Gas-Solid Flow. Manchester Institute of Mathematical Sciences, Manchester, 1-41.

[6] Giordano, M., Bonfiglioli, A. and Magi, V. (2006) A Parallel Finite Element Method for Two-phase Flows. European Conference on Computational Fluid Dynamics, Glasgow, 1-20.

[7] Kanarachos, S. and Flouros, M. (2014) Simulation of Air-Oil Mixture Flow in the Scanvenge Pipe of an Aero Engine Using Generalized Interface Momentum Exchange Model. WSEAS transaction on Fluid mechanics, 9, 144-153.

[8] Zhou, J.X., Shen, X., Yin, Y.J., Gou, Z. and Wang, H. (2015) Gas Liquid Two-phase Flow Modelling of Incompressible Fluid and Experimental Validation Studies in Vertical Centrifugal Casting. IOP Conference Series: Material science and Engineering, 84, 1-9. http://dx.doi.org/10.1088/1757-899X/84/1/012042

[9] Alias, A.Z.M., Koto, J. and Ahmed, Y.M. (2015) CFD Simulation for Stratified Oil-Water Two-Phase Flow in a Horizontal Pipe. International Society of Ocean, Mechanical and Aerospace Scientists and Engineers, 2, 1-6.

[10] Hu, H. and Cheng, Y.F. (2016) Modelling by Computational Fluid Dynamics Simulation of Pipeline Corrosion in CO2-Containing Oil-Water Two-phase Flow. Journal of Petroleum Science and Engineering, 146, 134-141.

[11] Nabil, T., EL-Sawaf, I. and EL-Nahhas, K. (2013) Computational Fluid Dynamics Simulation of the Solid-Liquid Slurry Flow in a Pipeline. Seventh International Water Technology Conference, 5-7 November 2013, Istanbul, Turkey.

[12] Rathod, R and Hebbal, O. (2013) CFD Simulation of Two-phase Flow and Heat Transfer in Helical Pipe. International Journal of Engineering Research and Technology, 2, 2150-2157.

[13] Kartushinsky, A., Siirde, A., Rudi, U. and Shablinsky, A. (2011) Mathematical Model of Two-phase Flows Loaded with Light and Heavy Particles to Analyze CFB Processes. Oil Shale, 28, 169-180. http://dx.doi.org/10.3176/oil.2011.1S.09

[14] Sanati, A. (2015) Numerical Simulation of Air-Water Two-phase Flow in Vertical Pipe Using k-ε Model. International Journal of Engineering and Technology, 4, 61-70.

https://doi.org/10.14419/ijet.v4i1.3965

[15] Yu, W.H. and Wiwatanapattaphee, B. (2006) Finite Element Method and Application. Mistercopy Publishing Company, Bangkok.

[16] Abed, E.M. and Auda, Z.A. (2014) Simulation and Experiment of Oil-Water Flow with Effect of Heat Transfer in Horizontal Pipe. International Journal of Mechanical Engineering and Applications, 2, 117-127.

http://dx.doi.org/10.11648/j.ijmea.20140206.16

[17] Shahriar, M., Monir, M.I. and Deb, U.K. (2016) Comparative Analysis of Hydrodynamics Behavior of Microalgae Suspensio Flow in Circular, Square and HexagonalShape Photo Bio-Reactors. America Journal of Computational Mathematics, 6, 320 -335. http://dx.doi.org/10.4236/ajcm.2016.64033

[18] Deb, U.K., Chayantrakom, K., Lenbury, Y. and Wiwatanapataphee, B. (2012) Numerical Simulation of Two-Phase Laminar Flow for CO2 and Microalgae Suspension in the HLTP. Latest Advances in System Science and Computational Intelligence, WSEAS Press, Haifa, 53-58.

[19] Cahn, J.W. and Hilliard, J.E. (1958) Free Energy of a Non-Uniform System in Interfacial Energy. Journal of Chemical Physics, 28, 258-267.

https://doi.org/10.1063/1.1744102

[20] Xu, X. (2007) Study on Oil-Water Two-Phase Flow in Horizontal Pipelines. Journal of Petrolium Science and Engineering, 59, 43-58.

https://doi.org/10.1016/j.petrol.2007.03.002