Investigation of the effects of pipe diameter of internal multiphase flow on pipe elbow vibration and resonance

Computational fluid dynamics modelling of internal two-phase flow induced transient forces at 90° elbows have been carried out to evaluate the effect of pipe diameter on the characteristics of multiphase flow induced vibration. Simulations of two-phase flows of slug, cap bubbly and churn induced vibration at a pipe elbow were carried out using the volume of fluid model for the two-phase flows and the k – ε model for turbulence. Modal analysis has been carried out to evaluate the risk of resonance. Results were compared across three geometrically similar pipes of different diameters. Simulation results showed that the behaviours of the flow induced forces at the pipe elbow as a function of gas velocity for internal diameters of 0.0525 and 0.2032 m are similar. However, the multiphase flow induced force characteristics are different in the 0.1016 m diameter (intermediate) pipe. It can be attributed to the transition behaviour of gas–liquid two-phase flows caused by Taylor instability in an intermediate sized pipe. The predicted root-mean-square flow induced forces as a function of Weber number were correlated with an existing empirical correlation for a wider range of pipe sizes and gas volume fractions between 40% and 80%. Furthermore, the pipe natural frequencies increase with the increase of gas volume fraction in smaller pipes and the resonance risk increases with the increase of pipe diameter.


Introduction
Internal multiphase flow induced vibration (MFIV) in pipe elbows poses significant fatigue failure risks in the piping network of oil and gas production, nuclear power plant and chemical process systems. 1,2 Current industry guideline has been developed considering single-phase flows, while the few existing MFIV models reported in literature are based on small-scale laboratory experiments which do not completely address the differing multiphase flow mechanisms between small and large pipes. In addition, literature shows that MFIV is two-phase flow regime dependent with the slug and churn flows having the most significant forcing functions. 1,3 Since flow regimes behaviours (especially slug and churn flow) and transition differ in small and large diameter pipes, the force correlations developed from experiments and modelling from small pipes can lead to costly flow assurance, pipeline integrity and safety problem if applied to large diameter pipes. Kataoka and Ishii 4 defined the critical diameter at which large slug bubble could no longer be sustained in a pipe as: where D H , σ, g and Δρ are the hydraulic diameter, surface tension, gravity and density difference between gas and liquid. Pipes having D * H ≥ 40 are termed large pipes in which cap bubbles rather than long slugs exist, and D * H < 18.6 are termed small pipes while 18.6 < D * H ≥ 40 represent intermediate size where transition behaviour involving the disintegration of long slugs into cap bubbles is present. Therefore, developing a database of dominant frequencies, magnitudes and RMS of force fluctuations due to two-phase flows in different sized pipes will be valuable to engineers at design stage of multiphase flow systems.
Experimental studies on MFIV in large pipes are few. The most relevant of such studies were conducted by Nennie and Belfroid 5 and Belfroid et al. 6 These experiments were also used to validate computational fluid School of Engineering, Robert Gordon University, Aberdeen, UK dynamics (CFD) studies that were conducted in a subsequent phase of the project. [7][8][9] These studies concluded that the fluctuations of liquid hold-up resemble the multiphase flow induced force fluctuations at pipe elbows. Belfroid et al. 10 summarised the results of all the test cases that were investigated by Nennie and Belfroid. 5 Their findings on the relationship between the nondimensional force fluctuation and Weber number compared fairly to the model by Riverin et al. 11 and data from Liu et al. 3 Most of the reported literature on FIV in pipe bends in large pipes are for nuclear power plant applications. 12 Yamano et al. 13 results prove that both pipe size and configuration have effects on MFIV characteristics. The cost and hazard of comprehensive life scale experiments of oil and gas and nuclear power plant pipework could be prohibitive. Hence recent FIV investigations were mainly focused on numerical methods. 14-16 Pontaza et al. 9 applied the volume of fluid (VOF) model for two-phase flows and large eddy simulation (LES) approach in their CFD investigation of MFIV and concluded that the standard deviation of flow induced forces agreed with the experimental measurements within ± 20%. A similar study was carried out by Emmerson et al. 7 using the CFD approach and their results showed that CFD was able to predict the magnitude of force fluctuation within 2-7%. Studies have applied U-RANS approach to capture the mechanism of FIV in large pipe bends due to pressure fluctuations, dominant in single-phase flows. 12,13,17 In general, the risk of FIV increases in multiphase flows. 18 Additionally, a complete assessment of the dynamic behaviour of a structure vibrating under loading includes the determination of its modal parameters. Kim and Srinil 19 investigated a full-scale M-shaped subsea jumper transporting slug flows and established the resonance risk because the predominant frequencies of the stress and displacement fluctuations were similar to the natural frequencies of the jumper. Finally, Wang et al. 20 observed in their numerical study that the first three orders of natural frequencies of a pipeline conveying gas-liquid two-phase flow, increased as the gas volume fraction at pipe inlet increased. Their findings also showed that high amplitude and frequency of pressure and displacement fluctuations observed within the slug flow regime could induce excessive stress and fatigue failure in the pipeline.
In addition, the effect of added mass due to internal multiphase flow is crucial in evaluating the natural frequencies of a pipe structure and should not be ignored in an MFIV assessment. However, this effect has not been quantified and compared for different pipe sizes so far. Therefore, the present study has implemented a validated CFD approach 21  (small) and the newly obtained data from intermediate and large pipes. Moreover, the effect of added mass due to internal multiphase flow on the natural frequencies of a pipe structure as a function of pipe diameter has been studied.

