Inlet conditions effect on bubble to slug flow transition in mini-channels

Chemical Engineering Science 102 (2013) 106–120

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Inlet conditions effect on bubble to slug flow transition in mini-channels Jurij Gregorc n, Iztok Žun Laboratory for Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, University of Ljubljana, Aškerčeva 6, 1000 Ljubljana, Slovenia

H I G H L I G H T S








Experimental investigation of bubble to slug flow transition in mini-channel. Porous media mixer and vs. cross-junction mixer. Numerical algorithm for inlet conditions in continuous mixing gas–liquid flow. Numerical prediction of bubble to slug flow transition. Bubble size distribution as a function of inlet conditions.

art ic l e i nf o

a b s t r a c t

Article history: Received 12 March 2013 Received in revised form 17 July 2013 Accepted 28 July 2013 Available online 7 August 2013

The paper deals with the impact of inlet conditions on bubble to slug flow transition in mini systems. A new experimental test loop with a glass mini-tube (D¼ 1.2 mm ID) has been constructed to assess the effects of inlet conditions on the two-phase flow pattern development in the spatial and temporal domains. The interchangeable inlet part of the test section allowed different geometrical combinations for the mixing of gas and liquid prior to it entering the mini-tube. Porous media mixer and cross-junction mixer were considered. High speed video recordings were taken of 70 combinations of flow rates, corresponding to superficial velocities ranging from 0.2 to 11 m/s and 0.25 to 3 m/s for air and water, respectively. The following discernible flow patterns are considered: bubbly, slug and semi-annular flow. No significant differences were found when comparing flow pattern maps for each mixer. Digital image post processing of high speed video recordings was used to estimate the equivalent diameter for every gas structure. A comparison of the distribution of bubbles with equivalent diameter revealed the inlet mixer's strong impact on bubble size and bubble distribution along the mini-tube and thus, the bubble to slug flow transition. The VOF method, implemented in ANSYS Fluent was used as the numerical tool to predict the flow patterns in a mini-tube of 1.2 mm ID. A novel approach to mimic continuous mixing is presented. By changing only the prescribed flow rates, different flow patterns can be simulated. Similar interfacial structures were obtained by numerical simulation and experiment for both mixers. Reasonable quantitative agreement was also achieved when analyzing bubble to slug flow transition. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Bubble to slug flow transition Flow pattern Interface CFD VOF Static mixer

1. Introduction Mini and micro fluidic systems have been extensively studied in recent years. The large interest in applying these solutions to real life applications is mainly due to the enhancement of heat and mass transfer, unit size reduction, process parallelization, continuous operation and many other factors. Micro chemical technology, biotechnology, microelectromechanical system (MEMS), micrototal analysis system (μTAS) and thermal management systems are only a few areas where mini/micro fluidic devices could

n

Corresponding author. Tel.: +386 1 47 71 404; fax: +386 1 47 71 447. E-mail address: [email protected] (J. Gregorc).

0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.07.047

be employed. Two-phase flow is frequently encountered in such devices. Primary flow patterns, e.g. bubbly flow, slug flow and annular flow can be observed on a macro and micro scale with the exception of stratified flow, which is generally suppressed as reported by Rahim et al. (2011). A significant difference can be found with regard to flow development. General flow pattern characteristics are largely decided in the inlet/mixing region. This fact has a large impact on flow pattern dependent processes. Different mixing geometries have been proposed over the years. Flow focusing, co-flowing and T-junction structures appear as unit operations and are frequently used in the flow of immiscible liquids, as reviewed by Zhao and Middelberg (2011). In publications referring to gas–liquid two-phase flow, similar geometrical arrangements have been used. T-junction mixers can be

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

differentiated based on the size of each part of the junction and the position of continuous phase introduction. Choi et al. (2011) used an equi-sized T-junction with axial introduction of deionized (DI) water (continuous phase) and side introduction of N2 (dispersed phase). Wang et al. (2011) on the other hand, introduced the dispersed phase in the axial direction and continuous phase from the side. The cross flow approach was used by Cubaud et al. (2005) with equal inlets and Triplett et al. (1999) with smaller inlets from the sides. These types of mixers are referred to as cross-component or cross-junction. Co-centric nozzles represent a very convenient way of producing annular type flows. Liquid is typically introduced at the channel walls, while gas flows through the nozzle that is positioned in the channel center. Many other mixing section constructions can be found. Mak et al. (2006) for example, used an axial introduction for air and porous wall introduction for liquid to create annular flow in a 5 mm pipe. In our previous work, a porous media mixer was presented in a combination with 1.2 mm hydraulic diameter semicircular channel (Zun et al., 2010). Many other variations can be found, however the listed examples cover the majority of the two-phase flow mixing devices used. Despite tremendous progress in computational fluid dynamics (CFD) and a significant increase in computational power, our ability to predict two-phase flow patterns over a wide range of flow rates remains far from complete. The multiscale nature of two-phase flow requires correct meshing and transient solving with small time steps, resulting in a very large computational effort to predict very short time samples. As a consequence, different simplifications have been used to reduce computational costs. Only interface tracking techniques are considered here as no form of two fluid model is capable of predicting interface between the phases. One of the most obvious simplifications is solving for 2D. This approach is widely used. The debate whether this approach is justified is still open among users, as reviewed by Talimi et al. (2012). According to Goel and Buwa (2009), care should be taken even when pipe flow is studied using axisymmetric formulation. The computational domain is reduced in many cases. Different values of minimum domain length are used but generally around 40D is needed for the flow to develop (Lakehal et al., 2008). Another way to reduce computational time is to reduce mesh density. Gupta et al. (2009) conducted an extensive study of mesh influence on the prediction of Taylor flow. As a size criterion they proposed a minimum of 5 cells in liquid film. In their study only equi-sided and right angled cells of equal size were considered (“perfect mesh”). Unfortunately, this meshing strategy is not applicable to many real engineering applications due to geometrical complexity and 3D geometry. In the current work, we present the impact of initial conditions on two-phase flow pattern development with an emphasis on

107

bubble to slug flow transition in a 1.2 mm ID horizontal mini-tube. A new methodology is presented showing how to deal with this issue numerically. Two inlet geometries are considered as the test cases: cross-junction mixer and porous media mixer. In Section 2 the impact of the initial conditions on two-phase flow patterns is discussed based on the experimental study. First, an experimental test section and procedure are given, followed by flow pattern analysis based on visual observation and the study of bubble to slug flow transition affected by the inlet mixer. Section 3 introduces a numerical prediction of different flow patterns. The commercial code ANSYS Fluent was used in 3D formulation to solve conservation equations and the VOF model was used to track the gas–liquid interface. Our main contribution in the numerical part is a novel approach with initial conditions that enable the modeling of different flow patterns. When searching for the flow pattern transition, indicators at different flow rates are the main goal. Finally, numerical predictions are compared and validated against experimental data in Section 4.

