INTRODUCTION TO

Electrochemistry Module
Introduction to the Electrochemistry Module
© 1998–2018 COMSOL
Protected by patents listed on [Link]/patents, and U.S. Patents 7,519,518; 7,596,474; 7,623,991; 8,457,932;
8,954,302; 9,098,106; 9,146,652; 9,323,503; 9,372,673; and 9,454,625. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/comsol-license-agreement) and may be used or copied only under the terms of the
license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Server, and LiveLink are either
registered trademarks or trademarks of COMSOL AB. All other trademarks are the property of their respective
owners, and COMSOL AB and its subsidiaries and products are not affiliated with, endorsed by, sponsored by, or
supported by those trademark owners. For a list of such trademark owners, see [Link]/trademarks.
Version: COMSOL 5.4

Contact Information
Visit the Contact COMSOL page at [Link]/contact to submit general inquiries, contact
Technical Support, or search for an address and phone number. You can also visit the Worldwide
Sales Offices page at [Link]/contact/offices for address and contact information.

If you need to contact Support, an online request form is located at the COMSOL Access page at
[Link]/support/case. Other useful links include:

• Support Center: [Link]/support

• Product Download: [Link]/product-download
• Product Updates: [Link]/support/updates
• COMSOL Blog: [Link]/blogs
• Discussion Forum: [Link]/community
• Events: [Link]/events
• COMSOL Video Gallery: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part No. CM024302

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
The Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
The Electrochemistry Module Physics Interfaces . . . . . . . . . . . . . 7
Physics Interface Guide by Space Dimension and Study Type . . . 10
Orange Battery. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Model Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Cyclic Voltammetry at a Macroelectrode in 1D . . . . . . . . . . . . 38
Model Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

|3
4|
Introduction

The Electrochemistry Module is intended for the modeling and simulation of

generic electrochemical cells.
The module defines components in 1D, 2D, and 3D geometries that describe the
electrochemical reactions, and other surface reactions, at the interface between a
metal electrode and an electrolyte, as well as the transport of ions and neutral
species in the electrolyte, including possible homogeneous reactions. The current
conduction within the metal electrode can also be modeled. The simulations can
be used to understand electrolytic processes as well as to design and optimize
electrochemical cells.
The descriptions in the Electrochemistry Module allow for the simulation of
systems at different scales and at different levels of detail, ranging from
rudimentary current distribution analysis of industrial cells in the range of meters
down to the detailed chemical and electrochemical phenomena within a pore of a
porous electrode, or at, and in the vicinity of, a microelectrode. In terms of size,
geometric complexity, and complexity in the described phenomena, the
Electrochemistry Module is able to handle these modeling extremes and anything
in between.

The Applications
The modeling and simulation capabilities of the Electrochemistry Module cover
processes such as electrode reactions, electrolyte transport and homogeneous
reactions. The capabilities also include, for example, transport in porous media
and heat transfer.
Figure 1 shows a typical simulation result from a tutorial example, which is
available in the module’s application library. The plot shows the current density
on and the concentration field around a wire mesh electrode. On the electrode, an
electrode reaction is consuming a reactant which is present in the surrounding
electrolyte. The simulation accounts for the concentration dependent electrode
kinetics, the flow of the electrolyte and the transport of the reactant in the
electrolyte. (You will explore how to model some of these phenomena in the first
modeling example of this introduction book.) The simulation predicts the current
distribution in the cell for different cell potentials. This type of simulation may be
relevant for any electrolytic manufacturing process. It is then important to

|5
estimate the optimum usage and distribution of catalytic material, which is often
based on expensive noble metals.

Figure 1: Current distribution and surrounding concentration field of a wire mesh electrode.

Possible applications include the study and design of chlor-alkali and chlorate
electrolysis, water electrolysis for hydrogen and oxygen production, waste water
treatment, desalination of seawater, electrophoretic separation of proteins,
fundamental electrochemical studies in electrocatalysis and electroanalysis, and
sensors for glucose, pH, hydrogen, and other gases. In the field of electroanalysis,
the Electrochemistry module offers tailor-made interfaces for cyclic voltammetry
and impedance analysis (Figure 2).

Figure 2: Simulated cyclic voltammograms (left) and Nyquist plots (right) in COMSOL..

6|
The Electrochemistry Module Physics Interfaces

The module has a number of physics interfaces, which describe the transport and
reaction of charged and neutral ions in aqueous electrolytes, surface reactions
(electrochemical and chemical) and the formation of adsorbed species or coatings
on metal surfaces. It also describes transport of current in electrolytes and metals,
transport of current in thin metal structures (treated as shells), fluid flow, and heat
transfer.
Figure 3 shows the available physics interfaces as displayed in the Model Wizard for
a 3D model.

Figure 3: The physics interfaces for the Electrochemistry Module as shown in the Model Wizard for a 3D
model.

|7
T HE C URRENT D ISTRIBUTION INTERFACES
The Electrochemistry interfaces include the generic Primary Current Distribution
( ), Secondary Current Distribution ( ), and Tertiary Current Distribution,
Nernst Planck ( ) interfaces.
The Primary Current Distribution interface neglects the activation losses for the
charge transfer reactions. It may only be used in the cases where the activation
losses are substantially smaller than the conductivity losses. The Secondary
Current Distribution interface accounts for the activation losses for the charge
transfer reactions.
The Primary and Secondary Current Distribution interfaces may be combined
with a Chemical Species Transport interface (described below) in order to
incorporate kinetics effects of active species in the electrolyte or adsorbed species
on an electrode surface.
The Tertiary Current Distribution, Nernst Planck ( ) interface accounts for the
transport of species through diffusion, migration, and convection and is therefore
able to describe the effects of variations in composition on the electrochemical
processes. The kinetics expressions for the electrochemical reactions account for
both activation and concentration overpotential.