Volume of fluid method
This study implements the VOF method to model the flow of air and water by solving a single set of momentum equations. In this method, an equation representing the conservation of the volume fraction α and one continuity equation for both phases are also solved. The VOF method tracks the volume fraction of each phase through the computational domain. 22,23 The volume fraction continuity equation of one phase given in equation (2) is used to track the phases: where subscript q represents each phase component, u is velocity vector, t is time and α represents the volume fraction of one of the phases. The volume fraction of the other phase is obtained from the relation The continuity and momentum equations are solved alongside the volume fraction equation. They are given as, continuity: Momentum: where ρ, μ and F are density, dynamic viscosity and external body forces, respectively. The properties appearing in the transport equation are calculated by summing the contributions of all phases in the multiphase flow according to the volume fractions of the phases. For example, equation (6) accounts for the contribution of each phase to the density of the two-phase flow.
ρ m is mixture density. The term F in the momentum equation is used to account for the surface tension force between liquid water and air. F is given as: Where, σ is the surface tension coefficient, and κ 1 is the surface curvature of the liquid droplet defined in terms of the divergence of the unit normal n 1 . κ 1 is given by The local gradients in the surface normal to the interface are employed to account for n 1 : s represents surface and n 1 is unit normal.

Turbulence model
The mixture k − ϵ turbulence model has been implemented in the present study. Launder and Spalding 24 gave the mixture turbulent kinetic energy k and energy dissipation rate ϵ as equations (10) and (11), respectively: the mixture density and velocities are given by: and μ t,m is the turbulent viscosity and is calculated from: and the production of turbulent kinetic energy, G k,m is computed from: The turbulent model constants are: The standard k-ϵ URANS turbulence model requires a modest computing resource to produce a reasonably accurate solution compared to the rest models and approaches of calculating turbulence. It is well established and performs within reasonable accuracy for many industrial flows. The unsteady k-ϵ model cannot predict the turbulence pressure fluctuations, however, unlike single-phase flow, the effect of the impact forces on bends due to density difference of the alternating slugs and large interface waves, is more significant than the effect of turbulence fluctuations in two-phase slug/churn FIV.

Computational geometry
Three different pipe sizes presented in Table 1 were studied, the pipe sizes were categorised according to the critical diameter criteria. The pipe geometry is shown in Figure 1. Base case pipe geometry represented in Table 1 as small is taken from Liu et al. 3

