Computational fluid dynamics modelling of multiphase flows in double elbow geometries

This study investigates the effects of elbow on the transition and development of multiphase flow using computational fluid dynamics modelling techniques. The Eulerian - Multifluid VOF model coupled with an Interfacial Area Transport Equation has been employed to simulate air-water two-phase flow in a pipe with two standard 90 degree elbows mounted in series. Turbulence effects were accounted for by the RNG k-ε model. The effects of separation distance on two-phase flow development have been studied for initial slug and churn flow regimes. Computational fluid dynamics simulation results of phase distribution and time series of void fraction fluctuations were obtained and they showed good agreement with available experimental data. The results show that for initial slug flow regime, there is no flow regime transformation upstream and downstream of the two elbows. While at initial churn flow regime, flow regime transformation occurs at different sections of the flow domain before and after the two elbows. It was noticed that irrespective of the flow regime, the amplitudes and frequencies of void fraction fluctuation become smaller as the fluid flows along the pipe. Changes in the separation distance between the two elbows have larger effects on the flow at churn flow regime.


Introduction
Two-phase gas-liquid flow occurs in many industries such as chemical, petroleum, food and drink, processing and manufacturing. Two-phase flow is characterised by differences in the interface between phases leading to different two-phase flow regimes. Gas dominant flow regimes such as pseudo slug, churn and annular are the most commonly observed flows in industry and many experimental and numerical techniques have been applied to study the features of these types of flows. [1][2][3][4][5][6][7][8][9][10][11] Many of these studies have focused on simple geometries where fluids flow through a vertical or horizontal pipe or through a single elbow. However due to the availability of space and economics of design, there is a need for using complex piping systems such as multiple bends pipe, and U-bend to transport fluids. The knowledge of two-phase flow phenomenon in bends and its development downstream of bends is important in many industrial applications. In addition to oil-gas processing platforms where two-phase flow undergoes several changes of directions due to bends under limited space, Zhao et al. 2 have shown two specific examples of two-phase flow at multiple bends (i) fire reboiler with serpentine tubing and (ii) a novel combined bend/T-junction phase separator.
Two-phase flow regimes in vertical and horizontal pipes are quite different. In vertical pipes, annular, churn, slug and bubbly flows are observed, while dispersed bubble, slug, stratified and annular flows are observed in horizontal pipes. 3 Flow regimes can change rapidly at the bends where the flow changes direction from vertical to horizontal or horizontal to vertical orientation.
The behaviour of two-phase flow upstream and downstream a sharp return elbow was studied using capacitance sensors by Kerpel et al. 12 with an elbow of 1 mm radius of curvature and 0.008 m internal diameter. R134a refrigerant was employed while mass flux and vapour liquid void fraction were varied between 200-400 kg/m 2 s and 0 to 1 respectively. To determine the effects of elbow, eight capacitance sensors were positioned at different locations upstream and downstream of the elbow to extract a series of void fraction data. In their experiment, flow regimes before the bend were slug, intermittent and annular flows. It was observed that in the pipe section downstream of the bend, no bend effect is evident for slug flow but there were ripples and disturbances at the phase interface. In intermittent and annular flows on the other hand, no visually detectable effects on the flow behaviour due to the elbow were evident. It was also concluded that the disturbance and ripple effects due to the elbow stretches out to at least 21.5 tube diameters downstream of the elbow. Abdulkadir et al. 13 also reported observations from the study of two-phase air-water churn-annular flow in a large diameter vertical 180 degrees return elbow using electrical conductance technique. The elbow has a diameter of 0.127 m and a curvature ratio (R/D) of 3. Superficial velocities of air were ranged from 3.5 to 16.1 m/s and those for water from 0.02 to 0.2 m/s. The flow patterns were identified to be within the slugchurn transition region using Probability Density Function (PDF) profiles of the time series of the mean film fractions. The flow pattern upstream of the elbow was confirmed to be churn flow while in the downstream, the liquid was observed to drain to the bottom of the elbow due to gravity effects, the gas was at the centre of the pipe leading to annular flow pattern. Asgharpour et al. 1 also performed experiments to investigate churn/annular and pseudoslug flow before and after pipe elbows in a 76.2 mm pipe. Wire mesh sensor (WMS) was used to determine the void fraction distributions in the upward-vertical orientation before the elbow and the horizontal orientation after the elbow for 11 m/s and 0.1 m/s gas and liquid velocities. They observed that the flow transitioned from annular -churn in the vertical section to a wavy stratified flow in the horizontal section. Vieira et al. 14 investigated the effect of 90 standard elbow on horizontal gas-liquid stratified and annular flow characteristics using dual wire-mesh sensors. The horizontal test section was a 0.0762 m ID, 18 m long pipe and it generated stratified-wavy and annular flows. Two 16 Â 16 wire-mesh configuration sensors were positioned at 0.9 m upstream and 0.6 m downstream of the elbow, and the experiments were conducted at superficial liquid velocities of 0.03 m/s and 0.2 m/s while superficial gas velocities ranged from 9 m/s to 34 m/s. The cross-sectional averaged void fraction time series showed similar structures in the upstream and downstream locations for stratified and annular flows. It was also reported that wave instability is the most notable feature of the stratified to slug transition in both upstream and downstream locations of the elbow, the cross-sectional time averaged void fraction values slightly increased after the 90 horizontal elbow and in the downstream section more liquid is in the perimeter of the pipe and waves are smaller as compared with the upstream conditions Zhao et al. 2 has reported an experimental study of two-phase flow regime transition in a double-bend pipe line. Zhao et al. 2 reported experimental investigations concerning gas-liquid two-phase flows in two 90 degree bend in series with fixed separation distance. The flow direction was vertical to horizontal then horizontal to horizontal configuration in their study. Phase distribution within the elbows and downstream of the elbows were studied via capacitance measurements. Their results showed that flow transformation occurred in two phase flows due to secondary flows as well as gravity effects at the bend. PDF analysis showed clear evidence that stratified plug, wavy and slug flows were developed in the horizontal section after the bend, while slug and churn flows occurred in the upstream of the elbow when the gas and liquid superficial velocities were 0.3 to 4 m/s and 0.21 to 0.91 m/s respectively. Further analyses showed that due to the bend, for a liquid superficial velocity of 0.38 m/s and the gas superficial velocity of 0.05 m/s, upstream bubbly flow transformed into stratified flow downstream of the elbow while at a liquid superficial velocity of 0.38 m/s and gas superficial velocity of 0.71 m/s the upstream slug flow transformed into plug flow after the elbow with large gas bubbles separated by a liquid layer. As the two-phase flow passed around the bend to the horizontal sections, the liquid phase drained to the bottom of the pipe due to the effect of gravity and transitioned into stratified flow which remained dominant across the horizontal sections. It was also observed that the minimum distance necessary for the development of the flow after the 90 bend is between 10 to 50 D from the bend depending on the flow rates.
The above literature review shows that there exists a gap in knowledge in the understanding of two-phase flow transition under multiple bends with different separation distances. These pipe configurations are very common in offshore oil and gas processing platforms, refineries, food processing plant and can be used in designing innovative bend phase separator. This paper aims to provide a fundamental understanding of two-phase flow transitions in multiple pipe bends.
The transient nature of gas-liquid two-phase flows still remains a key challenge in multiphase flow studies both experimentally and numerically and CFD can provide vital flow information and interactions between flow phases. Volume of Fluid (VOF) model is a widely used approach for simulating multiphase flows 4,10,[15][16][17][18][19] where there is distinct interface between phases such as stratified or slug flow regimes. In the present study, a Multifluid VOF model which is a Eulerian-Eulerian model with an added interface sharpening model from volume of fluid, has been employed in order to investigate the development of multiphase flows before and after the elbows in double elbow (vertical upward -horizontal -vertical downward) geometries for various separation distances.

