1 Introduction

Taylor flow is the main flow pattern in microchannels(diameter < 1mm)^[1-5]. Taylor flow consists of series of gas bubbles which are separated by liquid slugs. The length of gas bubble is larger than the channel diameter and a thin liquid film separates bubble from the channel wall. Taylor flow has advantage in intensifying heat and mass transfer because of the recirculation within liquid slugs and the weakening of axial mixing between liquid slugs^[6-10].

CFD simulation has been extensively used as and shown to be a tool in detailed study of flow characteristics of Taylor flow in microchannels.In earlier studies, some researchers^[11-13] assumed a bubble shape and set the flow conditions around the bubble using periodic boundary conditions to simulate Taylor flow. With the development of various interface capturing techniques, the volume of fluid (VOF) model^{[9, 14]}, level-set method^[15] and phase-field method^[16] were applied. The literatures reveal that the liquid film near the wall, which is an important fluid dynamic character in microchannels, could not be properly captured if the mesh resolutions for the CFD computing was low^{[14, 17-19]}. Some researchers were able to seize the liquid film, but at the expense of a veritable increase in computing resources^{[9, 17]}. Therefore, a trade-off between the fine resolution for the liquid film and the computing efficiency is an issue for the research of fluid dynamic of Taylor flow by CFD method.

The aim of the present paper is to numerically study the hydrodynamic characteristics of Taylor flow(i.e. liquid film thickness, gas hold-up, bubble velocity, pressure drop per unit length and slug length) in microchannels with a gradient mesh approach.Simulations are performed to propose proper rules for the gradient mesh approach so as to capture the liquid film. The simulations are validated by comparing the results with correlations and experimental data. Through analyzing the influence of liquid film on gas and liquid slug lengths, new correlations of slug length based on the previous correlations and our simulation data are also proposed.

2 Modeling of Taylor flow in a T-junction microchannel 2.1 Geometry model

As shown in Fig. 1, the gas and liquid phases are fed separately into the two inlets of the vertical mixing zone prior to entering the horizontal flow zone. The cross-sectional width of the channel is represented by d and the length of mixing zone and flow zone is 9d and 27d, respectively.

The assumptions adopted in the present work are as follows:Because the cross-sectional size of microchannel is small, the fluid flow is considered as laminar.Both the gas and liquid phases are incompressible since the pressure drop along the channel is relatively small.The physical properties of both phases are constant.As Hazel and Heil^[20] stated, the influence of gravity can be ignored when Bo < 1, which is the case in the present work.And the system is assumed to be isothermal.

2.2 Governing equations

In order to track the gas-liquid interface and capture the liquid film between the bubbles and the channel wall, the volume of fluid (VOF) method is adopted in our computation.The governing equations of the VOF formulations on two-phase flow are as follows:

Continuity equation:

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \vec v} \right) = 0 $

(1)

Momentum equation:

$ \begin{gathered} \frac{{\partial \left( {\rho \vec v} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \vec v\vec v} \right) = - \nabla p + \\ \nabla \cdot \left[ {\mu \left( {\nabla \vec v + {{\left( {\nabla \vec v} \right)}^T}} \right)} \right] + \rho \vec g + \vec F \\ \end{gathered} $

(2)

Volume fraction equation:

$ \frac{{\partial {\alpha _{\text{G}}}}}{{\partial t}} + \vec v \cdot \nabla {\alpha _{\text{G}}} = 0 $

(3)

In the present work, the continuum surface force (CSF) proposed by Brackbill is adopted as the surface tension model^[21].

2.3 Proposed gradient mesh approach

Different mesh approaches have been adopted in various numerical simulation cases by many researchers^{[9, 14, 17-19]}, including uniform coarse mesh, uniform fine mesh and non-uniform mesh.

However, it is still hardly achievable to predict the hydrodynamics of Taylor flow accurately with a high computing efficiency. In the present work, a gradient mesh approach is proposed to tackle this issue

The proposed gradient mesh is shown in Fig. 2. The meshed region for the computing is divided into two zones, the uniform mesh zone and near wall refined mesh zone. The size of the meshes in the uniform mesh zone is h. In the refined mesh zone, the size of the finest mesh nearest to the wall is φ and the size of the second row meshes near the wall is φ×k¹, the size of the third row is φ×k²and so forth. Thus, the width of the refined mesh zone is given by

$ {H_r} = \sum\limits_{j = 0}^i {\varphi \times {k^j}} $

(4)

wherei is the number of the last row of the refined mesh zone as shown in Fig. 2. According to the definition of the two zones, we also have

$ h = \varphi \times {k^i} $

(5)

There are five parameters in equations (4) and (5), and thus three of them need to be specified in order to define the gradient mesh system. In general, h, the mesh size in the uniform zone, can be determined by checking the mesh independency of the simulated results on the Taylor flow characters, like bubble length.

According to this rule, last raw width multiplied by coefficient k is equal to next raw width until the raw width φ×kⁱapproaches or equals to h.In the present work, i is set as eight to ensure that the liquid film is all in the grid refinement region. Through checking grid independence, the uniform mesh size (h) is finally set as 1/25 of the cross-sectional width (d) which is 0.1, 02, 0.3, 0.4, 0.5, 0.8 and 1.0 mm, respectively.

2.4 Numerical algorithm

In the present work, a finite volume method (FVM) based commercial CFD software, ANSYS FLUENT (Release 6.3.26, 2007), is used to simulate two-dimensional, double precision, unsteady, two-phase Taylor flow in microchannels.The boundary conditions are specified as follows: velocity inlet for gas and liquid feeds and pressure outlet for the microchannel outlet. The no-slip wall condition is applied to the walls. The implicit body force treatment for the body force formulation in the multiphase model panel is selected to improve solution convergence by accounting for the partial equilibrium of the pressure gradient and body forces in the momentum equations. Table 1 shows the properties of the gas and liquid phases used in the present study. To obtain more realistic looking rounded bubbles with low grid density, a contact angle of 0° is set for all cases. The pressure staggering options (PRESTO!) method is used for the discretization of the pressure interpolation. The momentum equations and the interface interpolation are discretized by second order upwind scheme and geometric reconstruction scheme respectively. The pressure-implicit with splitting of operators(PISO) scheme is chosen for pressure-velocity coupling and the initial value of water volume fraction is set as 1. In addition, the residual values, the time step and the maximum number of iterations per time step are properly adjusted to ensure convergence.

3 Results and discussion

The simulations are carried out on a workstation, four AMD Opteron CPUs (two cores and 2.2 GHz each) and 8 GB RAM. The air and water superficial velocities are determined based on the different flow regimes identified in microchannels^[1-5] and the value is mainly set as 0.1 m/s and 0.15 m/s, respectively.

3.1 Determine a mesh arrangement

In order to precisely capture the liquid film, series of simulations are conducted with different gradient mesh arrangements and the numerical results of the liquid film thickness are displayed in Table 2.According to the table, with the decrease of φ, the liquid film thickness becomes thinner. However, when φ is equal to or less than 1.5 μm, the liquid film thickness maintains a stable value. In summary, when φ is less than δ, the liquid film can be captured. If φ is equal to or less than 1/5 of δ, a stable film thickness can be obtained.Thus, based on the numerical accuracy and computational efficiency, φ and k is determined as 1.5 μm and 1.38, respectively.

3.2 Comparison of numerical methods

Square grids size of 5 μm×5 μm in the uniform mesh region and rectangular grids size of 2.4 μm×5.0 μm in the liquid film region, shown in Fig. 3a) were adopted in the study of Gupta^[9]. In the present research, the gradient mesh approach is applied, as shown in Fig. 3b). A single bubble and the adjacent liquid slug are defined as a unit cell.The computing time needed to acquire one unit cell with Gupta′s mesh approach would be considerable longer than that with our gradient mesh approach with the same computing resources, which is due to the reason that the grids number used in Gupta′s mesh approach is eleven times larger than that used in our gradient mesh approach.According to Fig. 3b), in the liquid film region, our simulation result turns out to be in good agreement with that from Gupta. However, it is hard to capture the gas-liquid interface accurately in the uniform mesh region. In curved microchannels, the finer meshes were used in the streamwise direction and the clearer gas-liquid interface could be seized^[18].Thus, the grid size is decreased to 10 μm along the direction of x-axis to investigate the effect of mesh size, as displayed in Fig. 3c). It is obvious that the numerical result of refined mesh resolution along the direction of x-axis is in good agreement with that from Gupta. However, the refined mesh resolution does not have significant influence on the hydrodynamic characteristics compared with the original mesh resolution. For instance, the bubble velocity is 0.260 m/s and 0.265 m/s in Fig. 3b) and Fig. 3c), respectively.As there is no chemical reaction in the microchannel, the mesh does not need to be refined in the streamwise direction.