2. Experiment The experimental setup used to study the effect of the inlet condition on two-phase flow is shown schematically in Fig. 1. The test section is composed of a glass mini-tube with D¼ 1.27 0.08 mm ID (D¼6 70.2 mm OD) and length equal to L¼ 240 mm which corresponds to L/D ¼200. Two different inlet mixers were used in the present study. The first was the porous media mixer as proposed by Zun et al. (2010) and the second was a conventional cross-junction type mixer. The porous media mixer (Fig. 2a) was composed of a typical Y-junction followed a mixing chamber made from Plexiglas. The mixing chamber was filled with glass spheres 0.5 mm in diameter. Air and water were introduced through Y-junction ports, respectively. In the case of the crossjunction mixer, two perpendicular channels with a diameter of 1.2 mm were precisely bored in a Plexiglas cube (Fig. 2b). Water was introduced in the axial direction, while equal flow rates of air were introduced via the side ports of the cross-junction (Fig. 2b). Each mixer could be precisely connected to the mini-tube forming one continuous channel. The support structure of the test section was made from optically transparent material to enable the best possible visual studies. The entire system was pressure driven to avoid eventual fluctuations. Air was taken from the laboratory main, filtered and set to constant pressure. Two identical sets of thermal mass flow controllers were used to measure and control air flow rates. Each set consisted of two EL-Flow series controllers; the first with a range of 0–3 l/h and the second with a range of 0–72 l/h, supplied by Bronkhorst High-Tech Instruments. The mass flow

Fig. 1. Schematics of experimental setup (dot line – air, solid line –water, dash-dot line – air–water two-phase flow).

108

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 2. 3D models of the mixer used: (a) porous media mixer (left), (b) cross-junction mixer (right).

controllers have an accuracy of 7(0.5%Rd+0.1%FS). The constant pressure tank was filled with clean tap water prior to measurements. Water flow rate was controlled by the Cori-Flow series coriollis mass flow controller with a flow rate range of 0–20 l/h, also supplied by Bronkhorst High-Tech Instruments. The Coriolis mass flow controller has an accuracy of 7(0.2%Rd+8 g/l). This arrangement of mass flow controllers enabled us to test different inlet conditions. After metering, both phases were introduced to the inlet/mixing part of the test section. Two-phase flow then entered the mini-tube where measurements to assess effects of the inlet/mixing part took place. Measurements were conducted over a wide range of air–water volume flow rate combinations. With water volume flow rates ranging from 1 to 12 l/h (corresponding to liquid superficial velocities from 0.25 m/s to 2.95 m/s) and gas volume flow rates ranging from 1 to 72 l/h (corresponding to gas superficial velocities from 0.25 m/s to 17.68 m/s) the following discernible flow regimes were achieved: bubbly flow, slug flow and semi-annular flow. 70 combinations of air–water volume flow rates were used to define flow regime transitions in the case of each inlet mixer. 2.1. Experimental procedure Different experimental series based on phase flow rate combinations were performed. Each measurement started with setting the appropriate flow rate of liquid, followed by gas flow rate. When both flow rates reached 99% of the respected set point value (typically 10–30 s) a relaxation period of 30 s started allowing two-phase flow to stabilize. The data acquisition period then commenced. Pressures were recorded for a 1 min time sequence and scanned at 1 kHz. Data acquisition was performed using a USB-6216 DAQ Card, supplied by National Instruments. To track the behavior of gas–liquid interface flow structures a SpeedCam Visario high speed video system by Weinberger was employed. Based on the studied phenomena and the selected volume flow rate, different sampling frequencies were chosen, ranging from 4000 fps to 10,000 fps. After the high speed video was recorded and the pressure acquisition time had elapsed, a certain amount of time was needed to transfer the high speed video data to the computer. During that time gas and liquid flow rates were stopped. The measurement procedure was developed prior to the measurement campaign. Doing so enabled us to achieve high reproducibility of two-phase flow as can be seen in Fig. 3. The experimental setup operation was fully automated and controlled by the measuring procedure developed in the LabView platform, also supplied by National Instruments.

Fig. 3. Bubble residence time distribution for two independent measurements with the same setup configuration, using the same flow rates, measured at 125D downstream from the mixer.

2.2. Flow patterns Two-phase flow patterns identification was done based on high-speed video recording at approximately 125D downstream from the mixer. Three major flow patterns have been recognized in the case of both mixers: bubbly flow, slug flow and semi-annular flow. Identified flow patterns resemble those found by Triplett et al. (1999) in a 1.097 mm glass mini-tube and Zun et al. (2010) in a 1.2 mm hydraulic diameter semi-circular channel made of Plexiglas. Bubbly flow developed when a low volume flow rate of gas and a high volume flow rate of liquid were introduced into the mixer. Bubbly flow is characterized by a large number of irregularly shaped bubbles that are smaller than the tube diameter. Fig. 4 represents a still image in the case of the porous media mixer (denoted by a) and the cross-junction mixer (denoted by b) at high bulk liquid Reynolds numbers. Evidently interfacial structures do not differ much. The bubbles exhibit rather irregular shapes, presumably due to the shear forces that are caused by rather high flow rates at moderate (mini) pipe diameter. A closer inspection of a time sequence by high speed video recording revealed larger bubble number density fluctuations in the case of the porous media mixer when compared with the cross-junction mixer at the

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

109

Fig. 4. Still images of bubbly flow: (a) porous media mixing, (b) cross-junction mixing (QL ¼ 12 l/h, QG ¼ 2 l/h in both cases, flow direction from left to right). ReL ¼3530, ReG ¼40.

Fig. 6. Still images of slug flow al lower flow rates: (a) porous media mixing, (b) cross-junction mixing (QL ¼ 1 l/h, QG ¼2 l/h in both cases, flow direction from left to right). ReL ¼ 294, ReG ¼ 40.

Fig. 5. Still images of bubbly flow: (a) porous media mixing, (b) cross-junction mixing (QL ¼ 3 l/h, QG ¼1 l/h in both cases, flow direction from left to right). ReL ¼882, ReG ¼20.

Fig. 7. Still images of slug flow at higher flow rates: (a) porous media mixing, (b) cross-junction mixing (QL ¼ 5 l/h, QG ¼6 l/h in both cases, flow direction from left to right). ReL ¼ 1471, ReG ¼121.

same flow rates. On the other hand, bubbly flow can be also discerned at lower liquid flow rates (Fig. 5). Here, very deterministic bubble shape and distribution along the mini-tube axis can be seen especially in the case of the cross-junction mixer. Similar shapes of the interfacial flow structures were also observed by Zhang et al. (2011) and Chen et al. (2002) for water and nitrogen in 0.5 mm micro-tube and 1.5 mm mini-tube, respectively. None of the flow rate combinations, however, which were tested in this work, produced bubbly flow with spherical bubbles. According to the literature, the existence of spherical bubbles does not depend so much upon the channel size (Sur and Liu, 2012;De Menech et al., 2008); it rather depends on mixer configuration (which phase is introduced at which inlet) and the bulk liquid Reynolds number. Zhang et al. (2011) and Sur and Liu (2012) observed spherical bubbles in micro-tube of approx. 0.3 mm ID for a Reynolds number of 600 and lower. Such flow conditions were not tested in this work. Decreasing the liquid volume flow rate eventually leads to slug flow. It is characterized by smooth elongated bubbles more than one tube diameter in length with a liquid film separating the gas structures from the wall. Fig. 6a and b represent slug flow at lower volume flow rates for the porous media mixer and cross-junction mixer, respectively. Even though the flow pattern can be classified as slug flow in both cases, it is evident that the elongated bubbles differ in length. In the case of the cross-junction mixing, the bubble length is shorter in comparison with the porous media mixer. This difference could play an important role in cases of heat and mass transfer or phase separation. At higher volume flow rates, slug flow is no longer affected by the initial mixer to such an extent. With both mixers elongated bubbles of varying lengths tend to form local train segments followed by liquid slug (Fig. 7). When increasing the gas volume flow rate beyond 10 l/h at a low liquid flow rate, annular-type flow emerges. Starting from slug–annular flow, the gas–liquid interface eventually develops into annular intermittent flow and in limited cases into stable annular flow (Fig. 8). Analyzing the time sequences revealed that liquid film was constantly being torn apart, especially in the case of the porous media mixer. Due to this behavior micro bubbles