Inlet and outlet boundary conditions
Based on the entrance length mechanism suggested by Taitel et al. 25 for slug-churn flow regime, with sufficiently long entrance length, multiphase flows introduced into a flow channel at individual phase velocities will eventually separate into the flow patterns corresponding to the superficial velocities. However, the concentric two inlets method described in Parsi et al. 26 has been utilised in the present study, where air and water are introduced into the pipe through two separate inlets. This method expedites the development of multiphase flow regimes. The central core in Figure 2 is used to introduce air by setting the air phase velocity while water was introduced through the outer annular region by setting the water phase velocity. Both are velocity inlets.
The phase velocities are calculated as:  and where v g , v l , v sg , v sl , A, A g and A l are specified gas inlet velocity, liquid inlet velocity, superficial gas velocity, superficial liquid velocity, cross-sectional area, gas area inlets and liquid area inlets, respectively. The turbulent intensity and hydraulic diameter given as I and D H , respectively, are the turbulence parameters which were set at the liquid and gas inlets. The turbulent intensity is calculated as: Re is Reynolds number and D H hydraulic diameter. Pressure outlet is specified at the geometry's single outlet. Since the outlet is at atmospheric condition, the gauge pressure is set as zero.

CFD solution methods
The choice of the time step size is conditioned by the Courant number-based stability criterion. Based on a preassigned maximum allowable Courant Number of 0.25 and the smallest time spent by fluids in cells (control volume) in the region near the free surface to empty out of the cell, the present numerical algorithm automatically refines the sub-time step for the volume fraction calculation. The sub-time step size is calculated as: where, C is the Courant number, V is the cell volume and U f are the outgoing fluxes in the cell. The explicit scheme is adopted for this VOF model formulation and the volume fraction cut-off and courant number are set as 1 × 10 −6 and 0.25, respectively. The water and air were assigned as the primary phase and secondary phase, respectively. The surface tension between water and air was set at a constant of 0.0728 N/m. In addition, the SIMPLE scheme was adopted for the pressure-velocity coupling, while PRESTO and second-order upwind were implemented for pressure and momentum discretization, respectively. The volume fraction equation was solved using the Geo-Reconstruct scheme and simulation time step was fixed at 0.00001 s.

Natural frequency model
The natural frequencies and mode shapes are important parameters in the design of a structure under dynamic loading. A modal analysis determines the free vibration characteristics by solving a equation of motion for an undamped system. 27 The equation is given in matrix form for the elements of the structure model ( Figure 3) as 28 : k, d m andd are stiffness, displacement, mass and acceleration, respectively. Assuming a linear system whose motion is harmonic under free vibration, the solution of Equation (19) above becomes an Eigenvalue problem of the form: Where ω i is the ith natural frequency of the pipe structure in radians per unit time and ϕ i is the displacement vector (eigenvector) representing the mode shape of the pipe structure at the ith natural frequency. M and K represent the global mass and stiffness matrices of the system. The multiphase flow in the pipe adds to the mass of the pipe depending on the global volumetric fractions of the gas and liquid flowing in the system. The natural frequency  of the pipe will change according to the added mass of each of the contained two-phase flow cases under study.
The mass M is given as: where The ANSYS Mechanical FEA tool employed in this study uses a special algorithm known as HSFLD242 3-D Hydrostatic Fluid elements to include contained fluid inside a solid. This model calculates the volume of the pipe inner channel, and the density is manually calculated and supplied as an input for each flow condition modelled. The density of fluid is given by: and M s is the mass of the pipe structure and M m represent the two-phase fluid mass, ρ m is the mixture density, V p is the volume of the pipe flow channel, β is the global volumetric gas fraction, ρ g and ρ f are the gas and liquid densities, respectively. Further details of the HSFLD242 3-D Hydrostatic Fluid Elements code can be found in ANSYS theory guide. The structural model geometry and properties are given in Table 2 and Supplemental  Table 3, respectively. The pipe thicknesses are selected according to the API 5L grade line pipes. The length, breadth and bend radius are as given in Table 1

