Soluble Lead-Acid Redox Flow Battery

In a redox flow battery electrochemical energy is stored as redox couples in the electrolyte, which is stored in tanks outside the electrochemical cell. During operation, electrolyte is pumped through the cell and, due to the electrochemical reactions, the individual concentrations of the active species in the electrolyte are changed.

The state of charge of the flow battery is determined by the electrolyte species concentrations, the total flowing electrolyte volume in the system (tank + pump + hoses + cell), and possibly also by the concentration of solid species on the electrodes. Depending on the cell chemistry the cell can have separated or combined anode and cathode compartments and electrolyte tanks.

Figure 1: Working principle of the soluble lead acid flow battery.

In the soluble lead acid flow battery one electrolyte solution is used. The active component in the electrolyte is the lead ion that reacts on the electrodes to form solid lead (negative electrode) or lead oxide (positive electrode). The electrode chemistry is similar to a traditional lead-acid battery, with the difference that solid lead sulfonate is not formed in the electrodes.

This example simulates a soluble lead-acid flow battery during an applied charge-discharge load cycle. The surface chemistry of the positive electrode is modeled by using two different lead oxides and two different positive electrode reactions in the model.

Model Definition

Cell Geometry and mesh

The electrochemical cell consist of two flat 10 cm square electrodes, placed in parallel with a 12 mm gap in between. The aspect ratio of the cell motivates modeling the cell in 2D. The cell geometry and mesh is shown in Figure 2.

Figure 2: Geometry and mesh of the electrochemical cell.

Due to the very high electrical conductivity of the electrodes, the potential gradients in the electrodes are neglected, and the electrodes are not included in the geometry.

To handle possible edge effects in the electrolyte, 1 mm regions are added at the inlet and outlet, outside the active electrode region.

A mapped rectangular mesh is used, and boundary meshing is used to resolve the steep gradients in the electrolyte close to the electrode surfaces.

Electrolyte Mass and current transport equations

The electrolyte is based on a mixture of lead methane sulfonate, methane sulfonic acid and water, which in this model is assumed to dissociate into an electrolyte consisting of Pb2+, H+, HSO4--ions dissolved in a bulk solution of zero-charged species (mainly water). Electroneutrality is assumed locally in the electrolyte. The combination of these assumptions allow for the use of Tertiary Current Distribution, Nernst-Planck interface for modeling the electrolyte transport.

The electric potential in the electrodes is assumed to be space independent. The negative electrode is grounded. On the positive electrode, an electrode potential is calculated in order to fulfill a current density condition defined by the load cycle (using the Electrode Surface boundary node).

A load cycle of 1 h charge, 20 s rest, 1 h discharge, 20 s is applied twice to the cell. During charge or discharge a constant current density corresponding to a mean current density in the cell 200 A/m2 is applied.

The species fluxes are defined on the electrode surfaces according to the electrode reactions below. An Inflow condition is used at the inlet with the inlet concentrations (

and

) taken from the tank model described below. An Outflow condition is set at the outlet. All other boundaries are isolated.

Negative Electrode Reaction

On the negative electrode the following electrode reaction occurs:

with the kinetics being described by a Butler-Volmer expression:

Where

+ is a rate constant and

is the concentration of lead ions in the electrolyte.

As reference electrode we use the negative electrode at reference conditions. The equilibrium potential for the negative electrode is assumed to follow the Nernst equation according to:

Positive Electrode Main Reaction

The positive electrode main reaction is:

with the kinetics being described by a Butler-Volmer expression:

where

is a rate constant,

is the electrolyte proton concentration and

is the proton reference concentration in the electrolyte at equilibrium.

The positive main reaction has the following equilibrium potential, described by the Nernst Equation:

(1)

Positive Electrode Side Reaction

Multiple types of lead oxides may form on the positive electrode. In this model the following side reaction is investigated:

where the electrode is kinetics is described by

(2)

