1 Introduction

A gas–liquid flow moving through one or several small channels (diameters are typically <1 mm) is often encountered in micro-reactors and in micro-scaled heat transfer equipment. An important design parameter for these types of devices is the pressure drop and not only for dimensioning of compressors and pumps and operation of the device. If the pressure drop is a significant fraction of the channel inlet pressure, then gas phase expansion along the channel length is not negligible. The values of all parameters depending on the volumetric gas flow rate, such as mass and heat transfer coefficients, are therefore also a function of the axial position in the channel. Therefore, the pressure profile needs to be known in order to properly determine the performance of the device.

Most experimental work available in the literature concerning pressure drop of gas–liquid flows in small channels covers a span of flow patterns. It should be noted that, in this context, small channels are defined as channels in which surface tension forces dominate gravitational forces (Angeli and Gavriilidis 2008). The pressure drop is usually described with a flow pattern independent correlation based on either a Lockhart and Martinelli (1949) or a homogeneous flow type of approach (Triplett et al. 1999; Kawahara et al. 2002; Chung and Kawaji 2004). The accuracy with which these correlations predict experimental results is usually not satisfactory and becomes worse when the correlation found by one author is used to predict the data of another. One reason for the failure of the Lockhart–Martinelli and homogeneous flow correlations at predicting the pressure drop of gas–liquid flows in small channels is the absence of the effect of surface tension on the pressure drop. Contrary to larger channels, for which these correlations were developed, the effect of surface tension on the pressure drop is not negligible compared to the influence of inertial and viscous forces. Chen et al. (2002) therefore developed a homogeneous flow type correlation for channels with a diameter smaller than 10 mm that includes surface tension effects in the form of Bond and Weber numbers. They collected 11 data sets from literature and compared their correlation and several other homogeneous flow and Lockhart–Martinelli type correlations to the data sets and none of them was capable of accurately predicting the pressure drops. The failure of the Lockhart–Martinelli and homogeneous flow type models in predicting pressure drops of gas–liquid flows in small channels, even when including surface tension effects, suggests that models are needed which take more details of the hydrodynamics, such as the spatial gas–liquid distribution and the velocity distribution of both phases, into account. This also means that there will not be one correlation or model capable of describing the pressure drop for all possible flow patterns, but that a detailed pressure drop model has to be developed for each individual flow pattern. In this respect, it should be noted that Lockhart and Martinelli themselves state that “slug flow, in which alternate slugs of liquid and gas move down the tube, is eliminated from consideration” in the development of their model (Lockhart and Martinelli 1949).

The focus in this work will henceforth be on obtaining a model for the pressure drop of gas–liquid Taylor flow in channels with a diameter smaller than 1 mm. The reason for considering gas–liquid Taylor flow in favor of other regimes is that it is one of the main flow regimes of interest for performing gas–liquid and gas–liquid solid reactions in micro-structured reactors and monoliths (Kreutzer et al. 2005a; Hessel et al. 2005; Günther and Jensen 2006; Angeli and Gavriilidis 2008). The gas–liquid Taylor flow regime consists of an alternating sequence of gas bubbles and liquid slugs moving through a small channel. The length of the gas bubbles is larger than the channel diameter and a thin liquid film separates the gas bubbles from the channel walls. Furthermore, the liquid in the slugs forms recirculation cells for capillary numbers (Ca b = μ 1 u b /σ) smaller than 0.5 (Thulasidas et al. 1997; Kolb and Cerro 1993). This further improves radial mass transfer rates compared to those of mass transfer driven solely by diffusion because of convection of liquid to and from the channel wall at, respectively, the back and front of the gas bubbles.

2 Previous work on pressure drop in gas–liquid Taylor flow

2.1 The pressure drop over a single gas bubble

The thickness of the liquid film surrounding a Taylor gas bubble is a defining parameter for describing the hydrodynamics of gas–liquid Taylor flow. It has been the subject of study since Fairbrother and Stubbs (1935) first showed that a single bubble in a small diameter tube moves faster than the average velocity in the tube. The liquid film around the gas bubble causes the bubble to move through a smaller cross-section than that of the channel. Fairbrother and Stubbs found the excess velocity of the bubble to be a function of the capillary number.

Bretherton (1961) expanded on this work and analytically derived an expression for the liquid film thickness in channels with a circular cross-section as a function of the capillary number. Bretherton’s analysis is valid for very small liquid film thickness d f and in absence of significant inertial and gravitational forces, i.e. Ca b → 0 and We b = ρ l D h u 2b /σ ≪ 1, and results in:

$$ \frac{{d_{\text{f}} }}{{D_{\text{c}} }} = 0.67Ca_{\text{b}}^{{\frac{ 2}{ 3}}} . $$
(1)

In the same work, Bretherton also presented an expression for the pressure drop over a single Taylor gas bubble moving through a liquid-filled channel with a circular cross-section:

$$ \Updelta P_{\text{b}} = 7.16\left( {3Ca_{\text{b}} } \right)^{{\frac{2}{3}}} \frac{\sigma }{{D_{\text{c}} }} $$
(2)

Several authors making use of finite element calculations have confirmed the results from Bretherton’s analysis for both the liquid film thickness and the pressure drop over a Taylor gas bubble for low capillary numbers (Ca b < 5 × 10−3). For larger capillary numbers, where the liquid film occupies a significant portion of the channel (d f > 0.02D c), deviations from Bretherton’s equations are observed, which was already indicated by Bretherton himself (Bretherton 1961). A few authors have also considered inertial effects in their finite element calculations, indeed finding deviations from Bretherton’s equations as the capillary number increases (Giavedoni and Saita 1997; Ratulowski and Chang 1989; Westborg and Hassager 1989; Heil 2001; Fujioka and Grotberg 2005). However, none of these authors suggests new expressions for predicting the liquid film thickness and the pressure drop over a Taylor gas bubble when the liquid film thickness occupies a significant portion of the channel and/or inertial effects are significant. Aussillous and Quéré (2000) have applied a scaling analysis to Bretherton’s result for the liquid film thickness to account for both effects. They, however, did not discuss the implications of their work for the pressure drop over a Taylor gas bubble. They also did not quantify their results other than the effect of non-negligible liquid film thickness, for which they fitted the result of their scaling analysis to Taylor’s data (Taylor 1961) to obtain:

$$ \frac{{d_{\text{f}} }}{{D_{\text{c}} }} = \frac{{0.67Ca_{\text{b}}^{{\frac{ 2}{ 3}}} }}{{1 + 3.34Ca_{\text{b}}^{{\frac{ 2}{ 3}}} }} $$
(3)

This result is consistent with Bretherton’s result, since Eq. 3 approaches Eq. 1 for Ca b → 0. They found the range of maximum Ca b, for which, Eq. 3 described their experimental results to scale with (Ca b /Re b)0.75 . Although the exact values of Ca b and Re b for which Eq. 3 can be applied thus depend on the physical properties of the liquid and the capillary radius, the equation was capable of describing all their experimental result for Ca b < 0.01 and Re b < 150.

2.2 Pressure drop model by Kreutzer et al. (2005b)

Bretherton’s (1961) work is rarely mentioned by authors who experimentally studied the pressure drop of gas–liquid flows in small channels, despite its success. On the other hand, the authors studying the problem using numerical methods do not compare their simulation results to experimental results, other than those of Bretherton (1961) and Taylor (1961) for single bubbles.

A notable exception to the former statement is the work of Kreutzer et al. (2005b), who performed both numerical calculations including inertial effects as well as experiments on gas–liquid Taylor flow and compared their results to Bretherton’s work. They identified two sources determining the pressure drop in a gas–liquid Taylor flow: (1) frictional pressure drop in the liquid slugs, and (2) the pressure drop over the gas bubbles due to effects near the caps of the bubbles, for which the limit for Ca b → 0 was described by Bretherton. The pressure drop due to frictional losses in the gas bubbles was not taken into account, because it is negligible when compared to the frictional losses in the liquid due to the large difference in viscosity of gas and liquid. In their experiments, the Bond number (Bo = ρ l gD 2c /σ) was in the order of 1 and gravity effects could therefore be neglected (Edvinsson and Irandoust 1996; Hazel and Heil 2002). Kreutzer et al. did not include gravity effects in their numerical calculations and used a no-shear boundary condition at the gas–liquid interface. Due to the absence of gravity effects and shear of the gas bubble, the liquid film surrounding the gas bubbles is stagnant and does not contribute to the pressure drop.

Kreutzer et al. consider the liquid flow in the slugs to be a fully developed Hagen–Poiseuille flow, which is disturbed by the gas bubbles, causing an excess pressure drop. The pressure drop for a fully developed Hagen–Poiseuille liquid flow in a pipe is given by:

$$ - \frac{{{\text{d}}P}}{{{\text{d}}z}} = f_{\text{Fanning}} \left( {\frac{1}{2}\rho U^{2} } \right)\frac{4}{{D_{c} }} $$
(4)

where the Fanning friction factor f Fanning is equal to 16/Re for pipes with a circular cross section and Re = ρD c U/μ. Kreutzer et al. account for effects near the bubble caps by adding an extra term to the friction factor, which is inversely proportional to the liquid slug length L s and further depends on the physical properties of the liquid and the channel diameter. Their modified friction factor for the flow in the liquid slugs f s is of the form:

$$ f_{\text{s}} = \frac{16}{{Re_{\text{gl}} }}\left( {1 + a\frac{{D_{\text{c}} }}{{L_{\text{s}} }}\left( {\frac{{Re_{\text{gl}} }}{{Ca_{\text{gl}} }}} \right)^{b} } \right) $$
(5)

in which a and b are constants. The Fanning friction factor f Fanning in Eq. 4 is then replaced with the liquid slug friction factor f s as defined in Eq. 5. The average velocity of the liquid in the slug is equal to the average cross-sectional velocity in the channel U g + U l. The equation for the pressure drop over the liquid slugs then becomes:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{s}} = f_{\text{s}} \left( {\frac{1}{2}\rho_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)^{2} } \right)\frac{4}{{D_{\text{c}} }} = \frac{{2f_{\text{s}} Re_{\text{gl}} \mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }} $$
(6)

where the Reynolds number is now given by Re gl = ρ l D c(U g + U l)/μ l. Footnote 1 It should be noted that f s Re gl in the right-hand term in Eq. 6 is a constant not depending on any velocity and equals 16 for a fully developed single phase Hagen–Poiseuille flow in a pipe with a circular cross-section.

Since only a part of the channel length is occupied by liquid slugs, the pressure drop over the liquid slugs must be multiplied with the fraction of channel length occupied by the slugs, L s/(L b + L s), to obtain the pressure drop over the channel (dP/dz)c:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{s}} \left( {\frac{{L_{\text{s}} }}{{L_{\text{b}} + L_{\text{s}} }}} \right) $$
(7)

Kreutzer et al. performed both numerical calculations at various flow conditions and experiments at varying gas and liquid velocities using several liquids. They fitted the constants a and b in Eq. 5 to both data sets and the results are given in Table 1. The pressure drop correlation fitted to their experiments predicts larger pressure drops for a given Taylor flow than the correlation fitted to their CFD simulations. They attribute the difference to Marangoni effects, where impurities in the gas and liquid may cause the gas–liquid surface to “harden” and the no-shear boundary condition at the gas–liquid interface should be replaced by a no-slip boundary condition. Kreutzer’s experiments were performed for approximately 3 × 10−3 < Ca gl < 4 × 10−2, 1.5 × 102 < Re gl < 1.4 × 103, and 1 < We gl < 15. The value of constant a in Eq. 5 is either 0.17 or 0.07 depending on whether impurities are present or not. For Reynolds numbers close to unity, they suggest using Eq. 8 for the liquid slug friction factor instead of Eq. 5, which they obtained from rewriting Eq. 2 and further using Eqs. 6 and 7.