Mesh sensitivity and validation studies
The present CFD approach was validated and reported in an earlier study by Hossain et al. 21 In the validation, void fraction and flow induced force signals obtained from CFD simulations were qualitatively and quantitively compared to the experimental report of Liu et al. 3 for the 0.0525 m pipe diameter. The simulation results conformed with the experimental results. Supplemental Table 4 shows  Table 1. Further mesh independency studies have been carried out for each of the pipe sizes. The pipe scaling was also adopted in reporting, the largest pipe scale is 1. Figure 5 shows a typical wall view of the meshes used in this study. Out of three different mesh sizes for each geometry, the meshes which were finally selected as a result of the mesh independence studies are 366912, 690688 and 269010 for the pipe diameters of 0.0525, 0.1016 and 0.2032 m, respectively, as shown in Figure 6(a)-(c). The geometries were created in ANSYS Design-Modeller while hexahedral meshes were created in ICEM CFD 18.0 meshing software. Supplemental Table 5 shows the refinement parameters that were adopted. The near wall cell size was predicted using the expression: and U τ is the friction or shear velocity and τ w is the wall shear stress. For internal pipe flows, friction coefficient, C f is given based on Blasius equation as: where Δy and U ∞ are the first cell height and the free stream velocity, respectively. y + values between 30 and 100 were chosen.
For mesh independence study, the slug flow condition was modelled as a homogenous mixture. Velocity profiles shown in Figure 7 Figure 8 and Supplemental Figure 9(b) represent the PDFs of the void fraction fluctuation extracted at a crosssection above 70D from the pipe inlet. Figure 8(b) shows that the two peaks in the PDFs of void fraction in the small pipe became almost equal while in the larger pipes the leading peak representing low void fraction of liquid  slug containing small bubbles was higher than the trailing peak which represents cap bubbles. These behaviours are consistent with reported experimental observations of slug and cap bubbly flow regimes characteristics. [30][31][32] The contour plot in Figure 8 shows that the liquid slug could still bridge the pipe cross-section. On the other hand, the contour plot in Supplemental Figure 9 shows that the liquid slugs could rarely bridge the pipe cross-section. The large liquid wave structures which are characteristic of churn flows are also visible in Supplemental Figure 9(a) and the observation is consistent with the findings in literature. [31][32][33][34] The PDFs in Supplemental Figure 9(b) all show peaks above 0.8 in the three pipe sizes. In summary, the CFD method was able to reproduce the basic physical behaviours of the slug, cap-bubbly and churn turbulent flow regimes expected in the pipes. Slug flow regimes which were observed in the vertical pipe section in small pipes transitioned into cap-bubbly flow in large pipes.

Two-phase flow induced force
The fluctuating forces acting on the elbow of the pipe were calculated using equations (29) to (31) based on a control volume around the elbow. The x is in the horizontal direction, while the y is in the vertical direction (see, Figure 1).

F x (t) =ṁ(t)V (t) + p(t)A
− at the exit plane of the bend (30) − at the inlet plane of the bend (31)    The momentum and pressure fluctuations arise from the alternating nature of liquid and gas phases in slug and churn flow through pipe bend control volume. The bend control volume analysis that was adopted in the present study to calculate force at the bend was based on momentum balance theory also implemented by Liu et al. 3 A control volume created at the bend starts at a plane positioned at the immediate upstream vicinity of the bend forming the inlet of the control volume. The control volume terminates at the immediate downstream region of the bend forming the control volume exit. Pressure and momentum are calculated on the inlet and outlet surfaces at each time step and the values depend on the fluid phase flowing through the surface at the time. Similar to the findings of Hossain et al. 21 for a small pipe of 0.0525 m I.D., Supplemental Figure 19(a) shows that the RMS of momentum flux fluctuation is significantly correlated to the RMS of force fluctuation in the x-direction for the intermediate and large pipes examined in the present study. This was not the case in the y-direction. Supplemental Figure 19