T HE E LECTROANALYSIS I NTERFACE
The Electroanalysis interface ( ) models mass transport of diluted species in
electrolytes using the diffusion-convection equation, solving for electroactive
species concentrations. The physics interface is applicable for electrolyte solutions
containing a large quantity of excess (supporting) electrolyte. Ohmic losses are
assumed to be negligible. The physics interface includes tailor-made functionality
for setting up cyclic voltammetry problems.

T HE C HEMICAL S PECIES T RANSPORT

The Transport of Diluted Species interface ( ) is available under the Chemical
Species Transport branch. This physics interface can be used to model transport of
species due to diffusion, migration and convection.
The Nernst-Planck-Poisson Equations ( ) interface can be used for investigation
of charge and ion distributions within the electrochemical double layer where
charge neutrality cannot be assumed. A requisite when using this interface is that
the double layer, which typically is in the range of tens nanometers, is fully
resolved in the mesh.
The Transport of Diluted Species in Porous Media interface ( ) is also available
and describes species transport between the fluid, solid, and gas phases in saturated
and variably saturated porous media. The interface for instance be used to model

8|
mass transport of non-charged particles within porous electrodes filled with
aqueous electrolyte.
The Electrophoretic Transport ( ) interface can be used to investigate the
transport of weak acids, bases, and ampholytes in aqueous solvents. The physics
interface is typically used to model various electrophoresis modes, such as zone
electrophoresis, isothachophoresis, isoelectric focusing, and moving boundary
electrophoresis, but is applicable to any aqueous system involving multiple
acid-base equilibria.
The Surface Reactions interface ( ) can be used model reactions and translateral
transport of surface (adsorbed) species.

P OROUS M EDIA F LOW

The Fluid Flow interfaces ( ) can be combined with the Electrochemistry
Module interfaces to model free and forced convection in electrochemical cells.
The Darcy’s Law interface ( ) is used to model fluid movement through
interstices in a porous medium where a homogenization of the porous and fluid
media into a single medium is done. Together with the continuity equation and
equation of state for the pore fluid (or gas) this physics interface can be used to
model low velocity flows, for which the pressure gradient is the major driving
force. The penetration of electrolyte through a porous electrode is a classic
example for the use of Darcy’s Law in electrochemical engineering.
Darcy’s law can be used in porous media where the fluid is mostly influenced by
the frictional resistance within the pores. Its use is within very low flows, or media
where the porosity is very small. Where the size of the interstices is larger, and the
fluid is also influenced by itself, the gradients in fluid velocity, pressure, and gravity
must be considered. This is done in the Brinkman Equations interface.
The Brinkman Equations interface ( ) is used to model compressible flow at
speeds of less than Mach 0.3, but control over the density and any of the mass
balances that are deployed must be maintained to help with this.
The Free and Porous Media Flow interface ( ) is useful for equipment that
contain domains where free flow is connected to porous media, such as concrete
structures immersed in water.

H EAT T RANSFER
The Heat Transfer interfaces ( ) have ready-made formulations for the
contribution of Joule heating, and other electrochemical heat sources, to the
thermal balance of electrochemical cells.
The features of the various Heat Transfer interfaces include Heat Transfer in
Fluids ( ), Heat Transfer in Solids ( ), and Heat Transfer in

|9
Porous Media ( ), and account for conductive and convective heat transfer.
These features interact seamlessly and can be used in combination in a single
model.
The detailed equations and assumptions that are defined by the physics interfaces
are formulated in the Electrochemistry Module User’s Guide and the COMSOL
Multiphysics Reference Manual.

Physics Interface Guide by Space Dimension and Study Type

The table lists the physics interfaces available with this module in addition to those
included with the COMSOL basic license.

PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE

DIMENSION

Chemical Species Transport

Surface Reactions sr all dimensions stationary (3D, 2D, and 2D

axisymmetric models only);
time dependent
Transport of Diluted tds all dimensions stationary; time dependent
Species
Transport of Diluted tds all dimensions stationary; time dependent
Species in Porous Media
Electrophoretic Transport el all dimensions stationary; stationary with
initialization; time
dependent; time dependent
with initialization
Nernst-Planck-Poisson tds+es all dimensions stationary; time dependent;
Equations stationary source sweep;
small-signal analysis,
frequency domain
Reacting Flow

Laminar Flow, Diluted — 3D, 2D, 2D stationary; time dependent

Species axisymmetric

10 |
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION

Electrochemistry

Primary Current cd all dimensions stationary; stationary with

Distribution initialization; time
dependent; time dependent
Secondary Current with initialization; AC
Distribution impedance, initial values; AC
impedance, stationary; AC
impedance, time dependent
Tertiary Current tcd all dimensions stationary; stationary with
Distribution, Nernst-Planck initialization; time
(Electroneutrality, dependent; time dependent
Water-Based with with initialization; AC
Electroneutrality, impedance, initial values; AC
Supporting Electrolyte) impedance, stationary; AC
impedance, time dependent
Electroanalysis elan all dimensions stationary; time dependent;
AC impedance, initial values;
AC impedance, stationary;
AC impedance, time
dependent; cyclic
voltammetry
Electrode, Shell els 3D, 2D, 2D stationary; time dependent
axisymmetric
Fluid Flow

