About this Journal Submit a Manuscript Table of Contents
Advances in Mechanical Engineering
Volume 2013 (2013), Article ID 987428, 8 pages
Research Article

A PLIC-VOF-Based Simulation of Water-Organic Slug Flow Characteristics in a T-Shaped Microchannel

1State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace, Xi'an Jiaotong University, No. 28, Xianning West Road, Xi'an 710049, China
2Department of Applied Chemistry, Okayama University of Science, Okayama, 700-0005, Japan
3Engineering Simulation and Aerospace Computing (ESAC), Northwestern Polytechnical University, P.O. Box 552, Xi'an, Shaanxi 710072, China

Received 18 February 2013; Revised 9 May 2013; Accepted 9 May 2013

Academic Editor: Cheng-Xian Lin

Copyright © 2013 Xian Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A water-organic slug flow in a T-shaped microchannel was numerically studied due to its importance in the microreactor system. Various factors affecting the flow mode were studied, for example, channel width, fluid viscosity, interfacial tension, and inlet velocity. The volume of fluid (VOF) method was used to track the liquid-liquid interface, and the piecewise-liner interface construction (PLIC) technique was adopted to get a sharp interface. The interfacial tension was simulated with continuum surface force (CSF), model and the wall adhesion boundary condition was taken into consideration. The results show that strong vortexes appear in both phases at the meeting sites of main and lateral channels where an organic slug is producing. Inlet velocity influences the slug length and flow mode greatly. The ratio between the slug lengths of two phases in the main channel is almost equal to the ratio between their inlet velocities. If the slug is produced, the interfacial tension and organic viscosity have less effect on the slug length for 200 µm microchannel. The slug producing rate is much higher in a narrow channel than that in a wide channel.

1. Introduction

Nowadays microreaction technology is widely used due to the excellent mass and heat transfer properties as well as uniform flow patterns and residence time distributions of the small microfabricated reactors. Some of these devices even are available commercially and being tested by a growing community of researchers [1]. The spectrum of applications includes gas and liquid flow as well as gas-liquid and liquid-liquid multiphase flow. Recently, immiscible liquid-liquid slug flows in microchannels absorb extensive attentions due to the segmentation of the two fluids, in which the chemical reaction, mixing, and diffusion of the solutes can be enhanced greatly. On the other hand, microdroplets can be produced in such a liquid-liquid system in the microchannel under some certain conditions, for example, small slug, large interfacial tension and small contact angle between fluid and solid wall. The efficiency and properties of the slug flow in microreactors are mainly influenced by the channel geometry, viscosity ratio between two fluids, interfacial tension, and the flux of two liquids [2]. The various affecting parameters mean that optimizing the slug-flow system in micro-eactors needs extensive experimental work. Therefore, numerical study on such a system is indispensable to provide a reasonable and economical designing process.

As we know, the variety and complexity of flow phenomena clearly pose major challenges to the modeling approach. Accurate simulation of a multiphase flow problems with a moving interface lies in the obtaining of a sharp interface, which is quite difficult and a great challenge for the CFD researchers. Several techniques were developed in the last 25 years and can be classified into three main categories according to physical and mathematical approaches: capturing (Lagrangian), tracking (Eulerian), and combined methods. Reference [3] presents the detailed review of these methods. Among the various approaches of multiphase flow simulation, volume-of-fluid (VOF) [4, 5], level set method [6], front-tracking method [7], and lattice Boltzmann method [8] are known as useful tools. Further new developments on capturing interface include particle-mesh method [9] and CIP method [10].

The VOF method is widely adopted by in-house codes and built-in commercial codes. It is a popular interface tracking algorithm which has proved to be a useful and robust tool since its development decades ago, in which a color function is set as 1 for the th fluid and 0 for the void, so that the interface locates in the computational cells in which . However, a sharp interface is hard to obtain when solving the color function equation due to the numerical diffusion, which is the key point for this method. Therefore, the reconstruction of the interface is required. Several interface reconstruction technique were developed. for example, donor acceptor, flux line-segment model for advection and interface reconstruction, simple line interface calculation, and piecewise linear interface calculation. The historical perspective and comparison of these interface reconstruction techniques were introduced in detail in [11].