Numerical modelling
Two-fluid Eulerian-Eulerian two-phase model has been utilised in the present simulations. In order to capture the distinct interface between liquid and gas phases especially in the slug and churn flow regime an additional interface capturing method has been used which is known as Multifluid-VOF model as known in Ansys Fluent. To account for different bubble sizes and their bubble breakage and coalescence within the flow without resolving to details of bubble size distribution such as in population balance model, a simplified modelling assumption of an Interfacial Area Concentration (IAC) transport equation was incorporated in this hybrid model.

Eulerian-Eulerian two-phase flow modelling
Continuity equation. The continuity equation for the phases are Where a i ; q i andṽ i are the volume fraction, density and velocity of the individual phases.

Momentum equation. The momentum equation for the phases are
Where P and s i are the pressure and stress-strain tensor of the individual phases,g is acceleration due to gravity andF ij is the interfacial force between the phases.
The general form of the interfacial forceF ij from equation (3) is expressed as where, K ij is the interphase momentum exchange coefficient between the phases. In the Eulerian-Eulerian two-phase flows, the second phase is treated as bubbles and the interphase exchange coefficient is given by where, X p is the interfacial area concentration, f is the drag co-efficient and s i is the particle relaxation time given by Schiller and Naumann 20 model is used to calculate drag coefficient f: Interfacial area concentration. The solution of interfacial area concentration allows the inclusion of bubble diameter distribution and coalescence/breakage effects within limited computational resources.
The transport equation for interfacial area concentration, IAC, is given as where X p is the interfacial area concentration and a g is the gas volume fraction. The first two terms on the right hand are the gas bubble expansion due to compressibility and mass transfer, m g is the mass transfer rate into the gas phase per unit mixture volume. S RC and S WE are coalescence sink terms due to random collision and wake entrainment respectively. S TI is the breakage source term due to turbulent impact. The source and sink terms are accounted for with the Hibiki-Ishii model. 21 The interfacial area concentration is coupled with the momentum equation via the interfacial force(s) in equation (5).
Multi-Fluid VOF. In slug and churn flow situations, there is a distinct and sharp interface between the liquid and gas phase. This sharp interface has been captured in the present study using a combination of volume of fluid method with Eulerian-Eulerian model explained earlier. At the interface of two phases based on the value of calculated volume fraction, a special interpolation treatment is applied to the cells that lie near the interface between the two phases to capture the shape interface using the High Resolution Interface Capturing (HRIC) scheme. 22