Porous Media and Subsurface Flow

Brinkman Equations br 3D, 2D, 2D stationary; time dependent

axisymmetric
Darcy’s Law dl all dimensions stationary; time dependent

Free and Porous Media fp 3D, 2D, 2D stationary; time dependent

Flow axisymmetric

| 11
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION

Heat Transfer

Heat Transfer in Porous ht all dimensions stationary; time dependent

Media

12 |
Orange Battery

This tutorial example serves as an introduction to electrochemistry modeling in

COMSOL Multiphysics and models the currents and the concentration of
dissolved metal ions in a battery (corrosion cell) made from an orange and two
metal nails.

Metal nails

Orange

Figure 4: Modeled geometry. Orange and two metal nails.

This type of battery is commonly used in chemistry class demonstrations. Instead

of an orange, lemons or potatoes can also be used.

Model Definition
The citric acid and various other ions in the orange serve as electrolyte, and using
nails of different metals as electrodes creates a galvanic potential over the cell.
In this example a zinc nail is used as one of the electrodes, giving rise to the
following electrode reaction:
2+ -
Zn ( s ) → Zn + 2e E eq, 0 = – 0.76 V

The other nail consists of copper, and here hydrogen evolution is assumed to take
place:
+ -
2H + 2e → H 2 ( g ) E eq, 0 = 0 V

| 13
Eeq,0 above denotes the equilibrium potentials at standard conditions versus a
standard hydrogen electrode (SHE). In the model, the equilibrium potentials are
corrected for the pH and zinc concentration of the orange pulp using the Nernst
equation.
The model for the currents in the orange and electrodes is set up using the
Secondary Current Distribution interface. The electrolyte current in the orange is
thereby solved for by Ohms law. The conductivity of the metal nails is so high that
the electrode domains are not included in the model, instead boundary conditions
on the nail surfaces are used to set the nail potentials. One nail is grounded and
the other one is set to a cell potential to comply with a total current condition.
This would correspond to a situation where the cell is controlled galvanostatically,
for instance by the use of a potentiostat.
Butler-Volmer type expressions, with concentration dependent exchange current
density for the zinc reaction, are used for the electrode kinetics on the surface of
the nails within the orange.
The initial values electrolyte potential is set to correspond to the potential of a cell
at open circuit (that is, no activation potential). Following the definition of the
overpotential:
η = φ s – φ l – E eq

the initial value becomes:

φ l, init = φ s – E eq – η = 0 – E eq, Zn – 0 = – E eq, Zn

In an extension of the model, the diffusion and migration of the dissolved zinc
ions in the orange from the zinc electrode reaction is modeled by the Transport
of Diluted Species interface in a time-dependent simulation. This assumes that the
zinc ion transport can be described by the Nernst-Planck equations (without an
electroneutrality condition due to a supporting electrolyte). In addition, the zinc
electrode kinetics are modified to be dependent on the zinc concentration, which
increases in the orange as more and more zinc is dissolved. The zinc concentration
is set to 0.001 mol/m3 at the start of the simulation. All boundaries except the
zinc electrode are insulated.

R ESULTS AND D ISCUSSION

Figure 5 shows the potential field in the orange. The potential decreases as the
current flows from the zinc electrode (left) to the upper electrode (right). The
main part of the cell voltage loss is due to ohmic losses in the electrolyte.

14 |
The performance of the battery could probably be increased by using an
electrolyte of higher conductivity (for example, a lemon instead of an orange) or
by decreasing the distance between the nails.

Figure 5: Potential field in the electrolyte at t=0.

Figure 6 shows a polarization plot as the total current of the battery increases from
0 to 1 mA. The large change in cell voltage seen at low currents is due to

| 15
overpotential losses at the zinc electrode. Increasing the area of the zinc electrode
would decrease this effect.

Figure 6: Polarization plot for the initial concentrations.

16 |
Figure 7 shows an isosurface for the 0.2 mol/m3 concentration level of zinc ions
after running the battery for five minutes.

Figure 7: 0.2 mol/m3 zinc concentration isosurface after five minutes.

| 17
Figure 8 shows how far the 0.2 mol/m3 isosurface level has reached after one hour.

Figure 8: 0.2 mol/m3 zinc concentration isosurface after one hour.

18 |
Figure 9 shows how the cell current evolves with time. Due to the increase of zinc
ions at the zinc nail electrode, the battery current decreases for a constant cell
current.

Figure 9: Cell current vs. time.

S UGGESTED E XERCISES AND E XTENSIONS OF THE M ODEL

Change the radius of the nails and the value of the electrolyte conductivity and
investigate how this affects the polarization plot.
The dissolved zinc ions may form a layer of zinc hydroxide on the zinc surface,
giving rise to an additional potential drop. You may use the Film Resistance
section on the Electrode Surface Boundary node to include this potential drop.
The value of the film resistance could for instance be a function of the zinc ion
concentration variable in the pulp. Alternatively you may add a Surface Reactions
interface to model the buildup of the surface concentration of zinc hydroxide, and
then couple the film resistance to the hydroxide concentration.

Model Wizard

Note: These instructions are for the user interface on Windows but apply, with
minor differences, also to Linux and Mac.