Fig. 8. Still images of annular flow: (a) porous media mixing, (b) cross-junction mixing (QL ¼1 l/h, QG ¼ 22 l/h in both cases, flow direction from left to right). ReL ¼ 294, ReG ¼ 443.

were often trapped in the liquid film, moving slowly along the mini-tube (Fig. 8b). The appearance of stable annular flow was rare. At a similar gas volume flow rate and higher liquid volume flow rate the disturbance of liquid film occurs more frequently. It is typically coupled with a breakup of the gas–liquid interface, forming churn like flow structures and is therefore named churn–annular flow. Trying to differentiate between the described annular based flow patterns did not supply us with any kind of justifiable criteria. Therefore, we decided to call the entire group of flow patterns semi-annular flow, following our previous studies (Zun et al., 2010). Fig. 9 shows flow pattern maps for the porous media mixer and cross-junction mixer constructed from data on interfacial structures obtained by visual analysis as described above. Three major regions can be noted: bubbly flow region, slug flow region and semi-annular flow region. Our data are compared with experimental transition lines from Triplett et al. (1999) in a 1.097 mm ID glass tube with cross-type mixer. They observed bubbly, slug, churn, slug–annular and annular flow patterns. Due to our decision to group annular-like flow patterns only bubble to slug flow transition (denoted as B–S in Fig. 9b) can be compared. Fair agreement can be seen in that region (Fig. 9b). When comparing the flow pattern maps, obtained by the porous media mixer and the cross-junction mixer, it may be concluded that the effect of the initial conditions on two-phase flow interfacial structures seems to

110

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 9. Flow pattern maps: (a) porous media mixer, (b) cross-junction mixer; B–S¼ bubble to slug, S–C¼ slug to churn, S–S/A ¼slug to slug-annular, S/A–A ¼ slug–annular to annular, A–C¼ annular to churn, S/A–C¼ slug–annular to churn.

be negligible. However, two important facts have to be taken into account 1. While the same flow patterns have been recognized with both mixers at the same flow rates, it was evident from the video analysis that actual gas–liquid interface flow structures differ depending upon the mixer used. The differences are in characteristics like bubble size and bubble frequency. The essence of the problem is clearly depicted in Fig. 6. While using the jL–jG plot for flow pattern maps offers the possibility for a fast comparison of data gained on different channel sizes or shapes, it clearly lacks information on other key aspects of two-phase flow e.g.: bubble size and frequency, liquid film thickness, void fraction, relative velocities of phases, etc. A more profound multiscale approach to construct flow pattern maps should be adopted to include this information, as proposed by Thome et al. (2013). 2. In both cases non distinctive (mixed) flow patterns emerged in the regions of flow pattern transitions. A subjective approach to flow pattern recognition in these cases is unreliable. Therefore, a more objective approach was employed based on a qualitative analysis of equivalent bubble diameter.

2.3. Detection of bubble to slug flow transition Bubble to slug flow transition has been the subject of numerous studies both in the past (Mandhane et al., 1974; Taitel and Dukler, 1976) and present. Most models predicting flow patterns and flow pattern transitions were developed based on theory and experimental data collected on pipe diameters ranging from 0.01 m to 0.15 m. In most cases, a substantial discrepancy was found when applying these models to mini-channels (Triplett et al., 1999; Ullmann and Brauner, 2007). In order to track bubble to slug flow transition as a function of the inlet mixer, a simple yet very rigorous method is needed. Considering that both numerical and experimental results have to be evaluated using the same approach we decided to estimate equivalent bubble diameter for every gas structure at a certain position in the mini-tube. The measurement position was 125D and 20D downstream of the mixer for experiment and numerical simulation, respectively. The values of equivalent bubble diameter were normalized by the mini-tube diameter and sorted into classes. A single value threshold was defined to declare when the flow was no longer bubbly. The threshold value used in this paper was set to unity (deq/Dch ¼1) meaning that as long as all gas

structures had an equivalent bubble diameter smaller than the tube diameter the flow was characterized as bubbly flow. This is in agreement with the bubbly flow definition at the beginning of Section 2.2.

2.3.1. Estimating equivalent bubble diameter High speed video recordings, in the case of experiments, and animations of two-phase flow in the case of numerical simulations, were used to collect data on equivalent bubble diameters. The three step digital image post processing approach was adopted to estimate equivalent bubble diameters from a 2D projection of every gas structure - Firstly, a sequence of filters and arithmetic operations was carried out. Starting with conversion to 8bit gray scale (Intensity image), followed by subtraction of background and low pass smoothing with a Gauss two dimensional filter to soften the image. Exponentiation was employed next with the intensity value of each pixel as a base. The value of the exponent was set in a way that after exponentiation none of the pixels exceed the maximum intensity of 255. - Secondly, edge enhancement was achieved using an erode two dimensional filter. In order to separate between gas and liquid, threshold intensity was defined by evaluating the image intensity gradient profile in a perpendicular direction to the mini-tube axis. - The last step was to measure the gas structure area as projected to face view and calculate the equivalent bubble diameter as a function of tube diameter. The equivalent bubble diameter is defined by (1), where deq,i represents equivalent bubble diameter and Ai is bubble projection area as recorded by high speed video. rffiffiffiffiffiffiffi 2 4Ai ð1Þ deq;i ¼ π The post processing and measuring routine was developed and performed using the commercial software ImagePro Analyzer 6.2 (2004). Due to memory size limitation only 1 s of high speed video could be acquired at once. To verify that the sample size was large enough to assure statistical relevance, we compared equivalent bubble diameter distributions for different sample lengths. Equivalent bubble diameter data were first extracted from one full length video recording. Next, the same video recording was split into two time sequences with equal duration. Equivalent bubble diameter

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 10. Equivalent bubble diameter distributions for different sample lengths at equal flow conditions.

