VOF Modeling and Analysis of the Segmented Flow in Y-Shaped Microchannels for Microreactor Systems

Abstract

Microscaled devices receive great attention in microreactor systems for producing high renewable energy due to higher surface-to-volume, higher transport rates (heat or/and mass transfer rates), and other advantages over conventional-size reactors. In this paper, the two-phase liquid-liquid flow in a microchannel with various Y-shaped junctions has been studied numerically. Two kinds of immiscible liquids were injected into a microchannel from the Y-shaped junctions to generate the segment flow mode. The segment length was studied. 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 simulated flow pattern presents consistence with our experimental one. The numerical results show that a segmented flow mode appears in the main channel. Under the same inlet velocities of two liquids, the segment lengths of the two liquids are the same and depend on the inclined angles of two lateral channels. The effect of inlet velocity is studied in a typical T-shaped microchannel. It is found that the ratio between the lengths of two liquids is almost equal to the ratio between their inlet velocities.

1. Introduction

To save energy or decrease energy consumption, the development of new advanced technologies is desired. Among them, the microreaction technology is considered an attractive solution way to replace the traditional energy resources like fossil energy resources. Nowadays, micro-reaction 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 [1]. The microreaction technology has other merits of perfect mixing, effective heat exchange, and small reaction volume. The representations in better performance of heat exchange are higher surface-to-volume ratio and better heat transfer magnitudes compared to normal equipment. The better mixing is attributed to the laminar flow where fast diffusion and thereby perfect mixing within seconds are achieved. Extensive reviews of advances in microreaction technology can be found in [2–6].

The spectrum of applications includes gas and liquid flow as well as gas-liquid and liquid-liquid multiphase flows. Recently, immiscible liquid-liquid slug flows in Y-shaped microchannels received extensive attention due to the segmentation of the two fluids, in which the chemical reaction, mixing, and diffusion of the solutes can be enhanced greatly [7–9]. On the other hand, microdroplets can be produced in such a liquid-liquid system in Y-shaped microchannels under some certain conditions. Such can be widely used in the areas of pharmacy, cosmetics, food, and various powder technologies. The efficiency of Y-shaped microchannels for various purposes depends on the segment lengths of both phases. For example, when the segment length is small, the enhanced reactions and small-sized particles can be obtained. Besides, there are various parameters affecting the segment length, such as the channel geometry, viscosity/density ratio between two fluids, interfacial tension, and the flux of two liquids. The various affecting parameters mean that optimizing the segment flow system in microchannels needs extensive experimental works, which would cause a large amount of human and financial resources. In this sense, numerical investigation (including simulation and prediction) on such flow system is indispensable to provide a reasonable and economical designing process for microreactor systems.

The variety and complexity of multiphase flow phenomena clearly pose major challenges to the modeling computational approaches. Tang et al. [10] presented the detailed review of the computational methods for multiphase flow simulation. Among the various approaches of multiphase flow simulation, volume of fluid (VOF) [11, 12], level set method [13], front-tracking method [14], and lattice boltzmann method [15] are known as useful tools for modeling multiphase flow phenomena. Further new developments on capturing interface include particle-mesh method [16] and constrained interpolation profile (CIP) method [17].

The VOF method is widely adopted by in-house codes and built-in in commercial codes. It is a popular interface tracking algorithm and has proved to be a useful and robust tool. The idea in VOF method is that 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 0<𝐶_𝑘<1. However, a sharp interface is hard to obtain when solving the color function equation due to the numerical diffusion. This is the key point for this method. Therefore, the reconstruction of the interface is required. Fortunately, several interface reconstruction techniques were developed successfully to fix the requirement, for example, Donor Acceptor, FLAIR (flux line-segment model for advection and interface reconstruction), SLIC (simple line interface calculation) and PLIC (piecewise linear interface calculation). The historical perspective and comparison of these interface reconstruction techniques were introduced in [18].

