PDF
Electroosmotic Flow in Porous Media
Introduction
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).
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.
Application Library path: Chemical_Reaction_Engineering_Module/Electrokinetic_Effects/electroosmotic_flow
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select AC/DC > Electric Fields and Currents > Electric Currents (ec).
3
Click Add.
4
In the Select Physics tree, select Mathematics > PDE Interfaces > General Form PDE (g).
5
Click Add.
6
In the Select Physics tree, select Chemical Species Transport > Transport of Diluted Species (tds).
7
Click Add.
8
Click  Study.
9
In the Select Study tree, select General Studies > Stationary.
10
Click  Done.
Global Definitions
Parameters 1
For your convenience, parameters and variables are made available in text files.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
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
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file electroosmotic_flow_variables.txt.
Now create the geometry. To simplify this step, insert the geometry sequence directly from the model’s Application Libraries file:
Geometry 1
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose mm.
4
In the Geometry toolbar, click Insert Sequence.
5
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
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Cathode in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object dif1, select Boundaries 4, 5, 11, and 12 only.
Anode
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Anode in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object dif1, select Boundaries 7, 8, 13, and 14 only.
Inlet
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Inlet in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object dif1, select Boundary 1 only.
Outlet
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Outlet in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object dif1, 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 in Fluids 1
1
In the Physics toolbar, click  Domains and choose Current Conservation in Fluids.
2
In the Settings window for Current Conservation in Fluids, locate the Domain Selection section.
3
From the Selection list, choose All domains.
4
Locate the Constitutive Relation Jc-E section. From the σ list, choose User defined. In the associated text field, type kappa0.
5
Locate the Constitutive Relation D-E section. From the εr list, choose User defined.
Electric Potential 1
1
In the Physics toolbar, click  Boundaries and choose Electric Potential.
2
In the Settings window for Electric Potential, locate the Boundary Selection section.
3
From the Selection list, choose Cathode.
Electric Potential 2
1
In the Physics toolbar, click  Boundaries and choose Electric Potential.
2
In the Settings window for Electric Potential, locate the Boundary Selection section.
3
From the Selection list, choose Anode.
4
Locate the Electric Potential 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 variable is pressure.
Electroosmotic Pressure
1
In the Model Builder window, under Component 1 (comp1) click General Form PDE (g).
2
In the Settings window for General Form PDE, type Electroosmotic Pressure in the Label text field.
3
Locate the Units section. Click  Select Dependent Variable Quantity.
4
In the Physical Quantity dialog, type pressure in the text field.
5
In the tree, select General > Pressure (Pa).
6
Click OK.
7
In the Settings window for General Form PDE, locate the Units section.
8
In the Source term quantity table, enter the following settings:
Source term quantity
Unit
Custom unit
1/s
9
Click to expand the Dependent Variables section. In the Field name (Pa) text field, type p.
10
In the Dependent variables (Pa) table, enter the following settings:
p
General Form PDE 1
1
In the Model Builder window, under Component 1 (comp1) > Electroosmotic Pressure (g) click General Form PDE 1.
2
In the Settings window for General Form PDE, locate the Conservative Flux section.
3
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.
4
Locate the Source Term section. In the f text field, type 0.
5
Locate the Damping or Mass Coefficient section. In the da text field, type 0.
Inlet - p=0
1
In the Physics toolbar, click  Boundaries and choose Dirichlet Boundary Condition.
As boundary conditions, set pressures at the system’s inlet and outlet.
2
In the Settings window for Dirichlet Boundary Condition, type Inlet - p=0 in the Label text field.
3
Locate the Boundary Selection section. From the Selection list, choose Inlet.
Outlet - p=p1
1
In the Physics toolbar, click  Boundaries and choose Dirichlet Boundary Condition.
2
In the Settings window for Dirichlet Boundary Condition, type Outlet - p=p1 in the Label text field.
3
Locate the Boundary Selection section. From the Selection list, choose Outlet.
4
Locate the Dirichlet Boundary Condition 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)
1
In the Model Builder window, under Component 1 (comp1) click Transport of Diluted Species (tds).
2
In the Settings window for Transport of Diluted Species, locate the Transport Mechanisms section.
3
Select the Migration in electric field checkbox.
Species Charges
1
In the Model Builder window, under Component 1 (comp1) > Transport of Diluted Species (tds) click Species Charges.
2
In the Settings window for Species Properties, locate the Charge section.
3
In the zc text field, type zn.
Make all necessary interface couplings in the Fluid node.
Fluid 1
1
In the Model Builder window, click Fluid 1.
2
In the Settings window for Fluid, locate the Model Input section.
3
From the T list, choose User defined. In the associated text field, type T.
4
Locate the Diffusion section. In the Dc text field, type D.
5
Locate the Convection section. Specify the u vector as
u_flow
x
v_flow
y
6
Locate the Migration in Electric Field section. From the V list, choose Electric potential (ec).
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
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
1
In the Physics toolbar, click  Boundaries and choose Flux.
2
Select Boundaries 1 and 10 only.
3
In the Settings window for Flux, locate the Inward Flux section.
4
Select the Species c checkbox.
5
In the J0,c text field, type -tds.nmflux_c.
Before solving, improve the Mesh.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Extremely fine.
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
1
In the Model Builder window, under Study 1 click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Transport of Diluted Species (tds).
Step 2: Time Dependent
1
In the Study toolbar, click  Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.1,2).
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 0.005.
6
Locate the Physics and Variables Selection section. In the Solve for column of the table, under Component 1 (comp1), clear the checkboxes for Electric Currents (ec) and Electroosmotic Pressure (g).
7
In the Study toolbar, click  Compute.
Results
Concentration (tds)
The surface plot in the third plot group shows the concentration at t = 2 s. Under Results > Export 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.
Arrow Surface 1
1
In the Model Builder window, expand the Concentration (tds) node.
2
Right-click Arrow Surface 1 and choose Disable.
Surface 1
1
In the Model Builder window, click Surface 1.
2
In the Settings window for Surface, locate the Data section.
3
From the Dataset list, choose Study 1/Solution 1 (sol1).
4
From the Time (s) list, choose 0.
5
Locate the Coloring and Style section. From the Color table list, choose TrafficLight.
6
Click the  Zoom Extents button in the Graphics toolbar.
7
In the Concentration (tds) toolbar, click  Plot.
8
Locate the Data section. From the Time (s) list, choose 0.6.
9
Click the  Zoom Extents button in the Graphics toolbar.
10
In the Concentration (tds) toolbar, click  Plot.
11
From the Time (s) list, choose 1.2.
12
Click the  Zoom Extents button in the Graphics toolbar.
13
In the Concentration (tds) toolbar, click  Plot.
14
From the Time (s) list, choose 1.8.
15
Click the  Zoom Extents button in the Graphics toolbar.
16
In the Concentration (tds) toolbar, click  Plot.
Proceed to study the effects of each of the two driving forces — the pressure gradient and the electric field — on the velocity field by reproducing the plots in Figure 2. First visualize the velocity field that results from the electric field alone.
Velocity Electroosmotic Term
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Velocity Electroosmotic Term in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Solution Store 1 (sol2).
4
Locate the Color Legend section. Select the Show units checkbox.
Velocity Magnitude
1
Right-click Velocity Electroosmotic Term and choose Surface.
2
In the Settings window for Surface, type Velocity Magnitude in the Label text field.
3
Click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Definitions > Variables > U_eo - Flow-velocity electroosmosis term, magnitude - m/s.
4
Locate the Expression section. From the Unit list, choose mm/s.
Velocity Field
1
In the Model Builder window, right-click Velocity Electroosmotic Term and choose Arrow Surface.
2
In the Settings window for Arrow Surface, type Velocity Field in the Label text field.
3
Locate the Expression section. In the x-component text field, type u_eo.
4
In the y-component text field, type v_eo.
5
Select the Description checkbox. In the associated text field, type Velocity Electroosmotic Term.
6
Click the  Zoom Extents button in the Graphics toolbar.
7
In the Velocity Electroosmotic Term toolbar, click  Plot.
Next, plot the velocity field that results from the pressure gradient.
Velocity Pressure Term
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Velocity Pressure Term in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Solution Store 1 (sol2).
4
Locate the Color Legend section. Select the Show units checkbox.
Velocity Magnitude
1
Right-click Velocity Pressure Term and choose Surface.
2
In the Settings window for Surface, type Velocity Magnitude in the Label text field.
3
Click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Definitions > Variables > U_p - Velocity pressure term, magnitude - m/s.
4
Locate the Expression section. From the Unit list, choose mm/s.
Velocity Field
1
In the Model Builder window, right-click Velocity Pressure Term and choose Arrow Surface.
2
In the Settings window for Arrow Surface, type Velocity Field in the Label text field.
3
Locate the Expression section. In the x-component text field, type u_p.
4
In the y-component text field, type v_p.
5
Select the Description checkbox. In the associated text field, type Velocity Pressure Term.
6
Click the  Zoom Extents button in the Graphics toolbar.
7
In the Velocity Pressure Term toolbar, click  Plot.
Finally, reproduce the cross-section plot of the concentration along the line y = 2.5 mm shown in Figure 4. First define a Cut Line 2D dataset.
Cut Line 2D 1
1
In the Results toolbar, click  Cut Line 2D.
2
In the Settings window for Cut Line 2D, locate the Line Data section.
3
In row Point 1, set x to -4 and y to 2.5.
4
In row Point 2, set x to 4 and y to 2.5.
Concentration Along x-axis @ y=2.5mm
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Concentration Along x-axis @ y=2.5mm in the Label text field.
3
Locate the Data section. From the Dataset list, choose Cut Line 2D 1.
4
From the Time selection list, choose From list.
5
In the Times (s) list, choose 0, 0.5, 1, 1.5, and 2.
6
Locate the Legend section. From the Layout list, choose Outside graph axis area.
Line Graph 1
1
Right-click Concentration Along x-axis @ y=2.5mm and choose Line Graph.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Transport of Diluted Species > Species c > c - Molar concentration, c - mol/m³.
3
Click Replace Expression in the upper-right corner of the x-Axis Data section. From the menu, choose Component 1 (comp1) > Geometry > Coordinate > x - x-coordinate.
4
Click to expand the Legends section. Select the Show legends checkbox.
5
From the Legends list, choose Automatic.
6
In the Concentration Along x-axis @ y=2.5mm toolbar, click  Plot.
Finish by renaming the second plot group.
Pressure
1
In the Model Builder window, under Results click Electroosmotic Pressure.
2
In the Settings window for 2D Plot Group, type Pressure in the Label text field.