extraction was taken from both sequences. A similar procedure was done for a quarter of the sequence that was randomly picked out of the initial video. The differences between half and full length samples did not exceed 3%. According to the test results (Fig. 10) half-length samples are representative. 2.4. Results and discussion In search of transition line it is well known that all flow patterns cannot be objectively discriminated. To elaborate such analysis one may use additional criteria like frequency domain characteristics shown in our recent paper (Zun et al., 2010), chaos theory (Wang et al., 2003), or neural network approach (Xie et al., 2004). In such cases we have used subjective judgment as an additional criterion to normalized equivalent bubble diameter distributions. Characterization of two-phase flow patterns based on a combination of subjective judgments and objective methods has also been frequently reported in the literature, like (Jones and Zuber, 1975, Mishima and Ishii, 1983, Matsui, 1984, Zhao and Bi, 2001, Zun et al., 2010). The aim of present studies is to search for mechanisms that govern the bubble to slug flow transition. In the number of laboratory data sets, the first one was taken at QL ¼ 12 l/h, QG ¼2 l/h (corresponding to ReL ¼3530, ReG ¼40). The bubbly flow regime obtained at such flow conditions is presented in Fig. 4. The liquid volume flow rate was then gradually decreased to QL ¼1 l/h, QG ¼2 l/h (corresponding to ReL ¼294, ReG ¼ 40) where slug flow was visually recognized (Fig. 6). All data sets obtained are marked by the semi-transparent zone in Fig. 9. We selected this zone at rather high liquid flow rates to enable sufficiently large data package at “transitional borderline” that demarcates bubbly from slug flow regime. In such a case, the interfacial structures are very deterministic. At low liquid flow rates, the bubble to slug flow changes are much more abrupt and generated by the instabilies in mixing devices that are beyond the scope of this paper. The paper aims to demonstrate numerical simulation capabilies of continuous gas and liquid flow streams which is believed to be the novelty in the literature of VOF modeling. Fig. 11 shows distributions of bubbles within a specific range of equivalent bubble diameters. Data for the porous media mixer and cross-junction mixer are presented using black square and gray circle symbols, respectively. A transition threshold between bubbly flow and slug flow regimes is defined by deq/Dch ¼ 1, hereafter denoted by the BTS (bubble to slug flow transition) line. The

111

number of classes is equal for both mixers and was determined using Sturges' criteria (Sturges, 1926). Every symbol has a horizontal error bar. It represents the smallest and the largest value of equivalent bubble diameter belonging to its class. Comparing relative frequencies between the two mixers at QL ¼12 l/h and QG ¼2 l/h (Fig. 11a), it can be noticed that the number of the smallest bubbles in the case of the porous media mixer is more than twice as high as in the case of the cross-junction mixer. Decreasing liquid volume flow rate to 8 l/h, while keeping gas volume flow rate at a constant value (Fig. 11b), produces bubbles that lay above the BTS line, in the case of the porous media mixer. Most of the gas entities, however, still resemble those typical for bubbly flow. This observation implies the onset of bubble to slug flow transition. In the case of cross-junction, bubbly flow pattern remains. A further decrease in liquid flow rate does not lead to any significant changes in terms of flow pattern transition (Fig. 11c). The majority of bubbles are still smaller than unity in both cases. A very deterministic formation of interfacial structures in the cross-mixer is evident, with over 65% of bubbles belonging to the same size class. Fig. 11d corresponds to the experimental run with 3 l/h of liquid and 2 l/h of gas flow rates used. A small percentage of bubbles lying above the BTS line indicate that the flow pattern is no longer pure bubbly in the case of both mixers. While very uniform flow remains in the case of the cross-junction mixer, a large scattering of equivalent bubble diameter data can be observed when the porous media mixer was used. In the case of 1 l/h of liquid and 2 l/h of gas flow rates the cross-junction mixer produces the majority of gas structures larger than the tube diameter. Only a very small number of small bubbles remained, as shown in Fig. 11e. On the other hand, bubbles represent more than one third of all gas structures in the case of the porous media mixer. From the graphs (Fig. 11) it is clearly noticeable that flow patterns are influenced by the type of mixer used. This fact can be confirmed by a direct comparison of experimental runs at equal flow rates, as well as by tracking the BTS transition, employing strict rules to discern between flow patterns. 2.4.1. Flow pattern vs. entrance length Residence time distribution stabilizes very quickly, while bubble shape requires a longer entrance length. Fig. 12 shows residence time distribution for elongated bubbles at different locations downstream. The resident time distributions at different positions downstream reveal that not much happens beyond 5D–10D in terms of coalescence and breakup. The flow pattern development length changes slightly depending on the flow regime, as reported by Talimi et al. (2012). According to our experimental observations and data available from the literature (Lakehal et al., 2008), the entrance length needed for bubble shape development is approx. 40D.

3. Numerical approach The commercial CFD code ANSYS FLUENT (Release 14.0, 2012) was chosen to predict two-phase flow patterns using the Volume of Fluid (VOF) numerical scheme. Assuming one fluid adiabatic flow, the continuity and momentum equations are given by (2) and (3), respectively. ∂ρ þ ∇  ðρvÞ ¼ 0 ∂t

ð2Þ

∂ ðρvÞ þ ∇  ðρvvÞ ¼ ∇p þ ∇τ þ ρg þ F; ∂t

ð3Þ

ρ represents density, v is velocity vector, p is static pressure, τ is stress tensor (given by (5)), ρg is gravitational force and F is representing external forces. Material properties (eg. density and viscosity) are defined as a function of the local phase fraction in

112

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 11. Relative frequency of specific normalized equivalent bubble diameter at different flow rate combinations.

Eq. (4). ρ ¼ αG ρG þ ð1αG ÞρL μ ¼ αG μG þ ð1αG ÞφL

ð4Þ

τ ¼ μ½ð∇v þ ∇vT Þ

ð5Þ

Gas–liquid interface is tracked via the volume fraction continuity equation (6) for gas (secondary) phase. Volume fraction of liquid (primary) phase is on the other hand determined based on Eq. (7).   1 ∂ ðαG ρG Þ þ ∇  ðαG ρG vÞ ¼ 0 ρG ∂t

ð6Þ

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

113

Fig. 13. Geometry model discretized using hexahedral cells.

Fig. 12. Residence time distribution of interfacial structures measured at different locations downstream.

αG þ αL ¼ 1

ð7Þ

Surface tension effects were accounted for by the Continuous surface model (CSF) described in Brackbill et al. (1992). The formulation is shown as F vol ¼ sL;G

ρκL ∇αL ; ð1=2ÞðρL þ ρG Þ

ð8Þ

here sL,G is the surface tension coefficient for chosen system of fluids, ∇αL is the normal to the interface, κL is interface curvature. 3.1. Flow solver settings ANSYS Fluent incorporates the Finite Volume framework to discretize and solve flow equations. In this work both phases were assumed to be incompressible at all times. Therefore, the pressure based solver was employed. Two-phase flow pattern formation and development were treated as transient phenomena. Temporal discretization was done using first order explicit and first order implicit scheme for the volume fraction equation and momentum equation, respectively. Pressure velocity coupling was done using the PISO method. Pressure discretization was done via the Presto! scheme. The second order upwind scheme was used to discretize the velocity field. Gradients were estimated with the Least Square cell based method. Gas–liquid interface was interpolated using the Geometric Reconstruction scheme. Variable time stepping was used to optimize simulation speed. Time step size was controlled by the Courant–Friedrichs–Levy condition with a Courant number of 0.25 for both the volume fraction and momentum equations. The convergence of numerical simulation was checked within every time step using the absolute convergence criteria of globally scaled residuals. The criteria of 1E  4 were used for all the variables. The maximum number of iterations per time step was set to 20. 3.2. Geometric model and spatial discretization A 3D geometry model was constructed to mimic the experimental setup shown in Fig. 1. The computational domain has a circular cross-section with a diameter of 1.2 mm. The length of the domain was shortened to L/D ¼40, based on the argument in Section 2.4. Spatial discretization was done using the ANSYS Gambit 2.6 pre-processor. The domain was split into 5 blocks enabling the construction of a fully structured hexahedral mesh