In the multiphase flow, interfacial tension force plays a crucial role. A popular model for the interfacial tension force is the continuum surface force (CSF) model, which was first developed by Brackbill et al. in 1992 [19]. In the 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 the CSF model, the wall adhesion/wetting effect can be conveniently simulated by a simple boundary condition on the interfacial tension force, which is an important phenomenon in a microsystem [20].

In this study, the influence of two factors on the segment length of a liquid-liquid flow system in Y-shaped microchannels is numerically studied. These are the inclined angles of Y-shaped junctions and the inlet flow rate of fluids. In the present simulation, 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.

2. Physical Model

Figure 1 shows the present model system. Two immiscible fluids were simultaneously injected into a Y-shaped microchannel from two lateral channels with the inlet velocities 𝑢₀₁ and 𝑢₀₂, respectively. Both phases were assumed to be incompressible and Newtonian at 1 atm and 25°C. The flow was treated as laminar one, which proves correct for the microflow under current situations [1]. The gravitational force was not taken into consideration due to its unimportance in the microflow.

Figure 1 

Computational model system.

The present computations are two-dimensional and steady, and the geometrical parameters of the computational domain (𝑥𝑙×𝑦𝑙) are also shown in Figure 1. The widths of the lateral and main channels ℎ are the same and equal to 200 μm. The lengths of the lateral channels and the main channel are 𝑙=5 h and 8𝑙=40 h, respectively. The inclined angles of two lateral channels (Y-shaped junction) are indicated as 𝜑₁ and 𝜑₂, respectively. An independent grid system of 800 × 220 cells is adopted for all simulations. The motivation of this study is to study the segmented flow at various inclined angles or inlet velocities of a typical Y-shaped microchannel with fixed inclined angles.

3. Governing Equations

The governing equations for the fluids are as follows.

Color function𝜕𝐶+𝜕()𝜕𝑡𝑢𝐶+𝜕()𝜕𝑥𝑣𝐶𝜕𝑦=0. (1)

Continuity equation𝜕𝑢+𝜕𝑥𝜕𝑣𝜕𝑦=0. (2)

Momentum equations𝜕𝑢𝜕𝑡+𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢1𝜕𝑦=−𝜌𝜕𝑝+𝜇𝜕𝑥𝜌𝜕²𝑢𝜕²𝑥+𝜕²𝑢𝜕²𝑦+1𝜌𝑓_sv𝑥,𝜕𝑣𝜕𝑡+𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣1𝜕𝑦=−𝜌𝜕𝑝+𝜇𝜕𝑦𝜌𝜕²𝑣𝜕²𝑥+𝜕²𝑣𝜕²𝑦+1𝜌𝑓_sv𝑦, (3) where 𝜌=𝜌₂+𝐶(𝜌₁−𝜌₂) and 𝜇=𝜇₂+𝐶(𝜇₁−𝜇₂). Subscripts 1 and 2 denote fluid 1 and 2, respectively. 𝐶 is the color function, which is set to 1 for fluid 1 and 0 for fluid 2, respectively. On the interface of two fluids, the value of 𝐶 is taken at (0, 1). In this study, the physical properties of the two fluids are assumed to be 𝜌₁=1000 kg/m³, 𝜇₁=0.001Pa·s and 𝜌₂=800 kg/m³, 𝜇₂=0.02Pa·s. The interfacial tension is set to be 𝜎=0.05 N/m. The magnitudes of inlet velocities are |𝑢₀₁|=|𝑢₀₂|=0.005 m/s, and their directions are parallel to the lateral channels (Figure 1). Accordingly, the corresponding typical dimensionless groups for microchannel are Re=𝜌₁𝑢₀₁ℎ/𝜇₁=1.0 and Ca=𝜇₁²/(𝜎𝜌₁ℎ)=10⁻⁴.