$$ f_{\text{s}} = \frac{16}{{Re_{\text{gl}} }} \left( {1 + \frac{{7.16D_{\text{c}} \left( {3Ca_{\text{b}} } \right)^{{\frac{2}{3}}} }}{{32L_{\text{s}} Ca_{\text{b}} }}} \right). $$
(8)

It should be noted that this result is based on Bretherton’s analysis and is, therefore, only valid for very thin liquid films and thus for Ca b < 5 × 10−3 and We b ≪ 1. Kreutzer et al. compared the pressure drop over a single Taylor gas bubble obtained from numerical calculations for Re gl = 1 to the values predicted using Eq. 2 and found a good agreement for 2 × 10−3 ≤ Ca gl ≤ 4 × 10−2. This suggests that the applicability of Eq. 8 extends to Ca gl ≤ 4 × 10−2, however, no experiments were performed in order to validate this.

Table 1 Values of the constants a and b in Eq. 5 in the model of Kreutzer et al. (2005b) from both their experiments and numerical calculations
Full size table

3 Motivation and scope of this work

Lockhart and Martinelli (1949) and homogeneous flow type models fail to accurately predict the pressure drop of gas–liquid Taylor flow in small channels, because surface tension effects and/or hydrodynamic details of the flow, such as the spatial gas–liquid distribution and the velocity distributions of both phases, are not taken into account. Kreutzer et al. (2005b) presented a semi-empirical model that successfully described their data obtained from experiments for Re b  > 150. Lower Reynolds numbers can easily be achieved in channels with a diameter in the order of 100 μm or smaller. For Reynolds numbers in the order of 1, Kreutzer et al. presented a model for the pressure drop of gas–liquid Taylor flows based on Bretherton’s (1961) results for the pressure drop over a single gas bubble. However, the model was neither compared to experimental results nor to numerical calculations and no suggestions are made for a pressure drop model for Reynolds numbers in the range 1 < Re b < 150.

Therefore, in this work, a model is developed for describing the pressure drop of gas–liquid Taylor flow with negligible liquid film velocities for Re b < 150 and Ca b < 0.01 using the approach of Kreutzer et al. The extension of Bretherton’s analysis by Aussillous and Quéré (2000) is used to obtain an expression for the pressure drop over a single gas bubble accounting for non-negligible liquid film thickness. A mass balance-based model for gas–liquid Taylor flow previously developed by the authors (Warnier et al. 2007) is used to describe the fraction of channel length occupied by liquid having a non-zero velocity causing the liquid frictional pressure drop. The mass balance based model is then also used to obtain the pressure drop per unit channel length from the liquid frictional pressure drop and the pressure drop over a single gas bubble. The pressure drop model is then compared to data obtained from experiments with nitrogen/water Taylor flow in a glass capillary with a circular cross-section and a diameter of 250 μm. The results from these experiments are also compared to model predictions using the model by Kreutzer et al. and the Lockhart–Martinelli–Chisholm (Lockhart and Martinelli 1949; Chisholm 1967) correlation (the model is described in Appendix A) to determine if the newly developed model truly is a significant improvement compared to currently available models and often used correlations. Even though gas–liquid Taylor flow was not considered in the development of the Lockhart–Martinelli–Chisholm model (Lockhart and Martinelli 1949), the model has been applied to the Taylor flow regime by various authors and is, therefore, considered in this work. The values of various hydrodynamic parameters are obtained from analysis of images of the flows. Since the pressure drop model is partly based on the mass balance-based model, the latter is briefly discussed first.

4 Mass balance-based model

Consider a gas–liquid Taylor flow moving through a channel with a cross-sectional area A. The gas bubbles have a velocity u b, a length L b and occupy a fraction of the cross-sectional area of the channel A b/A. The liquid slugs have a length L s. The flow is divided into unit cells consisting of one gas bubble, its surrounding liquid film, and one liquid slug and the unit cell length is L b + L s (Fig. 1).

Fig. 1
figure 1

Schematic of Taylor flow showing the definitions of the unit cell, gas bubble length L b and the liquid slug length L s. The lengths of the nose L nose and tail L tail sections of the gas bubble are also indicated

Full size image

Continuity requires that the overall average velocity through any cross-section of the channel perpendicular to the direction of flow is equal to the sum of the superficial gas U g and liquid U l velocities, which are based on the channel cross-section A. The flow through plane A 1 consists of the gas bubble moving with velocity u b through a cross-section A b, and the liquid film moving at an average velocity u f through a cross-section A − A b, giving:

$$ \frac{{A_{\text{b}} }}{A}u_{\text{b}} + \left( {1 - \frac{{A_{\text{b}} }}{A}} \right)u_{\text{f}} = U_{\text{g}} + U_{\text{l}} $$
(9)

The flow through plane A 2 consists solely of liquid, which occupies the whole cross-section of the channel A, and the average velocity of the liquid in the slug is therefore equal to U g + U l.

The superficial gas velocity U g is equal to the bubble volume times the bubble formation frequency F b, divided by the channel cross-sectional area A. If the bubble volume is taken to be the bubble length L b times its cross-sectional area A b, then it is overestimated because part of the volume A b(L nose + L tail) consists of liquid. This liquid volume depends on the bubble geometry, but can arbitrarily be written as the bubble cross-sectional area A b times a length δ. This length δ is then a correction on the bubble length to account for the overestimation of the gas bubble volume if it were taken to be equal to A b L b. In case of a bubble with a circular cross-section and hemispherical bubble caps, it can be shown that δ is 1/3 of the bubble diameter D b. The superficial gas velocity U g is then given by:

$$ U_{\text{g}} = \frac{{A_{\text{b}} }}{A}F_{\text{b}} \left( {L_{\text{b}} - \delta } \right). $$
(10)

The bubble formation frequency F b is equal to the number of unit cells passing a certain location in the channel per unit of time, which is the reciprocal of the time it takes for a unit cell to travel a distance equal to its own length: (L b + L s)/u b. The bubble formation frequency is therefore given by:

$$ F_{\text{b}} = \frac{{u_{\text{b}} }}{{L_{\text{b}} + L_{\text{s}} }} $$
(11)

If Eqs. 10 and 11 are substituted into Eq. 9, the following expression is obtained for the superficial liquid velocity U l:

$$ U_{\text{l}} = \frac{{A_{\text{b}} }}{A}F_{\text{b}} (L_{\text{s}} + \delta ) + \left( {1 - \frac{{A_{\text{b}} }}{A}} \right)u_{\text{f}} $$
(12)

The velocity of the liquid in the film surrounding the bubble u f is zero for horizontal Taylor flows as used in this work.

5 Pressure drop model

5.1 Model development

Aussillous and Quéré (2000) incorporated the effect of a non-negligible liquid film thickness in Bretherton’s (1961) analysis by replacing the term \( Ca_{\text{b}}^{{\frac{2}{3}}} \,{\text{with}}\,Ca_{\text{b}}^{{\frac{2}{3}}} /\left( {1 + 3.34Ca_{\text{b}}^{{\frac{2}{3}}} } \right). \) Bretherton’s result for the pressure drop over a single Taylor gas bubble can similarly be modified to account for a non-negligible liquid film thickness, giving:

$$ \Updelta P_{\text{b}} = 7.16\frac{{\sigma (3Ca_{\text{b}} )^{{\frac{2}{3}}} }}{{D_{\text{c}} \left( {1 + 3.34Ca_{\text{b}}^{{\frac{2}{3}}} } \right)}}. $$
(13)

The pressure drop over a unit cell, ΔP uc, is the sum of the frictional pressure drop of the liquid flowing in the slug, ΔP s , and the pressure drop over the Taylor gas bubble, ΔP b. The flow of liquid in the slug is assumed to be fully developed and is therefore given by the Hagen–Poiseuille equation. The deviations thereof near the bubble caps are exactly what are accounted for in Bretherton’s analysis and are thus incorporated by ΔP b. Continuity requires the average velocity in the liquid slug to be U g + U l. As shown in the derivation of the mass balance based model, to properly account for the total amount of liquid with an average velocity of U g + U l in one liquid slug, the liquid slug length is L s + δ. The frictional pressure loss in one liquid slug is then given by:

$$ \Updelta P_{\text{s}} = \left( {\frac{16}{{Re_{\text{gl}} }}} \right)\left( {\frac{4}{{D_{\text{c}} }}} \right)\left( {\frac{1}{2}\rho_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)^{2} } \right)(L_{\text{s}} + \delta ) = \frac{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }}(L_{\text{s}} + \delta ). $$
(14)

Note that the Reynolds number is based on the average velocity of the liquid U g + U l and is not equal to the bubble Reynolds number Re b, which is based on u b. The pressure drop over a unit cell is then given by:

$$ \Updelta P_{\text{uc}} = \Updelta P_{\text{s}} + \Updelta P_{\text{b}} = \frac{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }}(L_{\text{s}} + \delta ) + 7.16\frac{{\sigma \left( {3Ca_{\text{b}} } \right)^{{\frac{2}{3}}} }}{{D_{\text{c}} \left( {1 + 3.34Ca_{\text{b}}^{{\frac{2}{3}}} } \right)}}. $$
(15)

The length of a unit cell is L b + L s and thus the number of unit cells per unit length of channel is 1/(L b + L s). The pressure drop over a unit length of channel, −(dP/dz)c, is therefore given by:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = \frac{{\Updelta P_{\text{uc}} }}{{L_{\text{b}} + L_{\text{s}} }} = \frac{1}{{L_{\text{b}} + L_{\text{s}} }}\left( {\frac{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }}(L_{\text{s}} + \delta ) +7.16\frac{{\sigma \left( {3Ca_{\text{b}} } \right)^{{\frac{2}{3}}} }}{{D_{\text{c}} \left( {1 + 3.34Ca_{\text{b}}^{{\frac{2}{3}}} } \right)}}} \right) $$
(16)

Equation 16 can be rewritten to:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = \frac{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }}\frac{{\left( {L_{\text{s}} + \delta } \right)}}{{\left( {L_{\text{b}} + L_{\text{s}} } \right)}}\left( {1 + 7.16\frac{{D_{\text{c}} \sigma \left( {3Ca_{\text{b}} } \right)^{{\frac{2}{3}}} }}{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)(L_{\text{s}} + \delta )\left( {1 + 3.34Ca_{\text{b}}^{{\frac{2}{3}}} } \right)}}} \right) $$
(17)

This result can be interpreted as follows: the term 32 μ l(U g + U l)/D 2c is the pressure drop per unit length of liquid for a fully developed Hagen–Poiseuille liquid flow moving at the average velocity of U g + U l. The term (L s + δ)/(L b + L s) is the volumetric fraction of channel, in which liquid is moving at an average velocity of U g + U l. Therefore, these two terms combined are the pressure drop per unit channel length based on a fully developed Hagen–Poiseuille flow in the liquid slugs. The remaining term in Eq. 17 accounts for the additional pressure drop over the gas bubbles per unit channel length.

It can be shown from the mass balance-based model that (L s + δ)/(L b + L s) can be rewritten to U l/(U g + U l). The part of Eq. 17 describing the frictional pressure drop over a unit channel length due to fully developed liquid flow can then be rewritten as follows:

$$ \frac{{32\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}}{{D_{\text{c}}^{2} }}\frac{{\left( {L_{\text{s}} + \delta } \right)}}{{\left( {L_{\text{b}} + L_{\text{s}} } \right)}} = \frac{{32\mu_{\text{l}} U_{\text{l}} }}{{D_{\text{c}}^{2} }} $$
(18)

This result indicates that any change of the frictional pressure drop over a unit length of channel due to changing the liquids average velocity by varying the gas flow rate, is negated by the change of the volumetric channel fraction occupied by that liquid. Furthermore, the group σ/(μ l(U g + U l)) in the third term on the right hand side of Eq. 17 is the inverse of the capillary number based on the average velocity in the channel, Ca gl. This can be rewritten to the inverse of the capillary number based on the bubble velocity Ca b using Eq. 9 for a stagnant liquid film:

$$ \frac{\sigma }{{\mu_{\text{l}} \left( {U_{\text{g}} + U_{\text{l}} } \right)}} = \frac{A}{{A_{\text{b}} }}\frac{1}{{Ca_{\text{b}} }} $$
(19)

By substituting Eqs. 18 and 19 into Eq. 17, the following expression for the frictional pressure drop over a unit length of channel is obtained:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = \frac{{32\mu_{\text{l}} U_{\text{l}} }}{{D_{\text{c}}^{2} }}\left( {1 + \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\frac{{D_{\text{c}} }}{{\left( {L_{\text{s}} + \delta } \right)}}\left( {\frac{A}{{A_{\text{b}} }}} \right)\frac{1}{{\left( {Ca_{\text{b}}^{{\frac{1}{3}}} + 3.34Ca_{\text{b}} } \right)}}} \right) $$
(20)

The liquid slug friction factor f s is then given by:

$$ f_{\text{s}} = \left( {\frac{16}{{Re_{\text{gl}} }}} \right)\left( {1 + \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\frac{{D_{\text{c}} }}{{\left( {L_{\text{s}} + \delta } \right)}}\left( {\frac{A}{{A_{\text{b}} }}} \right)\frac{1}{{\left( {Ca_{\text{b}}^{{\frac{1}{3}}} + 3.34Ca_{\text{b}} } \right)}}} \right) $$
(21)

The cross-sectional area of the bubble A b directly depends on the liquid film thickness. The corrected liquid slug length L s + δ also depends on the bubble diameter and thus on the liquid film thickness. Both parameters can therefore be expressed as a function of the bubble velocity u b using Aussillous and Quéré’s result given in Eq. 3 and this can be substituted into Eqs. 20 and 21. At first glance, this would result in the pressure drop being a complicated function of the capillary number and thus of the bubble velocity. However, the term (L s + δ)−1(A/A b) in Eqs. 20 and 21 can be rewritten as a function of the superficial liquid velocity U l and the bubble frequency F b using Eq. 12 for a stagnant liquid film:

$$ \frac{1}{{\left( {L_{\text{s}} + \delta } \right)}}\frac{A}{{A_{\text{b}} }} = \frac{{f_{\text{b}} }}{{U_{\text{l}} }} $$
(22)

Although it is not clear at first glance, the combined term (L s + δ) −1(A/A b) in Eqs. 20 and 21 is not a function of bubble velocity and thus is not dependent on the pressure and this simplifies the calculation of the pressure profiles. Substituting this result into Eqs. 20 and 21 gives:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = \frac{{32\mu_{\text{l}} U_{\text{l}} }}{{D_{\text{c}}^{2} }}\left( {1 + \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\frac{{D_{\text{c}} f_{\text{b}} }}{{U_{\text{l}} }}\frac{1}{{\left( {Ca_{\text{b}}^{{\frac{1}{3}}} + 3.34Ca_{\text{b}} } \right)}}} \right) $$
(23)

and

$$ f_{\text{s}} = \frac{16}{{Re_{\text{gl}} }}\left( {1 + \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\frac{{D_{\text{c}} f_{\text{b}} }}{{U_{\text{l}} }}\frac{1}{{\left( {Ca_{\text{b}}^{{\frac{1}{3}}} + 3.34Ca_{\text{b}} } \right)}}} \right), $$
(24)

respectively. Note that, apart from the constant 3.34, none of the other constants in Eqs. 23 and 24 was obtained from a fitting procedure. The constant 3.34 was obtained by Aussillous and Quéré (2000) from fitting their equation to the data of Taylor (1961).

5.2 Comparison with the model of Kreutzer et al. (2005b)

Kreutzer et al. (2005b) found an expression for the liquid slug friction factor based on Bretherton’s (1961) analysis as given in Eq. 8. Bretherton’s analysis is valid for vanishing liquid film thickness and thus for Ca b → 0. The liquid slug friction factor obtained in this work as given in Eq. 21 can also be evaluated for these conditions. For Ca b → 0, the liquid film thickness approaches zero and the bubble cross-sectional area A b approaches the channel cross-sectional area A. Furthermore, Ca 1/3b  ≫ 3.34Ca b and Eq. 21 then becomes:

$$ f_{\text{s}} = \frac{16}{Re}\left( {1 + \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\frac{{D_{\text{c}} }}{{\left( {L_{\text{s}} + \delta } \right)}}\frac{1}{{Ca_{\text{b}}^{{\frac{1}{3}}} }}} \right) $$
(25)

The only difference between Eqs. 25 and 8 is that in Kreutzer et al.’s result, using Bretherton’s analysis, the liquid slug length is given by L s and in our work it is defined as L s + δ. It is not exactly clear how Kreutzer et al. have defined their slug length. In this work, it is defined as the shortest distance between the tail section of a gas bubble and the nose section of the trailing bubble. Kreutzer et al. take the volumetric fraction of the channel of liquid moving at an average velocity of U g + U l to be L s/(L b + L s). With our definition of L s, the amount of liquid surrounding the nose and tail sections of the bubbles having a non-zero velocity is not taken into account. As explained in the description of the mass balanced-based model, this volume is accounted for when L s is replaced with L s + δ. Applying this correction to Eq. 8 results in it being identical to Eq. 25 and the model developed in this work is indeed in accordance with Bretherton’s work as applied by Kreutzer et al.

The same correction on the liquid slug length is applied to the model of Kreutzer et al. in order to be able to properly compare their results with those obtained in this work. Equations 6 and 7 can be combined and (L s + δ)/(L b + L s) is substituted for L s/(L b + L s), following the same line of reasoning as in the derivation of Eq. 18, to give:

$$ - \left( {\frac{{{\text{d}}P}}{{{\text{d}}z}}} \right)_{\text{c}} = \frac{{2\mu_{\text{l}} U_{\text{l}} }}{{D_{\text{c}}^{2} }}f_{\text{s}}{Re_{\text{gl}}} $$
(26)