with 1,960,000 cells (Fig. 13). The chosen mesh type is in accordance with Rek and Žun’s (2013) best performance tests of full 3D VOF simulation on different mesh types. This is also in accordance with ANSYS Fluent documentation (2011), where the use of hexahedral or quadrilateral mesh is suggested when surface tension effects are important. In this work mesh refinement was done in the wall region. The ratio between the largest cell edge in the core of the domain and the smallest cell edge at the wall was equal to 5. The size of the wall cells in the radial direction was 5 μm. This yielded approx. 7 cells per liquid film. 3.3. Phases and material properties Simulations were done for water and air as the primary and secondary phases, respectively. Material properties used are given in Table 1. 3.4. Boundary and initial conditions The following boundary conditions were superimposed: Walls: no-slip, static contact angle of 601 (air–water on glass surface). Inlet: mass flux was prescribed (porous media mixer), uniform velocity (cross-junction mixer). Outlet: constant gauge pressure p¼ 0 Pa. Single phase flow simulations (water only) were run first to obtain the developed flow fields. These results were used to initiate the flow field for the two-phase flow simulation. In order to simulate a two-phase flow pattern, the initial conditions that enable continuous production of interfacial structures at the inlet and the evolution of gas–liquid interface further downstream must be defined. Generally speaking, at least two inlets are needed in order to form two-phase flow: one for each phase. The simplest mixer geometry available would therefore be a junction of three channels with two inlets and one outlet. Various types of mixing inlet parts can be used as reviewed in Section 1. When predicting the flow pattern numerically, geometric and operational parameters must be accounted for. An approach needs to be developed to mimic inlet mixing that produces relevant flow patterns at different flow rates. 3.4.1. The porous media mixer The porous media mixer shown in Fig. 14 was used as a test case. The same mixer has already been utilized in phase separation process as an inlet to mini manifold that has been reported by Zun et al. (2010). The mixer consists of two main components: the first is a typical Y-junction and the second is a mixing chamber filled with glass spheres with 0.5 mm diameter. The mixing chamber represents a miniaturized version of frequently used devices in macro systems in process engineering. At the current stage we did not test how the Y entrance to mixing chamber influences flow pattern formation. The focus was on the behavior of two-phase flow

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

114

Table 1 Material properties of fluids used in simulations. Fluid

Phase

Density [kg/m3]

Viscosity [kg/(m s)]

Surface tension coef. [N/m]

Water Air

Primary Secondary

998.2 1.225

1.003  10  3 1.7893  10  5

0.072

Fig. 14. Porous media mixer.

Table 2 Geometrical parameters of experimental mixer and numerical algorithm. Parameter

Experimental mixer (mm)

Numerical algorithm (mm)

〈dc 〉 dp dt;min

0.34 0.44 (average) 0.11

0.34 0.81 (maximum) 0.07

in the sphere packed chamber. Our experimental experience on two-phase flow in the trickle bed reactor (Žun et al., 1997) suggests that porous media filling behaves like multiple channels interlacing each other, which promotes interface breakup and coalescence. Different flow rates of gas and liquid also result in different twophase flow regimes in porous media, as noted by Žun et al, (1997). On the single passage scale this means that at certain conditions more liquid may be flowing, while at slightly different conditions only gas could be present. This behavior changes with time. Furthermore, the passage can exhibit a stall condition at which the flow velocity of either phase is zero. The latter phenomenon is referred to as pore activity. It could be a consequence of pressure– velocity fluctuations between nearby pores. Following trickle bed literature, average pore diameter 〈dp〉 as defined in equation (9) does not reflect the true passage size as it only represents large pore chambers, according to Ng (1986). Average channel diameter 〈dc〉, defined by (10) that is smaller than average pore diameter 〈dp〉 and larger than minimum pore throat dt,min, defined by (11), can be used instead. Nc denotes the number of passages per unit area defined by (12). ε denotes porosity of the porous media mixer that was measured in our work to be 0.41 when solid spheres with diameter ds ¼ 0.5 mm were randomly packed in the mixing chamber. 〈dp 〉 ¼

 ε 1=3 ds 1ε

sffiffiffiffiffiffiffiffi 4ε 〈dc 〉 ¼ πN c dt;min ¼

Nc ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 2 π 1 sin  ds π 3 2

6ð1εÞ πds

2

ð9Þ

ð10Þ

ð11Þ

ð12Þ

Note that Eqs. (9)–(12) are valid only if spherical filling material of equal size is used.

In the present numerical simulations, we did not simulate twophase flow in the porous media mixing chamber. Instead, we developed an inlet condition algorithm that accounts for the three major parameters of the porous media mixer: inlet geometry, gas flow rate and liquid flow rate. Multiple inlets of different size were created at the initial cross-section of the computation domain. Geometrical parameters, listed in Table 2, indicate general guidelines when forming inlets. The inlet cross-section of the computational domain was first split into small pieces, here referred to as “patches”. The split was done by “cutting” the inlet crosssection apart based on an 8  8 stencil with equal spacing in both directions (Fig. 15a and b). Stencil spacing was determined with a minimum throat diameter dt,min defined by (11). 60 inlet patches were formed (Fig. 15b) with a size ranging from minimal throat diameter which was in our case 0.07 mm to maximal patch size that was in our case 0.15 mm. A random merging of inlet patches was used to mimic gas phase flow penetrating through porous media thus forming various gas inlets as shown in Fig. 15c. At such formed inlets gas mass flux boundary conditions were superimposed. The prescribed values were computed based on experimental volume flow rates. In the first approximation, we neglected different slip conditions at gas–liquid interfaces generated in different inlet locations. This led to certain discrepancies comparing the integrated gas volume flow rates with those being measured. Due to the high costs of computational time any further iteration refinement was not considered since it would not have brought towards any significant changes in flow pattern predictions. Such inlets are limited by the experimental values of average channel diameter 〈dc〉 and overall number of passages Nc. Table 2 summarizes geometrical parameters in the case of the experimental porous media mixer and the numerical algorithm. Conversely, all the remaining patches were merged into a single inlet where the liquid phase was introduced. In the current implementation of the porous media mixing algorithm, both phases were fed into the flow domain as continuous streams. There were no abrupt changes of gas mass flux at the inlets that could promote bubble detachment from the inlet. This was done to ensure that the formation of interfacial structures was the result of the breakup phenomena in the initial part of the mini-tube. The pore activity and mechanism that would enable pressure–velocity fluctuations between the nearby pores have not been included in this work. Tests have disclosed a much slower convergence and pressure fluctuations resulting in unrealistic interfacial structures further downstream of the mixer. We are currently working on resolving this issue.

3.4.2. The cross-junction mixer On the other hand, the cross-junction mixer used in the experiments is schematically shown in Fig. 16a. When mapping real geometry to the numerical model, boundary conditions were not prescribed at the circular faces by the side arm inlets as marked in Fig. 16a. Instead, we generated two corresponding inlet regions along the mini-tube wall. Inlet regions are marked by gray shading in Fig. 16b. The construction of structured, body-fitted, hexahedral mesh on geometry that comprises two intersecting pipes of equal diameters usually results in poor mesh quality due to high cell skewness. A simplified situation for the cross-junction

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

115

Fig. 15. Schematics of inlet section splitting procedure.