3.3 Model validation

To confirm the validity of the model applied in the present study, comparisons of the hydrodynamic characteristics of Taylor flow between our simulation results and experimental data or correlation values are carried out, as discussed below.

3.3.1 Liquid film thickness

In the research of Han^[22], measurement of liquid film thickness in a micro parallel channel was conducted. Through analyzing liquid film theoretically and experimentally, a conclusion was drawn that the liquid film thickness in a parallel channel can be predicted well by the circular tube correlations^[22]. In their experiments, parallel channel was represented by the channel with a high aspect ratio (i.e. wide shallow microchannel).In the present work, series of 2D numerical simulations that simulate flow in an infinitely wide channel are performed in a plane of 500 μm high.Thus, the liquid film thickness of 2D simulation can be accurately predicted by the empirical correlations for circular tubes and can match well with the experimental data from experiments for circular tubes.

The correlation 2δ/d=1.34Ca^2/3 proposed by Bretherton^[23]is commonly used to predict the liquid film thickness and the correlation value under our simulation condition is 8.069 μm which is almost identical to our numerical prediction (8.098μm).

In the study of Han and Shikazono^[24], the liquid film thickness was measured in micro tube slug flow. Their experimental data are compared with our numerical results at the same bubble velocity in the present paper, as exhibited in Fig. 4.The simulation results are in good agreement with the experimental data. When bubble velocity is larger than 1.3 m/s, the simulation results and the experimental data are almost the same.

3.3.2 Gas hold-up

The gas hold-up can be estimated based on a unit cell. However, to avoid incidental error that may arise by relying on just one unit cell, the value of gas hold-up is calculated in the present work by dividing the volume of three gas bubbles by that of the three unit cells at different development time, as shown in Fig. 5.

Then, the average of such estimated values is compared with that predicted by the formula ${\varepsilon _{\text{G}}}=\left({\frac{{{A_{\text{B}}}}}{A}} \right){U_{\text{G}}}/\left({{U_{\text{G}}} + {U_{\text{L}}}} \right)$ proposed by Warnier^[25]in Table 3.In order to investigate the effect of liquid film on gas hold-up, simulations with uniform and gradient mesh approach are both performed. As shown in Table 3, since the liquid film can not be captured with uniform mesh approach, the gas hold-up is almost equal to U_G/(U_G+U_L). However, the gradient mesh approach is able to capture the liquid film and the error between the formula value and our simulation result is less than 1%.Thus, the gradient mesh approach is more accurate to predict the gas hold-up.

3.3.3 Bubble velocity

The velocity of the intersection point of central line and bubble front line is regarded as the bubble velocity.To calculate the bubble velocity, two different moments are selected as T₁ and T₂, as shown in Fig. 6. H is the distance that the first Taylor bubble passes from T₁ to T₂. The value is thus calculated by dividing the distance H by the time difference T₂-T₁.

Table 4 shows the simulation result and the value of formula U_B=U_G/ε_G used by Warnier^[25].According to the table, the simulation result differs within 1% from the formula value.

3.3.4 Pressure drop per unit length

Fig. 7 shows the pressure profile at the center line for a specific simulation condition, and the pressure drop of the microchannel is represented by that across a unit cell.

The expression for pressure drop of Taylor flow proposed by Kreutzer ^[26] is listed below:

$ \frac{{\Delta p}}{L} = {\varepsilon _{\text{L}}}f\left( {\frac{4}{d}} \right)\left( {\frac{1}{2}{\rho _{\text{L}}}{{\left( {{U_{\text{L}}} + {U_{\text{G}}}} \right)}^2}} \right) $

where ε_L, the dynamic liquid hold-up, is the fraction of channel length occupied by the liquid slugs. The friction factor f is given by

$ f = \frac{{16}}{{Re}}\left[ {1 + {\alpha _1}\frac{d}{{{L_L}}}{{\left( {\frac{{Re}}{{Ca}}} \right)}^{1/3}}} \right] $

where L_L is the liquid slug length, Re=d(U_L+U_G)ρ_L/μ_L is the Reynolds number and α₁ is a constant. According to Kreutzer ^[26], the value of α₁ was found to be 0.07 and 0.17 for numerical simulations and experimental data, respectively. As discussed in the literature^[26], such a large difference on the value of α₁ was attributed to the Marangoni effects generated by the existence of liquid impurity in the experiments.The pressure drop per unit length calculated by the correlation is 32.08 Pa/mm and 49.80 Pa/mm using the numerical and experimental value of α₁.Our simulation result is 23.12 Pa/mm while the prediction changes to 15.72 Pa/mm with uniform mesh approach. Thus, the simulation result with gradient mesh approach is closer to the correlation value than that with uniform mesh approach.The difference between these two predictions may result from the reason that the velocity gradient linear to shear force with uniform mesh approach is not captured as accurate as that with gradient mesh approach, while the pressure drop is dominated by shear force in microchannels.