Turbulence modelling
In the present study, turbulence is represented by mixture turbulence concept where turbulence is calculated using fluid mixture properties and was accounted for by using the RNG k-e model. 23 The RNG k-e transport equations are: Kinetic Energy, k G k is the generation of turbulence kinetic energy due to the mean velocity gradients. a k and a Ɛ are the inverse effective Prandtl numbers for k and Ɛ respectively. C 1Ɛ and C 2Ɛ are 1.42 and 1.68, l eff is the effective viscosity. The effective viscosity is given by Where, C l ¼ 0:085. RNG k-e is more accurate and reliable for a wider class of flows. They are also more responsive to the effects of rapid strain and streamline curvature and thus is suitable for flow through bends 23 and where, S is strain rate of fluid flow.

Solver control/numerical schemes
Ansys Fluent version 19.1 has been used for providing the numerical solutions of the governing equations of mass, momentum, phase volume fraction and turbulence. In this study, three dimensional (3-D) transient simulations were performed with the assumption that the liquid and gas phases were incompressible with no mass transfer between them. Momentum and turbulence equations are discretised using the upwind scheme and pressure and velocity equations are coupled with phase-coupled SIMPLE algorithm. Discretised algebraic equations use the under-relaxation factor of 0.3 for pressure, momentum and interfacial area concentration equations; 0.5 for volume fraction and 0.6 for turbulence kinetic energy and energy dissipation rate. The governing equations are solved in transient scheme with a time step of 0.001 s. Solutions are deemed converged when the residuals of all equations reach below 10 À6 . It took between 10 to 20 iterations to reach convergence in each time step.

Flow domain
The single elbow 3-D geometry in Figure 1(a) was modified to accommodate a second elbow as shown in Figure 1(b), and the normalized separation distance (L/D) between the two elbows mounted in series are 0, 5, 10, 15 and 20. The 3-D computational geometry in Figure 1(a) is similar to that of Parsi et al., 24 it consists of 3 m vertical and 1.9 m horizontal pipes, upstream and downstream of a standard 90degree elbow respectively and the flow of fluids is from upward vertical to horizontal. Pipe diameter and the radius of curvature of elbow are 0.0762 m and 1.5 respectively.