computational domain was no longer mesh dependent. Mesh density influence was first assessed by comparing the corresponding structural function obtained along the pipe axis. Due to the variable time step, structural function could not be evaluated at the same time in all cases. Simulation results at approximate time instances were used instead. For easier comparison all three data sets were slightly shifted so that the beginnings of the rightmost interface structure overlap (Fig. 17). A good agreement was found between the results on mesh 2 and mesh 3. This comparison was performed to assure that the distribution of interfacial structures along the computational domain center line is independent of the mesh. Next, bubble shapes were visually compared. Snapshots of numerical simulations were placed over the graph of structural function (Fig. 17). When comparing the bubble shapes, one has to consider the fact that at the end of the computational domain (L/Do20) the gas–liquid interface was still developing. Thus, only rough agreement was found. It was ascertained, that both mesh 2 and mesh 3 can predict flow patterns with a similar distribution of interfacial structures along the mini-channels. The mesh that was actually used has an average density between mesh 2 and mesh 3 and a further reduced aspect ratio (Table 3). 3.6. Simulation run times

Fig. 16. Initial condition geometry for cross-junction mixer: (a) real (experiment) geometry, (b) simplified geometry for numerical simulation.

mixing simulation is depicted in Fig. 16b. A uniform velocity profile was selected as a boundary condition type for both phases.

Each numerical simulation in this work represented about 0.03 s of real time and simulation required more than 1 month of computation time on the 24-core Supermicro HPC server. The large computational time is due to the small time steps that assured good temporal resolution and numerical stability. Consequently, achieving statistically relevant sample lengths proved to be a very demanding task to accomplish.

4. Results and discussion 3.5. Mesh independence study To assess mesh independence three different meshes were constructed with the same meshing strategy as described in Section 3.1.2. The computational domain length was further decreased to L/D ¼20 to reduce mesh size and thus, the computational time. Key mesh quality parameters, as well as average mesh densities are presented in Table 3. Simulation was run for each mesh using the same solver settings, models, boundary and initial conditions. Our aim was to find the mesh density at which the prediction of gas–liquid interfacial flow structures along the

Our main task was to predict relevant flow patterns by numerical simulation. The numerical algorithm for inlet conditions (as presented Section 3.2) can be successfully employed on different computational domain geometries and is capable of predicting relevant flow patterns over a wide range of inlet flow rates. Fig. 18 shows the three main flow patterns: bubbly flow (a), slug flow (b) and semi-annular flow (c) that are created at different air and water flow rates. Experimental snapshots in Fig. 18a–c were obtained from the mini header of a semi-circular cross-section. As shown in our publication on phase separation

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

116

Table 3 Key parameters for used meshes. Name

Average density (Cells/mm3)

Minimum aspect ratio

Maximum aspect ratio

Maximum skewness

Minimum cell volume (μm3)

Maximum cell volume (μm3)

Mesh 1 Mesh 2 Mesh 3 Used

25,792 30,951 46,702 36,105

3 3 3 2

24 18 21 16

0.4 0.4 0.4 0.4

18,585 15,487 10,133 13,275

56,168 46,806 31,111 40,120

Fig. 17. Structural function in last 7D of computational domain for different mesh densities overlaid by snapshots of interfacial structures represented by the isosurface of volume fraction 0.5.

Fig. 19. Flow pattern map for mini manifold header (Zun et al., 2010); larger symbols represent numerical simulations.

studies (Zun et al., 2010), header inlet conditions play a crucial role in determining flow patterns and consequently, phase separation. In order to successfully predict phase separation it is absolutely necessary to predict relevant flow patterns. Fig. 18d–f presents snapshots of numerically predicted interfacial structures (isosurfaces at volume fraction 0.5) using the porous media mixing algorithm with air and water flow rates that correspond to those used in the experiment. Fig. 19 presents the flow pattern map for the mini manifold header, based on experimental data. It can be seen that numerical simulations fall into the correct regions of the flow pattern map. Based on these results we decided to focus our attention on inlet conditions modeling which is the subject of the present paper. 4.1. Two-phase flow structures: numerical prediction vs. experiment

Fig. 18. Snapshots of interfacial structures in mini manifold header of semi-circular cross-section: bubbly flow (QL ¼ 12 l/h, QG ¼2 l/h for experiment (a) and QL ¼12 l/h, QG ¼ 2.2 l/h for numerical simulation (d)), slug flow (QL ¼1.2 l/h, QG ¼3 l/h for experiment (b) and QL ¼ 2 l/h, QG ¼ 3.5 l/h for numerical simulation (e)) and semiannular flow (QL ¼1.2 l/h, QG ¼ 50 l/h for experiment (c) and QL ¼ 1.4 l/h, QG ¼51 l/h for numerical simulation (f)). Numerical interface is represented by the iso-surface of volume fraction 0.5.

Four successful simulations of two phase flow in a 1.2 mm circular mini-channel are presented. Successful validation implies that we are being able to reflect understanding of two-phase flow pattern transition in a generalized, overall computational framework that could be eventually used for increasingly more complex predictive purposes. Ultimately, we can envision a gradual buildup in predictive capability of interfacial property changes, though this is just the first step in this direction at the moment. Fig. 20 presents an example of bubbly flow comparison between numerical simulation (a) and the experiment (b) at QL ¼12 l/h and QG ¼2 l/h, with the porous media mixer. Actual positions at which snapshots were taken are 20D downstream from the mixer for the numerical simulation and 125D downstream from the mixer for the experiment, respectively. Visual comparison reveals comparable interfacial structures. Smaller

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 20. Comparison of gas flow structures; simulation 1 (a) at QL ¼12 l/h, QG ¼2 l/h vs. experiment (b) at QL ¼ 12 l/h, QG ¼ 2 l/h. Numerical interface is represented by the iso-surface of volume fraction 0.5.

Fig. 21. Relative frequency of specific normalized equivalent bubble diameter class; simulation 1 vs. experiment; porous media mixing; skewed diamond symbol represents actual class width.

bubbles develop an irregular shape while larger ones (still deq oDch) develop a bullet like shape that is characteristic for elongated bubbles at higher Reynolds numbers (Kreutzer et al., 2005; Chen et al., 2009). Distributions of equivalent bubble diameters for the experiment and numerical simulation are plotted in Fig. 21. Here, experimental data are marked by black square symbols, while gray diamond symbols are used for numerically predicted data. It is apparent that in both cases bubble sizes remain smaller than the mini-tube diameter. Based on snapshots and equivalent bubble diameter comparison one may conclude that the flow pattern prediction is successful. Fig. 22 presents a comparison of simulated structures (a) at QL ¼ 8 l/h and QG ¼ 1.6 l/h with those obtained by the experiment (b) at QL ¼ 8 l/h and QG ¼ 2 l/h for the porous media mixer. While the bubble shape generation is again comparable, the distribution of equivalent bubble diameters (Fig. 23) indicates the onset of bubble size generation larger than the mini-tube diameter. Such bubbles do not appear during numerical simulation due to a slightly lower gas flow rate. Similar to the previous comparison, the number of the smallest bubbles is under predicted, whereas bubbles with 0.4 odeq/Do0.6 are over predicted. The degree of overall agreement between the prediction and measurement is quite good, even though flow pattern prediction is not entirely correct. Fig. 24 presents a comparison of simulated structures (a) at QL ¼6 l/h and QG ¼ 2 l/h with those obtained by the experiment (b) at QL ¼5 l/h and QG ¼2 l/h for the porous media mixer. Discrepancy in the liquid flow rate is not negligible in this case and has to be considered when analyzing the data. Due to large computational costs and time constrains we did not repeat the numerical

