Abstract

A Computational Fluid Dynamic (CFD) model is presented to analyze the flow dynamics of a two phase incompressible flow in the application sectors of different industries. The Finite Element method (FEM) which is based on the Galerkin approximation, has been implemented for this two phase flow model. Generally, two-phase flows can occur in different forms like gas-liquid, liquid-liquid and solid-liquid forms. The Oil and Water two-phase flow is an important phenomenon in petroleum industry for crude oil production and transportation. In our study, a laminar flow of liquid-liquid phase is considered to simulate the flow dynamics where the liquid phases are water and oil. The COMSOL Multiphysics software is used to perform the simulation including velocity profile, volume fraction, shear rate, pressure distributions and interfacial thicknesses at different times. A typical circular tube domain with radius 0.05 m and length 8 m is assumed for our simulation.

Keywords

Two-phase Flow, CFD, Galerkin Method, FEM, Simulation

Share and Cite:

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},$ (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{}\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=\bar{T}$ , where $\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³ and 2.518 m² 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$

$\left(r,v\right)=0\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=\bar{u}\text{in}{\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{d}\Omega$ (17)

The integral representation of (13) is presented as follows

${\displaystyle \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{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 }\bar{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,

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,\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},{\varphi }_k\right)}\end{array}$ $-b\left(\bar{T},{\varphi }_k\right)\}{v}_k=0$ (20)

${\displaystyle \underset{p=1}{\overset{M}{\sum }}\left(\nabla \cdot u,{\theta }_p\right)}{q}_p=0$ (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(\bar{T},{\varphi }_k\right)\end{array}$ (22)

$\left(\nabla \cdot u,{\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 {\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)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(\bar{T},{\varphi }_k\right)}\end{array}$ (24)

${\displaystyle \underset{k=1}{\overset{N}{\sum }}\left(\nabla \cdot {\varphi }_k,{\omega }_p\right)}{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(\bar{T},{\varphi }_k\right)$

$k,l=\left(1,2,\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.

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.

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

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.

Conflicts of Interest

The authors declare no conflicts of interest.

References