COMSOL 6.4 - Electroosmotic Micromixer

Electroosmotic Micromixer¹

Microlaboratories for biochemical applications often require rapid mixing of different fluid streams. At the microscale, flow is usually highly ordered laminar flow, and the lack of turbulence makes diffusion the primary mechanism for mixing. While diffusional mixing of small molecules (and therefore of rapidly diffusing species) can occur in a matter of seconds over distances of tens of micrometers, mixing of larger molecules such as peptides, proteins, and high molecular-weight nucleic acids can require equilibration times from minutes to hours over comparable distances. Such delays are impractically long for many chemical analyses. These problems have led to an intense search for more efficient mixers for microfluidic systems.

Most microscale mixing devices are either passive mixers that use geometrical stirring, or active mixers that use moving parts or external forces, such as pressure or electric field.

In a passive mixer, one way of increasing the mixing is by “shredding” two or several fluids into very thin alternating layers, which decreases the average diffusion length for the molecules between the different fluids. However, these mixers often require very long mixing channels because the different fluids often run in parallel. Another way of improving mixing efficiency is to use active mixers with moving parts that stir the fluids. At the microscale level moving parts in an active mixer are very fragile. One alternative is to use electroosmotic effects to achieve a mixing effect that is perpendicular to the main direction of the flow.

This model takes advantage of electroosmosis to mix fluids. The system applies a time-dependent electric field, and the resulting electroosmosis perturbs the parallel streamlines in the otherwise highly ordered laminar flow.

Model Definition

This example of a rather simple micromixer geometry (Figure 1) combines two fluids entering from different inlets into a single 10 μm wide channel. The fluids then enter a ring-shaped mixing chamber that has four microelectrodes placed on the outer wall at angular positions of 45, 135, −45, and −135 degrees, respectively. Assume that the aspect ratio (channel depth to width) is large enough that you can model the mixer using a 2D cross-sectional geometry. The material parameters relevant for the model are given in Table 1.

Figure 1: Geometry of the micromixer with four symmetric electrodes on the wall of the mixing chamber. This example does not model the two inlet channels. Here you assume a parabolic inflow at the beginning of the computational domain (the gray area).

The Navier–Stokes equations for incompressible flow describe the flow in the channels:

Here η denotes the dynamic viscosity (SI unit: kg/(m·s)), u is the velocity (SI unit: m/s), ρ equals the fluid density (SI unit: kg/m3), and p refers to the pressure (SI unit: Pa).

Because you do not model the two inlet channels, assume that the entrance channel starts at a position where the flow has a fully developed laminar profile. The mixed fluid flows freely out of the right end boundary, where you specify vanishing total stress components normal to the boundary:

When brought into contact with an electrolyte, most solid surfaces acquire a surface charge. In response to the spontaneously formed surface charge, a charged solution forms close to the liquid-solid interface. Known as an electric double layer, it forms because of the charged groups located on the surface that faces the solution. When the operator applies an electric field, the electric field generating the electroosmotic flow displaces the charged liquid in the electric double layer. This scheme imposes a force on the positively charged solution close to the wall surface, and the fluid starts to flow in the direction of the electric field. The velocity gradients perpendicular to the wall give rise to viscous transport in this direction. In the absence of other forces, the velocity profile eventually becomes almost uniform in the cross section perpendicular to the wall.

This model replaces the thin electric double layer with the Helmholtz–Smoluchowski relation between the electroosmotic velocity and the tangential component of the applied electric field:

In this equation, εw = ε0εr denotes the fluid’s electric permittivity (F/m), ζ0 represents the zeta potential at the channel wall (V), and V equals the potential (V). This equation applies on all boundaries except for the entrance and the outlet.

Assuming that there are no concentration gradients in the ions that carry the current, you can express the current balance in the channel with Ohm’s law and the balance equation for current density

where σ denotes conductivity (S/m) and the expression within parentheses represents the current density (A/m2).

The electric potentials on the four electrodes are sinusoidal in time with the same maximum value (V0 = 0.1 V) and the same frequency (8 Hz), but they alternate in polarity. The potentials on electrodes 1 and 3 are V0sin(2πft), whereas those on electrodes 2 and 4 are −V0sin(2πft) (see Figure 1).

Assume all other boundaries are insulated. The insulation boundary condition

sets the normal component of the electric field to zero.

At the upper half of the inlet (see Figure 1) the solute has a given concentration, c0; at the lower half the concentration is zero. Thus, assume that the concentration changes abruptly from zero to c0 at the middle of the inlet boundary. The mixed solution flows out from the right outlet by convection, and all other boundaries are assumed insulated.

Inside the mixer, the following convection–diffusion equation describes the concentration of the dissolved substances in the fluid:

Here c is the concentration, D represents the diffusion coefficient, R denotes the reaction rate, and u equals the flow velocity. In this model R = 0 because the concentration is not affected by any reactions.

Table 1: Model Input Data.
Parameter	Value	Description
ρ	1000 kg/m3	Density of the fluid
η	10-3 Pa·s	Dynamic viscosity of the fluid
U0	0.1 mm/s	Average velocity through the inlet
εr	80.2	Relative electric permittivity of the fluid
ζ	-0.1 V	Zeta potential on the wall-fluid boundary
σ	0.11845 S/m	Conductivity of the ionic solution
D	10-11 m2/s	Diffusion coefficient
c0	1 mol/m3	Initial concentration

Results and Discussion

Figure 2 shows a typical instantaneous streamline pattern. It reveals that electroosmotic recirculation of the fluid vigorously stirs the flow, typically in the form of two rotating vortices near the electrodes.