In the multiphase flow, interfacial tension force plays a crucial role. A popular model for the interfacial tension force is continuum surface force (CSF) model, which was developed by Brackbill in 1992 [12]. In CSF model, the interfacial tension was interpreted as a continuous, three-dimensional effect across an interface, rather than as a pressure boundary condition on the interface. It is widely used in the multiphase simulation due to its easy implementation and high accuracy. On the other hand, using CSF model, the wall adhesion/wetting effect can be easily simulated by a simple boundary condition on the interfacial tension force, which is an important phenomenon in a microsystem [13].

In the present study, using PLIC-VOF method combined with CSF interface tension model, a water-organic slug flow in a T-shaped microchannel was numerically studied including the effects of channel size, viscosity of organic phase, interfacial tension and inlet velocities of two fluids.

2. Model System

Figure 1 shows the present model system. Water and organic phases were simultaneously injected into the microchannel from the inlets of main channel and lateral channel, respectively. The inlet velocities were indicated as for water and for organic phase. Both phases were assumed to be incompressible, Newtonian and immiscible at 1 atm, 25°C. The flow was treated as laminar one, which proved correct for the microflow [13]. The gravitational force was not taken into consideration due to its less importance in the microflow. The width of the channel was indicated as , the length of the lateral channel was , and it is equal to . The length of the main channel was .

Figure 1: Computational model system.

In the present work, two-dimensional computation was considered. The computational domain is shown in Figure 1. Here, and . rectangular grids were adopted. Since the T-shaped flow region occupied only part of the computational domain , the velocities outside the T-shaped region were set to zero manually in the computation.

3. Computational Methods and Governing Equations

In the present computation, the immiscible interface of the two phases was tracked by volume-of fluid (VOF) method [4, 5] on Eulerian grids. The color function was set as 1 for water and 0 for organic phase. The interface was reconstructed by piecewise liner (PLIC) technique which was initially proposed by Youngs [14]. The surface tension was computed by continuum surface force (CSF) model [12].

3.1. Dimensional Governing Equations

The dimensional governing equations are as follows.

Color function:

Continuity equation:

Momentum equation:

where and . Subscripts 1 and 2 stand for water (fluid 1) and organic phase (fluid 2), respectively. In the present computation, the physical properties [kg/m3], [Pa·s] were adopted for the water phase and [kg/m3], for the organic phase. in (3) is the total external forces per volume on the fluids. In this work, only interfacial tension force was considered as an external force.

The flow field was solved with finite difference method on a staggered mesh system. The highly simplified marker and cell (HSMAC) algorithm [15] was adopted to solve the pressure field and to improve the velocity field.

3.2. Interfacial Tension-CSF Model

in (3) is the surface tension force per volume, and the direction of is pointing into the concave side of the interface. Assuming the center of the interface curve is in fluid 1 (water), for the CSF method, the surface tension is modeled as follows [12]: where is the interfacial tension coefficient [N/m] and is the local curvature of the interface curve [1/m]. To avoid the acceleration due to surface tension depending on the density, (4) is improved as Thus, when (5) is applied on (3), the last term, that is, acceleration due to surface tension, is independent on the density. The curvature is calculated from where, is the unit normal to the interface and . The normal vector , and the cell-centered is computed as follows: For example; Then, we get where and are the dimensions of a grid in and directions.

Cell-centered tension force can be rewritten as The required face-centered values , in (3) are obtained by interpolating from the two nearest cell-centered values.

The effect of wall adhesion at fluid interfaces in contact with rigid boundaries, which is very important in a microflow, can be estimated easily in CSF model by applying a boundary condition in (11) as follow: where is contact angle between fluids and solid wall. Contact angle is influenced by various factors, for example, the species of the fluids, the smoothness, and geometry of the solid wall. The static contact angle can be measured experimentally when the fluids are at rest. is the unit wall normal directly into the solid wall, and is the unit tangent normal pointing into fluid 1. For the points on the solid walls or near the solid walls (the case when the grid boundary does not overlap with the wall), the tangent component can be computed using color function reflected at the wall, and the wall component can be computed in terms of , that is, . In the present study, the contact angle was assumed to be 90° for the easiness of computation. Our work [16] shows the effect of contact angle on the slug shape.

3.3. PLIC Reconstruction of the Interface

To get a sharp interface, PLIC technique [14] was used to make reconstruction of the interface, in which a cell-centered line segment is used to simulate the interface and the flux of fluid is computed geometrically according to the interface location and slope. The line segment is described as where and were computed according to (9) and (10), respectively.

