Modelling Annular Micromixers
Modelling Annular Micromixers James P. Gleeson1 , Olivia M. Roche1 , Jonathan West2 , and Anne Gelb3 1 Applied Mathematics, University College Cork, Ireland. 2 National Microelectronics Research Centre, Lee Maltings, Cork, Ireland. 3 Department of Mathematics, Arizona State University, Tempe AZ 85287, U.S.A. Abstract Magnetohydrodynamic mixing of two fluids in an annular microchannel is modelled as a twodimensional, laminar, convection-diffusion problem, and examined using asymptotic analysis and numerical simulation. The time T required for mixing of a plug of solute depends on the Pecl´et number P e, and on the geometry of the annulus. Three scaling regimes are identified: purely diffusive, Taylordispersive, and convection-dominated; each with a characteristic power-law dependence of T upon P e. Consequences of these results for optimal micromixer design are discussed. Keywords: Laminar mixing, convection-diffusion, asymptotic analysis, microfluidics. AMS subject classifications: 76M45, 76R99.
1
-10
10
10
10
5
5
5
-5
5
10
-10
-5
5
10
-10
-5
5
-5
-5
-5
-10
-10
-10
10
Figure 1: Operation of an idealized micromixer at three times, neglecting diffusion.
1
Introduction
Recent advances in microfluidic and lab-on-a-chip technology have led to increased interest in laminar mixing of fluids [1, 2, 3] . Efficient mixing is vital for chemical reactions, but turbulence is absent at the low Reynolds numbers common in microscale devices, and molecular diffusion mixes on an unacceptably slow timescale. In this paper we discuss the mathematical modelling of an annular magnetohydrodynamic (MHD) micromixer, prototypes of which are under development at the Irish National Microelectronics Research Centre [4]. The device consists of an annular channel (see Figures 1 and 2), with inner and outer walls acting as electrodes, and an electromagnet underneath, which provides a vertical magnetic field. A radial electric field is imposed by applying a potential difference across the inner and outer electrodes, and the electric and magnetic fields produce an azimuthal Lorentz force which acts as a pumping mechanism for the fluid [5]. The idealized mixing action of this device is illustrated in Figure 1: in the absence of molecular diffusion, the initially separated fluids are convected through each other, increasing the interfacial length between them (linearly in time), and so promoting the mixing action of diffusion. In reality, the actions of convection and diffusion are felt simultaneously, and the goal of this paper is to examine their effect upon the efficiency of the mixer. The discussion here is limited to two dimensions, where the idealized limit of infinite depth has been taken. Previous studies of MHD pumping in an annulus have been motivated by liquid-metal flows and their stability [5, 6], but to our knowledge, this is the first investigation of the mixing effects in an annular geometry.
2
Notation and equations
The geometry of the annulus is shown in Figure 2. The radius of the center-line is R, and ρ represents the half-width of the channel, so the inner wall is located at r = R − ρ, and the outer wall at r = R + ρ. We describe the geometry using the nondimensional parameter γ=
ρ , R
(1)
which satisfies 0 < γ < 1. Note that the limit γ → 0 corresponds to a locally straight channel, and γ approaches 1 as the annulus becomes a punctured disk. The fluid velocity in the micromixer is found by solving the steady Navier-Stokes equations in the presence of the MHD body force [6]. In the two-dimensional case considered here, this reduces to an
2
R
r φ
ρ
Figure 2: Annular geometry, showing the center-line radius R, and the channel half-width ρ. ordinary differential equation for the azimuthal velocity v(r), all other velocity components being zero: d2 v 1 dv v α + − 2 =− , 2 dr r dr r r
(2)
where α represents the MHD forcing (specifically, α = −BI/4πhη, where I and B are the r.m.s. current and magnetic field strengths, h is the channel depth and η is the absolute viscosity of the fluid). The solution of (2) satisfying no-slip boundary conditions at the walls r = R − ρ and r = R + ρ is " µ ¶2 #−1 ω 1 (R2 − ρ2 )2 R−ρ v(r) = − ln × 8Rρr 4 16R2 ρ2 R+ρ · ¸ R−ρ r R+ρ × (R2 − ρ2 )2 ln + r2 (R − ρ)2 ln + r2 (R + ρ)2 ln , (3) R+ρ R−ρ r where the velocity is characterized by an average angular velocity ω, which is related to the MHD forcing parameter α by " µ ¶2 # 1 (R2 − ρ2 )2 R−ρ ω=α . (4) − ln 4 16R2 ρ2 R+ρ The profile (3) reduces to the parabolic Poiseuille profile for a straight channel when γ → 0: v(r) ≈
3ωR (r − R + ρ)(R + ρ − r) 2ρ2
(5)
with maximum at r = R:
3 ωR. (6) 2 The mixing of a plug of solute into a surrounding solvent is governed by the convection-diffusion equation (assuming both fluid phases have similar density, viscosity, etc.): vmax =
∂c + ∇ · (uc) − κ∇2 c = 0, ∂t
3
(7)
where u(x, t) is the fluid velocity vector, c(x, t) is the concentration of the solvent, and κ is the molecular diffusion coefficient. In the annular micromixer the convection-diffusion equation may be written in polar coordinates: µ ¶ ∂c v(r) ∂c κ ∂ ∂c κ ∂2c + − r − 2 2 = 0, (8) ∂t r ∂φ r ∂r ∂r r ∂φ with no-flux boundary conditions at the walls: ∂c (R − ρ, φ, t) ∂r ∂c (R + ρ, φ, t) ∂r
=
0
=
0.
(9)
In the following sections we examine solutions of (8) using analytical, asymptotic, and numerical methods. The dimensionless parameter used to compare the importance of convection and diffusion in (7) is the Pecl´et number, which we define for our system as Pe =
ωRρ . κ
(10)
Note that ωR is the characteristic linear velocity, while ρ is the smallest linear dimension in the system. Two natural groupings of parameters to give a dimensional time will be used in the sequel: the diffusion time R2 /κ, which measures the time for azimuthal mixing by diffusion in the absence of convection; and the convection time ω −1 , which is representative of the timescale of rotation of the fluid. The choice of appropriate time-scaling depends on whether diffusion or convection is dominant—we initially choose ω −1 , but will consider the alternative in section 6. A completely mixed solute-solvent system has the initial concentration c(x, 0) of solute spread evenly over the whole annulus. We define mixing efficiency using the timescale over which the concentration profile evolves to the uniform state. In order to measure the deviation from uniformity, we first introduce two averaging procedures, each operating on a function f (r, φ) defined on the annulus. Radial averaging is denoted by an overbar: Z R+ρ 1 f (r, φ)rdr, (11) f= 2Rρ R−ρ and angle brackets signify angular averaging: 1 2π
hf i =
Z
2π
f (r, φ)dφ.
(12)
0
Thus the average value of the concentration over the whole annulus is hci =
1 4πRρ
Z
R+ρ
Z
2π
c(r, φ, t)r dφ dr. R−ρ
(13)
0
In fact, it is straightforward to show from the convection-diffusion equation that hci is a constant—the total amount of solute is not changed over time, it is simply redistributed evenly over the annulus. In the following sections we adopt a simple initial concentration c(r, φ, 0) = 1 + cos φ,
4
(14)
for ease of asymptotic analysis. Later we show that the asymptotic results also apply to other initial conditions, for instance to the condition used in Figure 1: ½ 1 if 0 ≤ φ < π c(r, φ, 0) = . (15) 0 if π ≤ φ < 2π A mixing measure m(t) is a positive function of time characterizing the deviation of the concentration at time t from its uniformly mixed state hci. Define m(t) by D E 2 (c(r, φ, t) − hci) E, m(t) = D (16) 2 (c(r, φ, 0) − hci) so that m(0) = 1, and m(t) → 0 as t → ∞. The time TM for m(t) to decay from 1 to a specified value M is called the mixing time, and is defined by the condition m(TM ) = M.
(17)
From the point of view of experimentalists and design engineers, it is desirable to have simple formulas relating the mixing time TM to the Pecl´et number P e and the geometry ratio γ of the micromixer. We proceed to obtain asymptotic approximations to the solution of (8), and hence scaling laws for the mixing times TM .
3
Low Pecl´ et numbers
When the Pecl´et number (10) is sufficiently small, diffusive effects completely dominate convective motion. The mixing time by diffusion alone may be calculated by solving the diffusion equation µ ¶ ∂c κ ∂ ∂c κ ∂2c − r − 2 2 = 0, (18) ∂t r ∂r ∂r r ∂φ in the annulus, i.e., neglecting the convective term in (8). A series solution can be written as c = hci +
∞ X ∞ X
³ p ´ e−κλnj t Gn r λnj [αnj cos(nφ) + βnj sin(nφ)] ,
(19)
n=1 j=1
where the αnj and βnj are constants (determined from the initial condition) and λnm are eigenvalues, determined by the conditions dGn = 0 at r = R − ρ and r = R + ρ. (20) dr The eigenfunctions Gn are given in terms of Bessel functions as ³ ³ √ ´ √ ´ Jn−1 (R − ρ) λ − Jn+1 (R − ρ) λ √ √ √ ³ ³ (21) Gn (r λ) = Jn (r λ) − Yn (r λ) √ ´ √ ´ . Yn−1 (R − ρ) λ − Yn+1 (R − ρ) λ For the single-mode initial condition (14), only the n = 1 term of the double sum is present in (19). Consequently, the mixing measure (16) in this case is m(t) =
∞ X
e
¡ p ¢2 r λ1j ¡ p ¢2 . G1 r λ1j
−2κλ1j t G1
j=1
5
(22)
The form of m(t) at large times t is dominated by the first term (j = 1) in this sum. Indeed, in the limit γ → 0, we have m(t) ∼ exp(−2κλ11 t) as t → ∞, (23) with the eigenvalue given by λ11
· ¸ 1 1 2 4 = 2 1 + γ + O(γ ) for γ ¿ 1. R 3
The asymptotic mixing time TM as defined in (17) then follows from (23) and (24): µ ¶ 1 1 TM ∼ ln as M → 0, 2κλ11 M
(24)
(25)
and so the nondimensional mixing time is ωTM ≈
4
Pe ln 2γ
µ
1 M
¶· ¸ 1 1 − γ2 + . . . 3
for γ ¿ 1.
(26)
Intermediate Pecl´ et numbers
In this section we follow and adapt Taylor’s [7], [8] arguments to derive an effective dispersion under certain conditions on the Pecl´et number, and clearly identify the limits of validity. We note that Nunge et al. [9] have derived a Taylor dispersion coefficient in a curved channel (matching our equation (35)), but they do not investigate the values of P e where the analysis is valid. These bounding values of P e assume great importance in our later analysis, so we follow Taylor’s original formulation to understand when the approximations break down. Consider a reference frame rotating with angular velocity ω, i.e., with azimuthal angle measured by φ0 = φ − ωt. The convection-diffusion equation (8) in this rotating frame is µ ¶ µ ¶ ∂c v(r) ∂c κ ∂ ∂c κ ∂2c + −ω − r − 2 0 2 = 0, 0 ∂t r ∂φ r ∂r ∂r r ∂φ
(27)
(28)
with boundary conditions as before. We will work in this rotating frame for the remainder of this section, and so drop the prime on φ. Following Taylor, we assume the concentration is quasi-steady in the rotating ∂c frame, with angular gradient ∂φ independent of r. Moreover, the effect of azimuthal diffusion is neglected compared to radial diffusion. The resulting equation is readily integrated twice to yield Z Z 1 ∂c r 1 r2 c = c(R − ρ, φ) + [v(r1 ) − ωr1 ] dr1 dr2 . (29) κ ∂φ R−ρ r2 R−ρ We take the radial mean (as defined in equation (11)) to allow us to eliminate c(R − ρ, φ) and obtain c=c+ having used the notation
Z
r
g(r) = R−ρ
1 r2
1 ∂c [g(r) − g] , κ ∂φ
Z
(30)
r2
[v(r1 ) − ωr1 ] dr1 dr2 R−ρ
6
(31)
for brevity. The average flux of solvent through the line φ =constant (in the moving frame) is given by the radial average of the product of the concentration and angular velocity ω 0 (also relative to the moving frame): J
= = =
c ω0 ¸h Z R+ρ · i 1 v 1 ∂c c+ (g(r) − g) − ω r dr 2Rρ R−ρ κ ∂φ r Z R+ρ 1 ∂c 1 (v(r) − ωr) g(r)dr. κ ∂φ 2Rρ R−ρ
If we assume
∂c ∂c ≈ , ∂φ ∂φ
(32)
(33)
this implies that the radial mean concentration in the moving frame is governed by a one-dimensional diffusion equation ∂c ∂2c = D 2, (34) ∂t ∂φ where D is the Taylor dispersion coefficient for the annulus, defined by Z R+ρ 1 1 D = (v(r) − ωr) g(r) dr κ 2Rρ R−ρ Z R+ρ Z Z r 1 1 1 r2 = (v(r) − ωr) [v(r1 ) − ωr1 ] dr1 dr2 dr. κ 2Rρ R−ρ R−ρ r2 R−ρ This integral may be calculated for the velocity (3); after some algebra we find · µ ¶¸−2 ω 2 R2 1−γ D = 4γ 2 − (1 − γ 2 )2 ln × 24κγ 1+γ " µ ¶2 1+γ 1+γ 5 2 4 2 2 3 2 2 2 × 48γ (1 + γ ) + 120γ (1 − γ ) ln − 72γ (1 − γ ) (1 + γ ) ln 1−γ 1−γ µ ¶4 µ ¶5 # 1+γ 1+γ 2 4 2 2 4 2 4 −6γ(1 − γ ) (1 + γ ) ln + (1 − γ ) (3 + 10γ + 3γ ) ln . 1−γ 1−γ
(35)
(36)
In the γ ¿ 1 limit, this reduces to · D ≈ ωP e
¸ 2 346 3 γ+ γ ... . 105 1576
(37)
The solution of the one-dimensional diffusion equation (34) with initial condition from (14) is straightforward. Having found c, the concentration c and so the mixing measure m(t) are calculated from (30) and (16) respectively. The mixing measure has an exponential tail: m(t) ∼ exp(−2Dt) and so the mixing time TM is TM
1 ln ∼ 2D
µ
1 M
7
as t → ∞,
(38)
as M → 0
(39)
¶
in the Taylor dispersion regime. From (37), this yields the (nondimensional) asymptotic form µ ¶ · ¸ 1 1 105 18165 2 ωTM ≈ ln 1− γ + . . . for γ ¿ 1, γP e M 4 1576
(40)
when the assumptions made above are valid. The range of P e where the above approximate analysis holds will prove to be of great interest in later sections, so we now examine carefully the two basic assumptions that were made, and cast these into conditions on P e for the result (40) to hold. The first condition stems from the requirement that the azimuthal diffusion be negligible compared to ∂2c the radial diffusion; effectively this requires that the coefficient of ∂φ 2 in (28), actually its radial average, should be much less than the dispersion coefficient (35): κ ¿ Dr2 .
(41)
Consider this condition when γ ¿ 1, using the first term of (37): κ
¿
⇐⇒ P e2
À
2 γωP eR2 105 105 . 2
(42)
The resulting condition P e À 7.2 is analogous to that in a straight capillary, see for instance [8]. ∂c ∂c The second important approximation is the replacement of ∂φ by ∂φ in (33). This is valid if the ∂c coefficient of ∂φ in (30) is very small. The order of magnitude of this coefficient is g/κ, so the condition is Z R+ρ Z r Z 1 1 1 r2 r [v(r1 ) − ωr1 ] dr1 dr2 dr ¿ 1. (43) κ 2Rρ R−ρ R−ρ r2 R−ρ Again, this yields a simple condition if γ ¿ 1: 1 ωρ2 15 κ
¿
1
⇐⇒ P e
¿
15 . γ
(44)
Thus (42) and (44) provide rough bounds on P e for the lower and upper ranges of validity of the Taylor dispersion approximations.
5
Convection-dominated mixing: P e À 1
In this section we examine the solution of equation (8) when the molecular diffusion term is small in comparison to the convective term, i.e., the limit of large Pecl´et number. As the diffusion κ multiplies the terms with highest derivatives in (8), the limit κ → 0 yields a singularly-perturbed problem. Regular perturbation methods fail to give approximations which are uniformly valid in time for such problems, and so we introduce multiple time scales. Taking κ as a small parameter, the time scales are: t0 t1 t2 t3
= = = =
8
t κt3 κt2 κt.
(45)
As is standard in multiple scale expansions, these time scales are assumed to be effectively independent of each other, and so the time derivative in (8) may be replaced by ∂ ∂ ∂ ∂ ∂ = + 3t20 κ + 2t0 κ +κ , ∂t ∂t0 ∂t1 ∂t2 ∂t3 accurate to first order in κ. The corresponding expansion of the concentration is ¡ ¢ c = c0 + κc1 + O κ2 ,
(46)
(47)
where c0 and c1 are functions of all the time scales, as well as r and φ. Inserting the expansions (46) and (47) into equation (8) leads to the following equation for κ-independent terms: ∂c0 v(r) ∂c0 + = 0, (48) ∂t0 r ∂φ with solution
³ ³ v ´ v ´ c0 = 1 + A(t1 , t2 , t3 ) cos φ − t0 + B(t1 , t2 , t3 ) sin φ − t0 . r r Here A and B are functions satisfying the initial conditions A(0, 0, 0)
=
1
B(0, 0, 0) =
0.
In order to further specify A and B, it is necessary to consider the equation for terms of order κ: µ ¶ ∂c1 v(r) ∂c1 ∂c0 ∂c0 1 ∂ ∂c0 1 ∂ 2 c0 2 ∂c0 + = −3t0 − 2t0 − + r + 2 . ∂t0 r ∂φ ∂t1 ∂t2 ∂t3 r ∂r ∂r r ∂φ2
(49)
(50)
(51)
Note that the differential operator on the left hand side of this equation is the same as in (48), while the forcing terms on the right hand side all depend on c0 . Calculating the forcing terms using (49), and requiring that the solution c1 should not grow without bound as t0 → ∞ leads to the following equations for A and B: " # 2 ∂A 1 v2 2vv 0 v0 = − − 3 + 2 A ∂t1 3 r4 r r " # 2 ∂B 1 v2 2vv 0 v0 = − − 3 + 2 B ∂t1 3 r4 r r · 00 ¸ 0 ∂A 1 v v v = − − 2+ 3 B ∂t2 2 r r r · ¸ ∂B 1 v 00 v0 v = − 2+ 3 A ∂t2 2 r r r ∂A 1 = − 2A ∂t3 r ∂B 1 = − 2 B. (52) ∂t3 r
9
Here primes are used to denote differentiation with respect to r. The solutions of these equations satisfying the initial conditions (50) are: " Ã ! # · µ ¶ ¸ 2 1 v2 2vv 0 v0 1 v 00 v0 v 1 A = exp − − + t − t cos − − + t2 1 3 3 r4 r3 r2 r2 3 r r2 r3 " Ã ! # · µ ¶ ¸ 2 1 v2 2vv 0 v0 1 1 v 00 v0 v B = exp − (53) − 3 + 2 t1 − 2 t3 sin − − 2 + 3 t2 3 r4 r r r 3 r r r In the limit γ → 0, we approximate the velocity profile as in (5), and retain only the dominant terms in (53) to obtain " # µ ¶2 ¸ · 3 r−R 3 3 3 2 2 A = exp − ω t cos ω t γP e ρ 2γP e " # µ ¶2 · ¸ 3 3 r−R 3 3 B = exp − ω t sin ω 2 t2 . (54) γP e ρ 2γP e We have thus found a non-secular approximation to the concentration c for large P e and small γ. By combining the results (49) and (54), the corresponding mixing measure is calculated as # " µ ¶2 Z R+ρ 1 6 r−R 3 3 m(t) = exp − ω t r dr 2Rρ R−ρ γP e ρ µ 3 3¶ √ π 6ω t = F , (55) 2 γP e where F is defined as the monotonic function F (x) = x−1/2 erf(x1/2 ).
(56)
The nondimensional mixing time corresponding to a mixing measure value of M is therefore · ωTM ≈
γP e −1 F 6
µ
2M √ π
¶¸1/3 ,
(57)
and note in particular that this increases as P e1/3 when M and γ are fixed.
6
Numerical simulations
To check the asymptotic results derived in previous sections, and to extend to cases with initial conditions other than the simple form (14), we solve the convection-diffusion equation (8) numerically. A decomposition into a finite number N of Fourier modes in the azimuthal angle is employed: c(r, φ, t) =
N X 1 g0 (r, t) + gn (r, t) cos(nφ) + hn (r, t) sin(nφ), 2 n=1
(58)
and substitution into (8) yields partial differential equations for each gn and hn . Having solved the system of equations for the Fourier coefficients, the full concentration field may be constructed from (58), or the
10
0
0
10
10
Pe=8
−1
−1
10
10
Pe=46341 m(t)
m(t) Pe=2 −2
−2
10
10
Pe=16384
Pe=64
Pe=0.5
Pe=0.5
Pe=32
−3
10
−3
0
10
20
30
40
50
ωt
60
70
80
90
100
10
0
10
20
30
40
50
ωt
60
70
80
90
100
Figure 3: Mixing measure m(t) as a function of nondimensional time, calculated in numerical simulations with γ = 0.05 and various Pecl´et numbers as shown. Asymptotic predictions are also shown for: P e = 0.5 (dashed) using (23); P e = 64 (dotted) using (38); and P e = 16384 (dot-dash) using (55). mixing measures may be evaluated more directly from (16) by integrating over the angle and normalizing m(0) to unity: N ³g ´2 1 X 0 m(t) = − hci + g 2 + h2n . (59) 2 2 n=1 n Logarithmic plots of m(t) as a function of time at various P e values are shown in Figure 3 for γ = 0.05. For comparison we plot also the asymptotic forms of m(t), using (23) in the diffusive regime, (38) in the Taylor regime, and (55) in the convective regime. The asymptotic formulas fit the numerical results well, even at early times. Note that the straight-line asymptotes at very large times for the high-P e case are not reproduced by the multiple-time solution—this is apparently due to the neglect of the boundary conditions (9) in our analysis, and the subsequent importance of boundary layers at the inner and outer walls. From the numerical solution for m(t), it is straightforward to calculate the mixing times TM required for the measure to decay from its initial value of 1 to the value M . We choose three values of M for comparison with the asymptotic predictions: M = 0.3, 0.1, and 0.01, and investigate a wide range of Pecl´et numbers. For fixed values of γ and M , the asymptotic analysis of the previous sections predicts a mixing time proportional to P e in the diffusion-dominated regime, to P e−1 in the Taylor dispersion regime, and to P e1/3 when convection dominates. The numerical values of nondimensional times ωTM are plotted in Figure 4 along with the straight lines corresponding to the formulas (26), (40), and (57). Note the excellent agreement with predictions, except for the lowest value of M , when the multiple scale solution no longer accurately describes the time evolution of m(t). As much of our analysis has concentrated on the γ → 0 limit, it is worthwhile comparing numerics and asymptotics for a larger value of γ, noting that the maximum possible value of γ is 1. In Figure 5 we plot the results for γ = 0.2, and find that in general the correspondence between predicted and numerical values is excellent, but note that the limits on the Taylor regime now preclude the formation of a clear P e−1 scaling range. Another simplification in the analysis concerns the initial condition (14), which consists of only a single harmonic of the azimuthal variable. In order to examine the robustness of our predictions, we numerically
11
2
10
ωT
1
10
0
1
10
2
10
3
10
4
10
10
5
10
Pe
Figure 4: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.05. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles) 2
10
1
10 ωT
0
10
0
10
1
2
10
10
3
10
4
10
Pe
Figure 5: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.2. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles)
12
2
10
ωT
1
10
0
10
1
2
10
3
10
10
4
10
5
10
Pe
Figure 6: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.05, and an initial concentration given by (15) (slightly smoothed). Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles) solve the convection-diffusion equation for initial condition (15), replacing the discontinuous function by hyperbolic tangents to avoid Gibbs oscillations. The decay of the mixing measure is still dominated by the first angular harmonic, and so the mixing times are remarkably close to the asymptotic estimates, see Figure 6. We therefore expect the formulas to be useful for rather general initial distributions of the solute. The nondimensional time ωTM used in Figures 3 to 6 may be replaced by the alternative nondimensionalization mentioned in Section 2, i.e., the diffusion time κTM /R2 . The timescales are related by κTM γ = ωTM , (60) 2 R Pe and so the data of, for example, Figure 4 is easily recast in terms of the diffusion time, see Figure 7. The relevant mixing times in the diffusive, Taylor, and convective regimes are thus found by multiplying (26), (40), and (57) by γ/P e, yielding µ ¶· ¸ κTM 1 1 1 2 ≈ ln 1 − γ + ... , (61) R2 2 M 3 for P e ¿ 7.2, κTM 1 ≈ ln 2 R P e2 for 7.2 ¿ P e ¿ 15/γ, and
µ
1 M
¶
· ¸ 105 18165 2 1− γ + ... 4 1576
· µ ¶¸1/3 2M κTM − 23 43 1 −1 √ ≈ Pe γ F R2 6 π
13
(62)
(63)
0
10
−1
10
−2
10 κ T/R2
−3
10
−4
10
−5
10
0
10
1
10
2
3
10
10
4
10
5
10
Pe
Figure 7: Mixing times nondimensionalized by the diffusion time, κT /R2 , as a function of Pecl´et number, for γ = 0.05. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles)
14
for P e À 15/γ. Figure 7 is especially of interest to experimentalists working with a particular solute and solvent (so that κ is fixed) while varying the rate of rotation velocity ω of the micromixer to change P e. The mixing time is seen to decrease from the diffusion time through the Taylor regime (at a rate proportional to P e−2 ), and then continue to decrease at a slower rate beyond the Taylor regime. The slower rate corresponds to a P e−2/3 scaling when (55) is valid, i.e., for M ≥ 0.1, but is closer to P e−1/2 for very small values of M . This is again a consequence of the failure of the multiple-scale method to account for the very long time scales where boundary layer effects become non-negligible. From our analysis of the limits of validity of the Taylor dispersion description, see (44), we can estimate the end of the Taylor regime as 15γ −1 for small γ. Of interest to the micromixer designer is the influence of the geometry ratio γ: as the mixing time decreases as P e−2 until P e ≈ 15γ −1 , it seems advisable to decrease γ as much as possible in order to achieve faster mixing at lower velocities. However, this approach is overly simplistic, since the Pecl´et number also depends on γ—the linear velocity in the channel decreases with decreasing ρ, all other things being equal. A problem a great relevance to experimentalists is to find the optimal annular geometry to give lowest possible mixing time for a given species, and given values of the electric and magnetic pumping fields (so κ and α are fixed). Suppose the “footprint” R of the device is chosen, then we must find the channel half-width ρ to minimize the mixing time. Since P e may be shown to increase monotonically with ρ, we conclude from Figure 7 that the mixing time decreases as the channel width increases. Practical considerations such as the minimal inner radius for effective electrodes will also provide an upper bound on the acceptable geometry ratio γ.
7
Summary
Laminar mixing in a two-dimensional annulus has been examined by asymptotic analysis and numerical simulation of the convection-diffusion equation (8). The mixing measure defined in (16) has been shown to decay exponentially in the diffusion and Taylor regimes (see equations (23) and (38)), with time constants given in terms of the Pecl´et number and the ratio γ characterizing the annulus geometry. Corresponding mixing times TM are predicted to scale as P e0 and P e−2 relative to the azimuthal diffusion time, see equations (61) and (62). Figures 3 to 7 demonstrate the robustness of these predictions by comparison to numerical simulations. In the convection-dominated regime, multiple scale analysis yields the first approximation (55) to the mixing measure, with mixing time (63) showing a P e−2/3 dependence on the Pecl´et number. This is again found to agree well with numerical results, at least for M ≥ 0.1. Although the asymptotic formulas were derived in the limit γ → 0, we find they yield good predictions for larger values of γ also (Figure 5). An important limitation of the present analysis is the assumption that all motion is two-dimensional; it is known [10, 11] that the dispersion in a channel of finite depth may not equal the value calculated by assuming the depth to be infinite. Work incorporating three-dimensional effects may yield further insight into this important problem.
8
Acknowledgements
This work has been supported by funding from the Enterprise Ireland International Collaboration Fund, the National Microelectronics Research Centre, and by the Faculty of Arts Research Fund, University College Cork.
15
References [1] J. P. Gleeson and J. West, Magnetohydrodynamic Micromixing, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 318-321. [2] I. Meisel and P. Ehrhard, Simulation of electrically excited flows in microchannels for mixing application, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 62-65. [3] K. Mohseni, Mixing and impulse extremization in microscale vortex formation, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 392-395. [4] J. West et al., Application of magnetohydrodynamic actuation to continuous flow chemistry, Lab on a Chip, 2 (2002), pp. 224-230. [5] P. Tabeling and J. P. Chabrerie, Magnetohydrodynamic Taylor vortex flow under a transverse pressure gradient, Phys. Fluids, 24 (1981), pp. 406-412. [6] M-H. Chang and C-K. Chen, Hydromagnetic stability of current-induced flow in a small gap between concentric cylinders, J. Fluids Eng., 121 (1999), pp. 548-554. [7] G. I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proc. Roy. Soc. Lond. A, 219 (1953), pp. 186-203. [8] R. F. Probstein, Physiochemical Hydrodynamics, Wiley, New York, 1994. [9] R. J. Nunge, T.-S. Lin and W. N. Gill, Laminar dispersion in curved tubes and channels, J. Fluid Mech., 51 (1972), pp. 363-383. [10] H. Brenner and D. A. Edwards, Macrotransport Processes, Butterworth-Heinemann, 1993. [11] M. R. Doshi, P. M. Daiya and W. N. Gill, Three dimensional laminar dispersion in open and closed rectangular conduits, Chem. Eng. Sci., 33 (1978), pp. 795-804.
16
1
-10
10
10
10
5
5
5
-5
5
10
-10
-5
5
10
-10
-5
5
-5
-5
-5
-10
-10
-10
10
Figure 1: Operation of an idealized micromixer at three times, neglecting diffusion.
1
Introduction
Recent advances in microfluidic and lab-on-a-chip technology have led to increased interest in laminar mixing of fluids [1, 2, 3] . Efficient mixing is vital for chemical reactions, but turbulence is absent at the low Reynolds numbers common in microscale devices, and molecular diffusion mixes on an unacceptably slow timescale. In this paper we discuss the mathematical modelling of an annular magnetohydrodynamic (MHD) micromixer, prototypes of which are under development at the Irish National Microelectronics Research Centre [4]. The device consists of an annular channel (see Figures 1 and 2), with inner and outer walls acting as electrodes, and an electromagnet underneath, which provides a vertical magnetic field. A radial electric field is imposed by applying a potential difference across the inner and outer electrodes, and the electric and magnetic fields produce an azimuthal Lorentz force which acts as a pumping mechanism for the fluid [5]. The idealized mixing action of this device is illustrated in Figure 1: in the absence of molecular diffusion, the initially separated fluids are convected through each other, increasing the interfacial length between them (linearly in time), and so promoting the mixing action of diffusion. In reality, the actions of convection and diffusion are felt simultaneously, and the goal of this paper is to examine their effect upon the efficiency of the mixer. The discussion here is limited to two dimensions, where the idealized limit of infinite depth has been taken. Previous studies of MHD pumping in an annulus have been motivated by liquid-metal flows and their stability [5, 6], but to our knowledge, this is the first investigation of the mixing effects in an annular geometry.
2
Notation and equations
The geometry of the annulus is shown in Figure 2. The radius of the center-line is R, and ρ represents the half-width of the channel, so the inner wall is located at r = R − ρ, and the outer wall at r = R + ρ. We describe the geometry using the nondimensional parameter γ=
ρ , R
(1)
which satisfies 0 < γ < 1. Note that the limit γ → 0 corresponds to a locally straight channel, and γ approaches 1 as the annulus becomes a punctured disk. The fluid velocity in the micromixer is found by solving the steady Navier-Stokes equations in the presence of the MHD body force [6]. In the two-dimensional case considered here, this reduces to an
2
R
r φ
ρ
Figure 2: Annular geometry, showing the center-line radius R, and the channel half-width ρ. ordinary differential equation for the azimuthal velocity v(r), all other velocity components being zero: d2 v 1 dv v α + − 2 =− , 2 dr r dr r r
(2)
where α represents the MHD forcing (specifically, α = −BI/4πhη, where I and B are the r.m.s. current and magnetic field strengths, h is the channel depth and η is the absolute viscosity of the fluid). The solution of (2) satisfying no-slip boundary conditions at the walls r = R − ρ and r = R + ρ is " µ ¶2 #−1 ω 1 (R2 − ρ2 )2 R−ρ v(r) = − ln × 8Rρr 4 16R2 ρ2 R+ρ · ¸ R−ρ r R+ρ × (R2 − ρ2 )2 ln + r2 (R − ρ)2 ln + r2 (R + ρ)2 ln , (3) R+ρ R−ρ r where the velocity is characterized by an average angular velocity ω, which is related to the MHD forcing parameter α by " µ ¶2 # 1 (R2 − ρ2 )2 R−ρ ω=α . (4) − ln 4 16R2 ρ2 R+ρ The profile (3) reduces to the parabolic Poiseuille profile for a straight channel when γ → 0: v(r) ≈
3ωR (r − R + ρ)(R + ρ − r) 2ρ2
(5)
with maximum at r = R:
3 ωR. (6) 2 The mixing of a plug of solute into a surrounding solvent is governed by the convection-diffusion equation (assuming both fluid phases have similar density, viscosity, etc.): vmax =
∂c + ∇ · (uc) − κ∇2 c = 0, ∂t
3
(7)
where u(x, t) is the fluid velocity vector, c(x, t) is the concentration of the solvent, and κ is the molecular diffusion coefficient. In the annular micromixer the convection-diffusion equation may be written in polar coordinates: µ ¶ ∂c v(r) ∂c κ ∂ ∂c κ ∂2c + − r − 2 2 = 0, (8) ∂t r ∂φ r ∂r ∂r r ∂φ with no-flux boundary conditions at the walls: ∂c (R − ρ, φ, t) ∂r ∂c (R + ρ, φ, t) ∂r
=
0
=
0.
(9)
In the following sections we examine solutions of (8) using analytical, asymptotic, and numerical methods. The dimensionless parameter used to compare the importance of convection and diffusion in (7) is the Pecl´et number, which we define for our system as Pe =
ωRρ . κ
(10)
Note that ωR is the characteristic linear velocity, while ρ is the smallest linear dimension in the system. Two natural groupings of parameters to give a dimensional time will be used in the sequel: the diffusion time R2 /κ, which measures the time for azimuthal mixing by diffusion in the absence of convection; and the convection time ω −1 , which is representative of the timescale of rotation of the fluid. The choice of appropriate time-scaling depends on whether diffusion or convection is dominant—we initially choose ω −1 , but will consider the alternative in section 6. A completely mixed solute-solvent system has the initial concentration c(x, 0) of solute spread evenly over the whole annulus. We define mixing efficiency using the timescale over which the concentration profile evolves to the uniform state. In order to measure the deviation from uniformity, we first introduce two averaging procedures, each operating on a function f (r, φ) defined on the annulus. Radial averaging is denoted by an overbar: Z R+ρ 1 f (r, φ)rdr, (11) f= 2Rρ R−ρ and angle brackets signify angular averaging: 1 2π
hf i =
Z
2π
f (r, φ)dφ.
(12)
0
Thus the average value of the concentration over the whole annulus is hci =
1 4πRρ
Z
R+ρ
Z
2π
c(r, φ, t)r dφ dr. R−ρ
(13)
0
In fact, it is straightforward to show from the convection-diffusion equation that hci is a constant—the total amount of solute is not changed over time, it is simply redistributed evenly over the annulus. In the following sections we adopt a simple initial concentration c(r, φ, 0) = 1 + cos φ,
4
(14)
for ease of asymptotic analysis. Later we show that the asymptotic results also apply to other initial conditions, for instance to the condition used in Figure 1: ½ 1 if 0 ≤ φ < π c(r, φ, 0) = . (15) 0 if π ≤ φ < 2π A mixing measure m(t) is a positive function of time characterizing the deviation of the concentration at time t from its uniformly mixed state hci. Define m(t) by D E 2 (c(r, φ, t) − hci) E, m(t) = D (16) 2 (c(r, φ, 0) − hci) so that m(0) = 1, and m(t) → 0 as t → ∞. The time TM for m(t) to decay from 1 to a specified value M is called the mixing time, and is defined by the condition m(TM ) = M.
(17)
From the point of view of experimentalists and design engineers, it is desirable to have simple formulas relating the mixing time TM to the Pecl´et number P e and the geometry ratio γ of the micromixer. We proceed to obtain asymptotic approximations to the solution of (8), and hence scaling laws for the mixing times TM .
3
Low Pecl´ et numbers
When the Pecl´et number (10) is sufficiently small, diffusive effects completely dominate convective motion. The mixing time by diffusion alone may be calculated by solving the diffusion equation µ ¶ ∂c κ ∂ ∂c κ ∂2c − r − 2 2 = 0, (18) ∂t r ∂r ∂r r ∂φ in the annulus, i.e., neglecting the convective term in (8). A series solution can be written as c = hci +
∞ X ∞ X
³ p ´ e−κλnj t Gn r λnj [αnj cos(nφ) + βnj sin(nφ)] ,
(19)
n=1 j=1
where the αnj and βnj are constants (determined from the initial condition) and λnm are eigenvalues, determined by the conditions dGn = 0 at r = R − ρ and r = R + ρ. (20) dr The eigenfunctions Gn are given in terms of Bessel functions as ³ ³ √ ´ √ ´ Jn−1 (R − ρ) λ − Jn+1 (R − ρ) λ √ √ √ ³ ³ (21) Gn (r λ) = Jn (r λ) − Yn (r λ) √ ´ √ ´ . Yn−1 (R − ρ) λ − Yn+1 (R − ρ) λ For the single-mode initial condition (14), only the n = 1 term of the double sum is present in (19). Consequently, the mixing measure (16) in this case is m(t) =
∞ X
e
¡ p ¢2 r λ1j ¡ p ¢2 . G1 r λ1j
−2κλ1j t G1
j=1
5
(22)
The form of m(t) at large times t is dominated by the first term (j = 1) in this sum. Indeed, in the limit γ → 0, we have m(t) ∼ exp(−2κλ11 t) as t → ∞, (23) with the eigenvalue given by λ11
· ¸ 1 1 2 4 = 2 1 + γ + O(γ ) for γ ¿ 1. R 3
The asymptotic mixing time TM as defined in (17) then follows from (23) and (24): µ ¶ 1 1 TM ∼ ln as M → 0, 2κλ11 M
(24)
(25)
and so the nondimensional mixing time is ωTM ≈
4
Pe ln 2γ
µ
1 M
¶· ¸ 1 1 − γ2 + . . . 3
for γ ¿ 1.
(26)
Intermediate Pecl´ et numbers
In this section we follow and adapt Taylor’s [7], [8] arguments to derive an effective dispersion under certain conditions on the Pecl´et number, and clearly identify the limits of validity. We note that Nunge et al. [9] have derived a Taylor dispersion coefficient in a curved channel (matching our equation (35)), but they do not investigate the values of P e where the analysis is valid. These bounding values of P e assume great importance in our later analysis, so we follow Taylor’s original formulation to understand when the approximations break down. Consider a reference frame rotating with angular velocity ω, i.e., with azimuthal angle measured by φ0 = φ − ωt. The convection-diffusion equation (8) in this rotating frame is µ ¶ µ ¶ ∂c v(r) ∂c κ ∂ ∂c κ ∂2c + −ω − r − 2 0 2 = 0, 0 ∂t r ∂φ r ∂r ∂r r ∂φ
(27)
(28)
with boundary conditions as before. We will work in this rotating frame for the remainder of this section, and so drop the prime on φ. Following Taylor, we assume the concentration is quasi-steady in the rotating ∂c frame, with angular gradient ∂φ independent of r. Moreover, the effect of azimuthal diffusion is neglected compared to radial diffusion. The resulting equation is readily integrated twice to yield Z Z 1 ∂c r 1 r2 c = c(R − ρ, φ) + [v(r1 ) − ωr1 ] dr1 dr2 . (29) κ ∂φ R−ρ r2 R−ρ We take the radial mean (as defined in equation (11)) to allow us to eliminate c(R − ρ, φ) and obtain c=c+ having used the notation
Z
r
g(r) = R−ρ
1 r2
1 ∂c [g(r) − g] , κ ∂φ
Z
(30)
r2
[v(r1 ) − ωr1 ] dr1 dr2 R−ρ
6
(31)
for brevity. The average flux of solvent through the line φ =constant (in the moving frame) is given by the radial average of the product of the concentration and angular velocity ω 0 (also relative to the moving frame): J
= = =
c ω0 ¸h Z R+ρ · i 1 v 1 ∂c c+ (g(r) − g) − ω r dr 2Rρ R−ρ κ ∂φ r Z R+ρ 1 ∂c 1 (v(r) − ωr) g(r)dr. κ ∂φ 2Rρ R−ρ
If we assume
∂c ∂c ≈ , ∂φ ∂φ
(32)
(33)
this implies that the radial mean concentration in the moving frame is governed by a one-dimensional diffusion equation ∂c ∂2c = D 2, (34) ∂t ∂φ where D is the Taylor dispersion coefficient for the annulus, defined by Z R+ρ 1 1 D = (v(r) − ωr) g(r) dr κ 2Rρ R−ρ Z R+ρ Z Z r 1 1 1 r2 = (v(r) − ωr) [v(r1 ) − ωr1 ] dr1 dr2 dr. κ 2Rρ R−ρ R−ρ r2 R−ρ This integral may be calculated for the velocity (3); after some algebra we find · µ ¶¸−2 ω 2 R2 1−γ D = 4γ 2 − (1 − γ 2 )2 ln × 24κγ 1+γ " µ ¶2 1+γ 1+γ 5 2 4 2 2 3 2 2 2 × 48γ (1 + γ ) + 120γ (1 − γ ) ln − 72γ (1 − γ ) (1 + γ ) ln 1−γ 1−γ µ ¶4 µ ¶5 # 1+γ 1+γ 2 4 2 2 4 2 4 −6γ(1 − γ ) (1 + γ ) ln + (1 − γ ) (3 + 10γ + 3γ ) ln . 1−γ 1−γ
(35)
(36)
In the γ ¿ 1 limit, this reduces to · D ≈ ωP e
¸ 2 346 3 γ+ γ ... . 105 1576
(37)
The solution of the one-dimensional diffusion equation (34) with initial condition from (14) is straightforward. Having found c, the concentration c and so the mixing measure m(t) are calculated from (30) and (16) respectively. The mixing measure has an exponential tail: m(t) ∼ exp(−2Dt) and so the mixing time TM is TM
1 ln ∼ 2D
µ
1 M
7
as t → ∞,
(38)
as M → 0
(39)
¶
in the Taylor dispersion regime. From (37), this yields the (nondimensional) asymptotic form µ ¶ · ¸ 1 1 105 18165 2 ωTM ≈ ln 1− γ + . . . for γ ¿ 1, γP e M 4 1576
(40)
when the assumptions made above are valid. The range of P e where the above approximate analysis holds will prove to be of great interest in later sections, so we now examine carefully the two basic assumptions that were made, and cast these into conditions on P e for the result (40) to hold. The first condition stems from the requirement that the azimuthal diffusion be negligible compared to ∂2c the radial diffusion; effectively this requires that the coefficient of ∂φ 2 in (28), actually its radial average, should be much less than the dispersion coefficient (35): κ ¿ Dr2 .
(41)
Consider this condition when γ ¿ 1, using the first term of (37): κ
¿
⇐⇒ P e2
À
2 γωP eR2 105 105 . 2
(42)
The resulting condition P e À 7.2 is analogous to that in a straight capillary, see for instance [8]. ∂c ∂c The second important approximation is the replacement of ∂φ by ∂φ in (33). This is valid if the ∂c coefficient of ∂φ in (30) is very small. The order of magnitude of this coefficient is g/κ, so the condition is Z R+ρ Z r Z 1 1 1 r2 r [v(r1 ) − ωr1 ] dr1 dr2 dr ¿ 1. (43) κ 2Rρ R−ρ R−ρ r2 R−ρ Again, this yields a simple condition if γ ¿ 1: 1 ωρ2 15 κ
¿
1
⇐⇒ P e
¿
15 . γ
(44)
Thus (42) and (44) provide rough bounds on P e for the lower and upper ranges of validity of the Taylor dispersion approximations.
5
Convection-dominated mixing: P e À 1
In this section we examine the solution of equation (8) when the molecular diffusion term is small in comparison to the convective term, i.e., the limit of large Pecl´et number. As the diffusion κ multiplies the terms with highest derivatives in (8), the limit κ → 0 yields a singularly-perturbed problem. Regular perturbation methods fail to give approximations which are uniformly valid in time for such problems, and so we introduce multiple time scales. Taking κ as a small parameter, the time scales are: t0 t1 t2 t3
= = = =
8
t κt3 κt2 κt.
(45)
As is standard in multiple scale expansions, these time scales are assumed to be effectively independent of each other, and so the time derivative in (8) may be replaced by ∂ ∂ ∂ ∂ ∂ = + 3t20 κ + 2t0 κ +κ , ∂t ∂t0 ∂t1 ∂t2 ∂t3 accurate to first order in κ. The corresponding expansion of the concentration is ¡ ¢ c = c0 + κc1 + O κ2 ,
(46)
(47)
where c0 and c1 are functions of all the time scales, as well as r and φ. Inserting the expansions (46) and (47) into equation (8) leads to the following equation for κ-independent terms: ∂c0 v(r) ∂c0 + = 0, (48) ∂t0 r ∂φ with solution
³ ³ v ´ v ´ c0 = 1 + A(t1 , t2 , t3 ) cos φ − t0 + B(t1 , t2 , t3 ) sin φ − t0 . r r Here A and B are functions satisfying the initial conditions A(0, 0, 0)
=
1
B(0, 0, 0) =
0.
In order to further specify A and B, it is necessary to consider the equation for terms of order κ: µ ¶ ∂c1 v(r) ∂c1 ∂c0 ∂c0 1 ∂ ∂c0 1 ∂ 2 c0 2 ∂c0 + = −3t0 − 2t0 − + r + 2 . ∂t0 r ∂φ ∂t1 ∂t2 ∂t3 r ∂r ∂r r ∂φ2
(49)
(50)
(51)
Note that the differential operator on the left hand side of this equation is the same as in (48), while the forcing terms on the right hand side all depend on c0 . Calculating the forcing terms using (49), and requiring that the solution c1 should not grow without bound as t0 → ∞ leads to the following equations for A and B: " # 2 ∂A 1 v2 2vv 0 v0 = − − 3 + 2 A ∂t1 3 r4 r r " # 2 ∂B 1 v2 2vv 0 v0 = − − 3 + 2 B ∂t1 3 r4 r r · 00 ¸ 0 ∂A 1 v v v = − − 2+ 3 B ∂t2 2 r r r · ¸ ∂B 1 v 00 v0 v = − 2+ 3 A ∂t2 2 r r r ∂A 1 = − 2A ∂t3 r ∂B 1 = − 2 B. (52) ∂t3 r
9
Here primes are used to denote differentiation with respect to r. The solutions of these equations satisfying the initial conditions (50) are: " Ã ! # · µ ¶ ¸ 2 1 v2 2vv 0 v0 1 v 00 v0 v 1 A = exp − − + t − t cos − − + t2 1 3 3 r4 r3 r2 r2 3 r r2 r3 " Ã ! # · µ ¶ ¸ 2 1 v2 2vv 0 v0 1 1 v 00 v0 v B = exp − (53) − 3 + 2 t1 − 2 t3 sin − − 2 + 3 t2 3 r4 r r r 3 r r r In the limit γ → 0, we approximate the velocity profile as in (5), and retain only the dominant terms in (53) to obtain " # µ ¶2 ¸ · 3 r−R 3 3 3 2 2 A = exp − ω t cos ω t γP e ρ 2γP e " # µ ¶2 · ¸ 3 3 r−R 3 3 B = exp − ω t sin ω 2 t2 . (54) γP e ρ 2γP e We have thus found a non-secular approximation to the concentration c for large P e and small γ. By combining the results (49) and (54), the corresponding mixing measure is calculated as # " µ ¶2 Z R+ρ 1 6 r−R 3 3 m(t) = exp − ω t r dr 2Rρ R−ρ γP e ρ µ 3 3¶ √ π 6ω t = F , (55) 2 γP e where F is defined as the monotonic function F (x) = x−1/2 erf(x1/2 ).
(56)
The nondimensional mixing time corresponding to a mixing measure value of M is therefore · ωTM ≈
γP e −1 F 6
µ
2M √ π
¶¸1/3 ,
(57)
and note in particular that this increases as P e1/3 when M and γ are fixed.
6
Numerical simulations
To check the asymptotic results derived in previous sections, and to extend to cases with initial conditions other than the simple form (14), we solve the convection-diffusion equation (8) numerically. A decomposition into a finite number N of Fourier modes in the azimuthal angle is employed: c(r, φ, t) =
N X 1 g0 (r, t) + gn (r, t) cos(nφ) + hn (r, t) sin(nφ), 2 n=1
(58)
and substitution into (8) yields partial differential equations for each gn and hn . Having solved the system of equations for the Fourier coefficients, the full concentration field may be constructed from (58), or the
10
0
0
10
10
Pe=8
−1
−1
10
10
Pe=46341 m(t)
m(t) Pe=2 −2
−2
10
10
Pe=16384
Pe=64
Pe=0.5
Pe=0.5
Pe=32
−3
10
−3
0
10
20
30
40
50
ωt
60
70
80
90
100
10
0
10
20
30
40
50
ωt
60
70
80
90
100
Figure 3: Mixing measure m(t) as a function of nondimensional time, calculated in numerical simulations with γ = 0.05 and various Pecl´et numbers as shown. Asymptotic predictions are also shown for: P e = 0.5 (dashed) using (23); P e = 64 (dotted) using (38); and P e = 16384 (dot-dash) using (55). mixing measures may be evaluated more directly from (16) by integrating over the angle and normalizing m(0) to unity: N ³g ´2 1 X 0 m(t) = − hci + g 2 + h2n . (59) 2 2 n=1 n Logarithmic plots of m(t) as a function of time at various P e values are shown in Figure 3 for γ = 0.05. For comparison we plot also the asymptotic forms of m(t), using (23) in the diffusive regime, (38) in the Taylor regime, and (55) in the convective regime. The asymptotic formulas fit the numerical results well, even at early times. Note that the straight-line asymptotes at very large times for the high-P e case are not reproduced by the multiple-time solution—this is apparently due to the neglect of the boundary conditions (9) in our analysis, and the subsequent importance of boundary layers at the inner and outer walls. From the numerical solution for m(t), it is straightforward to calculate the mixing times TM required for the measure to decay from its initial value of 1 to the value M . We choose three values of M for comparison with the asymptotic predictions: M = 0.3, 0.1, and 0.01, and investigate a wide range of Pecl´et numbers. For fixed values of γ and M , the asymptotic analysis of the previous sections predicts a mixing time proportional to P e in the diffusion-dominated regime, to P e−1 in the Taylor dispersion regime, and to P e1/3 when convection dominates. The numerical values of nondimensional times ωTM are plotted in Figure 4 along with the straight lines corresponding to the formulas (26), (40), and (57). Note the excellent agreement with predictions, except for the lowest value of M , when the multiple scale solution no longer accurately describes the time evolution of m(t). As much of our analysis has concentrated on the γ → 0 limit, it is worthwhile comparing numerics and asymptotics for a larger value of γ, noting that the maximum possible value of γ is 1. In Figure 5 we plot the results for γ = 0.2, and find that in general the correspondence between predicted and numerical values is excellent, but note that the limits on the Taylor regime now preclude the formation of a clear P e−1 scaling range. Another simplification in the analysis concerns the initial condition (14), which consists of only a single harmonic of the azimuthal variable. In order to examine the robustness of our predictions, we numerically
11
2
10
ωT
1
10
0
1
10
2
10
3
10
4
10
10
5
10
Pe
Figure 4: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.05. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles) 2
10
1
10 ωT
0
10
0
10
1
2
10
10
3
10
4
10
Pe
Figure 5: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.2. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles)
12
2
10
ωT
1
10
0
10
1
2
10
3
10
10
4
10
5
10
Pe
Figure 6: Nondimensional mixing times as a function of Pecl´et number, for γ = 0.05, and an initial concentration given by (15) (slightly smoothed). Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles) solve the convection-diffusion equation for initial condition (15), replacing the discontinuous function by hyperbolic tangents to avoid Gibbs oscillations. The decay of the mixing measure is still dominated by the first angular harmonic, and so the mixing times are remarkably close to the asymptotic estimates, see Figure 6. We therefore expect the formulas to be useful for rather general initial distributions of the solute. The nondimensional time ωTM used in Figures 3 to 6 may be replaced by the alternative nondimensionalization mentioned in Section 2, i.e., the diffusion time κTM /R2 . The timescales are related by κTM γ = ωTM , (60) 2 R Pe and so the data of, for example, Figure 4 is easily recast in terms of the diffusion time, see Figure 7. The relevant mixing times in the diffusive, Taylor, and convective regimes are thus found by multiplying (26), (40), and (57) by γ/P e, yielding µ ¶· ¸ κTM 1 1 1 2 ≈ ln 1 − γ + ... , (61) R2 2 M 3 for P e ¿ 7.2, κTM 1 ≈ ln 2 R P e2 for 7.2 ¿ P e ¿ 15/γ, and
µ
1 M
¶
· ¸ 105 18165 2 1− γ + ... 4 1576
· µ ¶¸1/3 2M κTM − 23 43 1 −1 √ ≈ Pe γ F R2 6 π
13
(62)
(63)
0
10
−1
10
−2
10 κ T/R2
−3
10
−4
10
−5
10
0
10
1
10
2
3
10
10
4
10
5
10
Pe
Figure 7: Mixing times nondimensionalized by the diffusion time, κT /R2 , as a function of Pecl´et number, for γ = 0.05. Asymptotic results are shown as lines, and numerical results as symbols for values of the mixing measure: M = 0.3 (dashed line; squares), M = 0.1 (solid line, points), and M = 0.01 (dot-dash line, triangles)
14
for P e À 15/γ. Figure 7 is especially of interest to experimentalists working with a particular solute and solvent (so that κ is fixed) while varying the rate of rotation velocity ω of the micromixer to change P e. The mixing time is seen to decrease from the diffusion time through the Taylor regime (at a rate proportional to P e−2 ), and then continue to decrease at a slower rate beyond the Taylor regime. The slower rate corresponds to a P e−2/3 scaling when (55) is valid, i.e., for M ≥ 0.1, but is closer to P e−1/2 for very small values of M . This is again a consequence of the failure of the multiple-scale method to account for the very long time scales where boundary layer effects become non-negligible. From our analysis of the limits of validity of the Taylor dispersion description, see (44), we can estimate the end of the Taylor regime as 15γ −1 for small γ. Of interest to the micromixer designer is the influence of the geometry ratio γ: as the mixing time decreases as P e−2 until P e ≈ 15γ −1 , it seems advisable to decrease γ as much as possible in order to achieve faster mixing at lower velocities. However, this approach is overly simplistic, since the Pecl´et number also depends on γ—the linear velocity in the channel decreases with decreasing ρ, all other things being equal. A problem a great relevance to experimentalists is to find the optimal annular geometry to give lowest possible mixing time for a given species, and given values of the electric and magnetic pumping fields (so κ and α are fixed). Suppose the “footprint” R of the device is chosen, then we must find the channel half-width ρ to minimize the mixing time. Since P e may be shown to increase monotonically with ρ, we conclude from Figure 7 that the mixing time decreases as the channel width increases. Practical considerations such as the minimal inner radius for effective electrodes will also provide an upper bound on the acceptable geometry ratio γ.
7
Summary
Laminar mixing in a two-dimensional annulus has been examined by asymptotic analysis and numerical simulation of the convection-diffusion equation (8). The mixing measure defined in (16) has been shown to decay exponentially in the diffusion and Taylor regimes (see equations (23) and (38)), with time constants given in terms of the Pecl´et number and the ratio γ characterizing the annulus geometry. Corresponding mixing times TM are predicted to scale as P e0 and P e−2 relative to the azimuthal diffusion time, see equations (61) and (62). Figures 3 to 7 demonstrate the robustness of these predictions by comparison to numerical simulations. In the convection-dominated regime, multiple scale analysis yields the first approximation (55) to the mixing measure, with mixing time (63) showing a P e−2/3 dependence on the Pecl´et number. This is again found to agree well with numerical results, at least for M ≥ 0.1. Although the asymptotic formulas were derived in the limit γ → 0, we find they yield good predictions for larger values of γ also (Figure 5). An important limitation of the present analysis is the assumption that all motion is two-dimensional; it is known [10, 11] that the dispersion in a channel of finite depth may not equal the value calculated by assuming the depth to be infinite. Work incorporating three-dimensional effects may yield further insight into this important problem.
8
Acknowledgements
This work has been supported by funding from the Enterprise Ireland International Collaboration Fund, the National Microelectronics Research Centre, and by the Faculty of Arts Research Fund, University College Cork.
15
References [1] J. P. Gleeson and J. West, Magnetohydrodynamic Micromixing, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 318-321. [2] I. Meisel and P. Ehrhard, Simulation of electrically excited flows in microchannels for mixing application, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 62-65. [3] K. Mohseni, Mixing and impulse extremization in microscale vortex formation, in Proceedings of the Fifth International Conference on Modeling and Simulation of Microsystems 2002, Computational Publications, 2002, pp. 392-395. [4] J. West et al., Application of magnetohydrodynamic actuation to continuous flow chemistry, Lab on a Chip, 2 (2002), pp. 224-230. [5] P. Tabeling and J. P. Chabrerie, Magnetohydrodynamic Taylor vortex flow under a transverse pressure gradient, Phys. Fluids, 24 (1981), pp. 406-412. [6] M-H. Chang and C-K. Chen, Hydromagnetic stability of current-induced flow in a small gap between concentric cylinders, J. Fluids Eng., 121 (1999), pp. 548-554. [7] G. I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proc. Roy. Soc. Lond. A, 219 (1953), pp. 186-203. [8] R. F. Probstein, Physiochemical Hydrodynamics, Wiley, New York, 1994. [9] R. J. Nunge, T.-S. Lin and W. N. Gill, Laminar dispersion in curved tubes and channels, J. Fluid Mech., 51 (1972), pp. 363-383. [10] H. Brenner and D. A. Edwards, Macrotransport Processes, Butterworth-Heinemann, 1993. [11] M. R. Doshi, P. M. Daiya and W. N. Gill, Three dimensional laminar dispersion in open and closed rectangular conduits, Chem. Eng. Sci., 33 (1978), pp. 795-804.
16