| 19
1 To start the software, double-click the COMSOL icon on the desktop. When
the software opens, you can choose to use the Model Wizard to create a new
COMSOL model or Blank Model to create one manually. For this tutorial, click
the Model Wizard button.
If COMSOL is already open, you can start the Model Wizard by selecting
New from the File menu and then click Model Wizard .
The Model Wizard guides you through the first steps of setting up a model. The
next window lets you select the dimension of the modeling space.
2 In the Space Dimension window
click the 3D button .
3 In the Select Physics tree under
Electrochemistry>Primary and
Secondary Current Distribution,
click Secondary Current
Distribution .
4 Click Add then click the Study
button.
5 In the tree under General Studies,
click Stationary .
6 Click Done , when the Added
study and Added physics interfaces
lists look as to the right.

Geometry I

Start by drawing the geometry; one

sphere (the orange) and two cylinders
(the metal nails, zinc and copper).
Both shapes are available in the
Geometry toolbar.

20 |
Sphere I
1 On the Geometry toolbar, click Sphere .
2 In the Settings window for Sphere, locate the Size section and type 5e-2 in the
Radius text field.

Zinc nail
1 On the Geometry toolbar, select
Cylinder .
2 In the Settings window for Cylinder, type
Zinc nail in the Label text field.
3 Locate the Size and Shape section.
- In the Radius text field, type 2e-3.
- In the Height text field, type 5e-2.
4 Locate the Position section.
- In the x text field, type -2e-2.
- In the z text field, type 2e-2.
5 Locate the Selections of Resulting
Entities section. Select the Resulting
objects selection check box and choose
All levels in the Show in physics list. By
enabling this setting you can easily select
all boundaries of the nail later on when
setting up the physics.
Duplicate the cylinder and change the x
position to draw the second nail.

| 21
Copper nail
1 Right-click Component
1>Geometry 1>Zinc nail
and click Duplicate .
2 In the Settings window for
Cylinder, type Copper nail
in the Label text field.
3 In the Position section enter
2e-2 in the x text field.

4 Finally, to inspect in the

Graphics window that the
geometry becomes identical
to that in Figure 4, click
Build all on the
Geometry toolbar and
Zoom Extents on the
Graphics toolbar.

Global Definitions

Load the model parameters from a text file.

The location of the files used in this exercise may vary based on the installation.
For example, if the installation is on your hard drive, the file path might be similar
to C:\Program Files\COMSOL54\applications\.

Parameters
1 In the Model Builder under Global Definitions, click Parameters 1 .
Note: On Linux and Mac, the Home toolbar refers to the specific set of controls
near the top of the Desktop.
2 In the Settings window for Parameters click the Load from File button .

22 |
3 Browse to the application library folder, Electrochemistry Module>Tutorials,
and double-click the file orange_battery_parameters.txt.

Definitions

Add a manual explicit selection for the orange domain. This will facilitate if the
geometry needs to be changed in the future. Selections is generally a convenient
way to group different parts of the geometry together. Selections of the two nails
are already available under the Geometry node in the Model Builder Window.

Explicit I
1 On the Definitions toolbar, click Explicit .

2 In the Settings window for Explicit, type

Orange in the Label text field.
3 Select Domain 1 only.
Load some variable definitions from a file.
The location of the file should be the same
as the parameter file.

Variables I
1 In the Model Builder window, under
Component 1 right-click Definitions

| 23
 and select Variables . Alternatively, click the Definitions toolbar and select
Local Variables .
2 In the Settings window for Variables click the Load from File button .
3 Browse to the application’s Application Library folder and double-click the file
orange_battery_variables.txt.

Secondary Current Distribution

Now, start setting up the current distribution model. Change the selection of the
entire physics interface to the orange domain only.
1 In the Model Builder window, click Secondary Current Distribution.
2 In the Settings window for Secondary Current Distribution , locate the
Domain Selection section and choose Orange from the Selection list.

Electrolyte I
The selection is locked to all selected domains of the physics interface, which in
this case is the orange only.
1 In the Model Builder window, under Component 1>Secondary Current
Distribution click Electrolyte 1 . The ‘D’ in the upper left corner of a node
means it is a default node.
2 In the Settings window for Electrolyte,
locate the Electrolyte section. Here the
electrolyte conductivity can be set. From
the σl list choose User defined and type
sigma in the associated text field. sigma
is defined in the Parameters node.

Electrode Surface I
Use Electrode Surface nodes to define both a metal electrode potential and an
electrode electrolyte interface. Use a Bulter-Volmer expression for the zinc

24 |
electrode. The hydrogen kinetics are assumed to be very fast so that a linearized
Butler-Volmer expression is applicable. Start with the Zinc nail (electrode).
1 On the Physics toolbar, click Boundaries and choose Electrode Surface .

2 In the Settings window for Electrode Surface locate the Boundary Selection
section and select Zinc nail from the Boundary Selection.
3 Click the Transparency button on the Graphics toolbar to enable
transparency and inspect the active boundaries in the Graphics window.

Electrode Reaction I
1 In the Model Builder, expand the Electrode Surface 1 node and click Electrode
Reaction 1 .

| 25
2 In the Settings window for Electrode
Reaction, locate the Model Inputs
section and enter T in the T text field.
3 In the Equilibrium Potential section type
E_eq_Zn, i.e. the equilibrium potential
of zinc metal, in the Eeq text field.
4 In the Electrode Kinetics section the
reaction kinetics are set for the zinc
oxidation. From the Kinetics expression
type list, choose ‘Butler-Volmer’ kinetics
that fits this reaction well.
5 In the Exchange current density i0 text
field, type i0_Zn.
6 This reaction kinetics requires the anodic
and cathodic transfer coefficients to be
set. Do this with the entries alpha_a_Zn
and alpha_c_Zn in the αa and αc text
fields, respectively.