Figure 2 shows four typical illustrations of the moving interface reconstruction in a computational grid. The dark area is full of fluid 1 and the light one fluid 2. Based on the sign of , , 16 reconstruction images can be obtained. Assuming is the acute angle between the line segment and -coordinate, then cot . Assuming , the four basic types shown in Figure 2, for which , , can be decided by the color function and as follows:type (a): or type (b): ,type (c): ,type (d): or .

Figure 2: Sample illustrations of the moving interface reconstruction in a computational grid.

Other 12 types can also be decided by the same way. Besides, there are two special cases: and , that is, the interface is horizontal and vertical. They should be considered separately from the above 16 cases.

After deciding the type of the line segment, the flux of fluid 1 whose color function is equal to 1 can be calculated analytically. In the present computation, the out-flow flux in a grid was computed, and the in-flow fluxes was corrected by the out-flow fluxes in the four neighbor grids. This is because the in-flow flux should be computed according to the value of in the neighbor grids. Figure 3 shows the example of flux calculation. Assuming the velocities on the four sides are all positive, that is, , the out-flow fluxes on the left and bottom sides of the grid , are zero. The ones on the right and top sides , are the shaded areas and , respectively, which can be computed geometrically by and . The volume of fluid 1 in the grid (dark area) is for two dimensional computation. Then, the color function was computed as follows: where superscript * stands for the value before the current correction. It should be noted that, in the boundary girds, both the in-flow and out-flow fluxes at the boundary side should contribute to the flux change at the same time.

Figure 3: The sketch of flux calculation in a computational grid.

Figure 4 shows the performances of our PLIC-reconstruction method on Zalesak’s slotted cylinder test problem [17]. The solution region is and grids are . Time step . It shows clearly that the present PLIC reconstruction tracks the sharp boundary efficiently.

Figure 4: Performance of VOF-PLIC reconstruction on Zalesak’s slotted cylinder test problem using grids.

4. Result and Discussion

4.1. Transient Flow Mode and Velocity Field

Figure 5 shows the transient flow modes at μm,  mm/s,  N/m,  and   Pa·s. The light color fluid is water, and the dark one is the organic fluid. The initial condition is shown in the picture  s. For the multiphase flow computation, a reasonable initial condition for the color function is very important. In the present work, the same initial condition were set for all the computational cases. At  s, the organic fluid flows into the main channel and the first organic slug (dark one) starts to be produced near the crossing of the lateral and main channels. At  s, the second organic slug is produced and the same situation occurs as the first one at  s. Such a process repeats with time proceeding. A water-organic alternating slug flow pattern in the main channel is obtained. In this paper, symbols and are used to indicate the lengths of water and organic slugs shown in Figure 5 ( s).

Figure 5: The transient flow modes at μm,  mm/s,  N/m and  Pa·s.

Figure 6 shows the local fluid streak lines and velocity vectors at  s for , , , and . At , an old organic slug was just produced. The fluid streak lines are wavy and no vortex is founded. A complete image is refereed to picture  s in Figure 5. At  s, a new organic slug is growing. In this period, strong vortexes appear on the interface of two phases and inside the organic phase. As soon as the organic phase touches the bottom of the solid channel wall, the water-organic interface on the right-hand side becomes perpendicular to the wall at once due to , which happens at  s. In this period, the vortex appears mainly in the water phase. At  s, vortexes disappear and the streak lines become wavy again. At  s, the next organic slug starts to be produced and the old one flows downstream. The process from  s to 0.320 s shown in Figure 6 repeats in the main channel. An organic slug was produced in about 0.12 s under current conditions.

Figure 6: Local fluid streak lines and velocity vectors at  s for μm,  mm/s,  N/m and  Pa·s.
4.2. The Effects of Channel Width , Interfacial Tension, and Organic Viscosity

The channel width affects the length of slug to a large extent. The cases of , 200, 100 and 50 μm at , and are studied here. The results show that for the microchannels with , and 50 μm, the corresponding lengths of organic slugs (Figure 5) are 1970, 730, 353, and 177 μm, respectively. The corresponding producing rate of organic slug N , that is, the number of organic slug produced in a second, is about 3, 7, 21 and 39, respectively. It is concluded that with the decrease in the channel width , the organic slug length decreases, and the producing rate of organic slug increases greatly.

