Water–ethanol mixing in T-shaped microdevices

Chemical Engineering Science 95 (2013) 174–183

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Water–ethanol mixing in T-shaped microdevices Gianni Orsi, Mina Roudgar, Elisabetta Brunazzi, Chiara Galletti, Roberto Mauri n Department of Chemical Engineering, Industrial Chemistry and Materials Science, University of Pisa, Largo Lucio Lazzarino, 2,56126 Pisa, Italy

H I G H L I G H T S







Water–water (W–W) and water–ethanol (W–E) flows in a T-shaped micromixer are simulated. The model accounts for volume of mixing and fluidity of mixing of W–E mixtures. Water and ethanol streams are separated by a viscous layer. Engulfment occurs at larger Reynolds number for W–E (Re≅230) than for W–W (Re≅140).

art ic l e i nf o

a b s t r a c t

Article history: Received 26 July 2012 Received in revised form 6 March 2013 Accepted 11 March 2013 Available online 22 March 2013

Numerical simulations were performed to study the mixing dynamics of two miscible liquids within a Tshaped micromixer, comparing the case where the two inlet fluids are both water (W–W) with that where they consist of water and ethanol (W–E). We showed that at smaller Reynolds number, Reo 100, there is no vortex formation for either cases; therefore, mixing water and ethanol is slightly easier than mixing water with water because the residence time of the fluid occupying the interfacial region in the W–E case is larger than that of the W–W case. On the other hand, at larger Reynolds numbers, mixing water and ethanol may take considerably longer, as the onset of engulfment is retarded and occurs at larger Reynolds number. All these behaviors are due to the fact that, as a W–E mixture has a viscosity that is up to almost three times larger than that of water, at the confluence of the T-mixer the two streams are separated by a viscous layer of a W–E mixture, that hampers vortex formation and retards mixing. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Microfluidics Mixing Diffusion Water–ethanol mixture Simulation Multiphase flow

1. Introduction Mixing two different fluids in a micromixer is one of the most basic processes in microfluidics. Due to the small size of the device, pressure driven flows in simple channels (i.e. with smooth walls) are laminar and mostly uniaxial, so that confluencing liquids tend to flow side by side. On the other hand, as the Schmidt number (i.e. the ratio ν/D between kinematic viscosity and molecular diffusivity) is typically very large, the Peclet number is large, that is Pe¼ Ud/D⪢1, where U is the mean fluid velocity and d is the hydraulic diameter. Accordingly, in the absence of any transverse convection, complete mixing is reached after the distance L≅d Pe along the channel that can be prohibitively long. The easiest way to overcome this difficulty and enhance mixing is to induce transverse flows through clever geometries. The simplest one consists of a T-shaped micromixer, in which two inlet flows join the main mixing channel through Tshaped branches.

n

Corresponding author. Tel.: þ39 050 2217848; fax: þ39 050 2217866. E-mail address: [email protected] (R. Mauri).

0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.03.015

T-shaped micromixers have been investigated extensively in recent years (see, for example, Nguyen and Wu, 2005; Hessel et al., 2005; Kim et al., 2004; Aubin et al., 2005), as they are quite suitable to carry out fundamental studies to understand mixing at the microscale. As such, most of the investigators have confined themselves to the mixing of two identical liquids, assuming that in one of them a dilute solute is dissolved. Such systems were studied by Engler et al. (2004) both numerically and experimentally, finding that at the confluence of the two streams three flow regimes may occur, depending on the Reynolds number, namely stratified laminar flow, vortex flow and engulfment flow. As noted by Bothe et al. (2006, 2008), the appearance of an engulfment flow, with intertwinement of the inlet fluid streams, can achieve an efficient mixing by rolling up the initial planar contact area between the fluids and thus enhance mixing by reducing the diffusion lengths. Hoffmann et al. (2006) performed experimental investigations on T-shaped micromixers with rectangular microchannel in the stratified and engulfment flow regimes, reporting that the transition between them may occur even at rather low Reynolds numbers. In addition, as noted by Galletti et al. (2012), the morphology of the fluid flow at the confluence, and consequently the mixing efficiency as well, is strongly influenced

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

by the inlet flow conditions, as a symmetry break is observed when the length of the inlet channels is progressively decreased, so that the inlet flow at the entrance ceases to be fully developed. Surprisingly, much less investigation has been devoted to studying the mixing process between two different, albeit miscible, fluids from a fundamental perspective, despite two-fluid mixing being an essential process in many microfluidic devices. For example, various pharmaceutical, biomedical and biochemical processes involve the mixing of two fluids, including nano-drugs preparation by antisolvent methods, DNA purification, polymerase chain reaction (PCR), enzyme reaction, and protein folding, whose performance relies on the effective and rapid mixing of samples and reagents. So, for example, in the works by Wang et al. (2002), Jeon et al. (2005), Lin et al. (2007), Chung and Shih (2008), Ansari et al. (2010), Lin et al. (2011), where W–E mixtures were considered, the main objective was determining the best geometrical set up to optimize mixing. On the other hand, the present study focuses on the comparison between the case where the two inlet fluids are both water (W–W case) and that when they consist of water and ethanol (W–E case), using a simple T-type passive micromixer with Reynolds numbers ranging from 1 to 300. In order to do that, a commercial Computational Fluid Dynamic (CFD) code, Fluent 12.0 by Ansys Inc., was used to solve the three-dimensional equations of motion for a non-ideal W–E mixture that takes the non-ideality of the W–E mixture into account. In fact, the influence of the mixture non-ideality, above all its strong increase in viscosity, has not been accounted for in all previous works. At the end, our simulation results will be compared with the experimental findings by Wang et al. (2012), who, to the authors' knowledge, constitutes the only previous work on this subject.