Electrode Surface 2
Define the copper nail (electrode) in a similar way with an Electrode Surface.
1 On the Physics toolbar, Click Boundaries and choose Electrode Surface
2 In the Settings window for the Electrode Surface, locate the Boundary
Selection section. In the Selection list there, choose Copper nail.
3 This cell is under galvanic control.
Specify this in the Boundary Condition
section by selecting Total current in the
Boundary condition list and typing
-i_app in the Il,total text field that
appears. This parameter will also be used
in the study to perform a galvanic
polarization sweep over the cell. .

26 |
Electrode Reaction 1
1 In the Model Builder window, expand
the Electrode Surface 2 and click
Electrode Reaction 1.
2 In the Settings window for Electrode
Reaction, locate the Model Inputs
section and in the T text field, type T.
3 Locate the Equilibrium Potential
section. In the Eeq text field, type
E_eq_H2, i.e. the equilibrium potential
for hydrogen evolution since protons
reduce at the copper surface.
4 Type i0_H2 in the exchange current
density i0 text field available in the
Electrode Kinetics section. The
hydrogen kinetics are assumed to be very
fast so that a linearized Butler-Volmer
expression is applicable here.

Initial Values I
We are using non-linear kinetics in the model. Provide an initial value for the
electrolyte potential in order to reduce solver time and improve convergence. As
a rule of thumb one can often use the negative of the equilibrium potential of the
grounded electrode as initial value for the electrolyte potential.
1 In the Model Builder window, under Secondary Current Distribution click
Initial Values 1 .
2 In the Settings window for Initial Values, locate the Initial Values section.
3 As a rule of thumb one can often use the negative of the equilibrium potential
of the grounded electrode as initial value for the electrolyte potential.
Therefore, in the phil text field, type -E_eq_Zn0.

| 27
The interface node sequence in the Model Builder should match this figure:

Me sh 1

The default mesh settings can be used in this model and therefore no changes need
to be made here.

Study 1

The model is now ready for solving.

1 On the Study toolbar click Compute .

Results

Use a 3D isosurface plot for visualizing the potential field in the electrolyte.
1 On the Home toolbar, click Add
Plot Group and choose 3D Plot
Group .
2 In the Settings window for the 3D
Plot Group, type Potential
Isosurface in the Label text field.

Potential Isosurface
1 Right-click the Potential
Isosurface node and choose Isosurface .

28 |
2 In the Settings window for Isosurface,
locate the Levels section and in the Total
levels text field, type 25.
3 Click Zoom Extents and
Transparency on the Graphics
toolbar before clicking the Plot button
on the Potential Isosurface toolbar.
A plot identical to Figure 5 will now
appear in the Graphics window.

Study 1

Now use an Auxiliary Sweep available in Stationary study step to solve over a range
of cell currents in order to create a polarization plot.
1 In the Model Builder window, under Study 1 click Step 1: Stationary .
2 In the Settings window for Stationary, click to expand the Study extension
section and select the Auxiliary sweep check box.
3 Click the Add button under the table. In the Parameter name column, select
i_app (Applied Current) from the list that appears.
4 In the Parameter value column enter 0 10^range(-5,0.2,-3) and in the
Parameter unit column enter A.

5 On the Study toolbar click Compute .

Results

The results from this last run can be used to set up a polarization plot:

| 29
1 On the Home toolbar, click Add Plot Group and choose 1D Plot
Group .
2 In the Settings window for 1D Plot Group, type Polarization Plot in the
Label text field.

Polarization Plot
1 On the Polarization Plot toolbar, click Global .
2 In the Settings window for Global, locate the y-Axis Data section and click Add
Expression . Select Component 1>Secondary Current Distribution>Electric
potential (siec.phis_eebii2) from the list that appears.

3 On the Polarization plot toolbar, click Plot . This generates Figure 6 above.

Adding Mass Transport

Move on to extend the model to investigate the battery voltage over time at a
certain load current. Start by adding a physics interface to handle the mass
transport of zinc ions. The ions can be regarded as diluted in the orange interior
making the Transport of Diluted Species interface suitable for modeling the mass
transport.

30 |
1 On the Home toolbar, click Add
physics . This opens the Add Physics
window.
2 Here, locate the Transport of Diluted
Species (tds) under Chemical Species
Transport and click Add to
component .
3 To close the Add physics window, click
Add Physics on the Home toolbar.

Transport of Diluted Species

1 In the Model Builder window, under

Component 1 click Transport of Diluted
Species .

2 In the Settings window for Transport of

Diluted Species, locate the Domain
Selection section and in the Selection list
choose Orange.
3 The orange pulp is quiescent so
convection can be ignored in the model.
This is done in the Transport
Mechanisms section by clearing the
Convection check box. Additionally, due
to the presence of a lot of other ions in
the pulp, acting as supporting
electrolyte, we assume the potential
gradients to be small and hence also
ignore the effect of migrative transport
of the zinc ions.

Electrode Surface Coupling I

For this example we will use the default diffusion coefficient value in the Transport
Properties node so no settings are needed on the default domain node. The

| 31
following steps couple the electrochemical reaction currents to the ion flux at the
electrode surface.
1 On the Physics toolbar, click Boundaries and choose Electrode Surface
Coupling in which the coupling can be enabled.
2 In the Settings window for Electrode Surface Coupling, locate the Boundary
Selection section and in the Selection list select Zinc nail.

