MassTransportPDEComponent
MassTransferValue
MassConcentrationCondition
NDSolveValue

Chemistry
Partial Differential Equations
Finite Element Method

Numeric Cyclic Voltammetry

Solve a coupled reaction model

Open in Cloud
Download Notebook

Define a coupled reaction model:

masstransportmodel ︷ ∂ c A (t,x) ∂t +∇·(- d A ∇ c A (t,x))	=	masstransferboundary ︷  Γ x=0 k c c A (t,x)- k a c B (t,x)
masstransportmodel ︷ ∂ c B (t,x) ∂t +∇·(- d B ∇ c B (t,x))	=	masstransferboundary ︷  Γ x=0 k a c B (t,x)- k c c A (t,x)

The underlying reaction model is given by

+e⇌

and Butler–Volmer kinetics.

Set up the mass transport model variables

vars

In[1]:=

vars={{cA[t,x],cB[t,x]},t,{x}};

Set up a region

In[2]:=

Ω=Line[{{0},{1/40}}];

Specify material parameters:

In[3]:=

cAbulk=1;dA=10^-5;dB=10^-5;

Specify the electrochemical rate constants:

In[4]:=

k0=10^-5;α=1/2;β=1/2;ef0=0;rtbyf=25.7*10^-3;kc[t_]:=k0Exp[-α/rtbyf(e[t]-ef0)]ka[t_]:=k0Exp[β/rtbyf(e[t]-ef0)]

Specify the potential:

In[5]:=

ts=1;ν=-1;e1=1/2;e[t_]:=Piecewise[{{e1+νt,0≤t≤ts},{e1+2νts-νt,ts≤t≤2ts}}]

Specify mass transport model parameters, mass diffusivity

and

. The concentration of

at the far end is set to the bulk concentration and the concentration of

is set to 0:

In[6]:=

pars="DiffusionCoefficient"{{dA,0},{0,dB}},"BoundaryCondition1""AmbientConcentration"{ka[t]/kc[t]cB[t,x],kc[t]/ka[t]cA[t,x]},"MassTransferCoefficient"{kc[t],ka[t]},"BoundaryCondition2""MassConcentration"{cAbulk,0};

Set up the model:

In[7]:=

eqn={MassTransportPDEComponent[vars,pars]MassTransferValue[x0,vars,pars,"BoundaryCondition1"],MassConcentrationCondition[x1/40,vars,pars,"BoundaryCondition2"]}

Out[7]=



∇

{x}

·-

100000

.

∇

{x}

cA[t,x]+

(1,0)

[t,x],

∇

{x}

·-

100000

.

∇

{x}

cB[t,x]+

(1,0)

[t,x]NeumannValue0.00001

-19.4553

0.5-1.t	0.≤t≤1.
-1.5+t	1.≤t≤2.
0.	True

2.71828

-1.cA[t,x]+

38.9105

0.5-1.t	0.≤t≤1.
-1.5+t	1.≤t≤2.
0.	True

2.71828

cB[t,x],x0,NeumannValue0.00001

19.4553

0.5-1.t	0.≤t≤1.
-1.5+t	1.≤t≤2.
0.	True

2.71828

-38.9105

0.5-1.t	0.≤t≤1.
-1.5+t	1.≤t≤2.
0.	True

2.71828

cA[t,x]-1.cB[t,x],x0,DirichletCondition{cA[t,x]1,cB[t,x]0},x



Set up initial conditions such that only

is in bulk solution:

In[8]:=

ics={cA[0,x]cAbulk,cB[0,x]0};

Solve the model:

In[9]:=

tmax=2ts;{cAfun,cBfun}=NDSolveValue[{eqn,ics},{cA,cB},{t,0,tmax},{x}∈Ω]

Out[9]=

InterpolatingFunction

Domain: {{0.,2.},{0.,0.025}}

Output: scalar

,InterpolatingFunction

Domain: {{0.,2.},{0.,0.025}}

Output: scalar



Visualize the concentrations at various points in time:

In[10]:=

ManipulatePlot{cAfun[t,x],cBfun[t,x]},{x,0,1/20},

,

Rule[

]



Out[10]=

1.9

	cA
	cB

Visualize the cyclic voltammogram at various points in time:

In[11]:=

ManipulateParametricPlot{-(e[t]-ef0),kc[t](cAfun[t,0]-(ka[t]cBfun[t,0])/kc[t])},{t,0,tmax},

,

Rule[

]



Out[11]=

Related Symbols

MassTransportPDEComponent
MassTransferValue
MassConcentrationCondition
NDSolveValue

Publisher Information

Contributed by: Wolfram Staff

Wolfram Language Example Repository

Numeric Cyclic Voltammetry

See Also

Publisher Information