ρu⋅∇ϕ ¼ ∇:ðρD∇ϕÞ:

2.1. Governing equations Consider two fluids converging into a T junction. The two inlet streams have the same temperature, so that, as the heat of mixing has a negligible effect here (see below), we may assume that the process is isothermal. In general, the governing equations at steady state are: ∇⋅ðρuÞ ¼ 0,

ð1Þ

ρu⋅∇u þ ∇p ¼ ∇⋅½μð∇u þ ∇u þ Þ þ ρg,

ð2Þ

ð3Þ

Here, u denotes the velocity vector, ρ the fluid density, p the pressure, μ the viscosity, g the gravity acceleration, D the molecular diffusivity and ϕ the mass fraction of one of the two inlet fluids, say fluid A. If the two fluids are identical, we can imagine adding a very small amount of contaminant, i.e. a dye, to one of the fluids (which therefore continue to have the same physical properties) and therefore ϕ indicates the (normalized) dye mass fraction. In addition, both density ρ and viscosity μ are known functions of the mass fraction only, as temperature is constant, while their dependence on pressure can be easily neglected. As mentioned in the introduction, the characteristics of the velocity and concentration fields can be described through the Reynolds and Peclet numbers Re ¼

Ud ; ν

Pe ¼

Ud ¼ Re Sc; D

Sc ¼

ν , D

ð4Þ

where U is the mean velocity, while the characteristic fluid length d is assumed to be the hydraulic diameter, i.e., d¼

2WH , ðW þ HÞ

ð5Þ

where W and H are the mixing channel width and height, respectively (see Fig. 1). The most delicate part of this simulation concerns the dependence of both viscosity, μ, and density, ρ, of the W–E mixture on the mixture composition, ϕ. In fact, we know that ρ(ϕ) and μ(ϕ) are far from being linear functions, as it is assumed instead in the default setting of the Fluent 12.0 CFD code. For example, as density is the inverse of the specific volume, we know that, in general: 1 ϕ 1−ϕ ¼ þ þ Δvmix ðϕÞ, ρ ρA ρB

2. Simulation technique

175

ð6Þ

where ρΑ and ρB are the densities of pure fluids A and B, respectively, while Δvmix(ϕ) is the volume of mixing, that is the volume, per unit mass, that the mixture gains (or loses) upon mixing the two components with a given mass fraction ϕ. When Δvmix ¼0, as for example in regular mixtures (i.e. mixtures behaving like van der Waals fluids), volumes are conserved. That, however, does not apply to a W–E mixture, as it is well known that upon mixing 1 l of water with 1 l of ethanol, the volume of the resulting mixture is 1.9 l, corresponding to a 5% loss of volume (Perry and Green, 1997).

Fig. 1. Schematics of the T-micromixer.

176

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

Regarding viscosity, being it the resistance of a fluid against the diffusive transport of momentum, it is known (Kern, 1965) that for an ideal mixture, we have: m−1 ¼ϕmA−1 þ(1−ϕ)mB−1, where μΑ and μB are the viscosities of pure fluids A and B. Therefore, we can write 1 ϕ 1−ϕ ¼ þ þ Δf mix ðϕÞ, μ μA μB

ð7Þ

where Δfmix(ϕ) is a fluidity of mixing, taking into account the nonideality of the mixture. This term is particularly important in our case, since at 20 1C a mixture of ethanol and water with 0.3 oϕ o0.6 has a viscosity which is almost three times that of pure water (Simmonds, 1919). Similar behavior is observed for many aqueous mixtures of organic solvents, such as acetone, methanol, propanol and acetic acid (see Dizechi and Marschall, 1982; Wang et al., 2004; Laliberté, 2007). Finally, as mentioned above, we assume that the process is isothermal, therefore neglecting the temperature change due to the heat of mixing, Δhmix. Now, first of all, mixing water and ethanol is an endothermic process, with Δhmix≅−1.6  104 J/kg for a 50% and 50% mixture (Christensen et al., 1984). Therefore, the temperature of the mixture would actually decrease upon mixing, thus increasing viscosity and enhancing the effects that we study here. In addition, the temperature change due to the heat absorbed during the process is small enough to be negligible. In order to see that, assuming that water and ethanol mix completely, consider that the maximum change of temperature, ΔTmax, between inlet and outlet, occurs when the channel walls are _ Δhmix , equals the adiabatic. In that case, the heat generated, m _ ΔT, where m _ is the mass flow rate and c the heat convected, mc specific heat, yielding ΔTmax≅−4 1C. Therefore, considering that in our case the degree of mixing is at most 20%, ΔTmax reduces to about −1 1C. In this heat balance we have neglected heat conduction and transverse heat convection (i.e. towards the walls). The former, kL ΔT, where k is the thermal conductivity, is much smaller _ Δhmix , and can therefore be neglected, than the heat generated, m while the latter is very important in the vortex and engulfment regime, thus further reducing all temperature non-uniformities. So, we may conclude that the effect of the heat of mixing consists of a temperature decrease along the mixing channel of at most 1 1C, which amounts to a negligible change in the properties (i.e. density and viscosity) of the mixture.

