Electroosmotic Flow in Porous Media

This example treats the modeling of electroosmotic flow in porous media (Ref. 1, Ref. 2). The system consists of a compartment of sintered porous material and two electrodes that generate an electric field. The cell combines pressure-driven and electroosmotic flow.

The purpose of this example is to illustrate the modeling of electroosmosis and electrophoresis in porous media in COMSOL Multiphysics.

Model Definition

Figure 1 shows the system’s geometry. It is made up of a domain of porous material containing the electrolyte and two electrodes that generate a potential difference. The conductivity is very small, and the model assumes that the effects of electrochemical reactions at the electrode surfaces are negligible.

Figure 1: The modeled domain consists of an electrolyte contained in a porous structure and two electrodes that generate an electric field.

In the first part of the model, you solve the continuity equations for the flow velocity and the current density at steady state:

Here u denotes the velocity (SI unit: m/s) and i represents the current-density vector (SI unit: A/m2). The velocity includes two driving forces — a pressure term and an electroosmotic term:

In this equation, εp denotes the porosity, a is the average radius of the pores (SI unit: m), μ gives the fluid’s dynamic viscosity (SI unit: Pa·s), τ represents the tortuosity of the porous structure, εw is the fluid’s permittivity (SI unit: F/m), p gives the pressure (SI unit: Pa), ζ is the zeta potential (SI unit: V), and V equals the potential (SI unit: V). This equation is solved with a General Form PDE interface. The current density is modeled with the Electric Currents interface and solves the following:

where κ denotes the conductivity (SI unit: S/m).

At the solid walls, the normal velocity component vanishes:

At the inlet and the outlet, the pressure is fixed.

The boundary conditions for the current-density balance are insulating for all boundaries except the electrode surfaces, where the potential is fixed:

In the second model stage, you use the steady-state velocity and potential fields in a transient simulation of the concentration of a charged tracer species injected into the system, assuming that the tracer species does not influence the conductivity or the set potential in the porous structure. The mass-transport equation for the tracer is solved with the Transport of Diluted Species interface and reads:

(1)

where N is the flux vector given by the Nernst-Planck equation:

In this equation, D denotes the tracer’s diffusivity (SI unit: m2/s), c gives its concentration (SI unit: mol/m3), z represents the tracer’s charge number, and F is Faraday’s constant (SI unit: C/mol). The mobility, um (SI unit: mol·m2/(J·s)), is given by the Nernst-Einstein equation

where Rg = 8.314 J/(mol·K) is the gas constant and T (SI unit: K) is the temperature.

The boundary conditions for Equation 1, the mass-transport equation, are insulating except at the inlet and the outlet. There you use the Flux condition to set the diffusive and convective contributions to the flux through the boundaries to zero:

The initial concentration is given by a distribution that is uniform in the y direction and bell-shaped in the x direction:

Here ctop denotes the peak concentration, xm is the position of the peak along the x-axis, and pw equals the base width of the peak.

Results and Discussion

The upper plot in Figure 2 shows the flow distribution in the porous structure when only the electric field acts as a driving force. In the lower plot, it is instead the pressure gradient that drives the flow. In both cases, the velocity is large around the corners of the electrodes, where the field strength is large and the effect of the decrease in the geometry’s cross section is most pronounced. A comparison of the maximum flow-velocity values shows that the pressure gradient is the dominating driving force.

Figure 2: Velocity distributions in the cell with an applied electric field (top) and an applied pressure difference (bottom). The unit for the magnitude surface plot is mm/s.

Figure 3 shows the concentration field in the case where both pressure and electroosmotic forces are included. The figure shows that the migration of the charged tracer also influences the transport rate. In this case, the tracer is transported both by the movement of the flow, due to pressure and electroosmotic terms, and the electrophoretic effect on the charged tracer itself. At t = 0, the tracer is introduced as a bell-shaped vertical distribution near the inlet. In the subsequent simulation, this distribution is sheared and transported by diffusion, migration, and convection.

Figure 3: Concentration distribution in the domain at different times after injection. The effects of diffusion, convection, and migration deform the initial bell-shaped pulse.

Figure 4 shows the cross sections of the pulse along a line at y = 2.5 mm for the times t = 0, 0.5 s, 1.0 s, 1.5 s, and t = 2.0 s. Diffusion, mainly numerical diffusion, smears out the pulse, while migration and convection mainly translate and shear it, depending on the pressure and the electric field.

Figure 4: Cross-section plot of the concentration along a horizontal line at y = 2.5 mm for the times t = 0 (blue solid), 0.5 s, 1.0 s, 1.5 s, and 2.0 s.

Despite its geometrical simplicity, this example includes the components needed to model complex-shaped domains.

References

1. M-S. Chun, “Electrokinetic Flow Velocity in Charged Slit-like Microfluidic Channels with Linearized Poisson-Boltzmann Field,” Korean J. Chem. Eng., vol. 19, no. 5, pp. 729–734, 2002.