Reaction Coefficients I
1 In the Model Builder
window, expand the
Electrode Surface Coupling 1
node and click the Reaction
Coefficients 1 node.
2 In the Model Inputs section,
the Coupled reaction iloc list
setting should be that of the
Zinc nail. Therefore, choose
Local current density (siec/
eebiil/er1).
3 In the Stoichiometric
Coefficients section, enter In
the n text field type 2, and in
the νc text field type -1. These
stoichiometric coefficients refer to the stoichiometric number of the reacting
species when written as a reduction reaction.

Initial Values I
The initial zinc ion concentration at the start of the time-dependent simulation is
required.
1 In the Model Builder window, under Component 1 and Transport of Diluted
Species click Initial Values 1 .
2 Locate the Initial Values section and enter c_Zn20 in the c text field.

32 |
The node sequence of the interfaces in the Model Builder should now match this
figure:

Definitions

The zinc ion concentration is no longer constant. Modify the equilibrium

potential and exchange current density variables for the zinc electrode to be
dependent on the local concentration of zinc ions. This can be done by defining
the following variable.

Variables I
1 In the Model Builder window, under Component 1 and Definitions click
Variables 1 .
2 In the table, add the following settings:

Name Expression Unit Description

c_Zn2 c mol/m³ Zinc ion concentration
(used in above
expressions)

The name of this variable is 'c' by default in the Transport of Diluted Species
interface.

| 33
Add Study

Create a new time-dependent study for the

concentration simulation.
1 On the Home toolbar, click Add
Study to open the Add Study window.
2 In the Add Study window select Time
Dependent with Initialization under
Preset Studies for Selected Physics
Interfaces . This study type computes
the necessary initial values facilitating the
computations to start.
3 To close the Add Study window, click
Add Study on the Home toolbar.

Study 2

The Time Dependent with Initialization study consists of two steps.

Step I: Current Distribution Initialization

In Step I, use a secondary current distribution initialization step which works well
when concentration dependent electrode kinetics are modeled.
1 In the Model Builder, under Study
2 click Step 1: Current
Distribution Initialization .
2 From the Current distribution type
list available in the Study Settings section,
choose Secondary.

34 |
Step 2: Time Dependent
1 Now click Step 2: Time Dependent Under Study 2 to select the times that
are solved for.
2 In the Study Settings section enter range(0,60,3600) in the Times text field.
3 On the Home toolbar, click Compute .

Results

Plot the zinc concentration in the orange.

3D Plot Group 14
1 On the Home toolbar, click Add Plot Group and choose 3D Plot
Group .
2 In the Settings window for 3D Plot Group , type Concentration
Isosurface in the Label text field.
3 Locate the Data section. From the Data set list in the Data section, choose
Study 2/Solution 2. From the Time (s) list that appears, select 300.

Concentration Isosurface
1 Right-click Concentration Isosurface under Results and choose
Isosurface . Alternatively, this can be done by clicking Isosurface on the
Concentration Isosurface toolbar.
2 Click Replace Expression in the upper-right corner of the Expression
section in the Isosurface Settings window. From the menu, choose
Component 1>Transport of Diluted Species>c - Concentration -mol/m^3.
3 In the Levels section, from the Entry
method list, choose Levels and type 0.2
in the Levels text field.
4 After you have clicked the Zoom
Extents button on the Graphics
toolbar, click Plot on the
Concentration Isosurface toolbar. This
will create Figure 7 showing Zinc
concentration isosurface.

Export
An animation can be used to visualize the zinc ion isosurface during the simulated
time:

| 35
1 On the Results toolbar, click Animation and select
Player . This generates an Animation 1 node
under the Export node.

2 In the Settings window for Animation,

locate the Target section and select
Player from the Target list. Locate the
Scene section. From the Subject list,
choose Concentration Isosurface.
3 In the Animation Editing section,
choose From list in the Time selection
list.
4 At t=0 there is no concentration gradient
in the orange, so deselect the first
time-step. To do this, select all steps
except 0 in the Times (s) list.
5 Click Zoom Extents in the Graphics
toolbar followed by right-click of
Animation 1 under the Export
node and choose Play . The last frame
at 3600 s should match Figure 8.

36 |
F o r f u r th e r s i m u l a t i o n s . .

If you want to experiment with the model by changing parameter values and
simulate new polarization plots, you have to disable the Transport of Diluted
Species interface in the stationary study as follows:
1 In the Model Builder window under Study 1 click Step 1: Stationary .
2 Disable the Transport of Diluted Species interface in the table shown in the
Physics and Variables Selection section by clearing the check box in the ‘Solve
for’ column.

| 37
Cyclic Voltammetry at a Macroelectrode in 1D

Cyclic voltammetry is a common analytical technique for investigating

electrochemical systems. In this method, the potential difference between a
working electrode and a reference electrode is swept linearly in time from a start
potential to a vertex potential, and back again (see Figure 10). The resulting
current at the working electrode is recorded and is plotted against the applied
electrode potential in a voltammogram.

Figure 10: Potential of the working electrode during one voltammetric cycle. The potential is cycled
between the vertex potentials 0.4 V and –0.4 V. The scan rate is 1 mV/s.

Voltammetry is a valuable technique because information about both the

electrochemical reactivity and the transport properties of a system can be extracted
simultaneously. For quantitative interpretation of voltammetry, however, we must
use numerical methods to solve the physical equations that describe voltammetry.
Then, unknown physical quantities in the system can be inferred by ‘fitting’ to
experimental data.
This example demonstrates the use of a common approximation in which a large
electrode (macroelectrode) is assumed to have uniform transport behavior across
its surface, so only physics occurring normal to the surface need to be considered.
By simplifying the model to 1D, an efficient time-dependent analysis is possible.
In this model, a Parametric Sweep is used to compare voltammetry recorded at
different voltammetric scan rates.