117

Fig. 22. Comparison of gas flow structures; simulation 2 (a) at QL ¼8 l/h, QG ¼ 1.6 l/h vs. experiment (b) at QL ¼ 8 l/h, QG ¼2 l/h. Numerical interface is represented by the iso-surface of volume fraction 0.5.

Fig. 23. Relative frequency of specific normalized equivalent bubble diameter class; simulation 2 vs. experiment; porous media mixing; skewed diamond symbol represents actual class width.

Fig. 24. Comparison of gas flow structures; simulation 3 (a) at QL ¼6 l/h, QG ¼ 2 l/h vs. experiment (b) at QL ¼ 5 l/h, QG ¼2 l/h. Numerical interface is represented by the iso-surface of volume fraction 0.5.

run at QL ¼6 l/h. Liquid slugs between gas structures are larger in the numerically predicted data. This could be due to a larger liquid flow rate. It may still be concluded, however, that bubbles and bullet shape gas structures can be observed in both cases. In the case of elongated bubbles, stronger interfacial oscillations were observed in the experiment, especially at the bubble rear end. The distribution of equivalent bubble diameters (Fig. 25) shows the existence of a bubble to slug flow transition (BTS) flow pattern in both cases. The numerically predicted distribution is less dispersed. However, the distribution peak value is located in the same size class as experimentally measured and the overall distribution trend is quite close.

118

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

Fig. 25. Relative frequency of specific normalized equivalent bubble diameter class; simulation 3 vs. experiment; porous media mixing; skewed diamond symbol represents actual class width.

Fig. 26. Comparison of gas flow structures; simulation 4 (a) at QL ¼ 8 l/h, QG ¼ 1.79 l/ h vs. experiment (b) at QL ¼8 l/h, QG ¼2 l/h. Numerical interface is represented by the iso-surface of volume fraction 0.5.

For the case of the cross-junction mixing, only a visual comparison is made due to the much longer computational time required. The purpose of this test was only to demonstrate that the employed numerical methods can be used for different inlet configurations. Fig. 26 reveals that numerical prediction is comparable with the experimental snapshot when comparing bubble distribution along mini-tube as well as bubble sizes.

At the initial condition modeling stage crude assumptions were made in the cases of both mixers. These assumptions are rather extreme, as they distort the true physics that govern the mixing of the phases. Our intention is not detailed simulation of mixing itself; instead we are focusing on modeling/mimicking the mixing in a way that produces results comparable with experimental observations of two-phase flow patterns in downstream pipe region. The interfacial structures that were predicted by the numerical simulations are also affected by temporal and spatial evolution that differs from what is observed experimentally, even with the current state of the art numerics. While these crude assumptions distort the details in the mixing region, it can be argued that they do not drastically affect the flow further downstream. Despite the short simulation time, simplifications at geometry discretization and assumptions at the inlet conditions modeling stage, the porous media mixer and the cross-junction mixer still produce sound numerical prediction. In the majority of the presented cases one can notice discrepancies in shape when comparing interfacial structures gained by experiment and numerical simulation. The appearance of ribs or oscillations of the bubble’s rear end, which were commonly observed in the experiment, are not well predicted by present numerical simulations. This is due to the fact that very fine mesh is needed to accurately capture phenomena on the interface, leading to a very large computational domain. Additionally, complex geometry is usually considered in real engineering applications, which makes the construction of “perfect” orthogonal mesh, as proposed by software developers, a difficult task to accomplish. A full 3D approach has to be adopted to cover all geometrical and phenomenological aspects of the problem considered. Finally, the poor accuracy of curvature estimation by Young’s algorithm (Youngs, 1982) leads to the creation of parasitic currents across the gas–liquid interface (Lafaurie et al., 1994). Cummins et al. (2005) proposed employing a different algorithm for curvature estimation. Magnini et al. (2013) for example, substantially reduced parasitic currents by estimating curvature with the “Height Function” algorithm. Complete elimination of parasitic currents is, however, practically not possible, as shown by Mencinger and Žun (2007). Due to these setbacks we focused more on the prediction of flow patterns and disregarded the details of bubble shape as long as we were able to achieve a general resemblance of gas structures.

5. Conclusions 4.2. Discussion on numerical capabilities Numerical simulations applied to cross-junction mixer reveal comparable results to experimental observations of flow pattern formation. Distributions of equivalent bubble diameter are, in the case of numerical simulations, applying to porous media mixer, narrower and show that the number of the largest and smallest bubbles is under predicted. These differences are greatly emphasized with the onset of bubble to slug flow transition. Taking pore activity into account could improve the solution of numerical simulation by varying the mass flux at a specific inlet passage with respect to time. This would promote the formation of smaller bubbles. Pore activity would also affect the temporal fluctuation of bubble number density that was experimentally observed in the case of the porous media mixer. This behavior was not observed with the current implementation of the porous media mixer inlet condition. When discussing the distributions one has to take into account that the simulation time is too short to be statistically relevant.

The effects of inlet conditions on bubble to slug flow transition in mini geometries are considered in this work. An experimental investigation was done on a mini-tube of 1.2 mm ID and two different inlet mixers; the porous media mixer and the crossjunction mixer. Interfacial structures were analyzed using high speed video and digital image post processing. The commercial CFD code ANSYS Fluent with VOF model was used to predict experimentally observed two-phase flow patterns. Full 3D formulation was employed. Four successful simulations were considered at this stage; three with the porous media mixer inlet condition and one with the cross-junction mixer inlet condition. Validation against experimental data shows very promising results. The following conclusions can be drawn: 1. We have succeeded in developing a numerical algorithm of inlet conditions with continuous flow rates of both phases that produce flow patterns comparable with experimentally observed cases. 2. Two mixers were used as a benchmark test for inlet conditions: a porous media mixer and cross-junction mixer.

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

3. The inlet flow mixer affects the distribution of bubble size in both the temporal and spatial domains. 4. Inlet flow conditions affect the onset of the transition from bubbly to slug flow. 5. Comparison of distributions of equivalent bubble diameters obtained by numerical simulation and the experiment revealed reasonable agreement.

Nomenclature A D deq deq/Dch Dch ds dt,min 〈dc〉 〈dp〉 g ID j L Nc N/N0 OD p Q Re t v

bubble projection area (pixel) diameter (m) bubble projection equivalent diameter (pixel) normalized equivalent bubble diameter mini-channel projection diameter (pixel) diameter of solid spheres in packed mixing chamber (m) minimum throat diameter (m) average channel diameter (m) average pore diameter (m) acceleration of gravity (m/s2) inner diameter superficial phase velocity (m/s) length (m) number of passages pre unit area (m  2) relative frequency of events outer diameter pressure (Pa) volume flow rate (l/h) Reynolds number (¼ j ρ D/μ) time (s) velocity (m/s)

Greek symbols α ε κ μ ρ s τ

local phase fraction porosity of packed spheres curvature of interface (m-1) viscosity (Pa s) density (kg/m3) surface tension coefficient (N/m) viscous stress tensor (Pa)

Subscripts G L

gas liquid

Superscripts T

transpose