Table 1 shows the length of organic slug [μm] under various interfacial tensions and organic viscosities at and . It can be seen for two cases , and that the parallel flow is founded; that is, there is no slug produced in the mail channel and a slightly wavy interface appears at the horizontal centerline of the main channel. Except for the above two cases, the slug can be produced and the slug lengths have no large difference (maximum difference 130 μm). Therefore, as long as slugs can be produced in the water-organic system, the organic slug length is not affected by and so much for the channel of . It should be mentioned that the position where the organic slug is produced is different. For large and small , the organic slug is produced just at the crossing of the lateral and main channels. For small and large , this position lies in the middle or downstream of the main channel and away from the crossing; that is, there exists a parallel flow before an organic slug is produced.

Table 1: The length of organic slug (m) under various interfacial tensions (N/m) and organic viscosities (Pa·s). Other conditions are: m and 5 mm/s, = 90°, = 1000 kg/m3, =  Pa·s and  kg/m3.
4.3. The Effect of Inlet Velocity of Water

In this section, the effects of the inlet velocity of water phase on the organic slug length and water slug length are studied. Although it is reasonable to study the effect of inlet volume flow rate , due to the two-dimensional computation, the cross section area cannot be evaluated correctly. Since the lateral and main channel have the same-sized cross section, the inlet velocity ratio of two phases is equal to the ratio of their inlet volume flow rates, and we use the inlet velocity to study the effect of inlet volume flow rate. To show this effect clearly, the main channel length is increased to in this section.

Figure 7 shows the flow modes under various inlet velocities of water at , , , and . It can be seen that the inlet velocity of water affects flow mode and slug lengths greatly. With increasing, the organic slug length decreases and the spacing between neighboring organic slugs increases. Particular, at , there is no slug produced; instead of it, a parallel flow mode appears and the water phase occupies most of the channel. Figure 8 shows the slug lengths and [μm] against . Figure 9 shows the ratio of slug lengths of two phases () against the ratio of inlet velocities of two phases (). Other conditions in Figures 8 and 9 are same with these in Figure 7. It can be seen that organic slug length decreases with the increase in while water slug length increases (Figure 8). From Figure 9, we can see is almost linearly proportional to . That means that the ratio between the slug volumes of two phases in the main channel is almost equal to the ratio of their inlet volume flow rates.

Figure 7: Flow modes under various inlet velocities of water at  s, μm,  mm/s,  N/m and  Pa·s.
Figure 8: Slug lengths and [μm] against [mm/s]. Other conditions are: μm,  mm/s,  N/m and  Pa·s.
Figure 9: against . Other conditions are: μm,  mm/s,  N/m and Pa·s.

Figure 10 shows the producing rate of organic slug against [mm/s]. It can be seen that N increases rapidly with increasing when is in some range. For the current condition,  mm/s.

Figure 10: Producing rate of organic slug against [mm/s]. Other conditions are: μm,  mm/s,  N/m and  Pa·s.

5. Conclusion

The fluid-fluid flow plays an important role in micro-reactor for the enhancement of heat, mass, and momentum transfer. In the present work, we numerically studied water-organic flow in a T-shaped microchannel which is often met in the microreactor. Several important factors affecting the flow mode were studied, for example, the channel width, inlet velocity, fluid viscosity, and interfacial tension. These factors influence the flow mode, slug length, slug producing rate in the microchannel, especially the inlet velocity. It is also found that when a slug is producing, strong vortexes appear in both fluids. After a slug is produced, it moves along the channel, and there is no vortexes founded in the downstream of the main channel. The numerical results are expected to be a guide for various experimental works.


: Color function, [-]
: Flow flux at the bottom boundary of a grid, [m2]
: Flow flux at the left boundary of a grid, [m2]
: Flow flux at the right boundary of a grid, [m2]
: Flow flux at the top boundary of a grid, [m2]
: Total external force per volume on the fluid, [N/m3]
: The width of microchannel, [μm]
: The length of the lateral channel, [μm]
:The length of the slug, [μm]
: The ratio between organic slug length and water slug length, [-]
: Normal vector to the interface, [-]
Unit normal vector to the interface, [-]
: The producing rate of organic slug, [s−1]
: Time, [s]
: Time interval, [s]
: Velocity vector [m/s]
: Velocity component in direction, [m/s]
, [-]
: Velocity component in direction, [m/s]
: The length of solution domain, [m]
: The width of solution domain, [m].
:, [rad]
:, [rad]
: A constant term in the line segment equation, [-]
:Contact angle, [°]
: The curvature of the interface, [m−1]
: Viscosity of fluid, [Pa·s]
: Density of fluid, [kg/m3]
: Interfacial tension coefficient between water and organic fluid, [N/m].
: Indicators of cell-centered value
: Organic phase
: Wall tangent component
: Water phase or wall normal component
: Components in and directions
: Indicators for fluids 1 (water) and 2 (organic fluid).