Replacing L s with L s + δ in Eq. 5 then gives the following result for the slug friction factor f s:

$$ f_{\text{s}} = \frac{16}{{Re_{\text{gl}}}} \left( {1 + a\frac{{D_{\text{c}} }}{{L_{\text{s}} + \delta }}\left( {\frac{{Re_{\text{gl}} }}{{Ca_{\text{gl}} }}} \right)^{0.33} } \right) $$
(27)

These two equations will be used further throughout this work as the adapted model of Kreutzer et al.

6 Experimental

All experiments were carried out with nitrogen gas and demineralized water at a temperature of 20 ± 1°C. The gas flow was regulated by a mass flow controller (Bronkhorst F-200C) and the superficial gas velocity was varied between 0.06 and 0.33 m/s (these values are given at a temperature of 20°C and a pressure of 1 bar). The superficial liquid velocity was varied between 0.10 and 0.34 m/s using an Isco 100DM syringe pump. The ranges of Reynolds, capillary and Weber numbers used in this work are given in Table 2.

Table 2 Superficial gas U g and liquid U l velocities used in this work.
Full size table

The liquid flow was split into two equal flows, which were contacted with the nitrogen flow at an angle of 90° with respect to the direction of flow of the nitrogen in a stainless steel cross-mixer (Valco ZX.5) (Fig. 2). The inner diameter of the channels in the cross-mixer was 250 μm.

Fig. 2
figure 2

Schematic drawing of the cross-mixer. The liquid flow is split into two equal flows and contacted with the gas flow at an angle of 90°. All channels in the cross-mixer have a diameter of 250 µm

Full size image

A differential pressure sensor (Honeywell 26PC 05 KD, response time of 1 ms), measuring the pressure difference with ambient pressure, was connected at a distance of 13 cm downstream from the cross-mixer to ensure a fully developed Taylor flow passing the sensor. The flow then entered the glass capillary 3 cm downstream from the pressure sensor and exited the capillary under ambient conditions (T = 293 K and P exit = 1.03 × 105 Pa). The transparent, glass capillary had a length of 1 m and a circular cross-section with an inner diameter of 250 μm. All of the connectors and tubing had an inner diameter of 250 μm.

The recording time for each of the studied 45 combinations of gas and liquid velocities was 2 s and the recording of the pressure signal was synchronized with the recording of the images. A Redlake MotionPro CCD camera connected to a Zeiss Axiovert 200 MAT inverted microscope was used to record images of the Taylor flows at a frame rate of 2,500 frames per second. An exposure time of 20 μs was used, which was sufficiently short for preventing any significant motion blur. The length of one pixel represented 3.8 μm of channel length. A channel length of 3.04 mm was captured in the images of which the center point was located 29.5 cm downstream from the pressure sensor.

The gas bubbles in all the movies had a length shorter than the channel length captured in the images. The frame rate of 2,500 frames per second was sufficiently short to ensure that no gas bubble travelled a distance greater than the unit cell length between the recordings of two consecutive images so that the individual bubbles could be tracked. For every movie, the length of every individual bubble was averaged over all frames in which it was completely visible. These values were then averaged to obtain the average bubble length L b for that movie. The bubble frequency F b is the number of tracked bubbles in one movie divided by the measurement time. The average velocity of a single bubble was obtained by dividing the distance travelled by its centre of mass in the movie by the time that the bubble was present in that movie. The bubble velocity was first determined for every single bubble and then averaged over all bubbles in one movie to give the average bubble velocity u b. For each movie, the standard deviation of the bubble velocity was less than 0.5% of the average value, indicating that velocity fluctuations within the time of recording of one movie are negligible. The same procedure was followed for the liquid slugs for those movies in which the length of the liquid slugs was smaller than the channel length captured by the images. For those movies, in which the liquid slugs were not completely visible in one frame, Eq. 11 was used to calculate the slug length, since the bubble velocity, the bubble formation frequency, and the bubble length were known.

7 Results and discussion

In this section, for the sake of legibility, the adapted model of Kreutzer et al. (2005b) and that of Lockhart, Martinelli and Chisholm (Lockhart and Martinelli 1949; Chisholm 1967) will be referred to as, respectively, the “Kreutzer model” and the “LMC model”.

A summary of the values obtained for the various hydrodynamic parameters from image analysis is given in Table 3.

Table 3 Summary of the minimum and maximum values obtained for the various hydrodynamic parameters from image analysis
Full size table

The pressure difference with ambient conditions measured at the sensor location varied between 0.93 × 105 and 2.34 × 105 Pa. Since the pressure at the channel exit equals 1.03 × 105 Pa, the volumetric gas flow rate at the pressure sensor was between 0.53 and 0.31 times that at the channel exit, assuming ideal isothermal expansion. Expansion of the gas phase between the two locations is, therefore, significant and the superficial gas velocity and the gas bubble velocity are not constant over the channel length. This is especially important for the LMC model (see Appendix A) and the model developed in this work (see Eqs. 23, 24) since the pressure drop depends on, respectively, the superficial gas velocity and the gas bubble velocity. For these models, the pressure drop is a function of the axial position in the channel and the pressure is therefore not a linear function of the channel length. The pressure drop in the Kreutzer model, as described by Eqs. 26 and 27, does not directly depend on the superficial gas or on the bubble velocity. However, the correction on the liquid slug length δ is equal to D b/3 for gas bubbles with hemispherical caps, which is a valid assumption in this work, and therefore depends on the liquid film thickness. The liquid film thickness is a function of the bubble velocity as given in Eq. 3, and so Eq. 27 is a function of the bubble velocity.

In order to be able to properly compare experimental results and model predictions, the pressure difference between the sensor location and the channel exit is calculated and compared to the measurement of the differential pressure sensor. This procedure requires calculation of the pressure profiles rather than the pressure drop. The pressure profiles are calculated numerically using the assumptions that:

  • the gas phase obeys the ideal gas law,

  • gas expansion is isothermal,

  • absorption of the gas phase into the liquid has a negligible effect on the volumetric gas flow rate,

  • the liquid phase is incompressible.