38 |
Model Definition
The model contains a single 1D domain of length L, which is the maximum extent
of the diffusion layer over the duration of the voltammetry experiment. A
conservative setting for L is set to greatly exceed the mean diffusion layer
thickness:
L = 6 Dt max

Here, D is the diffusion coefficient of the reactant and tmax is the duration of the
cyclic voltammogram.

DOMAIN EQUATIONS
We assume the presence of a large quantity of supporting electrolyte. This is inert
salt that is added in electroanalytical experiments to increase the conductivity of
the electrolyte without otherwise interfering with the reaction chemistry. Under
these conditions, the resistance of the solution is sufficiently low that the electric
field is negligible, and we can assume φl = 0.
The Electroanalysis interface implements chemical transport equations for the
reactant and product species of the redox couple subject to this assumption. The
domain equation is the diffusion equation (also known as Fick’s second law) to
describe the chemical transport of the electroactive species A and B:
∂c i
------- = ∇ ⋅ ( D i ∇c i )
∂t

BOUNDARY EQUATIONS
At the bulk boundary (x = L), we assume a uniform concentration equal to the
bulk concentration for the reactant. The product has zero concentration here, as
in bulk.
At the electrode boundary (x = 0), the reactant species A oxidizes (loses one
electron) to form the product B. By convention, electrochemical reactions are
written in the reduction direction:
–
B+e ↔A

The stoichiometric coefficient is –1 for B, the “reactant” in the reduction

direction, and +1 for A, the “product” in the reduction direction. This
formulation is consistent even in examples such as this model where at certain
applied potentials, the reaction proceeds favorably to convert A to B. The number
of electrons transferred, n, equals one.

| 39
The current density for this reaction is given by the electroanalytical
Butler-Volmer equation for an oxidation reaction:
( n – α c )Fη – α c Fη
i loc = nFk 0  c A exp  ---------------------------- – c B exp  ----------------- 
  RT   RT  

in which k0 is the heterogeneous rate constant of the reaction, αc is the cathodic

transfer coefficient, and η is the overpotential at the working electrode. This
overpotential is the difference between the applied potential and the equilibrium
potential (formal reduction potential) of the redox couple of species A and B.
According to Faraday’s laws of electrolysis, the flux of the reactant and product
species are proportional to the current density drawn:
ν i i loc
– n ⋅ N i = -------------
nF

This is expressed in the Electrode Surface boundary condition in the

Electroanalysis Module.
The applied triangular waveform for the cyclic voltammetry study is specified in
the Electrode Surface boundary condition according to two vertex potentials—
forming a potential window between –0.4 V and +0.4 V, either side of the
equilibrium reduction potential—and a voltammetric scan rate, v (SI unit: V/s),
which is the rate at which the applied potential is changed.
In the 1D approximation, the total current is related to the current density simply
by multiplying by the electrode area A:
I el = i loc A

CYCLIC VOLTAMMETRY STUDY

In the cyclic voltammetry experiment, the potential applied to the working
electrode surface is varied linearly as a function of time. A Parametric Sweep is used
to compare the voltammetry recorded at different scan rates.

40 |
Results and Discussion
The shape of the cyclic voltammogram (Figure 11) shows the relation between
electrode kinetics and chemical species transport (diffusion).

Figure 11: Cyclic voltammetry recorded at a macroelectrode.

Initially, at reducing potentials, the oxidation reaction is not driven and negligible
current is drawn. As the potential is ramped conversion towards the reduction
potential of the redox couple, the oxidation reaction is accelerated and the current
increases. Once the oxidation reaction has consumed the reactant at the electrode
surface, the current becomes limited by the rate of transport of A towards the
working electrode. Therefore, a peak current is observed, and at higher potentials,
the voltammetric current falls off at a potential-independent rate; this region is
termed “diffusion-controlled” or “transport-controlled”.
On sweeping back towards more reducing potentials, the conversion of the
product B into the original reactant A gives a negative (cathodic, reduction)
current. Depletion of the reacting species B causes a negative peak current and
conversion thereafter proceeds at a diffusion-controlled rate.
The magnitude of the current on the forward peak, Ipf, is a common diagnostic
variable in voltammetry. For fast electrode kinetics and at a macroelectrode under
the 1D approximation, its value is given theoretically by the Randles–Ševcík equation (see
Ref. 1 and Ref. 2 for a detailed discussion and derivation):

nF
I pf = 0.446 nFAc -------- Dv
RT

| 41
where A is the electrode area, c is the bulk concentration of the reactant, and D is the
diffusion coefficient of the reactant.
The square-root relationship between peak current and scan rate is characteristic
of macroelectrode cyclic voltammetry under the above conditions.

R EFERENCES
1 R.G. Compton and C.E. Banks, Understanding Voltammetry, 2nd ed.,
London, 2011.
2 A.J. Bard and L.R. Faulkner, Electrochemical Methods, Fundamentals and
Applications, 2nd ed., Wiley, New York, 2001.

Model Wizard

1 If COMSOL is already open, you

