The meshless numerical simulation of Kelvin–Helmholtz instability during the wave growth of liquid–liquid slug flow
Introduction
Slug flows are often encountered in the practical applications during the hydrocarbons transportation pipeline system [1]. In general, it is intermittent and is formed by sequences of large gas bubbles followed by the liquid slug. It flows under the high frequencies in which large gas bubbles flow alternately with liquid slug [2], [3]. In the horizontal pipes, it is formed from a stratified flow due to the growth of hydrodynamic instabilities. Here small perturbations in the form of small waves naturally emerge. Those waves continue to grow, capturing the liquid that flows in front of them until the pipe cross-sectional area is fully blocked by the liquid, thus forming the slugs [4]. This wave growth mechanism is known as the Kelvin–Helmholtz instability (KHI) [5], [6].
Slug flow is one of the dangerous flow patterns in the petroleum industry. In order to ensure the safe operation of the pipeline, the knowledge of the slug flow mechanism is essential. The appearance of slug flow in the pipeline will increase the structural disadvantage that leads to crack and corrosion of the pipe [7], whereas high-frequency slug flows increase the rate of the pipeline corrosion [8]. Consequently, it will decrease the hydrocarbon production and increase the equipment maintenance cost [9]. From the viewpoint of multiphase-flow, the understanding of KHI is very important in pipeline design.
The researchers have developed the numerical methods to simulate the KHI in horizontal pipes or horizontal channel. Simmons et al. [10] explored the influence of KHI during the transition to slug flow in horizontal pipes. The basic mechanism of KHI is explained as the existence of the uniform shear stress between two different fluids. To simplify the calculation, the surface tension and viscosity were neglected, and, consequently, KHI is considered only as the fluid instability of a lighter fluid is placed on the top of a heavier one. San & Kara [11] calculated 2-D KHI in order to examine the performance of six Riemann flux formulas. Those formulas are the combination of the Rusanov scheme, Reo scheme, Harten–Lax–Van Leer scheme, first-order centered scheme, advection upstream splitting method, and Marquina scheme. Here, a single mode perturbation was used as a tested case. They concluded that the advection upstream splitting scheme results in the smallest numerical dissipation.
Patnaik et al. [12] developed a numerical simulation for KHI between two stratified fluids by using the finite difference method. They concluded that the minimum Richardson number required for the onset of KHI was at Reynolds number of 100. They claimed that their result is in a good agreement with that of the linear stability theory. On the other hand, Balabel et al. [13] introduced the level set method to solve the hydrodynamic instability. The numerical procedure was developed on the basis of the solution of Navier–Stokes equations and structural 2-D computational grid by using the control volume approach. The results showed that the onset of roll-up mechanism at the interface is largely affected by the considered flow regimes. Next, Lee et al. [14] examined the finite difference method to solve the KHI of multi-component fluids. It was modeled by both of modification of Navier–Stokes and the multi-component Cahn–Hilliard equations. They concluded that the linear growth rate for the KHI of two-component fluids was in good agreement with the analytical results. These finite difference and finite volume methods used the mesh to solve problems, and it can produce lot of numerical problems during the discretizing of complex domains with irregular boundaries [15]. Moreover, those methods were developed on the basis of the polynomial interpolations, whereas they are limited by the algebraic convergence rates [16].
Shadloo et al. [17] used the smoothed particle hydrodynamics (SPH) scheme to simulate the 2-D KHI problem of an incompressible two-phase immiscible fluid. They concluded that the numerical algorithm is able to handle the two-phase fluid flow under the variation of the density ratios. Yue et al. [18] studied also the growth of the KHI at the interface between two ideal gases using the SPH scheme. The linear growth rate in the numerical solution was found to be in a satisfactory agreement with that of the analytical solution, in which the average relative error was 13%. As the SPH is an approximate method, it is noticed that there are some disadvantages if it is not applied correctly in a specific problem. First, it is very sensitive to the number of particles used in the spatial approximation. Second, the precision of the SPH depends on the order of accuracy of the window function [19].
The above facts indicate that there is a strong disagreement from the previous works regarding the best solution to describe the physical phenomena around KHI. For this reason, a systematic numerical study to explain the involved phenomena around the transition to slug flow of liquid–liquid two-phase flow is needed. The aim of the present study is to carry out a numerical investigation of the KHI during the transition to slug flow. In the present numerical study, the radial basis function (RBF) method is used to simulate the phenomena using modification of the Navier–Stokes equations in terms of primitive variables and the Cahn–Hilliard equations in order to capture the evolution of the interface between two fluids.
An attractive and simple approach that is suitable for a wide range of the problems employs the capability of radial basis function (RBF) [20]. It has a good capability to solve the governing equations in the strong form, without any integration is required. In addition, it is also effective to solve the numerical terms with a low numerical diffusion [21]. The main challenge to implement the RBF method is the applicability to handle the large-scale domain. RBF method should be combined with the DDM, when handling large-scale domain. The basic idea of DDM is to divide one large domain into a number of smaller subdomains [22]. The DDM is able to handle large-scale problems, meanwhile it will decrease the calculation accuracy if more subdomains are employed [23].
Several numerical experiments were carried-out in order to examine the accuracy of RBF method combined with DDM. Mai-Duy et al. [22] used it to solve one and two-dimensional elliptic problems. The domain was divided into sixteen subdomains, and the results show a high level of accuracy and a fast rate of convergence. Therefore it presents more attractive than the one-domain method to solve the elliptic problems under a large computational domain. Zhou et al. [23] developed RBF method combine with DDM to solve the elliptic partial differential equations. They investigated the accuracy of the numerical solution during the increase of the numbers of subdomains. In case of 16 16 subdomains, they found that the root-mean square error (RMSE) was 0.005972, which represents a high accuracy of the numerical method.
In general, the RBF and SPH methods have several similarities. Those are as follows: (1) they are simple to implement, (2) they are a truly meshless method, and (3) they were developed on the basis of strong form of partial differential equations. Moreover, the RBF method has a better accuracy than the SPH method at the expense of more computational effort [24]. The need to invert the coefficient matrix makes the RBF more expensive than SPH method. However, the efficiency of the RBF method can be improved by using the domain decomposition method.
In the present numerical works, the combination of RBF and DDM was implemented, because it is a very efficient method and high accuracy with low numerical diffusion as suggested by [15] and [25]. The DDM is used to solve the multi small subdomain problems instead of single large domain problem. Next, the Cahn–Hilliard equation is used to capture the evolution of the interface, because it is a thermodynamically consistent equation in order to follow the interface evolution as reported by [26] and [27].
Section snippets
Radial Basis Function (RBF) method
The main advantage of RBF approximation is the ability to work for high dimensions of the arbitrary geometries. RBF is a function whose value depends only on the distance from some center point. By using the distance functions, it can be easily implemented to reconstruct a plane or surface using scattered data in two- and three-dimensional (2-D and 3-D) or higher dimensional spaces. It was initially developed for multivariate data and function interpolation. Hence, the meshfree nature has
Problem definition
For the simulation of KHI, it is considered that the 2-D motion of two immiscible fluids of different densities and velocities separated by an interface, in which the heavier fluid is at the bottom. In the numerical analysis, the Cartesian coordinates were selected. For simplicity, the x-dimension of the computational domain L (0 x L) is chosen to be equal to the domain height H (L H). The fluids move counter-currently at the uniform velocity U as shown in Fig. 2. In the calculation,
Results and discussions
In the next section, the meshless method on the basis of the RBF in the framework of the fractional step method is used to solve the KHI of two incompressible and immiscible fluids. The instability is described by the modification of 2-D incompressible Navier–Stokes equations as shown in Eqs. (31)–(32). Here the Cahn–Hilliard equation is used for capturing the interface of working fluids.
In the present numerical study, the point independence tests were carried-out firstly in order to determine
Conclusions
The interface evolution of the KHI in a rectangular cavity and complex domains was investigated numerically by using RBF and domain decomposition method. The phenomena around KHI were modeled by the modification of Navier–Stokes equations. Here the convective Cahn–Hilliard equations were applied to capture the interface between the fluids. The numerical method on the basis of the implicit Euler scheme for temporal discretization and RBF method for spatial discretization were also utilized. The
CRediT authorship contribution statement
Eko Prasetya Budiana: Conceptualization, Methodology, Investigation, Programming, Writing - original draft, Formal analysis. Pranowo: Investigation, Data checking, Formal analysis. Indarto: Resources, Writing - review & editing, Funding acquisition, Supervision. Deendarlianto: Analysis validation, Writing - review & editing, Resources, Funding acquisition, Supervision.
Acknowledgments
The work was carried out and partially funded by the Indonesian Ministry of Research and Technology/National Agency for Research and Innovation , and Indonesian Ministry of Education and Culture, under World Class University Program managed by Institut Teknologi Bandung. Next it was funded partially also within research programs “World Class Researcher 2019” and the Ph.D. scholarship for Eko Prasetya Budiana (first author) from the Directorate General of Higher Education, Ministry of Research,
References (36)
- et al.
Mechanism of slug formation in downwardly inclined pipes
Int. J. Multiph. Flow
(2000) - et al.
Prediction of the initiation of slugs with linear stability theory
Int. J. Multiph. Flow
(1986) - et al.
Carbon-steel corrosion in multiphase slug flow and CO2
Corros. Sci.
(2006) Investigation and prediction of slug frequency in gas/liquid horizontal pipe flow
J. Pet. Sci. Eng.
(2009)- et al.
New slug control strategies, tuning rules and experimental results
J. Process Control
(2005) - et al.
Transition from stratified intermittent flows in small angle upflows
Int. J. Multiph. Flows
(2001) - et al.
Two-dimensional Kelvin–Helmholtz instabilities of multi-component fluids
Eur. J. Mech. B Fluids
(2015) - et al.
Overlapping domain decomposition method by radial basis functions
Appl. Numer. Math.
(2003) - et al.
A comparison of efficiency and error convergence of multiquadric collocation method and finite element method
Eng. Anal. Bound. Elem.
(2003) - et al.
Improved multiquadric approximation for partial differential equations
J. Eng. Anal. Bound. Elem.
(1996)