Effect of flow pattern on natural frequency
Supplemental Figures 24(a) to (c) present the natural frequencies of the three similar pipe sizes. The added mass effect due to contained fluid is captured by defining the global volumetric gas fraction at the pipe inlet. Generally, as the gas fraction increases corresponding to increasing gas superficial velocities (Vsg) represented by the increasing trend of Vsg/Vsl in Supplemental Figure 24, the natural frequencies of the first three modes of vibration increases, approaching the highest value corresponding to gas volumetric fractions (β) of 0. Based on their experiment, Wang et al. 20 also reported that increasing gas volume fraction at inlet had similar effect on the first three natural frequencies of an acrylic material pipeline riser of 0.0514 m I.D. and thickness of 0.0058 m conveying gas-liquid two-phase flow. However, the present study shows that the increase in natural frequencies with gas superficial velocities and volume fraction were most significant in the smallest pipe of ¼ scale than the ½ and 1.0 scale sized pipes. Frequencies as low as 1 Hz were recorded in the 1.0 scale pipe while the lowest frequency in the ½ scale pipe is ∼2 Hz.

Resonance risk assessment
Risk of resonance in the three pipe sizes, subjected to identical internal multiphase flow conditions were examined. The 1.0 scale pipe presented the highest resonance occurrence risk. Supplemental Figure 25(a) to (c) presents the comparison of the first three natural frequencies of the pipe bends with fixed supports (representing modes 1-3) as a function of gas superficial velocity and the dominant frequencies of force fluctuation. In Supplemental Figure 23(c), the dominant frequencies of force fluctuation were close to the first natural frequency of the pipe structure with and without contained fluid at gas superficial velocities of 1.7 and 2.765 m/s. Generally, the risks of resonance were most likely within slug/cap bubbly and churn turbulent flow regimes up to V sg of 5 m/s. In summary, although the risk of resonance in the ¼ scale (Supplemental Figure 23(a)) is insignificant, further studies might show frequencies of pipe deformation or displacement coinciding with natural frequencies or high frequency stress circles which could also reduce the fatigue life of the structure.

Conclusion
A numerical investigation to evaluate the effect of diameter on MFIV characteristics and the effect of internal two-phase flow on the resonance behaviour of three different sizes of geometrically similar pipe bends have been carried out. The significant findings in this study are summarised as follows: fluctuations with Riverin et al. 11 correlation was observed for the ¼ scale pipe where the constant C is 10. However, when the larger pipe sizes were considered, the correlation using C = 10 could not predict the non-dimensional RMS of forces well. A better prediction was obtained using C = 20. 3. In all cases, the presence of internal two-phase fluid changed the natural frequencies of pipe. This change was more significant for small pipes. Increasing the gas volume fraction increases the natural frequencies of the pipe, especially in the small pipe. 4. However, the risk of resonance was higher in the large pipe since the observed low natural frequencies were closely matched with the low peak frequencies of force fluctuations.

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.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Commonwealth Scholarship Commission in the UK for Nkemjika Mirian Asiegbu's PhD study. Greek letters α volume fraction of one phasẽ α A g area averaged void fraction ϵ energy dissipation rate θ w contact angle κ l surface curvature of the liquid droplet κ g surface curvature of the gas bubble κ Von Karman constant given as 0.4187 μ dynamic viscosity, kg/m s μ f dynamic viscosity of liquid, kg/ms ρ density, kg/m 3 ρ f liquid density, kg/m 3 ρ m mixture density, kg/m 3 ρ g gas density, kg/m 3 ρ l water density, kg/m 3 Δρ density difference between gas and liquid, kg/m 3 σ surface tension, N/m ∅ arbitrary flow property τ w wall shear stress, N/m 2