Boundary conditions
Velocity-inlet and pressure-outlet boundary conditions have been set at the inlet and outlet respectively (see Figure 2). Standard log-law of the wall was used at the wall. Turbulence intensity of 5% was set at the inlet. To aid faster flow development, the pipe inlet surface was split into two as shown in Figure 3 similar to CFD simulation of Parsi et al. 25 Gas was introduced into the domain via the middle of the inlet (red patch) while the liquid was introduced circumstantially (blue patch). The liquid and gas phases were introduced into the domain based on the velocities obtained from equations (13) and (14), respectively. The whole domain was initially filled with the liquid phase at zero velocity.
Where, V sg and V sl are superficial gas and liquid velocities, A p is the area of the pipe and A g and A l are inlet areas of gas and liquid.

Mesh/grid generation
A structured grid was generated across the flow domain as shown in Figure 4. Three different grids were considered in conducting a mesh sensitivity analysis in this study comprises of 80000 (mesh 1), 105000 (mesh 2) and 181608 (mesh 3) cells. Figure 5 shows the cross-sectional slices of the three grids in the perpendicular direction of the flow.

Results and discussion
Mesh Sensitivity study  Figure 7, time averaged void fraction and its standard deviation have been compared with the experimental data of Parsi et al. 24 Table 1 shows the mean void fraction and standard deviation of the times series of cross-sectional averaged void fraction. Meshes 1 and 2 produced similar averages while Mesh 3 produced an average higher than the experimental data. The standard deviation however shows that Mesh 2 produced fluctuation data that are much closer to the experimental data of Parsi et al. 24 than Meshes 1 and 3. Mesh 2 was therefore selected for all other simulations.

Single elbow validation
The predicted area-averaged volume fraction at 1 m upstream of the bend has been compared with both previously reported experimental and simulation data of Parsi et al. 24  Furthermore, Table 2 shows the mean and standard deviation of the time series of cross-sectional averaged void fraction and its standard deviation of the present study compared with available numerical and experimental data. It is evident from Table 2