2. S. Yao, D. Huber, J. Mikkelsen, and J.G. Santiago, “A Large Flowrate Electroosmotic Pump with Micron Pores,” Proc. IMECE, 2001 ASME International Mechanical Engineering Congress and Exposition, New York, November 2001.

Notes About the COMSOL Implementation

The COMSOL Multiphysics implementation is straightforward, and the only aspect to remember is to solve the steady-state problem first and then the time-dependent problem.

Chemical_Reaction_Engineering_Module/Electrokinetic_Effects/electroosmotic_flow

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 .

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

For your convenience, parameters and variables are made available in text files.

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file electroosmotic_flow_parameters.txt.

In these expressions, the predefined physical constants epsilon0_const and R_const refer to the permittivity of vacuum and the gas constant, respectively.

Definitions

Variables 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file electroosmotic_flow_variables.txt.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Now create the geometry. To simplify this step, insert the geometry sequence directly from the model’s Application Libraries file:

In the toolbar, click .

Browse to the model’s Application Libraries folder and double-click the file electroosmotic_flow.mph.

Rectangle 1 (r1)

Make all the necessary geometry selections.

Cathode

In the toolbar, click

and choose .

In the window for , type Cathode in the text field.

Locate the section. From the list, choose .

On the object , select Boundaries 4, 5, 11, and 12 only.

Anode

In the toolbar, click

and choose .

In the window for , type Anode in the text field.

Locate the section. From the list, choose .

On the object , select Boundaries 7, 8, 13, and 14 only.

Inlet

In the toolbar, click

and choose .

In the window for , type Inlet in the text field.

Locate the section. From the list, choose .

On the object , select Boundary 1 only.

Outlet

In the toolbar, click

and choose .

In the window for , type Outlet in the text field.

Locate the section. From the list, choose .

On the object , select Boundary 10 only.

In the Electric Currents interface set the electrolyte properties to compute the current density within the device.

Electric Currents (ec)

Current Conservation 1

In the window, under click .

In the window for , locate the section.

From the σ list, choose . In the associated text field, type kappa0.

Locate the section. From the εr list, choose .

Electric Potential 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Electric Potential 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the V0 text field, type V_anode.

Continue with the General Form PDE interface and set up a convective velocity depending on pressure and electroosmosis. The dependent value is the pressure.

Electroosmotic Pressure

In the window, under click .

In the window for , type Electroosmotic Pressure in the text field.

Locate the section. Click

In the dialog box, type pressure in the text field.

Click

In the tree, select .

Click .

In the window for , locate the section.

In the table, enter the following settings:


Source term quantity	Unit
Custom unit	1/s

Click to expand the section. In the text field, type p.

In the table, enter the following settings:


p

General Form PDE 1

In the window, under click .

In the window for , locate the section.

Specify the Γ vector as


u_flow	x
v_flow	y

To get the correct equation, finish by setting the source and mass terms to zero.

Locate the section. In the f text field, type 0.

Locate the section. In the da text field, type 0.

Inlet - p=0

In the toolbar, click

and choose .

As boundary conditions, set pressures at the system’s inlet and outlet.

In the window for , type Inlet - p=0 in the text field.

Locate the section. From the list, choose .

Outlet - p=p1

In the toolbar, click

and choose .

In the window for , type Outlet - p=p1 in the text field.

Locate the section. From the list, choose .

Locate the section. In the r text field, type p1.

Last, fill in the mass transport properties of the charged tracer species.

Transport of Diluted Species (tds)

In the window, under click .

In the window for , locate the section.

Select the check box.

Make all necessary interface couplings in the Transport Properties node.

Transport Properties 1

In the window, under click .

In the window for , locate the section.

From the T list, choose . In the associated text field, type T.

Locate the section. In the Dc text field, type D.

Locate the section. In the zc text field, type zn.

Locate the section. Specify the u vector as


u_flow	x
v_flow	y

Locate the section. From the V list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

In the c text field, type c_init.

Set the sum of the convective and diffusive fluxes at the inlet and outlet boundaries to zero with a Flux condition equal to the migration flux (electrophoretic flux=tds.nmflux_c).

Flux 1

In the toolbar, click

and choose .

Select Boundaries 1 and 10 only.

In the window for , locate the section.

Select the check box.

In the J0,c text field, type -tds.nmflux_c.

Before solving, improve the .

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Solve the model in two steps: First, compute a stationary (steady-state) solution for the velocity and current density distributions. Then, use this solution to compute the time-dependent tracer concentration distribution.

Study 1

Step 1: Stationary

In the window, under click .

In the window for , locate the section.

In the table, clear the check box for .

Time Dependent

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type range(0,0.1,2).

From the list, choose .

In the text field, type 0.005.

Locate the section. In the table, clear the check boxes for and .

In the toolbar, click

Results

Concentration (tds)

The surface plot in the third plot group shows the concentration at t = 2 s. Under in the model tree, you can add an Animation node to see how the distribution evolves with time. Alternatively, to reproduce the plot in Figure 3, showing the concentration distribution at four different times, follow the steps below.

Streamline 1

In the window, expand the node.

Right-click and choose .

Surface 1

In the window, click .

In the window for , locate the section.