The flow field was solved with finite difference method on a staggered mesh system, and the HSMAC (highly simplified marker and cell) algorithm [21] is employed to solve the pressure field and improve the velocity field. For HSMAC, only the local pressure value 𝑝^𝑛+1_𝑖,𝑗 (or pressure correction 𝑝^′_𝑖,𝑗) is regarded as unknown in the Poisson equation of pressure and is solved by Newton-Raphson method. The velocities in the neighbor grids are corrected by 𝑝^′_𝑖,𝑗 until the continuity equation is satisfied.

Since 𝐶 is not a continuous variable, (1) cannot be solved as usual. To get a sharp interface, the PLIC technique [22] is adopted, 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 can be described as 𝑛_𝑥𝑥+𝑛_𝑦𝑦+𝜂=0. (4)𝑛_𝑥, 𝑛_𝑦, and 𝜂 can be computed according to 𝐶. Based on the values of 𝑛_𝑥, 𝑛_𝑦, and 𝜂, the flux of fluid 1, whose color function is equal to 1, can be calculated analytically. In the present computation, the outflow flux in a grid is computed, and the in-flow flux is corrected by the out-flow fluxes of the four neighboring grids. This is because the in-flow flux of a grid should be computed according to the values of 𝐶 in its neighboring grids. Our previous work [23, 24] shows the detailed implementation and performance of our PLIC-reconstruction method.

𝑓_sv in (3) is the surface tension force per volume. Its direction is pointing into the concave side of the interface. Assuming the center of the interface curve is in fluid 1, for the CSF method, the surface tension is modeled according to [15]: 𝑓_sv=𝜎𝜅∇𝜌𝜌₁−𝜌₂=𝜎𝜅∇𝐶, (5) where 𝜎 is the interfacial tension coefficient and 𝜅 is the local curvature of the interface curve. To avoid the acceleration due to surface tension depending on the density, (5) is improved as follow: 𝑓_sv=𝜌𝜌₁+𝜌₂/2𝜎𝜅∇𝐶 (6) Thus, when (6) is applied on (3), the last term (i.e., acceleration due to surface tension) is independent of the density. The curvature 𝜅 is calculated from 1𝜅=−∇⋅𝐧=|𝐧|𝐧|𝐧|⋅∇|𝐧|−(∇⋅𝐧), (7) where 𝐧 is the unit normal to the interface and 𝐧=𝐧/|𝐧|. The normal vector 𝐧=∇𝐂. The detailed computational equations of every term in (7) are shown in [20].

The effect of wall adhesion at fluid interfaces in contact with rigid boundaries, which is very important in a microflow, can be estimated with CSF model by applying a boundary condition in (7) as follows: 𝐧=cos𝜃_ca𝐧_𝑤+sin𝜃_ca𝐧_𝑡, (8) where 𝜃_ca is the contact angle between fluid and solid wall. It is set to be 90° in the present computation. 𝐧_𝑤 is the unit wall normal directed into the solid wall and 𝐧_𝑡 is the unit tangent normal pointing into fluid 1.

4. Results and Discussion

4.1. Comparison with Experiment

The validity of our computational code was proved by comparing with our experimental result. Figure 2 is the experimental photograph. The red liquid is cyclohexane (with red color dye) and the green one is water (with blue color dye). The channel width is 200 μm, the inlet velocities are 𝑢₀₁=𝑢₀₂=0.025 m/s, and 𝜑₁=𝜑₂=30^∘. The interfacial tension is 0.025 N/m. The segment lengths of cyclohexane and water are about 0.67 and 0.7 mm, respectively. It should be noted that the exact segment length is very hard to be measured because the real interfaces shown in Figure 2 are not as clear as the computational one, which is due to the irregular contact angle caused by the complicated solid wall condition. Figure 3 shows the corresponding numerical result of our code under the same conditions with experiment. The blue and green liquids in Figure 3 are organic and aquatic phases, respectively. The computational segment lengths of cyclohexane and water are 0.68 and 0.71 mm. The images of slug generation near the Y-shaped junction, which is indicated in the circle in Figure 2, are quite similar for experimental and numerical results. Therefore, from both qualitative and quantitative points of view, the numerical and experimental results agree well.