3.3.5 Slug length

Table 5 shows the comparison between the experimental bubble lengths from Santos and Kawaji^[27]and the simulation results obtained under the identical condition.Even though the structure of the inlet T-junction used by Santos and Kawaji^[27]is different from that in the present work, there is not much difference in the slug length between the channels with these two inlets, which was reported by Qian and Lawal^[14]. However, Table 5 reveals a large discrepancy between numerical predictions and experimental data. This discrepancy is likely due to the flow disturbances produced in the experiments that is not considered in simulations using the incompressible model and the fixed velocity boundary condition^[27]. Still, our simulations support the trend of the slug lengths’ variation with the increase of U_L obtained from the experimental data.According to the table, our simulation results are much closer to the experimental data than those from Santos and Kawaji^[27].

4 Newly developed correlations of slug length

Mass transfer and pressure drop in microchannels are highly dependent on slug length. Some researchers^[28-30]had proposed mass transfer coefficients that were related to slug length and some researchers^{[25, 30-32]}had measured slug length using various methods (MRI method, conductivity method and photographic method). Laborie^[33] measured bubble and liquid slug length and realized that slug lengths were correlated with Re and Eo. Qian and Lawal^[14]performed the numerical simulations of Taylor flow in T-junction microchannels with varying cross-sectional width (1~3 mm), using commercial CFD soft ware and VOF model. Qian and Lawal proposed the slug length correlationsand illustrated the decisive influence of gas hold-up and the limited effects of Re and Ca.Kumar^[18] and Sobieszuk^[34]conducted experimental studies of the bubble and liquid slug length in Taylor flow. The measurements agreed well with the values obtained from the correlation of Qian and Lawal^[14]. Sobieszu kimproved and simplified Qian and Lawal′s correlation and derived a new correlation only concerning ε_G and We.

According to the formula used by Warnier^[25], gas hold-up is dependent on the liquid film thickness.Moreover, the dimensionless slug length is mainly depended on the gas hold-up in the study of Qian and Lawal^[14].Therefore, liquid film thickness in directly influences slug length.On the basis of dimensional analysis method and the numerical simulation results, Qian and Lawal^[14] proposed the correlations of gas and liquid slug length in a T-junction microchannel. However, the liquid film was not captured in the numerical simulations performed by Qian and Lawal, which resulted in a 20%—30% difference between the slug length ratio and the volume flow rate ratio, causing the numerical slug lengths to be in consistent with the experimental data. Thus, new correlations of slug length in a T-junction microchannel are further developed based on the correlations of Qian and Lawal and the simulation results in the present work. The influences of gravity, viscosity, density and contact angle are not considered.Dimensional analysis method (the details can be found in the research of Qian and Lawal^[14]) is used and the dimensionless slug length has the following dependence:

$ {L^*} = {L^*}\left( {Re,Ca,{\varepsilon _{\text{G}}}} \right) $

(6)

Through analyzing 28 sets of numerical data, the following correlations of the dimensionless slug length are determined:

$ {L_{\text{G}}}/d = 6.8\varepsilon _{\text{G}}^{1.1858}{\left( {1 - {\varepsilon _{\text{G}}}} \right)^{0.1414}}R{e^{ - 0.0908}}C{a^{ - 0.1303}} $

(7)

$ \begin{gathered} {L_L}/d = 6.8\varepsilon _{\text{G}}^{1.1858}{\left( {1 - {\varepsilon _{\text{G}}}} \right)^{0.1414}}\left( {1 - } \right. \\ \left. {{\varepsilon _{\text{G}}} - 1.34C{a^{2/3}}} \right)R{e^{ - 0.0908}}C{a^{ - 0.1303}} \\ \end{gathered} $

(8)