recommendation by Hussong et al. (2009), regarding how to accurately describe the velocity fields. A second order discretization scheme was used to solve all the governing equations; using higher-order scheme did not lead to any significant change, thereby ruling out the importance of numerical diffusion. Simulations were typically considered converging when the normalized residuals for velocities were stationary with iterations and fell below 1  10–6, although smaller residuals were required in specific cases, especially near the engulfment (Galletti et al., 2012). The steadiness of the solution with iterations was also assessed by checking the velocity and concentration field in the outlet section of the mixing channels at different iterations. In addition, a grid independence study was performed on the velocity field, showing that, repeating the numerical simulation with smaller cubical elements (i.e. with H/ 64 edge), the simulations performed at maximum flow rate, Re¼ 300, give results (in particular, the degree of mixing) that do not differ from the original simulation. A quite similar analysis was also conducted by Roudgar et al. (2012), using a slightly different geometric setup. As mentioned above, here we discretize the concentration field using the same mesh size as the velocity field. Clearly, since in our

2.2. Geometry and numerical model The geometric setting of our simulation consists of the T-shaped microdevice shown in Fig. 1; it is positioned horizontally, so that gravity acts along the z-axis. The mixing channel has a rectangular cross section with a 2:1 aspect ratio (i.e. width W ¼2H and depth H, so that d ¼4H/3), with length L ¼5H, while the inlet channels are identical, with square cross section (i.e. width W'¼ H and depth H) and length L'¼ L¼ 5H. In particular, the dimensions indicated in all our figures are referred to the particular case with H ¼100 μm, as this corresponds to the geometric setup used by many investigators (see, for example, the numerical work by Bothe et al. (2006, 2011), and the experimental work by Hoffmann et al. (2006)). In addition, flow rates were varied so that Re in the mixing channel ranges from 1 to 300, while the Schmidt number is quite large, Sc ¼ Oð103 Þ. Simulations were performed using the commercial software Ansys 13, with the FLUENT fluid dynamics package, where length and velocity were scaled in terms of d and U, i.e. yn ¼ y/d and un ¼u/ U. The grid consists of cubic elements with H/34 edge, corresponding, approximately, to the viscosity-dependent typical dimension μ/ρU. This leads to a mesh of 34  68 elements in each cross section of the mixing channel, thus in agreement with the

Fig. 2. (a) Density vs. ethanol mass fraction and (b) viscosity vs. ethanol mass fraction. Dotted lines show the mixture density and viscosity when a linear dependence on mass fraction is assumed.

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

case Sc ¼ Oð103 Þ, so that Pe⪢Re, a rigorous simulation of mass diffusion seems to require a grid resolution that is much finer than that of the velocity field (see the extensive discussion in Ottino (1994), about this issue). This, however, would lead to unacceptably long computation times, which can be resolved only by using clever numerical schemes, such as the one used by Bothe et al. (2011). Using, as we do here, the same mesh size, effectively amounts to studying mass diffusion on a coarse-grained lengthscale. Clearly, this approximation is unacceptable when mixing is a diffusion-driven process, as it happens when, in the W–W case at Re⪡1, the two inlet steams flow side by side, so that a concentration boundary layer of thickness δ∼d Pe−1/2 forms at the interface. However, in our case (even at Re¼1, as we will see), convection at the confluence occurs also along a transverse direction, i.e., perpendicularly to the mean flow, as each inlet flow penetrates into the other by a typical distance, R (as we saw, R≅d in the engulfment regime). This is particularly true in the W–E case, where a transverse flow is induced by the term (dm/dϕ)∇ϕ  ∇u appearing in Eq. (2), so that purely diffusive fluxes can be neglected compared to the convective fluxes induced by transverse fluid motion. Therefore, we may conclude that an approximately correct concentration field can be determined by assum-

177

ing that the fluid particles are passive tracers and that we can follow their motion in the coarse-grained scale of the velocity field. Clearly, the explanation that we have given above to justify the use of a coarse-grained concentration field is heuristic; in the Appendix, a multiple scale analysis is carried over to justify it rigorously. There, we see that the leading-order term of Eq. (3) at the microscale (i.e. over distances of O(δ)) is: u⋅∇φ ¼ 0, showing that the concentration remains unchanged when we move along the fluid streamlines. Now, in the W–W case, as viscosity is uniform, the Navier–Stokes Eqs. (1) and (2) can be solved independently of the concentration field, and therefore the concentration gradient will be locally perpendicular to the velocity field, as it happens when the two inlet streams flow side by side. In the W–E case, instead, Eqs. (2) and (3) are intrinsically coupled to each other, and therefore u⋅∇φ ¼ 0 implies that ϕ is uniform within the microscale, therefore justifying the fact that u and ϕ vary on the same, coarse-grained scale. In any case, it should be stressed that our simulations allow only a qualitative analysis of the concentration fields and the mixing quality, while further studies (and experiments, above all) are required to better quantify them.

Fig. 3. Ethanol mass fraction contour plots for W–E systems: (a) at Re¼0.1 both close to the bottom wall and at the centerline of the mixing channel; at the outlet cross-section when (b) Re¼ 0.1; (c) Re¼1 and (d) Re¼10.

178

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

2.3. Boundary and thermodynamic conditions The boundary conditions consist of no-slip velocity and nomass-flux at the channel walls, a constant ambient pressure at the exit, while at the entrance a given flow profile was imposed, assuming fully developed fluid conditions. In this way, we avoid the undesired complications connected to the inlet fluid conditions, that we have analyzed in a separated work (Galletti et al., 2012). We also checked that the length of the mixing channel did not influence the flow pattern at the confluence, by comparing our case, with L ¼5H, with that for L ¼15H in the most critical case, namely the W–E mixture at Re¼300, verifying that the concentration profile and the degree of mixing at a 4H distance from the confluence is the same in the two cases. The fully developed velocity profile in a closed rectangular conduit can be easily obtained by solving the Navier–Stokes equations with no-slip boundary conditions at the walls and a constant axial pressure gradient G. For our purposes, the most convenient form of this solution is (Happel and Brenner, 1973; Chatwin and Sullivan, 1982): uðx,zÞ ¼ −


 G 4GX 2 1 x z kπ z ∑ sinðkπ Þ Coshðkπ Þ−Tanhð ÞSinhðkπ Þ , xðX−xÞ− 2μ X X 2η X μπ 3 k odd k3

ð8Þ where u is the fluid velocity along the y-axis, X and Z are the sizes of the conduit, while η ¼X/Z is the aspect ratio (η ¼1 in our case). From this expression, we can derive the pressure gradient G as a function of the mean velocity U, finding: "  #−1 12μU 192 1 kπ G¼− 1− η ∑ Tanh , ð9Þ 2η π 5 k odd k5 X2 The values of density and viscosity of the pure fluids were set equal to ρΑ ¼998 kg m−3 and μΑ ¼10–3 kg m−1 s−1 for water (fluid A) and ρΒ ¼ 789 kg m−3 and μΒ ¼ 1.2  10–3 kg m−1 s−1 for ethanol (fluid B), respectively, which correspond to their respective values at a 20 1C temperature. In addition, the density and the viscosity of the mixture were evaluated through Eqs. (6) and (7), where the volume of mixing and the fluidity of mixing are fitted from the experimental measurements by Dizechi and Marschall (1982), obtaining: Δvmix ðϕÞ ¼

1 a0 ϕð1−ϕÞ, ρ

ð10Þ

1 6 ∑ bn ϕn , μn¼0

ð11Þ

the cross section, while φ is its mean mass fraction within the cross section at position x, i.e., Z 1 1 N ϕðxÞ ¼ ϕðrÞdy dz ¼ ð13Þ ∑ ϕ: S S Ni¼1 i Clearly, since the velocity field is not uniform within a cross section, φ will vary along the channel. At steady state, it tends asymptotically (at large x, when the concentration field is uniform) to the ratio between the mass flow rate of fluid A to the total mass flow rate R _A ϕðrÞρðrÞ uðrÞ dy dz m ∑ϕi ρi ui , ð14Þ ¼ ϕb ¼ ¼ SR _ tot Nρu m S ρðrÞ uðrÞ dy dz where ui and ρi are the axial velocity and the density of the mixture at point i within the cross section, while ρ and u are the mean density and the mean velocity, with: Z Z 1 1 _ tot ¼ ρuS, ρðxÞ ¼ ρðrÞdy dz, uðxÞ ¼ ρðrÞuðrÞdy dz, m S S SρðxÞ S ð14aÞ Note that, as the volume of the mixture is not necessarily conserved, the total volumetric flow rate may change as we move along the channel. This is the “cup mixing average”, or “bulk”, mass fraction (Middleman, 1998), that is the mass fraction that we would measure when we collect the efflux from the channel in a cup and mix the contents to yield a uniform composition. In order to show the difference between the two averages, consider two immiscible fluid species, say A and B, having the same density and such that at the inlet they occupy the right and _ A SA and left side of a 2D Poiseuille channel flow. Denoting by m _ B SB the mass flow rate and the cross sectional area of fluid A and m _ tot ¼ m _ A þm _ B and S¼ SA þ SB), of fluid B, respectively (with m

and Δf mix ðϕÞ ¼

where a0 ¼0.145, while bn are appropriate fitting coefficients. In Fig. 2 both density and viscosity of a W–E mixture at 20 1C are represented as functions of the ethanol mass fraction. It is worth noting that the default setting of the Fluent 12.0 CFD code, i.e. a linear dependence of both ρ and μ on ϕ, appears to work reasonably well for ρ (this is obvious a casual coincidence), while it is completely wrong for μ. 2.4. Degree of mixing In most of the previous studies, the mixing quality was evaluated through the mean square deviation of the normalized mass fraction, which is defined at any cross section (in the plane yz) normal to the flow direction x as Z 1 1 N s2 ðxÞ ¼ ½ϕðrÞ−ϕðxÞ2 dy dz ¼ ð12Þ ∑ ðϕ −ϕÞ2 , S S Ni¼1 i where r¼ (x,y,z), N is the number of sampling points within the cross section of area S, ϕi is the mass fraction of fluid A at point i in

Fig. 4. Velocity contour plots at the outer cross section for W–E (top) and W–W (bottom) mixtures when Re¼ 10.

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

_ A =m _ B ¼ ðSA =SB Þ2 . Therefore, we see that, as the we easily find that m minority fluid species has a lower mean velocity, it occupies more space within the channel than its share based on mass flow rates, so that its mean mass fraction is less than its bulk mass fraction. Based on the above considerations, we will use a definition of mixing efficiency based on material fluxes, instead of concentrations, defining a volumetric flow variance as R ðϕðrÞ−ϕb Þ2 ρðrÞuðrÞdy dz 1 N R s2b ðxÞ ¼ S ¼ ð15Þ ∑ ðϕi −ϕb Þ2 ρi ui , Nρu ρðrÞuðrÞdy dz i¼1 S _ A ¼ ϕb ρu is the prescribed mass flow rate of species A per where m unit area. Therefore, we can define the degree of mixing as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s ð16Þ δm ¼ 1− b , smax ¼ ϕb ð1−ϕb Þ, smax where smax is the maximum value of the variance, that is achieved when the two streams are completely segregated. This definition is similar, although not identical, to the one used by Glasgow and Aubry (2003) and Galletti et al. (2012). Obviously, we expect that δm increases monotonically as we move down the mixing channel with x, tending asymptotically to 1 as the two fluids mix completely.

179

3. Results and discussion In this section we present the numerical simulation of the mixing process between two fluids in the T-type passive micromixer shown in Fig. 1, with Reynolds numbers ranging from 1 to 300. The W–W case, where the two inlet fluids are both water, is compared with W–E case, when they consist of water and ethanol, imposing the same volumetric flow rates at the two inlets, i.e. using the same mean inlet velocity U, and assuming that inlet velocity profiles are fully developed. In the same way, the Reynolds number is defined using the mean inlet velocity, U, the hydraulic diameter d of the mixing channel and the kinematic viscosity of water. First of all, it should be noted that, as the T-mixer is placed horizontally, in the W–W case a vertical interface region is stable and, in laminar flow with no mass diffusion, could remain unchanged as we move down the channel; on the other hand, in the W–E case this horizontally segregated morphology is unstable, due to the water larger density, and, in fact, as we move along the mixing channel, we see in Fig. 3 that the water stream will progressively move down and occupy the lower part of the pipe, while the ethanol stream will move up. A simple dimensional

Fig. 5. Ethanol mass fraction contour plots at Re¼ 100 (a) for a W–E system both close to the bottom wall and at the centerline of the mixing channel; at the outlet cross-section for (b) W–E, and (c) W–W systems. Streamlines of W–E and W–W systems are also shown in (d) and (e), respectively, seen from the upper end (see Fig. 1).

180

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

analysis reveals that the interface will start deforming from its initial vertical position when the vertical, transversal velocity u will balance the viscous velocity, μ⧸ρd. Therefore, considering that U≅uL/d from continuity, where U is the longitudinal velocity, we find: L/d≅Re¼ρUd/μ. Since in our case L/d¼ 2.5, that explains why, at the exit of the mixing channel, when Re¼0.1 the water stream has already occupied most of the lower part of the cross section, while when Re¼ 1 this process is just starting to occur, and when Re¼10 the W–E interface at the exit is still vertical (see the concentration contour plots in Fig. 3). The viscosity effects on the velocity field are well illustrated in Fig. 4, where the contour plots of the longitudinal velocity at the outer cross section of the mixing channel are shown for both W–E and W–W mixtures at the same Re¼10. As expected, we see that in the W–E case the velocity profile is smoother and the velocity at the interface is 15–20% smaller than in the W–W case. Then, at Re¼ 100, the heavier liquid, i.e. water, tends to occupy the center of the cross section, where the velocity is larger (see Fig. 5a and b), due to its larger inertia, as was clearly shown by many investigators (see for example, Joseph et al., 1984). As a comparison, in Fig. 5c we see that, in the W–W case at the same Re¼100, the contour plots exhibit mirror symmetry, although a secondary flow occurs in the form of a double vortex pair, due to the instabilities induced by the centrifugal forces at the confluence [see Kockmann et al. (2005)]; the corresponding streamlines are shown in Fig. 5d and e. Finally, in the engulfment regime occurring beyond a certain critical Reynolds number, fluid elements reach the opposite side of the mixing channel, thus largely increasing the degree of mixing. Using a micromixer having a 200  100 μm2 mixing channel (i.e. the same aspect ratio as ours), Bothe et al. (2006, 2008) found that in the W–W case engulfment occurs at Re ¼146; similar results were obtained by Hoffmann et al. (2006) and Engler et al. (2004), who indicated the engulfment to occur at 135o Reo 150 and Re¼147, respectively, and observed similar engulfment flow patterns even at larger Reynolds number. As our results are very similar, we only show a typical example, namely the engulfment contour plot at Re¼200 for the W–W case (see Fig. 6a).

On the other hand, in the W–E case, the onset of engulfment is retarded, occurring at larger Reynolds numbers compared to the W–W case. In fact, as shown in Fig. 6b, when Re¼ 200 the concentration plot is quite similar to that at Re¼100 (and therefore quite different than that of the W–W case at the same Re). Then, as we progressively increase the inlet flow rates, at Re¼ 230 we see the appearance of an engulfment flow (Fig. 6c), with the concentration fields remaining basically unchanged as we move up to Re¼300 (Fig. 6d). This behavior is well explained in Fig. 7, where the viscosity contour plots at the exit cross section are represented for different Re. Considering that the viscosity of a W–E mixture is quite larger than that of both components, we see that as Re increases from 0.1 to 10 the width of the mixing region narrows; then, as Re increases to 100 and 200, we see the more viscous fluid tending to occupy the center of the cross section until, at Re¼ 230, engulfment appears and the extent of the mixing region increases. From these viscosity plots we can conclude that when we try to mix water and ethanol, vortex formation is hampered by the more viscous mixture that emerges from the contact between the two streams at the confluence, thus retarding the onset of engulfment. This result agrees well with the experimental data (see Fig. 6b in Wang et al., 2012), showing that at Re¼200, in a T-mixer with W¼Win ¼300 μm and H ¼200 μm, the degree of mixing of the W–W system is very high, indicating the presence of an engulfment flow, while the W–E system exhibits a much lower degree of mixing, that is typical of vortex flows. Our results are summarized in Fig. 8, where the degree of mixing is plotted as a function of the Reynolds number for both W–W and W–E cases. As Re varies between 0.1 and 1, the degree of mixing is consistently higher in the W–E than in the W–W case, due to the longer residence time resulting from slower fluid velocity. Then, for 1 oReo10, the degree of mixing is about the same in the two cases, as the effect of the larger residence time for the W–E mixture is compensated by that of its smoother flow fields. Finally, at the onset of the engulfment flow, the degree of mixing greatly increases; that happens at Re≅140 in the W–W case and at Re≅230 in the W–E case.

Fig. 6. Ethanol mass fraction contour plots at the outlet cross-section: (a) for a W–W system at Re¼ 200; for a W–E system; (b) at Re¼200; (c) at Re¼230; and (d) at Re¼ 300.

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

181

Fig. 7. Viscosity contour plots for W–E systems when: (a) Re¼ 0.1; (b) Re¼1; (c) Re¼10; (d) Re¼ 100; (e) Re¼ 200; (f) Re¼230, and (g) Re¼ 300.

It should also be noted that, if we assume viscosity to be a linear function of the mass fraction, which is the default choice of most CFD programs, we would find that in the 1 oReo10 range the degree of mixing is overestimated by 10–20%. That could explain, at least qualitatively, why, using quite different geometries, Lin et al. (2007) and Chung and Shih (2008) found experimental values of the degree of mixing that were smaller than those predicted by simulations in which viscosity is assumed to depend linearly on mass fraction.

4. Conclusions Numerical simulations were performed to study the flow fields and mixing characteristics of liquid flows converging in a T-shaped

micromixer, when the two inlet fluids are either both water (W–W case) or water and ethanol (W–E case). We showed that, at smaller Reynolds numbers, Reo100, mixing is controlled by transverse diffusion, and therefore by the residence time of the fluid mixture occupying the interfacial region. Accordingly, as ethanol is slightly more viscous, and therefore slower, than water, the degree of mixing in the W–E case is slightly larger than that of the W–W case. On the other hand, at larger Reynolds numbers, mixing water and ethanol may take considerably longer, as the onset of engulfment is retarded and occurs at larger Reynolds numbers, namely it increases from Re≅140 in the W–W case to Re≅230 in the W–E case. In fact, at the confluence between the water and the ethanol streams, the W–E mixture is up to three times more viscous than water alone; this viscous layer separating the two streams will hamper any vortex formation, thus opposing mixing.

182

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

and concentration fields, can be expanded as ~ ¼ ∇x þ ε−1 ∇y , ∇

ϕðx,yÞ ¼ ∑ εn ϕn ðx,yÞ, n

~ uðx,yÞ ¼ ∑ εn un ðx,yÞ: n

ðA4Þ At leading order, we obtain ∇y ⋅u0 ¼ 0,

ðA5Þ

þ ~ y 2 u0 , ¼ 0, μ′ð∇ ~ y ϕ0 Þ⋅ð∇y u0 þ ∇y u0 Þ þ μ∇

ðA6Þ

u0 ⋅∇y ϕ0 ¼ 0,

ðA7Þ

The solution of this problem is u0 ¼ u0 ðxÞ and ϕ0 ¼ ϕ0 ðxÞ,

Fig. 8. Degree of mixing as a function of the Reynolds number for W–W and W–E systems.

Acknowledgment This work was supported in part by the Italian Ministry of Education and Research (MIUR), Grant PRIN, no. 2009-3JPM5Z.

ðA8Þ

indicating that the velocity and the concentration fields, at leading order, are both coarse-grained quantities. We see that when the viscosity is constant, i.e. μ~ ¼ 1, there is another possible solution, namely u0 ¼ u0 ðxÞ and ϕ0 ¼ ϕ0 ðx,yÞ, with u0 :∇y ϕ ¼ 0. This reflects the fact that, in this case, the Navier–Stokes Eqs. (A1) and (A2) are decoupled from the convection-diffusion equation (A3) and therefore it may occur that a solute diffuses in a direction that is the perpendicular to the fluid flow. When the viscosity of the fluid mixture depends on its concentration, though, this solution is no more possible. At the next order, Eqs. (A1)–(A3) give ∇y ⋅u1 þ ∇x ⋅u0 ¼ 0,

ðA9Þ

∇y P~ ¼ μ∇ ~ 2y u1 ,

ðA10Þ

Re u0 ⋅ð∇x φ0 þ ∇y φ1 Þ ¼ 0,

ðA11Þ

Appendix. Multiple scale analysis Consider a binary mixture whose viscosity, μ, depends on the composition, φ, of one of its two components, i.e. μ ¼ μ0 μðϕÞ, ~ where μðϕÞ ~ is an O(1), non-dimensional viscosity. Assume that the velocity field of the mixture varies over a viscosity-dependent distance, L ¼ μ0 =ρU, where U is a typical (e.g., mean) velocity and ρ is the density. For sake of simplicity, fluid motion due to changes of the fluid density, ρ, is neglected, which means assuming that the Archimedes number (i.e. the ratio between gravitational and viscous forces) is small, i.e., Ar ¼ ðΔρÞρgL4 =μ0 U⪡1, so that ρ can be assumed constant. With these assumptions, rewrite Eqs. (1)–(3) in terms of the non-dimensional quantities, r~ ¼ r~ =L and u~ ¼ u=U, obtaining,

where (A8) have been considered. Now, applying the solvability condition, that is taking the integral over all possible position in the y-space, Eqs. (A9) and (A11) become: ∇x ⋅u0 ¼ 0,

ðA12Þ

and Re u0 ⋅∇x ϕ0 ¼ 0,

ðA13Þ

while (A10) is identically satisfied. Finally, at the next order, after repeating a similar analysis, Eq. (A2) gives: þ Re u0 ⋅∇x u0 þ ∇x P~ ¼ μ′ð∇ ~ ~ 2x u0 : x ϕ0 Þ⋅ð∇x u0 þ ∇x u0 Þ þ μ∇

ðA14Þ

~ u~ ¼ 0, ∇⋅

ðA1Þ

Finally, when we go back to dimensional quantities, Eqs. (A12)– (A14) yield the following governing equations:

~ u~ þ ∇ ~ P~ ¼ μ′ð ~ ~ þ ∇u ~ þ Þ þ μ~ ∇ ~ 2 u, ~ ∇ ~ Re u⋅ ~ ∇ϕÞ⋅ð ∇u

ðA2Þ

∇⋅u ¼ 0,

ðA15Þ

ρu⋅∇u þ ∇P ¼ ∇⋅½μð∇u þ ∇u þ Þ,

ðA16Þ

u⋅∇ϕ ¼ 0:

ðA17Þ

~ ¼∇ ~ ϕ, ~ ∇ϕ Re Sc u⋅ 2

ðA3Þ

~ ¼ L∇ and P~ ¼ P=ðμ0 U=LÞ, P indicating the dynamic preswhere ∇ sure. Here, Re ¼ ρUd=μ0 and Sc ¼ μ0 =ρD denote, respectively, the Reynolds number and the Schmidt number, where D is the constant molecular diffusivity, while μ′ ~ ¼ dμ=dϕ ~ is an O(1) quantity, indicating the dependence of the fluid viscosity on composition. Now we apply the method of homogenization (Bensoussan et al., 1978; Mauri, 1995, 2003) to derive the coarse-grained effective transport equations in the limit when Sc ¼ 1=ε2 , with ε⪡1. We assume that each quantity depends separately on the macroscopic position vector, x, and on the stretched coordinate, ~ εy, so that ϕ ¼ ϕðx,yÞ and u~ ¼ uðx,yÞ. Note that, while the x-scale corresponds to the viscosity-dependent distance L, the y-scale corresponds to a microscopic distance l, with l/L≅Sc−1/2; in turbulent flow, L and l would be the Kolmogorov and the Batchelor scale. In addition, the gradient operator, together with the velocity

Therefore, we see that, at leading order, the diffusive term in the convection-diffusion equation (3) can be neglected, as it contributes only with O(Sc−1/2)-terms. The same conclusion is reached when we study the transport of colloidal particles in random velocity fields at large Pe: again, at leading order, the concentration field is constant within the microscale, so that molecular diffusivity does not play any role. Accordingly, the average solute transport is governed by an effective, eddy, diffusivity, Dn, that does not depend on the molecular diffusivity, D, and equals instead the time integral of the velocity autocorrelation function (Mauri (2003), and references therein). For example, when transverse convection is due to the presence of suspended spheres of radius R, as in fixed beds, then we find an effective diffusivity, Dn ¼UR/3 (Koch and Brady, 1985; Mauri, 1995), that does not depend on D. These results show

G. Orsi et al. / Chemical Engineering Science 95 (2013) 174–183

that, just as in the binary mixture case that we have studied here, if we follow a solute particle in a random velocity field neglecting its molecular diffusion altogether, we would capture all its average transport properties. References Ansari, M.A., Kim, K.Y., Kim, S.M., 2010. Numerical study of the effect on mixing of the position of fluid stream interfaces in a rectangular microchannel. Microsyst. Technol. 16, 1757–1763. Aubin, J., Fletcher, D.F., Xuereb, C., 2005. Design of micromixers using CFD modeling. Chem. Eng. Sci. 60, 2503–2516. Bensoussan, A., Lions, J.L., Papanicolaou, G., 1978. Asymptotic Analysis for Periodic Structures. North-Holland, Amsterdam. Bothe, D., Stemich, C., Warnecke, H.J., 2006. Fluid mixing in a T-shaped micromixer. Chem. Eng. Sci. 61, 2950–2958. Bothe, D., Stemich, C., Warnecke, H.J., 2008. Computation of scales and quality of mixing in a T-shaped microreactor. Comput. Chem. Eng. 32, 108–114. Bothe, D., Lojewski, A., Warnecke, H.J., 2011. Fully resolved numerical simulation of reactive mixing in a T-shaped micromixer using parabolized species equations. Chem. Eng. Sci. 66, 6424–6440. Chatwin, P.C., Sullivan, P.J., 1982. The effect of aspect ratio on the longitudinal diffusivity in rectangular channels. J. Fluid Mech. 120, 347–358. Christensen, C., Gmehling, J., Rassmussen, P., Weidlich, U., Holderbaum, T., 1984. Heats of Mixing Data Collection, vol. III., Part I Chemistry Data Series. Dechema, Frankfurt. Chung, C.K., Shih, T.R., 2008. Effect of geometry on fluid mixing of the rhombic micromixers. Microfluid. Nanofluid. 4, 419–425. Dizechi, M., Marschall, E., 1982. Viscosity of some binary and ternary liquid mixtures. J. Chem. Eng. Data 27, 358363. Engler, M., Kockmann, N., Kiefer, T., Woias, P., 2004. Numerical and experimental investigations on liquid mixing in static micromixers. Chem. Eng. J. 101, 315–322. Galletti, C., Roudgar, M., Brunazzi, E., Mauri, R., 2012. Effect of inlet conditions on the engulfment pattern in a T-shaped micro-mixer. Chem. Eng. J. 185–186, 300–313. Glasgow, I., Aubry, N., 2003. Enhancement of microfluidic mixing using time pulsing. Lab Chip 3, 114–120. Happel, J., Brenner, H., 1973. Low Reynolds Number Hydrodynamics. Noordhoff, Amsterdam. Hessel, V., Lowe, H., Schonfeld, F., 2005. Micromixers—a review on passive and active mixing principles. Chem. Eng. Sci. 60, 2479–2501. Hoffmann, M., Schluter, M., Rabiger, N., 2006. Experimental investigation of liquid– liquid mixing in T-shaped micro-mixers using mu-LIF and mu-PIV. Chem. Eng. Sci. 61, 2968–2976.

183

Hussong, J., Lindken, R., Pourquie, M., Westerweel, J., 2009. Numerical study on the flow physics of a T-shaped micro mixer. In: Ellero, M., et al. (Eds.), IUTAM Symposium on Advances in Micro- and Nano-fluidics. Springer, Heidelberg, pp. 191–205. Jeon, M.K., Kim, J.H., Noh, J., Kim, S.H., Park, H.G., Woo, S.I., 2005. Design and characterization of a passive recycle micromixer. J. Micromech. Microeng. 15, 346–350. Joseph, D.D., Nguyen, K., Beavers, G.S., 1984. Non-uniqueness and stability of the configuration of flow of immiscible fluids with different viscosities. J. Fluid Mech. 141, 319–345. Kern, D.Q., 1965. Process Heat Transfer. McGraw-Hill, New York, p. 161. Kim, D.S., Lee, S.W., Kwon, T.H., Lee, S.S., 2004. A barrier-embedded chaotic micromixer. J. Micromech. Microeng. 14, 798–805. Koch, D.L., Brady, J.F., 1985. Dispersion in fixed beds. J. Fluid Mech. 154, 399–427. Kockmann, N., Engler, M., Haller, D., Woias, P., 2005. Fluid dynamics and transfer processes in bended microchannels. Heat Transfer Eng. 26, 71–78. Laliberté, M., 2007. Model for calculating the viscosity of aqueous solutions. J. Chem. Eng. Data 52, 321–335. Lin, Y.C., Chung, Y.C., Wu, C.Y., 2007. Mixing enhancement of the passive microfluidic mixer with J-shaped baffles in the tee channel. Biomed. Microdevices 9, 215–221. Lin, Y., Yu, X., Wang, Z., Tu, S.T., Wang, Z., 2011. Design and evaluation of an easily fabricated micromixer with three-dimensional periodic perturbation. Chem. Eng. J. 171, 291–300. Mauri, R., 1995. Heat and mass transport in random velocity fields with application to dispersion in porous media. J. Eng. Math. 29, 77–89. Mauri, R., 2003. Heat and mass transport in nonhomogeneous random velocity fields. Phys. Rev. E 68, 066306. Middleman, S., 1998. An Introduction to Mass and Heat Transfer. Wiley, New York, pp. 217–218. Nguyen, N.T., Wu, Z.G., 2005. Micromixers—a review. J. Micromech. Microeng. 15, 1–16. Ottino, J.M., 1994. Mixing and chemical reactions: a tutorial. Chem. Eng. Sci. 49, 4005–4027. Perry, R.H., Green, D.W., 1997. Perry's Chemical Engineers' Handbook, 7th ed McGraw-Hill, New York (Tables 2–111). Roudgar, M., Brunazzi, E., Galletti, C., Mauri, R., 2012. Numerical study of split T-micromixers. Chem. Eng. Technol. 35, 1291–1299. Simmonds, C., 1919. Alcohol, Its Production, Properties, Chemistry, and Industrial Applications. Macmillan, New York. Wang, H., Iovenitti, P., Harvey, E., Masood, S., 2002. Optimizing layout of obstacles for enhanced mixing in microchannels. Smart Mater. Struct. 11, 662–667. Wang, P., Andenko, A., Young, R.D., 2004. Modeling viscosities of concentrated and mixed-solvent electrolyte systems. Fluid Phase Equilibr. 226, 71–82. Wang, W., Zhao, S., Shao, T., Jin, Y., Cheng, Y., 2012. Visualization of micro-scale mixing in miscible liquids using l-LIF technique and drug nano-particle preparation in T-shaped micro-channels. Chem. Eng. J. 192, 252–261.