Figure 2 

Experimental photograph for the segment flow in a Y-shaped microchannel at ℎ=200 μm, |𝑢01|=|𝑢02|=0.025 m/s and 𝜑1=𝜑2=30∘. Red: cyclohexane; green: water.

Figure 3 

Numerical simulation result for the segment flow in a Y-shaped microchannel at ℎ=200 μm, |𝑢01|=|𝑢02|=0.025 m/s, and 𝜑1=𝜑2=30∘. Blue: cyclohexane; green: water.

4.2. The Effect of Inclined Angles of Y-Shaped Junction

Figure 4 shows three typical images of simulated segment flow modes in the Y-shaped microchannels. Table 1 shows the relative segment length at various inclined angles of Y-shaped junction. Here, ℎ is the channel width (Figure 4) and 𝑤=(𝑤₁+𝑤₂)/2. 𝑤₁ and 𝑤₂ are the segment lengths of liquid 1 and 2, respectively (Figure 4). The two values should be same theoretically since the inlet fluxes are the same. We used average value 𝑤 as the result in Table 1 for the sake of a slight difference between 𝑤₁ and 𝑤₂ caused by some unavoidable computational errors.

Table 1 

Relative segment length 𝑤/ℎ at various inclined angles of Y-shaped junction at 𝜌₁=1000 kg/m³, 𝜇₁=0.001 Pa·s, 𝜌₂=800 kg/m³, 𝜇₂=0.02 Pa·s, 𝜎=0.05 N/m, and |𝑢₀₁|=|𝑢₀₂|=0.005 m/s.

(a) 𝜑1=𝜑2=45∘

(b) 𝜑1=𝜑2=150∘

(c) 𝜑1=0∘,𝜑2=45∘

(d) 𝜑1=0∘,𝜑2=90∘

Figure 4 

Simulated segment flow mode in a Y-shaped microchannel at various inclined angles 𝜑1 and 𝜑2.

It can be seen from Table 1, under current conditions, when 𝜑₂ is equal to 15°, the segment length is largest for both cases (𝜑₁=0^∘and 𝜑₁=𝜑₂). When 𝜑₂ is equal to 90°, although the segment length is not the smallest, it keeps to a small value for both cases (𝜑₁=0^∘ and 𝜑₁=𝜑₂). For the cases of 𝜑₁=0^∘ (Table 1 column 𝜑₁=0^∘), the value of 𝜑₂ affects the segment length greatly. When 𝜑₂<60^∘, the segment length decreases with the increase in 𝜑₂. When 𝜑₂=60^∘, the segment length reaches its smallest value. When 𝜑₂>60^∘, the segment length increases with 𝜑₂. However, when 𝜑₁=𝜑₂ (Table 1 column 𝜑₁=𝜑₂), the value of the inclined angle does not affect the slug length so much. Except the cases of 𝜑₁=𝜑₂=15^∘ and 60^∘, the segment lengths for other cases are almost the same.

Thus, the segment length can be controlled by adjusting the inclined angles of Y-shaped junction. T-shaped junction tends to generate more interfaces of two liquids (𝜃_ca=90^∘) or get the small slugs and droplets (𝜃_ca≠90^∘), which is very important for the mixing, chemical reaction, and making micro/nanosized particles using microreactor.

4.3. The Effect of Inlet Flow Rate

