EULERIAN-EULERIAN MODELING OF DISPERSE TWO-PHASE FLOW IN A GAS-LIQUID CYLINDRICAL CYCLONE Miguel A. Reyes-Gutiérrez Departamento de Termodinámica Universidad Simón Bolívar José Colmenares Gerencia de Exploración y Producción PDVSA-INTEVEP Luis R. Rojas-Solórzano Departamento de Conversión de Energía Universidad Simón Bolívar Juan C. Marín-Moreno CEMFA Universidad Simón Bolívar ABSTRACT This work presents a three-dimensional CFD study of a two-phase flow field in a Gas-Liquid Cylindrical Cyclone (GLCC) using CFX4.3TM, a commercial code based on the finite volume method. The numerical analysis was made for air-water mixtures at near atmospheric conditions, while both liquid and gas flow rates were changed. The two-phase flow behavior is modeled using an Eulerian-Eulerian approach, considering both phases as an interpenetrating continuum. This method computes the inter-phase phenomena by including a source term in the momentum equation to consider the drag between the liquid and gas phases. The gas phase is modeled as a bimodal bubble size distribution to allow for the presence of free- and entrapment gas, simultaneously. The results (interface vortex shape and liquid angular velocity) show a reasonable match with experimental data. The CFD technique here proposed, demonstrates to satisfactorily reproduce angular velocities of the phases and their spatial distribution inside the GLCC. The experiments showed gas volume fractions smaller than 150 ppm along the liquid exit, however it was not possible to reproduce numerically these small volume fractions, since they had approximately the same order of magnitude of the error in the numerical method. INTRODUCTION During the last two decades, large efforts have been dedicated to research aiming to the development of multiphase flow compact separators. Today, compact separators are widely used in the industry and, particularly, in the oil production process, where expensive, heavy and bulky static separators have been traditionally used. The Gas-Liquid Cylindrical Cyclone (GLCC), patented by the University of Tulsa in 1994, is one of these compact devices, developed to separate gasliquid streams. The GLCC consists on a tangential pipeline connected to a vertical cylindrical body. The incoming multiphase stream runs into the cylindrical body through the tangential pipeline, which is also inclined 30˚ downward. The gas-liquid mixture enters to the separator and a first separation Antonio J. Meléndez-Ramírez Departamento de Termodinámica Universidad Simón Bolívar stage may occur; if the incoming flow velocity is too low, the gravity dominates the inertia and the mixture falls down, promoting the stratification or static-like separation of phases. While, if the flow enters the GLCC body with a mid-high velocity, then a swirling motion promotes an inertia-dominated initial separation process. The centripetal/centrifugal and buoyancy forces drive the gas toward the body centerline and upward to the top, while impelling the liquid toward the wall and bottom of the cylindrical body (see Fig. 1). A nozzle enhances the swirling motion of the mixture at the entrance. However, the excessive inflow velocity may decrease the separator performance, since it also promotes the turbulent mixing of the phases. The GLCC efficiency is determined by the quality of the outlet streams, i.e., if only gas flows through the top exit and only liquid flows through the bottom exit, then the separator efficiency is 100%. The gas carry under (GCU) and liquid carry over (LCO) are also practical indexes to compute the GLCC efficiency of separation for each exit stream. During regular operation, inside the GLCC it is established a bubble, slug, churn or annular flow pattern above the injection section, while a vortex sets in underneath. The gas phase accommodates in the vortex core, whereas the liquid phase occupies its periphery [1]. At the present time, it is estimated that there are around 600 GLCC units installed worldwide. Nevertheless, it is not still widely used because of the lack of models to accurately predict its performance. Current models are hardly close yet to predict the complex transport phenomena occurring inside the separator. 1 Copyright © 2004 by ASME Figure 1: Gas-Liquid Cylindrical Cyclone Recent research aiming to predict the GLCC performance has been based on CFD modeling, since this approach allows a deeper understanding of the detailed hydrodynamics inside the separator. Numerical results have leaded to the development of mechanistic models. The accurate prediction of the vortex inside the separator, which affects the GCU, is the main focus of most of the precedent research. In 1997, Erdal et al. [2] performed CFD simulations, using CFXTM, to reproduce the tangential velocity measured in a reduced-scale GLCC by Farchi [3]. In their work, Erdal et al. [2], performed single- and twophase flow simulations using the k-ε turbulence model. Singlephase simulations were carried out to compare the flow fields computed from the almost inexpensive axis-symmetrical simulation and the memory/time consuming three-dimensional simulation. From those simulations, they found insignificant differences and therefore performed the two-phase simulations using an axis-symmetrical geometry. However, when comparing to experiments, Erdal et al. [2] found similar tendencies, but over-prediction of the decay of the swirling strength. More recently, Motta et al. [4] compared the prediction capabilities of two CFD codes: CFXTM and TUCFDTM. Both codes are based on the finite volume method, but CFXTM uses SIMPLEC, while TU-CFD uses SIMPLER, as their pressure-computation algorithm. Furthermore, CFXTM results were obtained using the standard k-ε and the multi-fluid models, while TU-CFD computations employed the zero equation/mixing length and drift-flux models, for turbulence and inter-phase interaction, respectively. For the comparison, Motta et al. [4] used an axis-symmetric model to reproduce experimental data from Farchi [3]. Both numerical codes performed quite similar, but CFXTM results were closer to experimental data. Erdal et al. [5] compared multiphase models available in CFXTM (homogeneous and multi-fluid) using an axis-symmetrical model of a laboratory-scale GLCC. The k-ε model was employed in that study. Numerical results were compared with visualization experiments. The authors concluded that both models might predict the free surface shape with accuracy, noticing, in general, a good agreement between experimental and numerical results. Erdal et al. [6] studied the trajectories of small bubbles underneath the GLCC injection section and the influence of them on the GCU. Single-phase, two-phase and bubble trajectory simulations were carried out by them using an axissymmetrical computational model of a laboratory-scale GLCC. In their work, the numerical results were validated with visualization experiments. Results showed that the bubble trajectories and GCU are significantly affected by turbulence dispersion. Finally, Erdal et al. [7] carried out local velocity measurements in laboratory model of the GLCC using Laser Doppler Velocimetry (LDV). The tests were conducted with single-phase liquid flow only. Experimental data were compared with flow fields obtained from three-dimensional simulations. Additionally, the performance of two turbulence models was studied. Numerical results qualitatively agreed with experimental data, but failed to accurately predict the vortex wavelength and velocity profiles. As it is shown, most of the previous numerical works on the GLCC have been limited to simplified models to tackle the high computational cost of a full 3-D simulation. Nevertheless, the complexity of the transport phenomena inside the GLCC limits the obtained results from properly predicting the efficiency of this device. The main objective of the present study is to further the CFD analyses of the GLCC by performing a more accurate multiphase and three-dimensional simulation of the flow field using the finite volume method including the multi-fluid model and the Eulerian-Eulerian approach for continuous-disperse two-phase flow. Two-phase flow experiments on a laboratory-scale GLCC are also performed and the data is compared with computational results. FLOW FIELD SIMULATION Single-Phase Flow Simulation The flow field simulations were performed using CFX 4.3TM. This code is based on the finite volume method and uses SIMPLEC as pressure-coupling algorithm. The k-ε model has been used to compute the turbulent transport. Figure 2 depicts the main dimensions of the separator under study. The mesh refinement analysis was performed by comparison of the angular velocity along a vertical line located at 3/8D from the central axis. Meshes of 37000, 63000 and 122000 elements were compared. A mean difference of 1.2% and a maximum difference of 1.9% were encountered between the 63000- and 122000-element mesh results, normalized by the inlet angular velocity. According to those results, the intermediate mesh with 63000 elements was adopted for the rest of the study. Details of the selected mesh are displayed on Fig. 3. After the mesh sensitivity analysis was carried out, the computed results, in terms of angular velocities, were compared with experimental data. Simulations were performed with turbulence k-ε and Low Reynolds k-ε models. No significant differences in the angular velocities along the axis of the separator were found between results from both models. Higher order upwind differencing scheme was used for the convective term in the momentum equation. Further details about the numerical scheme may be found in [8] y [9]. Lower order difference schemes, such as 2 Copyright © 2004 by ASME upwind and hybrid, over-predicted the vortex dissipation along the separator body. The experimental angular velocity of the gas-liquid mixture was measured using a paddle-wheel meter at several locations below the nozzle section. Reported data corresponds to an arithmetic mean calculated over a large enough period of time. The computed numerical angular velocity was determined by averaging the local angular velocities at 4 coplanar points, shown in Fig. 4, with respect to the separator axis. Despite the vortex axis does not lie along the separator axis, the paddlewheel center of rotation does, making this methodology suitable for the calculation of the angular velocity. The comparison between computational and experimental mean angular velocities is shown in Fig. 5, demonstrating a reasonable agreement. However, the computed swirling decay appears slightly deeper than in the experiments. Indeed, the average relative error in Fig. 5 is 13%, while the maximum relative error is 27%. Figure 4: Coplanar sampling points for determining the computational mean angular velocity Angular velocity vs (z) position 30 Q1 Sim angular velocity (rps) Q2 Sim Q3 Sim Q1 exp 20 Q2 exp Q3 exp 10 0 0 Figure 2: Dimensions of the GLCC Model 0.2 0.4 0.6 0.8 1 1.2 z position from botton (m) Figure 5: Comparison of computational and experimental mean angular velocity for single-phase flow. Q1=0.00078 m3/s, Q2=0.00157 m3/s and Q3=0.0023 m3/s (a (b Figure 6 shows a three-dimensional view of the computed zero axial-velocity surface, also named the capturing surface, for the liquid single-phase simulation. Inside the volume bounded by this surface, the axial velocity is upward, while outside, it is downward. Mechanistic models assume that any bubble capable to reach this zone can be separated. This result demonstrates that there is an important difference between the three-dimensional performance of the separator and the axissymmetric assumption used in mechanistic models. Figure 3: (a) Surface-mesh at separator mid-section highlighting element agglomeration around the injection point; (b) Typical transversal section of volume-mesh in separator 3 Copyright © 2004 by ASME        r   U       T       r      r S t           (4) The subindexes  and ß denote the different phases present in the domain. The inter-phase momentum transfer is only possible between the disperse phases and the continuous phase; only one continuous phase may be defined. The disperse phases might be present in form of bubbles, drops or spherical solids. In general, the coefficient of inter-phase momentum transport is defined as follows[8]: d   c Figure 6: Computed capturing surface for single-phase flow. Q3=0.00078 m3/s Two-Phase Flow Model Two-phase flow simulations were carried out using CFX 4.3TM multi-fluid model. This model considers each phase as an interpenetrating continuum. This means that, each phase is present in each control volume. The volume fraction of a phase represents the fraction of the control volume that is occupied by that phase. Two types of phase may be defined: continuous phase and disperse phase. There is one solution field for each phase separately. Transported quantities interact by means of inter-phase transfer terms. The multi-fluid model is implemented using the Inter-Phase Slip Algorithm (IPSA) of Spalding. For each phase the continuity (1), momentum (2), k (3) and ε (4) transport equations are solved [8]:   r   U    0    d  U   U   r B  p    c   k       24  5.48 Re 0.573  0.36 Re (6) To calculate the turbulent viscosity and the source terms of the turbulence parameters equations, the following customary expressions are used [8]: T  C  k2  Sk  P  G     S   1        r   U  k     T  k     r  k   r S k t CD  (2) NP (5) Ihme correlation is used for the calculation of the bubble drag coefficient (eqn. 6) along with Gidaspow volume fraction correction (power coefficient of 1.6). This correction prevents viscous drag misestimating at control volumes where disperse fraction is closer to 1 [8]. (1)  r U  T    r  U   U    U   U   t 3 CD r  U   U 4 d  k C1 P  C3 maxG ,0  C2   (7) (8) (9) The equations of continuity, momentum and turbulent properties are solved in iterative way through the Semi-Implicit Method for Pressure Linked Equation Corrected (SIMPLEC) of Patankar and Spalding [8]. (3) Two-Phase Flow Simulation The experiments were carried out with air-water mixtures. The mean pressure inside the separator was held at 17.7 psia during the tests. In the simulation, water was assumed to be the continuum phase and two bubble diameters were considered for air, as the disperse phase. The bubble diameter of the disperse phase corresponding to free air was 500 m; while, the bubble diameter of the disperse phase corresponding to entrapment air was 150 m. The free air bubble diameter was chosen to be as small as possible yet to be totally separated (i.e., not carried under) at all the simulated conditions. At the same time, the 4 Copyright © 2004 by ASME entrapment air bubble diameter was chosen as the equilibrium diameter for the smallest flow rate explored. The equilibrium diameter is obtained by balancing buoyant and drag forces. Boundary Conditions Inlet: A homogeneous inflow tree-phase mixture (free gas, entrapment gas and liquid) was considered. The liquid volume fraction (rl) was taken according to the values for each experiment (rl=ql/qt). The entrapment gas volume fraction (reg) was taken equal to the liquid volume fraction (rl), under the assumption that it must be the largest gas volume fraction possible within the liquid stream once the free gas has been separated. The free-gas volume fraction ‘rfg’ was obtained as ‘rg- reg’. The free-gas bubble diameter was chosen relatively larger than the entrapment gas bubble diameter in order to ensure its buoyant release from the liquid continuous phase. Liquid flow rates (ql) of: 0.00073, 0.0011 y 0.00147 m3/s; and gas flow rates (qg) of: 0.0023, 0.0047 y 0.007 sm3/s were prescribed in the simulations. Gas exit: A reference pressure of 0 Pa was set. Liquid exit: During the experiments, the liquid level inside the separator was maintained at 1.15 m above the separator bottom regulating control valves (i.e., 0.1 m under the nozzle). For the liquid exit, an appropriate pressure difference, with respect to the gas exit, was used in order to guarantee the proper location of the vortex free surface, thus avoiding the necessity of adjustments in flow rate and phase volume fractions at these outlets. The results show a maximum deviation around 17.6% between simulations and experiments with an average of 11% over the total sample. Figure 7: Comparison of simulated and experimental free surface for qg = 0.0047 sm3/s. Left: ql = 0.00073 m3/s; Right: ql = 0.0011 m3/s Results Figures 7 and 8 illustrate the vortex free surface for several liquid flow rates. Pictures of experiments shown in Figs. 7 and 8, correspond to the average vortex size between the minimum and maximum sizes observed in each experiment. The minimum and maximum vortex sizes for every experiment were obtained from 15-second movies, each of which depicted approximately between 15 and 35 vortex-length fluctuations, depending upon the flow conditions. The typical difference between the minimum and maximum vortex length in an experiment was about 1.5 times the separator diameter. Tables 1 and 2 show the results obtained for the flow angular velocity at 0.35 m above the injection point at different gas and liquid flow rates. Table 1: Mean angular velocity for different gas flow rates and ql = 0.0011 m3/s Angular velocity (rps) qg (Scms) 0.0024 0.0047 0.0071 Exp 9.64 10.87 11.63 Sim 9.76 12.46 13.68 Error % 1.2 14.6 17.6 Table 2: Mean angular velocity for different liquid flow rates and qg = 0.0047 sm3/s Angular velocity (rps) ql (Scms) 0.00074 0.00110 0.00147 Exp 8.15 10.87 14.24 Sim 9.48 12.46 14.86 Error % 16.3 14.6 4.3 Figure 8: Comparison of simulated and experimental free surface for qg = 0.0047 sm3/s and ql = 0.00143 m3/s The simulated free surface, presented in Figs. 7 and 8, is obtained by choosing the surface at which the volume fraction rg = rl = 0.5. A satisfactory match between experimental and numerical results is shown; in fact, vortex shape and length are quite similar. 5 Copyright © 2004 by ASME visualization of the corresponding experiment and indicates that the proposed numerical model reproduces appreciably well the performance of the experimental model. It may be also observed that the entrapment gas is removed from the liquid within the vortex, while free-gas, originally separated from the stream right before the entrance as mentioned, remains out of the liquid within the separator body. Table 3, presents simulation and experimental results for different phase volume fractions at separator outlets for liquid and gas flow rates of ql = 0.0011 m3/s and qg = 0.0047 sm3/s, respectively. Table 3: Liquid and Gas volume fractions for two-phase flow with ql = 0.0011 m3/s and qg = 0.0047 sm3/s Liquid Leg Gas Leg Figure 9: 3D view of the capturing surface for qg = 0.0047 sm3/s, ql = 0.0011 m3/s. Figure 9 shows the capturing surface for qg = 0.0047 sm3/s and ql = 0.0011 m3/s. This image scale was reduced 50% vertically for visual proposes. As in the single-phase flow (Fig. 6), this surface presents a notorious 3D shape. It disagrees with the axis-symmetric and mechanistic models assumption as it was previously commented. Liquid Entrapment gas Free gas Total gas Figure 10: Computed volume fractions at mid-longitudinal plane of GLCC for qg = 0.0047 sm3/s and ql = 0.0011 m3/s Figure 10 depicts the volume fraction plot at the mid longitudinal plane of the separator. This result agrees with the rg rl rg rl Exp 0.00005 0.99995 1 0 Sim 0.007 0.993 0.9923 0.0077 As noticed in Table 3, the magnitude of measurements of the gas void fraction at the bottom stream (GCU) and at the top stream holdup (LCO) are smaller than the error associated to the entire modeling process, which considers the contribution of all possible error sources. In consequence, these values are not reproduced, regardless of the bubble diameter. CONCLUDING REMARKS Three-dimensional simulation of the flow field inside a Gas-Liquid Cylindrical Cyclone (GLCC) separator have been carried out. Both single-phase (water) and two-phase (airwater) flows have been simulated using the finite volume method. The obtained numerical results, in terms of vortex free surface, agree satisfactorily with images of validation experiments. Also, the mean angular velocities calculated from simulations show a reasonable match with experimental data. Since the experiments here reported were performed using air-water mixtures, the gas carry-under volume fraction resulted too small and therefore, impossible to capture numerically using Eulerian-Eulerian models. These models, despite of being largely reliable for engineering purposes, are still limited to errors larger than the very small amount of gas carry-under here explored. The computed capturing surface resulted to be a helicoidal cone, unlike the paraboloid traditionally assumed in mechanistic models. Although it is beyond the scope of this research, it is possible to think of the assumed paraboloidcapturing surface as a potential source of errors in current mechanistic models. Further studies, oriented by the results from this investigation could shed some light on this respect. NOMENCLATURE GLCC = gas liquid cylindrical cyclonic =dissipation of turbulent kinetic energy k = turbulent kinetic energy U = velocity vector 6 Copyright © 2004 by ASME Sk = term of generation of turbulent kinetic energy r = volume fraction. G = Production of turbulent kinetic energy due to the external forces. NP = Number phases. q = volumetric flowrate, m3/s.  = fluid density, kg/m3  = fluid viscosity, Pa.s C = constant empiric. t = time, s m3/s = cubic meters per second. sm3/s = cubic meters per second at standards conditions. Scms= cubic meters per second at standards conditions. rps= revolutions per second. SUBSCRIPTS α ,β = Phases p = phase D = Drag 1,2,3 = Index to constant different. eg = entrapment gas. fg = free gas. g = gas. l = liquid. FEDSM98-5206, 1998 ASME Fluid Engineering Division Summer Meeting, June 1998. [6] Erdal, F., Shirazi, S., Mantilla, I., Shoham, O., 1998, “CFD Study of Bubble Carry-Under in Gas-Liquid Cylindrical Cyclone Separators” SPE 49309, 1998 SPE Annual Technical Conference and Exhibition, September 1998. [7] Erdal, F., Shirazi, S., 2001, “Local Velocities Measurement and Computational Fluid Dynamics (CFD) Simulations of Swirling Flow in a Cylindrical Cyclone Separator” ETCE 2001-17101, Engineering Technology Conference on Energy, February 2001. [8] AEA Technology, “CFX 4.3 Solver Manual”, 1997. [9] Thomson, C. and Wilkes, N., 1982, “Experiments with Higher-Order Finite Difference Formulae”, AERE-R 10493. SUPERSCRIPTS T = Transpose Tensor d = Drag. ACKNOWLEDGMENTS The authors would like to thank FONACIT-Venezuela for the financial support to this project. Antonio MelendezRamirez wishes to thanks “Decanato de Investigación y Desarrollo” (DID) at the Universidad Simón Bolívar for supporting his M. Sc. Studies. REFERENCES [1] Shoham, O., Kouba, G., 1998, ''The state of the art of GasLiquid Cylindrical Cyclone Separator,'' Journal of Petroleum Technology, Vol. 50, No. 7, July 1998, pp. 58-65. [2] Erdal, F., Shirazi, S., Shoham, O., Kouba, G, 1997, ''CFD Simulation of Single-Phase and Two-Phase Flow in Gas-Liquid Cylindrical Cyclone Separators,'' SPE 36645, SPE 71st Annual Meeting, SPEJ, Vol. 2, December 1997, pp 436-446. [3] Farchi, D., “A Study of Mixer and Separator for Two-Phase Flow in M. H. D. Energy Conversion Systems” M. S. Thesis (in Hebrew), Ben-Gurion University, Israel, 1990. [4] Motta, B., Erdal, F., Shirazi, S., Shoham, O., Rhyne, L., 1997, ''Simulation of Single-Phase and Two-Phase Flow in Gas-Liquid Cylindrical Cyclone Separators,'' FEDSM97-3554, 1997 ASME Fluid Engineering Division Summer Meeting, June 1997. [5] Erdal, F., Mantilla, I., Shirazi, S., Shoham, O., 1998, “Simulation of Free Interface Shape and Complex Two Phase Flow Behavior in a Gas-Liquid Cylindrical Cyclone Separator” 7 Copyright © 2004 by ASME