Acknowledgements The work presented in this paper was supported by the Slovenian Ministry of Higher Education and Science under Contract no. P2-0162. References ANSYS Fluent 14.0 Documentation, ANSYS Inc., Cotonsburg, PA, 2011. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. Journal of Computational Physics 100, 335–354. Chen, W.L., Twu, M.C., Pan, C., 2002. Gas–liquid two-phase flow in micro-channels. International Journal of Multiphase Flow 28, 1235–1247.

119

Chen, Y., Kulenovic, R., Mertz, R., 2009. Numerical study on the formation of Taylor bubbles in capillary tubes. International Journal of Thermal Sciences 48, 234–242. Choi, C., Yu, D.I., Kim, M., 2011. Surface wettability effect on flow pattern and pressure drop in adiabatic two-phase flows in rectangular microchannels with T-junction mixer. Experimental Thermal and Fluid Science 35, 1086–1096. Cubaud, T., Tatineni, M., Zhong, X., Ho, C.-M., 2005. Bubble dispenser in microfluidic devices. Physical Review E 72, 037302. Cummins, S.J., Francois, M.M., Kothe, D.B., 2005. Estimating curvature from volume fractions. Computers & Structures 83, 425–434. De Menech, M., Garstecki, P., Jousse, F., Stone, H.A., 2008. Transition from squeezing to dripping in a microfluidic T-shaped junction. Journal of Fluid Mechanics 595, 141–151. Goel, D., Buwa, V.V., 2009. Numerical simulations of bubble formation and rise in microchannels. Industrial & Engineering Chemistry Research 48, 8109–8120. Gupta, R., Fletcher, D.F., Haynes, B.S., 2009. On the CFD modelling of Taylor flow in microchannels. Chemical Engineering Science 64, 2941–2950. ImagePro Analyzer: 6.2, User Manual, Media Cybernetics Inc, Bethesda, MD, 2004. Jones, O.C., Zuber, N., 1975. The interrelation between void fraction fluctuations and flow patterns in two-phase flow. International Journal of Multiphase Flow 2, 273–306. Kreutzer, M.T., Kapteijn, F., Moulijn, J.A., Heiszwolf, J.J., 2005. Multiphase monolith reactors: chemical reaction engineering of segmented flow in microchannels. Chemical Engineering Science 60, 5895–5916. Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., Zanetti, G., 1994. Modelling merging and fragmentation in multiphase flows with SURFER. Journal of Computational Physics 113, 134–147. Lakehal, D., Larrignon, G., Narayanan, C., 2008. Computational heat transfer and two-phase flow topology in miniature tubes. Microfluidics and Nanofluidics 4, 261–271. Magnini, M., Pulvirenti, B., Thome, J.R., 2013. Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel. International Journal of Heat and Mass Transfer 59, 451–471. Mak, C.Y., Omebere-Iyari, N.K., Azzopardi, B.J., 2006. The split of vertical two-phase flow at a small diameter T-junction. Chemical Engineering Science 61, 6261–6272. Mandhane, J.M., Gregory, G.A., Aziz, K., 1974. A flow pattern map for gas–liquid flow in horizontal pipes. International Journal of Multiphase Flow 1, 537–553. Matsui, G., 1984. Identification of flow regimes in vertical gas liquid two-pahse flow using differential pressure fluctuations. International Journal of Multiphase Flow 10, 711–720. Mencinger, J., Žun, I., 2007. On the finite volume discretization of discontinuous body force field on collocated grid: Application to VOF method. Journal of Computational Physics 221, 524–538. Mishima, K., Ishii, M., 1983. Flow regime transition criteria for upward two-phase flow in vertical tubes. International Journal of Heat and Mass Transfer 27, 723–737. Ng, K.M., 1986. A model for flow regime transitions in cocurrent down-flow tricklebed reactors. AIChE Journal 32, 115–122. Rahim, E., Revellin, R., Thome, J., Bar-Cohen, A., 2011. Characterization and prediction of two-phase flow regimes in miniature tubes. International Journal of Multiphase Flow 37, 12–23. Rek, Z., Žun, I., 2013. CFD based mini vs. micro system delineation in elongated bubble flow regime. International Journal of Multiphase Flow. Sturges, H.A., 1926. The choice of a class interval. Journal of the American Statistical Association 21, 65–66. Sur, A., Liu, D., 2012. Adiabatic air–water two-phase flow in circular microchannels. International Journal of Thermal Sciences 53, 18–34. Taitel, Y., Dukler, A.E., 1976. A model for predicting flow regime transitions in horizontal and near horizontal gas–liquid flow. AIChE Journal 22, 47–55. Talimi, V., Muzychka, Y.S., Kocabiyik, S., 2012. A review on numerical studies of slug flow hydrodynamics and heat transfer in microtubes and microchannels. International Journal of Multiphase Flow 39, 88–104. Thome, J.R., Bar-Cohen, A., Revellin, R., Zun, I., 2013. Unified mechanistic multiscale mapping of two-phase flow patterns in microchannels. Experimental Thermal and Fluid Science 44, 1–22. Triplett, K.A., Ghiaasiaan, S.M., Abdel-Khalik, S.I., Sadowski, D.L., 1999. Gas–liquid two-phase flow in microchannels Part I: two-phase flow patterns. International Journal of Multiphase Flow 25, 377–394. Ullmann, A., Brauner, N., 2007. The prediction of flow pattern maps in minichannels. Multiphase Science and Technology 19, 49–73. Wang, S.F., Mosdorf, R., Shoji, M., 2003. Nonlinear anlysis on fluctuation feature of two-phase flow through a T-junction. International Journal of Heat and Mass Transfer 46, 1519–1528. Wang, S.F., He, K., Huang, J.Z., 2011. Phase splitting of a slug-annular flow at a horizontal micro-T-junction. International Journal of Heat and Mass Transfer 54, 589–596. Xie, T, Ghiaasiaan, S.M., Karrila, S, 2004. Artifcial neural network approach for flow regime classification in gas–liquid-fiber flows based on frequency domain analysis of pressure signals. Chemical Engineering Science 59, 2241–2251. Youngs, D.L., 1982. Time-dependent multi-material flow with large fluid distortion. In: Morton, K.W., Baines, M.J. (Eds.), Numerical Methods for Fluid Dynamics, pp. 273–285. Zhang, T., Cao, B., Fan, Y., Gonthier, Y., Luo, L., Wang, S., 2011. Gas–liquid flow in circular microchannel. Part I: Influence of liquid physical properties and

120

J. Gregorc, I. Žun / Chemical Engineering Science 102 (2013) 106–120

channel diameter on flow patterns. Chemical Engineering Science vol. 66, 5791–5803. Zhao, C.-X., Middelberg, A.P.J., 2011. Two-phase microfluidic flows. Chemical Engineering Science 66, 1394–1411. Zhao, T.S., Bi, Q.C., 2001. Co-current air-water two-phase flow patterns in vertical triangular microchannels. International Journal of Multiphase Flow 27, 765–782.

Zun, I., Gregorc, J., Perpar, M., 2010. Two-phase flow maldistribution in a minimanifold system. Multiphase Science and Technology 22, 95–114. Žun, I., Perpar, M., Filipič, B., 1997. Probe measurements of local carrying gas fraction in trickle-bed reactor. Experimental Thermal and Fluid Science 15, 163–173.