Multiphase flows in double elbow geometries
Two different gas velocities in the slug/churn multiphase flow region were studied at a constant liquid velocity. Superficial gas velocities were 0.9 and 10.3 m/s while superficial liquid velocity was kept constant at 0.3 m/s. The normalised separation distance (L/D) between the upstream and downstream elbows (Elbows 1 and 2) was varied from 0 to 20 at intervals of 5. Supplementary Figure 9 shows the position of four (4) monitoring surfaces for which   the data considered in this study were extracted. For L/D of 0, Locations 1 and 2 are not available.
Flow visualisation. Supplementary Figure 10(a) and (b) show the contour plots of the cases 1 (V sg ¼ 10.3 m/s; V sl ¼ 0.3 m/s) and 2 (V sg ¼ 0.9 m/s; V sl ¼ 0.3 m/s) studied after 10 s of simulation run time when the normalized separation distance is 0, the figures show that the flow patterns in the vertical upstream pipe before Elbow 1 are characterised by churn and slug flows respectively. Churn flow has the attributes of both slug and annular flows, however, there are no clear boundaries between the liquid and gas phases but there is the presence of waves and discontinuous gas cores. Some authors consider this type of flow to be fundamentally annular flow in nature with large disturbance waves carried by the gas flow 28,29 and just as in this study, huge waves in churn flow was also observed by Parsi et al. 4,24 and Riva and Col. 10 Unlike the churn flow, there are clear and defined boundaries between the liquid and gas phases in Case 2 (Supplementary Figure 10(b)). The small   bubbles between the two Taylor bubbles in Supplementary Figure 10(b) were also observed by previous researchers. [28][29][30] As the normalised separation distance increases, the multiphase flow in the horizontal section after Elbow 1 (upstream elbow) was analysed to study the flow development and transition towards and after Elbow 2 (downstream elbow). Supplementary  Figure 11(a) shows the contour plots of the averaged void fraction across a symmetric plane of the flow at 10 s when superficial gas velocity (V sg ) is 10.3 m/s and superficial liquid velocity (V sl ) is 0.3 m/s. It is observed that in the churn flow, after the first elbow, gas and liquid separates out into wavy stratified flow. As the normalized separation distance increases, a slight but notable difference is observed in the flow development in the horizontal section (between the elbows), when L/D is 0 it is observed that the liquid phase flows through the elbows more as a thin film pushed towards the extrados of Elbow 2 (downstream elbow), this is due to the effects of secondary flow and centrifugal force in the elbow. When the normalized separation distance (L/D) increased to 5, the presence of liquid phase becomes more pronounced in the horizontal section between Elbows 1 and 2 and the effects of secondary flow and centrifugal forces diminished as the separation distance between elbows increased. At this stage, gravity effects on the flow become more significant with the fluid with the higher density (liquid phase) flowing at the bottom of the pipe and the lighter gas phase moving at the top of the liquid phase. At the normalised separation distance of 10, 15 and 20, after an equivalent normalised distance of about L/D of 15, the liquid phase is observed to be settled and flowing at an almost uniform level at the bottom of the horizontal pipe section. After Elbow 2 (downstream elbow), secondary flow effects on the flow become significant again, the liquid and gas phases are seen to flow at either end of the downstream vertical pipe section. Supplementary Figure 11(b) shows instantaneous flow development across the 4 monitoring surfaces in the flow domains across all L/Ds, it can be seen that in the vertical section before Elbow 1, the flow behaves more like annular flow as stated with the presence of a gas core while at the outlet of upstream elbow (Elbow 1) as the flow navigates into the horizontal section, the phases separate into a stratified flow due to gravity effects. As shown in Figure 6, the flow regime in horizontal pipe should be slug, the secondary flow created at the elbow enhanced flow separation leading to stratified flow in the horizontal section after the bend.
Supplementary Figure 12(a) shows the contour plots of the averaged void fraction of the flow in the symmetrical plane across the flow domains at 10 s for the case of gas superficial velocity of 0.9 m/s and liquid superficial velocity of 0.3 m/s. Although the contour plots show the flow to be stratified in the horizontal section, however, the liquid content at different sections of the pipe occupies more than 50% of the flow domain and in some sections the whole domain, thereby separating the gas phase into pockets. These features are indicative of slug flows in horizontal pipes where the liquid slug at the bottom of the pipe bridges the gas phase at different sections to form gas pockets known as the Taylor Bubbles. An increase in the normalized separation distance aided better development and separation of the flow, it is seen that after a separation distance (L/D) of 5, the presence of the liquid slug and elongated Taylor bubble becomes obvious. Supplementary  Figure 12(b) shows the instantaneous flow development across the four (4) monitoring surfaces in the flow domain across all L/Ds, the gas core is indicative of the Taylor bubble passing through at the particular instance on the monitoring surface before Elbow 1.  Quantitative analysis of two-phase flow behaviours due to double bend. The two-phase flow development due to double bends has been analysed through area-average void fraction along the length of the pipe as well as probability density function and power spectral density function. Statistical analysis of the behaviour of void fraction is one method that can also be used to further ascertain flow patterns in double bend pipe arrangement. According to Costigan and Whalley 28 and Lowe and Rezkallah 29 every two-phase flow pattern exhibit a specific probability density function (PDF) of void fraction. Churn flow is characterized by a PDF curve with a single broad peak at high void fractions and a broad tail at the low void fractions. The single peak at high void fraction is indicative of the flow's proximity to annular flow and the broad tail represents the passage of unstable slugs. According to Lowe and Rezkallah, 29 a typical PDF curve of churn flow is between an average void fraction of 0.7 and 0.9. In slug flow on the other hand, the PDF of time series averaged void fraction has two peaks, one at the low void fractions and another at the higher ones. These peaks represent the periodic passage of the two specific features of slug flows; the Taylor bubble at high void fractions and liquid slug at low void fractions. In this study, the Ksdensity function of MATLAB was used to generate the Probability Density Function (PDF) profiles of the times series of average void fractions. Power spectral density (PSD) of the time domain signals of the crosssectional averaged void fractions is another important parameter to ascertain two-phase flow characteristics. Power Spectral Density (PSD) shows the strength of the variations as a function of frequency. PSD is used to determine the dominant frequencies and amplitudes of signals in a time series data. Averaged spectral coefficients that are independent of time are produced with the Fourier transform in the ANSYS Post-Processing Software, CFD-Post, this is useful to identify the dominant frequency in a signal in various flow regimes. The time domain signals of the crosssectional averaged void fractions are converted into a frequency domain from which magnitude and PSD of the dominant frequencies are identified. When the PSD displays more than one peak, the frequency with the highest peak is said to be the dominant frequency. 31 PSD has also been used to identify different flow regimes by researchers such as Liu et al., 26 Hanafizadeh et al., 31 Ye and Guo, 32 Franca et al., 33 and Bouyahiaoui et al.. 34 Supplementary Figure 13 shows the time series of cross-sectional averaged void fraction, PDF and PSD when L/D is 0 and flow condition is Case 1 (V sg ¼ 10.3 m/s and V sl ¼ 0.3 m/s). Due to the absence of a separation distance when L/D is 0, the data shown are for the monitoring surfaces at Locations 1 and 4 only. Supplementary Figure 13(a) shows the time series of cross -sectional averaged void fraction at Locations 1 and 4, The amplitude of void fraction fluctuations diminishes at Location 4 (after the downstream elbow), compared to that at Location 1. In Supplementary Figure 13(b), the flow is further characterised by a PDF with a broad tail and single peak at Location 1 while at Location 4 the tail has become narrower compared to the former, this indicates the transition of the flow from churn flow before Elbow 1 to wavy-annular flow after Elbow 2. Furthermore, despite the transition in flow pattern, the PSD at both locations shows that the flow is dominated by broad band fluctuations (see Supplementary  Figure 13  Supplementary Figure 18 shows the time series of cross-sectional averaged void fraction, PDF and PSD when L/D is 0 and flow condition is that of Case 2 (V sg ¼ 0.9 m/s and V sl ¼ 0.3 m/s). In Supplementary  Figure 18 Supplementary Figure 21(a) shows that the twophase flows remain slug from throughout the flow domain. PDF plots in the supplementary Supplementary Figure 21(b) shows that the flow is dominated by two peaks indicating slug flows however as the separation distance between the elbows increases, the liquid body becomes more aerated with the lower peak shifting towards the higher void fraction. PSD profiles in Supplementary Figure 21(c) shows that there is a slight drop in the dominant frequency in the horizontal section between the elbows compared to the upstream and downstream vertical sections.
In Supplementary Figure 22(a) it is seen that when the separation distance increased to 20, the flow pattern along the pipe length remains slug flow. In Supplementary Figure 22(b), the PDF plots show typical representation of slug flow with double peak in the profile, but with various level of aeration in liquid body. At this separation distance, a slight drop in dominant frequency is observed as the fluid flow along the pipe as shown in Supplementary  Figure 22(c).
In summary, the initial churn flows in vertical pipe transition into stratified wavy flows, while the initial slug flows remain slug flows in horizontal section after the first elbow. The increase of separation distances allows separated stratified flow to develop with the reduction of frequencies of void fraction fluctuations in stratified flow, while in slug flow regimes increased aeration was observed in the liquid with the increase of separation and a reduction of predominant frequency of slug flows.

Conclusions
Computational Fluid Dynamics (CFD) has been used to simulate air-water two-phase flow with the use of the Eulerian-Multifluid VOF with RNG k À turbulence model in a pipe with two standard 90 degree elbows mounted in series in upward verticalhorizontal-downward vertical configuration. The non-dimensional separation distance between the elbows was varied between 0 and 20 to ascertain the effects of bend on the two-phase flow transition and development under churn and slug flow regimes.
Simulated results show that initial churn flow in upward vertical section transitioned into a wavy stratified flow in the horizontal section between the two elbows, and wavy annular flow in the vertical pipe after the second elbow. While the flow pattern remained the same in slug flow after the first elbow, the Taylor bubbles are more elongated in the horizontal section between the first and second elbows, and they drop in size after the second elbow as the normalised separation distance (L/D) between the elbows increases. On the other hand, the initial slug flow in upward vertical section remains slug flow in the horizontal and downward vertical sections after the bend, though the liquid body is more aerated as the separation distance increases. These findings have practical implications in designing boiling/condensing heat exchanger, novel combined bend/T-junction separator and oil and gas pipelines in offshore processing platform.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.