The boundary condition for all models is that the pressure at the exit of the channel equals 1.03 × 105 Pa. The superficial gas velocity at the exit conditions is, therefore, also known. For the LMC model the calculations are then straightforward, since the pressure drop depends on the local superficial gas velocity U g (Appendix A), which can easily be related to the pressure using the ideal gas law. Calculating the pressure profile from the model presented in this work as well as from the Kreutzer model requires an additional step since the pressure drop depends on the bubble velocity u b and not directly on U g. The bubble velocity is related to the superficial gas velocity and the cross-sectional area of the gas bubble A b as described in Eq. 9 and thus depends on the liquid film thickness. The liquid film thickness depends on the bubble velocity as described by Aussillous and Quéré (Eq. 3). Therefore, the bubble velocity and the liquid film thickness are calculated by an iteration of Eqs. 9 and 3 for every axial position in the channel while calculating the pressure profile.

Kreutzer et al. (2005b) found two values for the constant a in Eq. 27: 0.17 and 0.07 based on, respectively, their experimental and numerical data. They explain the difference to be caused by impurities in the gas and/or liquid phases during the experiments. The experimental results obtained in this work are compared to the Kreutzer model for both values of a and the results are given in Fig. 3. The Kreutzer model was developed for Re gl > 150 and indeed fails to describe the experimental data for Re gl < 150 using either value for a. The model consistently over predicts the experimental results for a = 0.17 and under predicts them for a = 0.07. The accuracy of the model predictions is expected to increase when Re gl approaches and exceeds 150. The opposite trend is observed for a = 0.17, e.g. for increasing Re gl the deviation of the model predictions from the experimental values becomes larger. The expected trend is observed for a = 0.07, but the scatter is large. Therefore, the constant a was fitted to the data in this work in order to improve the performance of the model. The best agreement between the Kreutzer model and the experimental data is obtained for a = 0.10. Henceforth, this value of parameter a is used when comparing the Kreutzer model to the experimental data and the other two models.

Fig. 3
figure 3

The pressure difference between the channel exit and the sensor location, as calculated with the model by Kreutzer et al. (2005b) divided by the measured pressure difference is plotted as a function of the Reynolds number Re gl for three values of the constant a. The capillary number Ca gl varied between 2.3 × 10−3 and 8.8 × 10−3

Full size image

Figure 4 shows the calculated pressure difference between the sensor location and the channel exit versus the measured pressure difference for all three models. Figure 5 shows the same data, but then plotted as the ratio of the model predictions and the measured data versus Re gl in order to better visualize the differences between model and measured values.

Fig. 4
figure 4

The pressure difference between the channel exit and the sensor location as calculated with the three models versus the measured pressure difference

Full size image
Fig. 5
figure 5

The ratio between the pressure difference between the channel exit and the sensor location as calculated with the three models and the measured pressure difference versus Re gl. The capillary number Ca gl varied between 2.3 × 10−3 and 8.8 × 10−3

Full size image

The LMC model under predicts the measured data at the smallest measured pressure differences, while it over predicts the data at the largest measured pressure differences. With increasing measured pressure difference, the accuracy of the model first increases until the model predicts larger values than the measured data. The accuracy of the model then decreases with increasing measured pressure difference. The LMC model predicts values of −14 to +10% of the measured value. The Kreutzer model is more accurate than the LMC model and its predictions vary between −15 and +4% of the measured values. Contrary to the predictions of the LMC model, the accuracy of the predictions by the Kreutzer model increases with increasing Re gl, as expected. The model developed in this work predicts the pressure drop within a range of −4 to +3% of the experimental data, which is a significantly smaller range than those of both the Kreutzer and LMC models.

The inability of the LMC model to accurately describe the experimental results obtained in this work confirms that both surface tension and hydrodynamic properties of the two phase flow, such as the spatial gas–liquid distribution and the velocity distributions of both phases, have to be considered when developing a model for the pressure drop of gas–liquid Taylor flows in small channels. The Kreutzer model does include surface tension effects and considers hydrodynamic properties in the form of the liquid slug length L s, but it does not include the local velocities of the phases. Furthermore, it is a semi-empirical model and apart from the two values of the fitting parameter a reported by Kreutzer et al. themselves, a different value for a was found in this work. Kreutzer et al. state that they expect impurities in the gas and/or liquid phases in their experiments to be the cause of finding a different value for a than the value obtained from their numerical calculations. The value found for a in this work is in between the limiting values obtained by Kreutzer et al. If impurities are indeed the cause of these differences, then the value of a obtained in this work implies some concentration of impurities present in our experiments as well. These effects should then also be of significance on the accuracy of the model developed in this work. However, without considering possible impurities, the model developed in this work describes the experimental data within an error of ±4%. This suggests that no significant concentration of impurities was present during our experiments and Marangoni effects therefore cannot explain the different value of the fitting parameter a obtained in this work. It should be noted that the introduction of the group a(Re/Ca)0.33 in the Kreutzer model is based on their experimental observations, but physical interpretation of the group a(Re/Ca)0.33 is lacking. The parameter a may very well be a function of other hydrodynamic parameters causing the differences in values observed for a. In order to understand the nature of the fitting parameter a, the model developed in this work is compared in more detail with the Kreutzer model.

An expression for a is obtained from combining the expression for the liquid slug friction factor f s obtained by Kreutzer et al. as given in Eq. 27 with the expression for f s derived in this work as given in Eq. 21:

$$ a = \frac{{7.16 \times 3^{{\frac{2}{3}}} }}{32}\left( {\frac{A}{{A_{\text{b}} }}} \right)\left( {\frac{{Ca_{\text{gl}} }}{{Re_{\text{gl}} }}} \right)^{{\frac{1}{3}}} \left( {Ca_{\text{b}}^{{\frac{1}{3}}} + 3.34Ca_{\text{b}} } \right)^{ - 1} $$
(28)