where the overpotential, η, is the same as for the positive electrode main reaction (Equation 1). (The deviation of the equilibrium potential of the side reaction versus the positive main reaction equilibrium potential is controlled by the rate parameters.)

In Equation 2

and

are rate constants, and

and

are the surface concentration of the lead oxides (mol/m2).

Tank model

The electrolyte flowing out from the cell flows into the tank, undergoes mixing, and is then led into the cell again on the inlet side.

Assuming good mixing in the tank the inlet concentrations,

and

, are governed by the following ODEs:

Where V is the total volume of flowing electrolyte in the tank, and L is the height of the electrodes. (

and

denote the molar fluxes of the respective electrolyte species in the normal direction to the boundary).

The two ODEs are modeled using an ODEs and DAEs interface.

Fluid flow Equations

The fluid is led into the cell at a velocity Vin of 2.3 cm/s. The relevant Reynolds number for the flow between the plates is:

where the parameter values for water are used for the density ρ, 1000 kg/m3, and viscosity μ, 10-3 Pa·s. We can assume that the flow is in the laminar regime (Re<2000), and hence the Laminar Flow interface is used to model the fluid flow.

Vin is applied at the inlet, a pressure condition is applied to the outlet, and no slip conditions are applied to the electrode surfaces and channel walls. The induced convection at the electrode surfaces due to the electrochemical reactions is assumed to be negligible. In this way the flow model is stationary and only solved for once. The convective flow is used as a model input to the Tertiary Current Distribution, Nernst-Planck interface.

Surface concentrations on the positive electrode

Two different lead oxides, PbO and PbO2, may be formed on the positive electrodes due to the electrochemical reactions. The surface concentrations of these two species,

and

(SI unit: mol/m2), are modeled using the Dissolving-Depositing Species section of the Tertiary Current Distribution, Nernst-Planck interface.

Results and Discussion

Figure 3 shows the flow field and pressure drop for the cell. The parabolic velocity profile is expected for this rectangular geometry (Poiseuille flow).

Figure 3: Velocity field and pressure.

Figure 4 shows the cell voltage during the load cycling. The first charge cycle voltages differs from the second.

Figure 4: Cell potential versus time.

Figure 5 shows the average surface concentrations of PbO and PbO2 at the positive electrode during the load cycle. The build-up of PbO2 starts during the first charge cycle, whereas there is only small amounts of PbO formed until the beginning of the discharge cycle. The presence of PbO alters the kinetics of the positive electrode during the second charge cycle. Figure 7 shows the difference in local current densities between the different parts of the load cycle. The modified kinetics on the positive electrode impacts the overpotentials, which in turn explains the difference in cell voltages during the first and second charge cycles in Figure 4.

Figure 5: Average surface concentrations of the two different lead oxides on the positive electrode.

Figure 6 shows the inlet concentrations of lead ions and protons from the tank model during the load cycling. The lead ion concentrations is lower and the proton concentration is higher at the beginning of the second charge cycle, compared to the initial values. The reasons for these variations are due to the lead ion consumption to form the lead oxide layer, and the proton release from water molecules in the same process.

Figure 6: Inlet lead ion and proton concentrations versus time.

Figure 7: Local current densities of the electrode reactions. The PbO2 and PbO reactions occurring on the positive electrode, Pb on the negative.

Figure 8 and Figure 9 depict the Pb2+ concentration distribution in the electrolyte at the end of the first charge and discharge step, respectively. Large gradients are present in the boundary layer close to the electrode surfaces.

Figure 8: Lead ion concentration in the electrolyte at the end of the first charging cycle.

Figure 9: Lead ion concentration in the electrolyte at the end of the discharge cycle.

Figure 10 and Figure 11 plot the proton concentration distributions in the electrolyte at the end of the first charge and discharge step. Also for this species large gradients are present in the boundary layer close to the electrode surfaces.

Figure 10: Hydrogen ion concentration in the electrolyte at the end of the first charging cycle.