can start the Model Wizard by
selecting New from the File
menu and then click Model Wizard
.
2 In the Space Dimension Window
click the 1D button .
3 In the Select Physics tree under
Electrochemistry, click
Electroanalysis (elan) .
4 Click Add.
5 This model will solve for the two
concentrations of a redox couple,
change the default concentration
variable names to cA and cB in the
Concentrations table.
6 Click the Study button.
7 In the tree under Preset Studies for
Selected Physics Interfaces, select
Cyclic Voltammetry .
8 Click Done when the Added
study and Added physics interfaces
lists look as to the right.

42 |
Global Definitions

Load the model parameters from a text file.

The location of the files used in this exercise may vary based on the installation.
For example, if the installation is on your hard drive, the file path might be similar
to C:\Program Files\COMSOL54\applications\.

Parameters
1 In the Model Builder under Global Definitions, click Parameters 1 .
Note: On Linux and Mac, the Home toolbar refers to the specific set of controls
near the top of the Desktop.
2 In the Settings window for Parameters click the Load from File button .
3 Browse to the application library folder, Electrochemistry Module>Tutorials,
and double-click the file cyclic_voltammetry_1d_parameters.txt.
The figure shows what the parameter list should look like after loading the text
file. For any model, you can create your own parameter lists and store these in text
files to reuse them in your models. To create your own text file, click the Save to
File button .

G e o me t r y I

Build the model geometry as a single interval, where the left boundary will later
be defined as the electrode surface, and the right boundary will be the boundary
towards the bulk.

| 43
Interval 1
1 On the Geometry toolbar, click
Interval .
In the Settings window for Interval type L in
the Right endpoint text field. This sets the
interval length to L. The value of L is defined
in the Parameters node. Note that L will vary
with the Voltammetric scan rate parameter v,
which is also defined in the Parameters node.
This is done to ensure the length of the
electrolyte domain is always greater than
diffusion length.

2 On the Geometry toolbar, click Build all . The completed geometry is now
built and appears in the Graphics window:

Ele c tr oana lysis

Start defining the physics by setting the diffusion coefficients for the two species
of the redox couple.

44 |
Transport Properties 1
1 In the Model Builder window, under Component
1>Electroanalysis, Click Transport Properties
. The ‘D’ in the upper left corner of a node
means it is a default node.
2 Locate the Diffusion section. In the DcA and DcB text fields, type DA and DB,
respectively.

Concentration 1
Set the boundary to the right to bulk concentration values.
1 On the Physics toolbar, click Boundaries and choose Concentration .

2 Select Boundary 2 only.

3 In the Concentration section of the Concentration Settings window, select the
Species cA and Species cB check boxes.
4 In the c0,cA text field, type c_bulk.

Electrode Surface 1
Set up the electrode surface and the cyclic voltammetry settings on the left
boundary.

| 45
1 On the Physics toolbar, click Boundaries
and click Electrode Surface .
2 Select Boundary 1 only.
3 In the Boundary section of the Electrode
Surface Settings window, choose Cyclic
voltammetry from the Boundary
condition list.
4 In the Linear sweep rate, type v.
5 The voltage will be cycled between the
vertex potentials. When the start
potential is not specified, the sweep will
start at Vertex potential 2. Enter the
following:
- E_vertex1, in the Vertex potential 1 text
field.
- E_vertex2, in the Vertex potential 2 text
field.

Electrode Reaction 1
Specify the electrode reaction as an Electroanalytical Butler-Volmer reaction,
which is concentration dependent as defined by the stoichiometric coefficients in
the reaction.

46 |
1 In the Model Builder window, expand
the Electrode Surface 1 node and click
Electrode Reaction 1.
2 Under Model Inputs, type T in the
Temperature text field.
3 Locate the Electrode Kinetics section
and in the k0 text field, type k0.
4 Locate the Stoichiometric Coefficients
section and in the νcA and νcB text
fields, type 1 and -1, respectively.

Double Layer Capacitance

1 In the Model Builder window,
right-click Electrode Surface 1 and
choose the boundary condition Double
Layer Capacitance.
2 Click the Double Layer Capacitance
node and type Cdl in the Cdl text field.

Initial Values 1
Continue with specifying the initial concentration values. This will set the
concentration values when the simulation starts at t=0.
1 In the Model Builder window under Component 1>Electroanalysis, click Initial
Values 1 .
2 In the Settings window for Initial Values, locate the Initial Values section.
3 In the cA text field, type c_bulk-cB0*(1-x/L).
4 In the cB text field, type cB0*(1-x/L).
The node sequence in the Model Builder under the Electroanalysis interface
should match this figure:

| 47
Me sh 1

The default mesh settings can be used in this model and therefore no changes need
to be made here.

Study 1

The problem is now ready for solving. Use a parametric sweep to solve the
problem for various sweep rates.

Parametric Sweep
1 On the Study toolbar, click Parametric Sweep .
2 Under Study Settings in Parametric Sweep Settings window, click the Add
button under the table. In the Parameter name column, select
v (Voltammetric scan rate) from the list that appears.
3 In the Parameter value column enter (10^range(-3,1,0))[V/s].

4 On the Study toolbar click Compute .

Results

A number of plots are created by default. The first default plot group shows
voltammograms that are automatically created together with the Cyclic
Voltammetry feature in the Electrode Surface node.

48 |
Cyclic Voltammograms
1 On the Cyclic Voltammograms (elan) toolbar, click Plot to see the simulated
voltammograms in the Graphics window.

2 If the Legend box covers the graph, go to the Settings window for the 1D Plot
Group denoted Cyclic Voltammograms (elan) and click to expand the
Legend section. Switch for instance to Upper left in the Position list. Doing so
and clicking Plot will reproduce Figure 11.

| 49
50 |