The parameter a thus depends on the liquid properties and is velocity dependent. A/A b is a function of the liquid film thickness and thus of the bubble velocity following the equation of Aussillous and Quéré (Eq. 3). The values of a can be calculated for all experiments performed in this work using the bubble velocities obtained from image analysis and the local superficial gas velocities at the imaging location. The latter is known from the pressure profiles and the ideal gas law. Note that these values of a are estimates, since a is a function of the axial position in the channel and only the values at the imaging location are used. The result is shown in Fig. 6. When fitting the Kreutzer model to a data set, one value of a is obtained. It was already shown that the Kreutzer model described our experimental data best for a = 0.10. The data displayed in Fig. 6 are averaged in order to be able to make the comparison with the one value obtained in the fitting procedure. When doing so, a = 0.11 is obtained and this is, as expected, close to the value resulting from fitting the model to the experimental data. This result confirms that parameter a as described by Eq. 28 is a lumped parameter containing all velocity effects not explicitly taken into account in the Kreutzer model. The value found for a upon fitting the model to a data set depends on the exact combination of velocities covered in the experiments and the liquid properties and this can explain the differences between the value of a found in this work and the values found by Kreutzer et al. Furthermore, the variation of parameter a with varying bubble velocity decreases with increasing bubble velocity and thus with increasing Reynolds and capillary numbers. Although the capillary numbers used in this work are within the range of those used by Kreutzer et al., the Reynolds numbers used in their experiments are larger than those in this work and the influence of the bubble velocity on the value of parameter a is less in their experiments. This may be the reason why Kreutzer et al. observed no or little influence of the bubble velocity on the value of parameter a. However, for Re gl < 150 the effect of the bubble velocity on the additional pressure drop caused by the presence of the gas bubbles has to be taken into account.

Fig. 6
figure 6

Parameter a in the model of Kreutzer et al. (2005b) is calculated using Eq. 28 and plotted as a function of Re gl. The capillary number Ca gl varied between 2.3 × 10−3 and 8.8 × 10−3

Full size image

8 Conclusions

In this paper, a model has been developed to describe the pressure drop of gas–liquid Taylor flows in small round channels (diameter typically smaller than 1 mm). The model takes two sources of pressure drop into account: (1) frictional pressure drop caused by laminar flow in the liquid slugs, and (ii) an additional pressure drop over a single gas bubble due to the bubble disturbing the otherwise parabolic velocity profile in the liquid slugs. The pressure drop over a single gas bubble is described by the classical theory of Bretherton (1961), of which, the applicability is extended up to a capillary number Ca gl of 0.01 and Reynolds numbers Re gl in the order of 102 using the scaling analysis of Aussillous and Quéré (2000). The pressure drop per unit channel length is obtained by combining the contributions of both sources using a mass balance based Taylor flow model previously developed by the authors (Warnier et al. 2007). The model shows that (1) the additional pressure drop caused by the presence of the gas bubbles depends on the bubble frequency and the bubble velocity, and (2) its contribution to the overall pressure drop relative to that of the frictional pressure drop in the liquid phase decreases with increasing gas bubble velocity.

The model is compared to experimental data obtained for nitrogen–water Taylor flow in a glass channel with a circular cross-section and an inner diameter of 250 μm. The Reynolds number Re gl is varied between 41 and 159 and the capillary number Ca gl varied between 2.3 × 10−3 and 8.8 × 10−3. The model developed in this work describes the experimental data within an accuracy of ±4% without fitting any of the constants to the data obtained in this work.

Both the model and the experimental data are compared to the Lockhart–Martinelli–Chisholm model (Lockhart and Martinelli 1949; Chisholm 1967) and the model of Kreutzer et al. (2005b). The model of Kreutzer et al. is a semi-empirical model and is fitted to the experimental data obtained in this work and the fitted model predicts the measured pressure drop with increasing accuracy as the Reynolds number increases to 150. The Lockhart–Martinelli–Chisholm model was not developed for gas–liquid Taylor flow, which was explicitly indicated by the authors themselves (Lockhart and Martinelli 1949). Their model is, therefore, not expected to accurately describe the pressure drop in gas–liquid Taylor flow, yet it has been applied to this flow regime by other authors and was therefore, included in this work. Contrary to the other two models, it takes neither surface tension nor any hydrodynamic properties of gas–liquid Taylor flow into account and, as expected, the values of the pressure drop predicted by the Lockhart–Martinelli–Chisholm model deviate more from the measured values than the predictions of the other two models.

The model of Kreutzer et al. has been specifically developed for gas–liquid Taylor flow for Reynolds numbers larger than 150 and considers the same two sources of pressure drop as the model developed in this work. Furthermore, both models similarly take the spatial gas–liquid distribution into account in the form of, respectively, the liquid slug length and the gas bubble frequency. The main difference between the two models is that the dependence of the pressure drop on the gas bubble velocity is not incorporated in the model of Kreutzer et al. while it is included in the model developed in this work. The decreasing accuracy of the model of Kreutzer et al. with decreasing Reynolds number, and thus with decreasing gas bubble velocity, is shown to be due to the model not taking the effect of the gas bubble velocity on the pressure drop into account. This result, therefore, shows that other than the spatial distribution of the gas and liquid phases, the velocity distribution of the gas and liquid phases needs to be included when describing the pressure drop of gas–liquid Taylor flow in small channels.

The model developed in this work is valid for gas–liquid Taylor flow in round channels for Ca b < 0.01, Re b < 150 and negligible velocity of the liquid in the film surrounding the gas bubbles. When these conditions are met, the expression for the pressure drop over a single Taylor gas bubble is known from combining the analysis of Bretherton and that of Aussillous and Quéré, as shown in this work. Unfortunately, such an expression is not available when either the channel is not axisymetric, the liquid flow is not laminar, inertial effects are significant, or the liquid film surrounding the Taylor gas bubbles is not stagnant. Efforts to describe the pressure drop of gas–liquid Taylor taking one or more of these effects into account have been made, such as the semi-empirical model of Kreutzer et al., which includes inertial effects, but the understanding is still far from complete.