Figure 11: Hydrogen ion concentration at the end of the discharge cycle.

References

1. R. Willis, J. Collins, D. Stratton-Campbell, C. Low, D. Pletcher, and F. Walsh, “Developments in the Soluble Lead-acid Flow Battery,” J. Appl. Electrochem, vol. 40, pp. 955–965, 2010.

2. A. Shah, R. Wills, and F. Walsh, “A Mathematical Model for the Soluble Lead-Acid Flow Battery,” J. Electrochemical Society, vol. 157, pp. A589–A599, 2010.

3. J. Newman and K. Thomas-Alyea, Electrochemical Systems, p. 284, Table 11.1, John Wiley & Sons, 2004.

Battery_Design_Module/Flow_Batteries/pb_flow_battery

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 12[mm].

In the text field, type 10[cm].

Click

Global Definitions

Parameters 1

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 pb_flow_battery_parameters.txt.

Definitions

Inlet

In the toolbar, click

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

Locate the section. From the list, choose .

Select Boundary 2 only.

Outlet

In the toolbar, click

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

Locate the section. From the list, choose .

Select Boundary 3 only.

Positive electrode

In the toolbar, click

In the window for , type Positive electrode in the text field.

Locate the section. From the list, choose .

Select Boundary 4 only.

Negative electrode

In the toolbar, click

In the window for , type Negative electrode in the text field.

Locate the section. From the list, choose .

Select Boundary 1 only.

Materials

Use the material parameter values for water from the model library.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Laminar Flow (spf)

Inlet 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the Uav text field, type U_in.

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the check box.

Mesh 1

The rectangular geometry makes a mapped mesh suitable for this problem.

Edge 1

In the toolbar, click

Select Boundaries 1 and 4 only.

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type 30.

In the text field, type 5.

Edge 2

In the toolbar, click

Select Boundaries 2 and 3 only.

Distribution 1

Right-click and choose .

Mapped 1

In the toolbar, click

Boundary Layers 1

In the toolbar, click

In the window for , click to expand the section.

Clear the check box.

Boundary Layer Properties

In the window, click .

Select Boundaries 1 and 4 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 5e-5.

In the window, right-click and choose .

In the window for , click

Study 1

In the toolbar, click

Results

Arrow Surface 1

In the window, right-click and choose .

In the window for , locate the section.

Find the subsection. In the text field, type 5.

Locate the section. From the list, choose .

In the toolbar, click

Flow

In the window, under click .

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

Component 1 (comp1)

Now add the electrochemistry to the model, start by adding the appropriate physics interface.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click to expand the section. Click

In the table, enter the following settings:


cPbII
cH
cHSO4

Find the subsection. In the table, clear the check box for .

Click in the window toolbar.

In the tree, select .

In the table, clear the check box for .

Click in the window toolbar.

In the toolbar, click

to close the window.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Definitions

Rectangle 1 (rect1)

Use a number of rectangle functions to set up the charge/discharge cycle that will be a function of time.

In the toolbar, click

and choose .

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

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

In the text field, type t_charge.

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

Rectangle 2 (rect2)

In the toolbar, click

and choose .

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

Locate the section. In the text field, type t_charge+t_rest.

In the text field, type t_charge+t_rest+t_discharge.

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

Rectangle 3 (rect3)

In the toolbar, click

and choose .

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

Locate the section. In the text field, type t_charge+t_rest+t_discharge+t_rest.

In the text field, type 2*t_charge+t_rest+t_discharge+t_rest.

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

Variables 1

Now the defined analytical functions can be used to set up variables on the negative and positive electrodes. Load the variables from a text file.

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 pb_flow_battery_variables.txt.

(The tcd.phisext variable is the electric potential on the boundary.)

Tertiary Current Distribution, Nernst-Planck (tcd)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Electrolyte 1

In the window, under click .

In the window for , locate the section.

From the u list, choose .

Locate the section. In the DcPbII text field, type D_PbII.

In the DcH text field, type D_H.