The authors would like to acknowledge the financial support for this work provided by the National Natural Science Foundation of China (no. 11242010) and by the project of Cooperation of Innovative Technology and Advanced Research in Evolutional Area by Okayama Prefecture Industrial Promotion Foundation from the Ministry of Education, Culture, Sports, Science and Technology, Japan.


  1. W. Ehrfeld, V. Hessel, and H. Lowe, Microreactors, New Technology for Modern Chemistry, WILEY-VCH, Weinheim, Germany, 2000.
  2. N. Harries, J. R. Burns, D. A. Barrow, and C. Ramshaw, “A numerical model for segmented flow in a microreactor,” International Journal of Heat and Mass Transfer, vol. 46, no. 17, pp. 3313–3322, 2003. View at Publisher · View at Google Scholar · View at Scopus
  3. H. Tang, L. C. Wrobel, and Z. Fan, “Tracking of immiscible interfaces in multiple-material mixing processes,” Computational Materials Science, vol. 29, no. 1, pp. 103–118, 2004. View at Publisher · View at Google Scholar · View at Scopus
  4. C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981. View at Scopus
  5. B. D. Nichols, C. W. Hirt, and R. S. Hotchkiss, “SOLA-VOF: a solution algorithm for transient fluid flow with multiple free boundaries,” Los Alamos Scientific Laboratory Report LA-8355, 1980.
  6. M. Sussman, “A level set approach for computing solutions to incompressible two-phase flow,” Journal of Computational Physics, vol. 114, no. 1, pp. 146–159, 1994. View at Publisher · View at Google Scholar · View at Scopus
  7. G. Tryggvason, B. Bunner, A. Esmaeeli et al., “A front-tracking method for the computations of multiphase flow,” Journal of Computational Physics, vol. 169, no. 2, pp. 708–759, 2001. View at Publisher · View at Google Scholar · View at Scopus
  8. M. R. Swift, W. R. Osborn, and J. M. Yeomans, “Lattice Boltzmann simulation of nonideal fluids,” Physical Review Letters, vol. 75, no. 5, pp. 830–833, 1995. View at Publisher · View at Google Scholar · View at Scopus
  9. J. Liu, S. Koshizuka, and Y. Oka, “A hybrid particle-mesh method for viscous, incompressible, multiphase flows,” Journal of Computational Physics, vol. 202, no. 1, pp. 65–93, 2005. View at Publisher · View at Google Scholar · View at Scopus
  10. T. Yabe, F. Xiao, and T. Utsumi, “The constrained interpolation profile method for multiphase analysis,” Journal of Computational Physics, vol. 169, no. 2, pp. 556–593, 2001. View at Publisher · View at Google Scholar · View at Scopus
  11. W. J. Rider and D. B. Kothe, “Reconstructing volume tracking,” Journal of Computational Physics, vol. 141, no. 2, pp. 112–152, 1998. View at Publisher · View at Google Scholar · View at Scopus
  12. J. U. Brackbill, “A continuum method for modeling surface tension,” Journal of Computational Physics, vol. 100, no. 2, pp. 335–354, 1992. View at Publisher · View at Google Scholar · View at Scopus
  13. V. Hessel, S. Hardt, and H. Lowe, Chemical Micro Process Engineering, Fundamentals, Modeling and Reactions, WILEY-VCH, Weinheim, Germany, 2004.
  14. D. L. Youngs, “Time-dependent multi-material flow with large fluid distortion,” in Numerical Methods for Fluid Dynamics, K. W. Morton and M. J. Baines, Eds., pp. 273–285, Academic Press, New York, NY, USA, 1982.
  15. C. W. Hirt, B. D. Nichols, and N. C. Romero, “SOLA: a numerical solution algorithm for transient fluid flow -Addendum-,” Los Alamos Scientific Laboratory Report LA-5852, 1975.
  16. X. Wang, H. Hirano, and N. Okamoto, “Numerical investigation on the liquid-liquid two phase flow in a Y-shaped microchannel,” The Australian and New Zealand Industrial and Applied Mathematics Journal, vol. 48, pp. C963–C976, 2008.
  17. S. T. Zalesak, “Fully multidimensional flux-corrected transport algorithms for fluids,” Journal of Computational Physics, vol. 31, no. 3, pp. 335–362, 1979. View at Scopus