In this section, the effects of the inlet velocities of two phases (𝑢₀₁ and 𝑢₀₂) on segment lengths (𝑤₁ and 𝑤₂) are studied with a typical Y-shaped case of 𝜑₁=0^∘ and 𝜑₂=90^∘ (hereafter called the T-shaped microchannel). Figure 5 shows the simulated segment flow modes at |𝑢₀₁|=0.001 and 0.02 m/s. |𝑢₀₂| is not changed and kept to 0.005 m/s. It can be seen that the segment length of phase 1, 𝑤₁, increases with its inlet velocity |𝑢₀₁|, while 𝑤₂ decreases. The segment length can be controlled easily by the inlet flux of the liquids. Figure 6 shows the plot of 𝑤 against 𝑢 for the T-shaped microchannel. Here, 𝑤 is the ratio of the segment lengths of two phases, 𝑤₁/𝑤₂, and 𝑢 is the ratio of their inlet velocities |𝑢₀₁|/|𝑢₀₂|, which is equal to their inlet volume flow rates. Here, different values of 𝑢 are obtained by different |𝑢₀₁|. From Figure 6, it can be seen that 𝑤 is almost linearly proportional to 𝑢. This means the ratio between the segment volumes of two phases is almost equal to the ratio of their inlet volume flow rates. It should be noted that under current conditions, when |𝑢₀₁|≥0.04 m/s, the segment flow mode cannot be obtained. Instead, a parallel flow mode appears in the T-shaped microchannel (𝜑₁=0^∘ and 𝜑₂=90^∘).

(a) |𝑢01|=0.001 m/s

(b) |𝑢01|=0.02 m/s

Figure 5 

Simulated segment flow mode in a T-shaped microchannel at |𝑢02|=0.005 m/s, and (a) |𝑢01|=0.001 m/s; (b) |𝑢01|=0.02 m/s.

Figure 6 

The plot of 𝑤 against 𝑢 at |𝑢02|=0.005 m/s for the T-shaped microchannel of 𝜑1=0∘ and 𝜑2=90∘.

5. Conclusion

Microchannels are gradually used in microreactor systems with microreaction technology, where a typical application includes liquid-liquid multiphase flows. In the present work, a two-phase liquid-liquid segmented flow in a microchannel with various Y-shaped junctions has been studied numerically by means of volume of fluid (VOF) approach. Using liner interface construction (PLIC) technique, a sharp interface between two phases was simulated. The influence of the inclined angle of Y-shaped junction and the inlet velocity on the segmented flow length has been observed.

A segment flow mode is found in the Y-shaped microchannel when the contact angle is set to 90°. Under current situations, the segment length is largest at 𝜑₂=15^∘, and is kept to a small value at 𝜑₂=90^∘ for both cases (𝜑₁=0^∘ and 𝜑₁=𝜑₂). This means using T-shaped junction, small slug, and more interfaces of two fluids can be obtained. The ratio between the segment lengths of two liquids is almost equal to the ratio of their inlet velocities for the case of Y-shaped microchannel (𝜑₁=0^∘ and 𝜑₂=90^∘).

Nomenclature

𝐶:	Color function (−)
𝑓sv:	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 (−)
𝑡:	Time (s)
Δ𝑡:	Time interval (s)
𝑢:	Velocity vector (m/s)
𝑢:	Velocity component in 𝑥 direction (m/s)
𝑢:	𝑢=|𝑢02|/|𝑢01|(−)
𝑣:	Velocity component in 𝑦 direction (m/s)
𝑥𝑙:	The length of solution domain (m)
𝑦𝑙:	The width of solution domain (m).

Greeks

𝜂:	A constant term in the line segment equation (−)
𝜃ca:	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)
𝜑:	The angle between lateral channel and the extension central line of main channel in its opposite direction (°).

Subscripts

𝑖,𝑗:	Indicators of cell-centered value
𝑜:	Orangic phase
𝑡:	Wall tangent component
𝑤:	Water phase, or wall nomal component
𝑥,𝑦:	Components in 𝑥 and 𝑦 directions
1,2:	Indicators for fluids 1 (water) and 2 (organic fluid).

Acknowledgments

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