In the DcHSO4 text field, type D_HSO4.

Locate the section. In the zcPbII text field, type 2.

In the zcH text field, type 1.

In the zcHSO4 text field, type -1.

Inflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the c0,cPbII text field, type cPbII_in.

In the c0,cH text field, type cH_in.

Locate the section. From the list, choose .

Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Negative Electrode

In the toolbar, click

and choose .

In the window for , type Negative Electrode in the text field.

Locate the section. From the list, choose .

Electrode Reaction 1

In the window, expand the node, then click .

In the window for , locate the section.

In the n text field, type 2.

In the νcPbII text field, type -1.

Locate the section. From the list, choose .

In the i0,ref(T) text field, type i0ref_neg.

In the table, enter the following settings:


Electrolyte species	γ? (1)
cPbII	1

In the αa text field, type 1.

In the αc text field, type 1.

Positive Electrode

In the toolbar, click

and choose .

In the window for , type Positive Electrode in the text field.

Locate the section. From the list, choose .

In the il,average text field, type i_cycle.

In the φs,ext,init text field, type E0_pos.

Electrode Reaction 1

In the window, expand the node, then click .

In the window for , locate the section.

In the n text field, type 2.

In the νcPbII text field, type 1.

In the νcH text field, type -4.

Locate the section. In the Eeq,ref(T) text field, type Eeq_pos.

Click to expand the section. In the table, enter the following settings:


Electrolyte species	Reference concentrations (mol/m^3)
cH	0.5[M]

Locate the section. From the list, choose .

In the i0,ref(T) text field, type i0ref_pos.

In the table, enter the following settings:


Electrolyte species	γ? (1)
cPbII	1
cH	1

Positive Electrode

In the window, click .

Electrode Reaction 2

In the toolbar, click

and choose .

Main Reaction

In the window, under click .

In the window for , type Main Reaction in the text field.

Side Reaction

In the window, under click .

In the window for , type Side Reaction in the text field.

Locate the section. In the n text field, type 2.

In the νcH text field, type -2.

Locate the section. From the Eeq list, choose . Locate the section. From the iloc,expr list, choose . In the associated text field, type i_PbO.

Positive Electrode

In the window, click .

In the window for , click to expand the section.

Click

In the table, enter the following settings:


Species	Density (kg/m^3)	Molar mass (kg/mol)
PbO2	9.38[g/cm^3]	0.2392
PbO	9.53[g/cm^3]	0.2232

Main Reaction

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Species	Stoichiometric coefficient (1)
PbO2	-1
PbO	0

Side Reaction

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Species	Stoichiometric coefficient (1)
PbO2	-1
PbO	1

Initial Values 1

In the window, click .

In the window for , locate the section.

In the cPbII text field, type c0_PbII.

In the cH text field, type c0_H.

In the phil text field, type -E0_neg.

Global Definitions

Default Model Inputs

Set up the temperature value used in the entire model.

In the window, under click .

In the window for , locate the section.

In the tree, select .

Find the subsection. In the text field, type T.

Definitions

The tank model is based on two ODEs and the integrals of the ion fluxes over the electrode boundaries.

Integration 1 (intop1)

In the toolbar, click

and choose .

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

Locate the section. From the list, choose .

From the list, choose .

Integration 2 (intop2)

In the toolbar, click

and choose .

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

Locate the section. From the list, choose .

From the list, choose .

Global ODEs and DAEs (ge)

Global Equations 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (u_t0) (1/s)	Description
cPbII_in	cPbII_int-L/V*(comp1.int_outlet(comp1.tcd.tflux_cPbIIy)-comp1.int_inlet(comp1.tcd.tflux_cPbIIy))	c0_PbII	0	Pb ion inlet concentration
cH_in	cH_int-L/V*(comp1.int_outlet(comp1.tcd.tflux_cHy)-comp1.int_inlet(comp1.tcd.tflux_cHy))	c0_H	0	Proton inlet concentration