The fundamental processes of effective mixing involve a combination of repeated stretching and folding of fluid elements in combination with diffusion at small scales. As the system applies the AC field (Figure 3), the resulting electroosmotic flow perturbs the laminar pressure-driven flow such that it pushes the combined stream pattern up and down at the beginning of the mixing chamber, causing extensive folding and stretching of material lines.

Figure 2: Fluid streamlines in an electroosmotic micromixer at t = 0.0375 s.

Figure 3: Electric potential lines for an electroosmotic micromixer. The contour lines show the shape when the device uses maximal potentials (±V0).

The following plots further exemplify how the mixer operates. Figure 4 shows the concentration at steady state when the electric field is not applied. The flow is laminar and the diffusion coefficient is very small, so the two fluids are well separated also at the outlet. When the alternating electric field is applied, the mixing increases considerably owing to the alternating swirls in the flow. Figure 5 depicts the system at the instant when the electric field and the electroosmotic velocity have their largest magnitudes during the cycle (that is, when |sin ωt| = 1). From the plot you can estimate that the concentration at the output fluctuates with the same frequency as the electric field. Thus, this mixer should be further improved to get a steadier output.

Figure 4: Steady-state solution in the absence of an electric field.

Figure 5: Time-dependent solution at the time when the alternating electric field has its largest magnitude.

This example demonstrates a rather simple and effective use of electrokinetic forces for mixing. The scheme is easy to implement, and you can easily control both the amplitude and the frequency. At low Reynolds numbers the inertial forces are small, which makes it possible to calculate stationary streamlines patterns using the parametric solver to control amplitude.

Notes About the COMSOL Implementation

Cummings and others (Ref. 3) have shown that in order to use the Helmholtz–Smoluchowski equation at the fluid-solid boundaries, the electric field must be at least quasistatic to neglect transient effects. In other words, the time scale of the unsteady electric field must be much larger than that of the transient flow. Y.T. Zhang and others (Ref. 1) estimated that the time scale of the transient effect in the modeled micromixer (with a channel width of 10 microns) is roughly 0.0127 s. In this simulation, the frequency of the applied electric potential is 8 Hz, which corresponds to a time scale of the electric field 10 times larger than that of the flow.

Because you can model the time-dependent electric field as a product of a stationary electric field and a time-dependent phase factor (sinωt), it is possible to reduce the simulation time and memory requirements by dividing the solution into two stages. In the first, calculate the amplitude of the electric potential field and the initial state for the time-dependent flow model using a stationary solver. In the second stage, you deactivate the Electric Currents interface and calculate the transient solution for the Laminar Flow and the Transport of Diluted Species interfaces. You obtain the tangential electric field components used in the electroosmotic velocity boundary condition by multiplying the stationary DC solution by sin(ωt). This approach is permissible because there is only a unidirectional dependence between the electric field and the fluid fields.

References

1. H. Chen, Y.T. Zhang, I. Mezic, C.D. Meinhart, and L. Petzold, “Numerical Simulation of an Electroosmotic Micromixer,” Proc. Microfluidics 2003 (ASME IMECE), 2003.

2. Y.T. Zhang, H. Chen, I. Mezic, C.D. Meinhart, L. Petzold, and N.C. MacDonald, “SOI Processing of a Ring Electrokinetic Chaotic Micromixer,” Proc. NSTI Nanotechnology Conference and Trade Show (Nanotech 2004), vol. 1, pp. 292–295, 2004.

3. E. Cummings, S. Griffiths, R. Nilson, and P. Paul, “Conditions for Similitude Between the Fluid Velocity and the Electric Field in Electroosmotic Flow,” Anal. Chem., vol. 72, pp. 2526–2532, 2000.

Microfluidics_Module/Micromixers/electroosmotic_mixer

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
U0	0.1[mm/s]	1E-4 m/s	Mean inflow velocity
sigma_w	0.11845[S/m]	0.11845 S/m	Conductivity of the ionic solution
eps_r	80.2	80.2	Relative permittivity of the fluid
zeta	-0.1[V]	-0.1 V	Zeta potential
V0	0.1[V]	0.1 V	Maximum value of the AC potential
omega	2pi[rad]8[Hz]	50.265 Hz	Angular frequency of the AC potential
t	0[s]	0 s	Start time
D	1e-11[m^2/s]	1E-11 m²/s	Diffusion coefficient of the solution
c0	1[mol/m^3]	1 mol/m³	Initial concentration

You need the constant t (used in the scalar expressions below) when first solving the model using a stationary solver. In the time-dependent simulation, the internal time variable, t, overwrites this constant (the red color is just a warning signaling that t is an internal variable).

Now define a smoothed step function that you will later use to impose a step in the concentration in the middle of the channel entrance.

Step 1 (step1)

In the toolbar, click

and choose > .

In the window for , click to expand the section.

In the text field, type 0.1e-6.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 80.

In the text field, type 10.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 15.

To add vertices for the electrode endpoints on the outer boundary, use the Partition Edges.

Partition Edges 1 (pare1)

In the toolbar, click

and choose .

On the object , select Boundaries 1–4 only.

In the window for , locate the section.

In the table, enter the following settings:


Relative arc length parameters
0.47
0.53

Circle 2 (c2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

Compose 1 (co1)

In the toolbar, click

and choose .

Click in the window and then press Ctrl+A to select all objects.

In the window for , locate the section.

In the text field, type (r1+pare1)-c2.

Clear the checkbox.

Form Union (fin)

In the window, click .

The model geometry is now complete.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1e3[kg/m^3]	kg/m³	Basic
Dynamic viscosity	mu	1e-3[Pa*s]	Pa·s	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_w	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	eps_r	1	Basic