Fig. 8 illustratesthe comparison between our numerical predictions and the values obtained from our newly developed correlation.Eq.(7) is obtained through correlating 28 sets of numerical simulation data within a channel width of 1 mm and other parameter values are in the range: 0. < ε_G < 0.6, 10 < Re < 300, 0.001 < Ca < 0.01. Eq.(8) is obtained by multiplying Eq.(7) by ${L_L}/{L_{\text{G}}}={\varepsilon _{\text{G}}}/\left({1 -{\varepsilon _{\text{G}}} -1.34C{a^{\frac{2}{3}}}} \right)$, for ${\varepsilon _{\text{G}}}=\left({\frac{{{A_{\text{B}}}}}{A}} \right){U_{\text{G}}}/\left({{U_{\text{G}}} + {U_{\text{L}}}} \right)$, $\frac{{{A_{\text{B}}}}}{A}=\frac{{d -2\delta }}{d}=1 -1.34C{a^{\frac{2}{3}}}$ and ${U_{\text{G}}}/\left({{U_{\text{G}}} + {U_{\text{L}}}} \right) \approx {L_L}/\left({{L_{\text{G}}} + {L_L}} \right)$ under our simulation condition. The average error between the numerical predictions and the correlation values is less than 10%.Fig. 9 reveals the comparison among the experimental bubble lengths from Santos and Kawaji^[27], the values of the correlation of Qian and Lawal and those calculated from Eq.(7). As shown in Fig. 9, when U_G is equal to 0.317 m/s, the prediction of slug length by our correlation is more accurate than that conducted by Qian and Lawal. However, due to the limited application range of our correlation, when U_G is 0.04 m/s or 0.079 m/s, only part of the correlation results (ε_G>0.3) are closer to the experimental data than those computed from the correlation of Qian and Lawal.The comparison demonstrates that the newly developed correlations of slug length can accurately predict gas and liquid slug length in the specific range.

5 Conclusions

In the present study, a gradient mesh approach is applied to simulate Taylor flow in a T-junction microchannel using the VOF method. Series of numerical simulations are performed to generate proper rules for the approach and determine an effective mesh arrangement to capture the liquid film steadily and the comparison of numerical methods is also carried out. The validity of our model is confirmed through comparing the numerical predictions with experimental data and correlation values.Finally, through analyzing 28 sets of numerical simulation data, new correlations of slug length in a T-junction microchannel are further developed.The conclusions of this research are as follows:

1) The gradient mesh approach with our newly proposed rules can capture the liquid film steadily. The mesh arrangement is able to make a good balance between the computing accuracy and efficiency for predicting the hydrodynamic characteristics of Taylor flow.

2) When the finest mesh size (φ) is less than the liquid film thickness(δ), the liquid film can be captured. If the finest mesh size (φ) is less than 1/5 of the film thickness(δ), an accurate and stable film thickness can be achieved.

3) The newly derived correlationsof slug length are able to accurately predict the gas and liquid slug length in the microchannel under some certain conditions.

Notation:

A——channel cross-sectional area, m²;

A_B——Taylor bubble cross-sectional area, m²;

Bo——Bond number (=ρgd²/4σ);

Ca——Capillary number (=μ_L(U_G+U_L)/σ);

d——channel cross-sectional width, m;

Eo——Eotvos number (=(ρ_L-ρ_G)d²g/σ);

f——friction factor;

F——volume force, N/m³;

g——gravitational acceleration, N/kg;

H——distance, m;

h——uniform mesh size, m;

k——the ratio of rows in gradient mesh;

L——unit cell length, m;

L^*——dimensionless slug length;

L_G——Taylor bubble length, m;

L_L——liquid slug length, m;

p——static pressure, Pa;

Re——Reynolds number (=ρ_L(U_G+U_L)d/μ_L);

t——time, s;

T₁——moment one, s;

T₂——moment two, s;

U_B——bubble velocity from correlation, m/s;

U_G——gas superficial velocity, m/s;

U_L——liquid superficial velocity, m/s;

We——Weber number (=ρ_LdU_B²/σ);

Δp——pressure drop, Pa.

Greek letters:

α——volume fraction of one phase;

α₁——empirical constant;

ε——hold-up of one phase;

μ——viscosity of one phase, mixture viscosity (=(1-α_G)μ_L+α_Gμ_G), Pa·s;

ρ——density of one phase, mixture density (=(1-α_G)ρ_L+α_Gρ_G), kg/m³;

δ——film thickness, m;

λ——Interface thickness, μm;

σ——surface tension, N/m;

υ——velocity vector, m/s;

φ——the finest mesh size in gradient mesh, m.

Subscripts and superscripts:

B——Taylor bubble;

G——gas phase;

L——liquid phase.