УДК 536.25
On Optimization of Mixing Process of Liquids in Microchannels
Andrey V. Minakov*
Siberian Federal University, Svobodny 79, Krasnoyarsk, 660041, Russia
Valery Ya. Rudyak
Novosibirsk State University of Architecture and Civil Engineering,, Leningradskaya 113, Novosibirsk, 630008,
Russia
Andrey A. Gavrilov Alexander A. Dekterev
Institute of Thermophysics SB RAS, Lavrentyev 1, Novosibirsk, 630090,
Russia
Received 10.12.2009, received in revised form 10.02.2010, accepted 20.03.2010 A method for simulating fluid flows in microchannels is proposed. The method is tested using available experimental data obtained in micro-PIV studies of microchannel flows. Flow regimes in Y-, T- and S-types micromixers are studied. Passive and active mixers are considered. The dependence of the mixing efficiency on the Peclet number and the Reynolds number is examined, and the possibility of using hydrophobic and ultrahydrophobic coatings is analyzed.
Keywords: Microchannels, micromixers, hydrodynamics, mixing.
Introduction
Fluid mixing in microchannels is an extremely important problem in various applications of microflows. In macroscopic flows, mixing usually occurs in a turbulent flow regime. Microflows, however, are usually laminar and mixing under standard conditions is carried out by means of diffusion processes. Because of the extremely low values of the molecular diffusion coefficient D (D ~ 10-9 ^ 10-11 m2/sec) this channel of fluid flow mixing is extremely inefficient. To increase the mixing rate, it is necessary to use special devices - micromixers. Micromixers are therefore key devices in various MEMS applications.
There are passive and active methods of increasing the mixing rate. In the first case, the geometrical shape of the channels is varied; various types of inserts are used, etc. In the second case, various external fields (acoustic, electric, magnetic) are applied, and the fluid flow rate is varied. The purpose of the present paper is to develop a method for simulating microflows based on solutions of hydrodynamic equations and to study the mixing efficiency for both mixing methods. Several types of mixers are considered. Flows are simulated using the software package for solving the Navier-Stokes equations developed by Rudyak et al. [1] (2008) and Minakov
* [email protected] © Siberian Federal University. All rights reserved
[2] et al. (2008). This package has been tested on available data obtained using the micro-PIV technique. Both stationary and quasistationary flows with tunable velocity profiles are considered. To study the effect of hydrophobic or even ultrahydrophobic coatings on mixing and pressure drop the slip boundary conditions are used.
1. Numerical Simulation Method
It is known that, under normal conditions, flows of both liquids and not very rarefied gases can be described by methods of continuum mechanics. For microchannels, however, the situation changes drastically. The term microchannel usually refers to a channel in which one of the characteristic sizes h (for example, the height of a flat channel or the diameter of a cylindrical channel) is in the range from 1 to approximately 300 pm. Under these conditions, liquid and gas flows should be described differently, as a rule. In the present paper we consider the flows of the homogeneous liquid in channels with the smallest size more than or equal to 50 pm. Such flows are described within the framework of the continuum model using the Navier-Stokes equations with no-slip or slip boundary conditions. A multicomponent flow is simulated using the one-velocity hydrodynamics. Incompressible fluid flow is described by the Navier-Stokes equations
dp + V■ (p v) = 0, (1)
dpV + V- (Pvv) = -Vp + VT, (2)
where p is the fluid density, p is the pressure, v is the fluid velocity, and T is the viscous stress tensor, whose components are defined as follows:
T / duj dui
ij P \ dxi dxj
where p is the viscosity of the mixture, ui are the components of the velocity vector.
The mixture viscosity is computed on the basis on a simple mass fraction average of the pure species viscosities pi:
P = 53 fiPi, i
here fi is the mass fractions of the mixture components.
The density of the mixture is expressed in terms of the mass fractions fi of the mixture components and the partial densities pi
-1
p ^ ' " '
Efi/pi)
The transport equation for the components of the mixture in mass fractions is written as
^p/ + V • (pfiv) = V • (p AV/i), (3)
where Dj is the diffusion coefficient of the ¿-component.
On a stream input the Dirichlet boundary conditions are used. The three components of the velocity vector and the concentration of mixture component are considered set: v = vin, fi = fin, where vin is the value of a vector of velocity on an input, fin is the value of a mass fractions components on an input.
<9$
onditions are set:
where n is the vector of an external normal, $ are scalar quantities.
On output boundary for all considered quantities Neumann's conditions are set: -7— = 0,
on
For hydrophilic surface on the solid wall non-slip conditions are set, normal and tangential components of the velocity on walls are equal to zero: vn = 0, vt = 0.
dvt
For hydrophobic surface on the solid wall slip conditions are set: vn = 0, vt = b—— ,
dn wall
where b is the slip length.
For mass fractions of mixture components on all solid walls Neumann's conditions are set:
dfi
—— = 0, here n is normal vector to a wall. on
The difference analogue of the convection-diffusion equations (1)-(3) is found using a finite-volume method [3] for structured multiblock grids. In this case, the resulting scheme is automatically conservative. The method consists of partitioning the calculated domain into control volumes with the subsequent integration of the initial conservation equations over each control volume to obtain finite-difference relations.
The convective terms of the transport equations (2) and (3) are approximated by the second-order upwind QUICK [4] and TVD schemes [5], respectively. A second-order implicit scheme is use to approximate the unsteady terms of the hydrodynamic equations.
Diffusion fluxes and source terms are approximated by finite-volume analogs of the central-difference relations of second-order accuracy. The coupling between the velocity and pressure fields, which provides satisfaction of the continuity equation (1), is achieved by the SIMPLEC algorithm on aligned grids [6]. Oscillations of the pressure fields are prevented using the Rhie-Chow method involving the introduction of a pressure smoothing term into the equations for pressure correction [3]. The difference equations obtained by discretization of the initial system of differential equations are solved by an iterative method using an algebraic multigrid solver.
The algorithm described here has been employed to solve a wide range of problems of external and internal flow [1, 2, 7]. The applicability of this algorithm was investigated earlier [8] (Gavrilov et al., 2009). Because below the Y-type mixers will be analyzed mainly only one example of comparison of numerical and experimental data is here considered namely for this type mixer. Kim et al. [9] studied mixing of water and a glycerol solution in this type micromixer. Inlet channels were each 150 ^m wide, and the mixer was 300 ^m wide. The heights of the channels were identical and equal to 50 ^m. In experiment [9], water was mixed with 20 %, 50 %, and 60 % glycerol solutions, whose density and viscosity coefficient increase successively. The calculations of mixing in the Y-type micromixer for 20 % and 60% glycerol solutions shown in Fig. 1. In all cases, the calculation is in good agreement with the experimental data. Fig. 2 gives a comparison of the longitudinal flow velocity profile in section C for 20 % and 60% solutions. Here the solid curve corresponds to calculated data and filled circles to experimental data. The calculated and experimental data are in a good agreement with each other.
a) 20 % glycerol solutions b) 60 % glycerol solutions
Fig. 1. Mixing of water and a glycerol solution in the Y-type micromixer
Fig. 2. Comparison of calculated and experimental values of longitudinal flow velocity in the section C (see Fig. 1)
2. Mixing in an Y-type Mixer
An Y-type micromixer is one of the simplest models. Mixing in this micromixer has been studied extensively during the recent decade. However, systematic data on the effect of mixing on flow parameters are not available. For this reason, we start the investigation of micromixers from this model.
So, we consider the mixing of two identical fluids in the micromixer presented in Fig. 1. Let, for definiteness, the first fluid enter channel B, and the second fluid enter channel A. If the flow rates in both inlet channels are identical, the flow Reynolds numbers in channels A and B are identical too and equal to Rew = pUqw/p, where p is the viscosity coefficient of fluid. The mixing time is determined only by the fluid diffusion coefficient D and is of order Tm ~ w2/D. If the channel has length L, the time required for a fluid particle to pass through the channel is of order tl ~ L/Uq ~ (pLwh)/Q. The mixing efficiency in the channel is determined by the ratio of these two times
TL
w2UQ w „ v
-Q--Rew —,
DL L wD1
(4)
where v = p/pis the kinematic viscosity coefficient. Therefore, for the given channel and given fluids, the mixing time is proportional to the Reynolds number. If one uses water as the fluids being mixed, then, according to (4), the mixing in a channel of length L = 1cm is efficient if
m
Tm < 10Rewtl, which is only valid for very small Reynolds numbers. From this, in particular, it follows that, starting from some values of the Reynolds number, the mixing efficiency changes little.
From Eq. (4), it is easy to obtain another estimate
Tm W
— ^ — Pe, tl L
(5)
where Pe is the Peclet number Pe = wUq/D (or the diffusion Peclet number). For a given Peclet number, the mixing length Lm ~ UQTm ~ Pew is the smaller the smaller the channel width, and this length increases with increasing Peclet number. It should be noted, however, that the dependence of the mixing efficiency on the Peclet number is nonlinear. It is given in Fig. 3 for Y-mixer considered (Fig. 1). Here the mixing efficiency M is determined by the formula
M
W
1 - 2 j(f - 0.5)dy
• 100%,
(6)
where W = 300^m is the width of the channel in which mixing occurs (Fig. 1) and f is the concentration of mixture in the section at 4000 ^m from the site of flow merging. This dependence is well described by the formula
M
113.2
1 +0.037 Pe0-7'
In Fig. 3, it is shown by the solid curve and the points correspond to calculations.
Fig. 3. Mixing efficiency in a Y-type mixer versus Peclet number
Next, we studied the effect of the angles a and p. Both angles were varied from 0o to 90o(see Fig. 1, in last case we have so named T-type mixer). The degree of opening of channels A and B was found to have little effect on the mixing efficiency.
Since diffusion is the most important factor of the mixing process in microchannels, the channel length should be sufficiently large to provide flow mixing. Naturally, this leads to a significant pressure drop due to friction on the channel walls. On the other hand, it is obvious
that such losses are much lower in the case of a flow slip rather than no-slip conditions on the channel walls. At the macroscopic level of description, the boundary conditions are always specified as no-slip conditions. Slip boundary conditions are only employed to describe flows of fairly rarefied gas. This is due to the fact that, in a gas, the characteristic wall slip length b is of the order of the Knudsen number b ~ Kn (here slip length b is determined by the following relation v = b(dv/dy), where v is the flow velocity along the x axis). In macroscopic fluid flows under usual conditions, the slip length varies from a few nanometers to several tens of nanometers. This slip can be ignored. In microflows, however, such a slip can be significant. In addition, in microflows, the slip length can reach a few hundred nanometers. This is related to a change in the short-range order of fluid near surfaces, the possibility of gas release on the channel walls, etc.
Walls slip leads to a significant decrease of the friction on the channel walls and, hence, to a reduction in the pressure drop. To increase the slip length, it is common to use various hydrophobic or even ultrahydrophobic coatings (see, for example [10, 11]). In this case, the slip length can reach tens of micrometers. Systematic calculations performed for a Y-type micromixer have shown that as the slip length increases to 5 pm, the pressure loss decreases by almost half. It should be emphasized that, in the presence of wall slip, the mixing efficiency changes only slightly.
Table 1. Parameters of the mixing liquids
p, k/m3 p, Pa sec Dx10 9 m2/sec
Acetone 800 0,00029 4,56
Glycerol 1260 1,48 0,0083
Ethanol 790 0,0012 1,24
Isopropyl alcohol 786 0,0029 0,38
Till now we discussed the mixing of the same fluids. In practice the mixing the different liquids is interesting. Certainly our algorithm permits to study these situations too. As an example we present the mixing of water with different liquids (see table 1) in Y-type mixers. The water is supplied through the entry A (Fig. 1). Because the physical characteristics of these liquids are very different the mixing efficiency is essentially differed too. In particular, the parameter M for acetone, glycerol, ethanol and isopropyl alcohol is equal to 36%, 3%, 18% and 9%, respectively. In these cases the pressure drop between exits A and B is different and it is changes from 239.5 kPa for the glycerol till 165 Pa for acetone.
3. Mixing in a T-type Mixer with Inserts
The relations (4), (5) give an estimate of the mixing time and length and, hence, the efficiency of micromixers. From general considerations, it is clear that the mixing time can be reduced considerably by multiple splitting of the mixing flow. This principle underlies the design of lamination mixers (see, for example [12]). A similar idea is to place a number of inserts in mixing flow to change the flow structure and thus to accelerate mixing. It is clear that, in a laminar flow, symmetric arrangement of inserts is of little effectiveness. In addition, the characteristic size of the inserts should be comparable to the sizes of the channel. Moreover, the mixer is easy to
optimize for this parameter. Series of such calculations for various inserts were performed in the present paper. As an example, Fig. 4 shows mixing of two fluids in a mixer with asymmetrically arranged three and seven rectangular inserts and without inserts. The mixing channel is 100 ^m wide, the inlet channels are 50 ^m wide, and the channel height is 50 ^m. The Reynolds number Re = 2, and the Peclet number Pe= 5 • 103. All inserts are identical and have a cross-sectional area of 50 ^mx20 ^m. So far as with increasing number of inserts length of channel also increase, it is reasonable to analyze the specific mixing efficiency and the specific pressure loss. This data is presented in the Table 2.
Table 2. Specific mixing efficiency and specific pressure loss in a mixer with rectangular inserts
N M %/>m. Ap, Pa/^m.
0 0.0094 0.37
3 0.026 0.92
5 0.035 1.44
7i 0.045 1.92
72 0.085 7.38
Fig. 4. Comparison of mixing of two fluids in a mixer of T-type with asymmetrically arranged inserts
The mixing efficiency increases with increasing number of inserts, but the pressure loss also increases. It can be reduced by using hydrophobic coatings. Thus, for a slip length b =1 ^m, the pressure loss can be reduced by 10 %. The mixing efficiency can be increased if the inserts connect with the wall. Such situation presents on the last picture in Fig. 4. Here the length of the insert is equal to 70 ^m (last row in Table 2).
Mixing efficiency depends on the distance between inserts. In the general case this dependence is not monotonous. As an example in Fig. 5 the dependence of mixing efficiency (Fig. 5a) and pressure drop (Fig. 5b) from distance between inserts for micromixer with five inserts is presented. The inserts size is a 50 ^mx20 ^m. The Reynolds number Re = 2, and the Peclet number Pe = 5 • 103. So far as with increasing of distance between inserts length of channel also
increase, here reasonable to analyze specific mixing efficiency and specific pressure loss. On Fig. 5 the efficiency of mixture and pressure loss divided into length of the channel are given. The distance between inserts is divided into a thickness of an insert a = 20pm. As Fig. 5a shows, specific mixing efficiency reaches a maximum at a distance between inserts equal 40 microns.
Fig. 5. Specific mixing efficiency (a) and specific pressure loss (b) versus distance between inserts
For the given type of a micromixer dependence of mixing efficiency on Reynolds number Rew = pUqw/p has been investigated. The number of inserts equaled 5, the distance between inserts has been chosen to be 20 microns. The Peclet number has been fixed and equaled 5000. The Reynolds number varied by means of the mass flow rate. Dependence of mixing efficiency on Reynolds number at the fixed value of Peclet number Pe = 5000 is presented on Fig. 6a. Apparently, mixing efficiency at a creeping flow (Re<2) practically does not depend from Reynolds number. At higher values of the Reynolds number vortex zones are formed behind inserts, that considerably improves mixing. However, regimes of flow with Re>10 in microchannels meet seldom. Thus, it is possible to make a conclusion, that the variation of the Reynolds number does not lead to mixing substantial improvement in microchannels.
Re Pe
a) b)
Fig. 6. Mixing efficiency in a T-type mixer versus Reynolds (a) and Peclet number (b)
For the chosen configuration of T-type micromixer dependence of mixing efficiency on the Peclet number is presented in Fig. 6b (Re = 2). It is clear that changes of the Peclet number lead to considerable change of mixing efficiency. With growth of the Peclet number a diffusion influence decreases and mixing efficiency mixer strongly falls.
4. Mixing in a S-type Mixer
One more widespread kind of a micromixer is an S-type mixer. In Fig. 7 the typical geometry of such a mixer is presented. For numerical simulation a mixer with the following geometrical
Fig. 7. Mixing of two fluids (a) and pressure on the wall (b) in a mixer of S-type
sizes has been chosen: the width of the mixing channel 100^m, the height of channel 50^m, the length of the channel perpendicular axes of the mixer 200 ^m, the length of the channel parallel axes of the mixer 300 ^m. For a given micromixer the dependence of the mixing efficiency of two Newtonian liquids on the number of sections of the mixing channel has been investigated. In Fig. 8 mixing (Fig. 8a) and pressure drop (Fig. 8b) in the channel with four curvilinear sections is shown. All calculations were carried at the following parameters: Re = 2, Pe = 5000. So far as at addition of curvilinear sections in the channel its length varies, here as for a micromixer T-type we will analyze specific mixing efficiency and specific pressure loss.
The analysis shows that specific mixing efficiency attains a maximum for a mixer with one curvilinear section. The growth of mixing efficiency with the increase of quantity of section for a given type of mixer appears to be insignificant. Therefore the specific mixing efficiency decreases as the quantity of sections increase. Comparison of specific mixing efficiency for S- and T-mixers (Fig 5a, Table 2 and Fig. 8a) shows that for a T-mixer the given value on the average in some time above. However, specific pressure loss of T-type mixers exceeds specific pressure loss of S-type mixers in many times (Fig. 5b, Table 2 and Fig. 8b).
Conclusion
The numerical algorithm for simulation of processes of hydrodynamics and diffusion of liquids in microchannels is developed and adapted. The present algorithm considers unsteady and three-dimensional liquid flow in microchannels, adapts for real geometry of any complicated object, it is easily generalized on two-phase and not isothermal problems. Testing of the developed numerical algorithm for the solving of some hydrodynamics problems in microchannels has shown
m
Fig. 8. Specific mixing efficiency (a) and specific pressure loss (b) versus sections number of channel
its high efficiency. Numerical solutions of all considered test problems are in good agreement with experimental data. It shows accuracy of numerical algorithm and confirms validity of application of the equations of Navier-Stokes for simulation of flow fluids in microchannels.
By the numerical simulation the process of mixing of liquids in microchannels has been investigated. The most widespread forms of channels applied to mixing improvement have been thus considered. Research of processes mixing of liquids in the micromixers Y, T, S - type has shown, what even insignificant optimization of the form of channels allows to considerably increase the mixing efficiency. In certain cases this increase attains hundreds of percents, at insignificant growth of resistance. For a reduction of resistance of microchannels the influence of hydrophobic and ultrahydrophobic coverings of walls of the channel has been studied. Essential falling of resistance is shown at almost invariable characteristics of mixing.
In the future it is planned to essentially expand the area of applicability of the mathematical model and numerical algorithm, because in practice in microchannels the greatest interest is represented by two-phase, not isothermal flow, and processes of chemical reaction.
This work was supported in part by the Russian Foundation for Basic Researches (grant No. 07-08-00164) and by the grant of the President of the Russian Federation for Support of Leading Scientific Schools (project No. NSh-454.2008.1).
References
[1] V.Ya. Rudyak, A.V. Minakov, A.A. Gavrilov, A.ADekterev, Application of new numerical algorithm of solving the Navier-Stokes equations for modeling the work of a viscometer of the physical pendulum type, Thermophysics & Aeromechanics, 15(2008), 333-345.
[2] A.V. Minakov, A.A. Gavrilov, A.A. Dekterev, Numerical algorithm for solving spatial problems of hydrodynamics with moving solids and a free surface, Siberian J. Industrial Math., 11(2008), no. 4, 94-104 (in Russian).
[3] J.H. Ferziger, M. Peric, Computational methods for fluid dynamics, Springer Verlag, Berlin, 2002.
[4] B.P. Leonard, A stable and accurate convective modeling procedure based on quadratic upstream interpolation, Camp. Math. Appl. Mech. Eng., 19(1979), 59-98.
[5] F.S.Lien, M.A. Leschziner, Upstream monotonic interpolation for scalar transport with application to complex turbulent flows, Int. J. Num. Math. Fluids, 19(1994), 527-534.
[6] S.V. Patankar, Numerical heat transfer and fluid flow, Washington, DC, Hemisphere, 1980.
[7] A.V. Minakov, A.A. Gavrilov, A.A. Dekterev, V.Ya. Rudyak, Numerical simulation of operation process of turning-fork vibroacoustic viscosity sensor, J. Siberian Federal University., Mathematics & Physics, 2(2009), no. 4, 457-469 (in Russian).
[8] A.A. Gavrilov, A.A. Dekterev, A.V. Minakov, V.Ya. Rudyak, Micromixers simulation, Proc. of 2nd Micro and nano Flows Conf., West London, Brunel Univ., 2009, MNF32-1-MNF32-8.
[9] W.J. Kim, Y.Z. Liu, H.J. Sung, MicroPIV in two-fluid flow, Proc 5th Int. Symp. on Particle Image Velocimetry, Busan. Lauga E., Brenner M.P. Stone H.A. 2005, 8.
[10] J.Ou, B.Perot, P. Rothstein, Laminar drag reduction in microchannels using ultrahy-drophobic surfaces, Phys. Fluids., (2004), no. 16, 4635-4643.
[11] Microfluidics: the no-slip boundary condition. Handbook of experimental fluid dynamics, Eds. Foss J., Tropea C., Yarin A.Ch., 2006, 150.
[12] R. Karnik, Microfluidic mixing encyclopedia of microfluidics and nanofluidics, Ed. Li D., Springer, 2008, 1177-1186.
Об оптимизации процесса смешения жидкостей в микроканалах
Андрей В. Минаков Валерий Я. Рудяк Андрей А.Гаврилов Александр А. Дектерев
В работе предложена математическая модель и численный алгоритм для моделирования течений жидкости в микроканалах. Приведены результаты тестирования численной методики и сопоставления полученных результатов с экспериментальными данными микро-PIV измерений. Изучены режимы течения и закономерности смешения жидкостей в микроканалах Y-, T-, S-типов. Для данных видов каналов получены зависимости эффективности смешения от чисел Рейнольдса и Пекле. Проведена оптимизация процесса смешения в микроканале Т типа. Показана возможность использования гидрофобных и ультрагидрофобных покрытий для уменьшения сопротивления каналов.
Ключевые слова: микроканалы, микромиксеры, гидродинамика, смешение.