# Mass Transport Model

I've asked similar questions before about Mathematica's Mass Transport model. My aim is to model these systems and show how they change by manipulating various parameters.

This time it's the following system.

Edit:

The reaction that the system is modeling and the equilibrium constants are given below (My apologies for not uploading them from the start but my question was predominantly about those boundary conditions):

End of Edit

The system above should yield a voltammogram like this:

I tried implementing the model using the following code (excluding the plotting of results).

Needs["NDSolveFEM"]
ClearAll["Global*"]
(*Experimental Parameters*)
k1 := 0; k2 := 0 (*10^8*);
ef0AB := 0; ef0BC = -0.4;
α1 := 0.5; α2 := 0.5;
k10 := 1; k20 := 1;
ar := 1; cAbulk := 10^-3;
dA := 10^-5; dB = 10^-5; dC := 10^-5;
rtbyf := 25.7 10^-3(*volt*);
f := 96485.33;

ts := 1; tmax = 2 ts; ν := -1; e1 := -0.3; ef0 := 0;
e[t_] := Piecewise[{{e1 + ν t,
0 <= t <= ts}, {e1 + 2 ν ts - ν t, ts <= t <= 2 ts}}]

large = 6 Sqrt[dA tmax];

i[t_, x_] := f*ar ( D[2 dA *cA[t, x] + dB cB[t, x]]) /. x -> 0

vars = {{cA[t, x], cB[t, x], cC[t, x]}, t, {x}};
pars = <|
"DiffusionCoefficient" -> {{dA, 0, 0}, {0, dB, 0}, {0, 0, dC}},
"MassReactionRate" -> {{Subscript[k, 2] cC[t, x], 0, 0}, {0,
2 Subscript[k, 1] cB[t, x], 0}, {0, 0,
Subscript[k, 2] cA[t, x]}},
"MassSource" -> {{Subscript[k, 1] cB[t, x]^2}, {2 Subscript[k, 2]
cA[t, x] cC[t, x]}, {Subscript[k, 1] cB[t, x]^2}},

"BoundaryConditionMassFlux" ->
<|"MassFlux" -> {D[-dB cB[t, x] - dC cC[t, x], x] ,
D[-dA cA[t, x] - dC cC[t, x], x],
D[-dA cA[t, x] - dB cB[t, x], x]}|>,

"BoundaryConditionConcentration" ->
<|"MassConcentration" -> {cB[t, x] Exp[rtbyf^-1 (e[t] - ef0AB)],
cA[t, x] /Exp[rtbyf^-1 (e[t] - ef0AB)],
cA[t, x] /Exp[rtbyf^-1 (e[t] - ef0BC)]}|>,

"BoundaryConditionInf" -> <|"MassConcentration" -> {cAbulk, 0, 0}|>|>;

ops = MassTransportPDEComponent[vars, pars];

ics = {cA[0, x] == cAbulk, cB[0, x] == 0, cC[0, x] == 0};
Γflux =
MassFluxValue[x == 0, vars, pars, "BoundaryConditionMassFlux"];
Γcond =
MassConcentrationCondition[x == 0, vars, pars,
"BoundaryConditionConcentration"];
Γcondinf =
MassConcentrationCondition[x == large, vars, pars,
"BoundaryConditionInf"];

{cAfun, cBfun, cCfun} =
NDSolveValue[{ops == Γflux, Γcond, \
Γcondinf, ics}, {cA, cB, cC}, {t, 0, tmax}, {x, 0,
large}, Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> large/1000}}}];


I get two errors; one of which says:

NDSolveValue::fembcdepderiv: Derivatives of dependent variables in boundary conditions are not supported with the Finite Element Method in this version of NDSolve.


The other one says that the lists are not the same shape which again has me confused because NDSolveValue should return a list with three elements.

I tried to test it with a different model by removing the derivatives but then it returned similar errors with DirichletCondition. So I think I'm doing something wrong here.

Thank you to everyone in advance.

• 1. "NDSolveValue should return a list with three elements." No, NDSolveValue already fails after the first warning, what's returned is an unevaluated NDSolveValue[…]. Just execute the NDSolveValue[…] separately and observe. 2. As indicated by the description of NDSolveValue::fembcdepderiv, you can not have derivatives in NeumannValue. (Please observe what's inside Γflux.) Do notice MassTransportPDEComponent, etc. are no more than generator of PDEs and b.c.s, in other words, if NeumannValue isn't able to do something, MassTransportPDEComponent, etc. won't help either. – xzczd Mar 29 at 5:34
• 3. As I've said for several times, if the domain is always regular, consider using the old good TensorProductGrid instead of FiniteElement. – xzczd Mar 29 at 5:43
• Can you share where this model comes from? I see a few issue to be tackled. Everything seems clear, except the last BC. There you have a relation of Dirichlet values - that might work if Dirichlet cross coupling were implemented but even then I am not sure this is what is needed. Additionally you have a constraint that the sum of Neumann values is 0. This might be doable with an integral constraint. All in all the FEM version is not doable from an NDSolve level right now. It might be doable with the low level FEM code but would requite considerable work. Do you want to go there? – user21 Apr 2 at 6:55
• Yes I'll be happy to share them. My next target was to solve problems involving second order reactions. The models come from Understanding Voltammetry by R.G. Compton. This one is page 145 (chapter 4). I don't want anything too complicated just yet. But if you can point to some basic literature and examples regarding it, that would be helpful. – Walser Apr 4 at 10:57
• Thanks for the reference. Unfortunately. I do not have any literature to point you to, but this post has an integral constraint. But again, I am very unsure if this would work. This would need some time to tackle. – user21 Apr 5 at 6:41

OK. Let me extend my comments to an answer. First of all, I'd like to point out why OP's attempt doesn't work:

1. "NDSolveValue should return a list with three elements." No, NDSolveValue already fails after the first warning, what's returned is an unevaluated NDSolveValue[…]. If you still feel confused, execute the NDSolveValue[…] separately and observe.

2. As indicated by the description of NDSolveValue::fembcdepderiv, you can not have derivatives in NeumannValue. (Please observe what's inside Γflux.) Do notice MassTransportPDEComponent, etc. are no more than generator of PDEs and b.c.s, in other words, if NeumannValue isn't able to do something, MassTransportPDEComponent, etc. won't help either.

3. Units of parameters should be converted to SI units.

4. e1 should be 0.3. Notice the label of horizonal axis is $$\color{red}{-}\text E/\text V$$.

5. There's a typo in the b.c. at $$x=0$$, $$\frac{[\text A]_{x=0}}{[\text C]_{x=0}}$$ should be $$\frac{[\color{red}{\text B}]_{x=0}}{[\text C]_{x=0}}$$.

6. Another typo: $$-\text {I}/\text {mA}$$ should be $$\text {I}/\text {mA}$$.

Then, how to solve? It's not obvious to me how one can set the b.c. at $$x=0$$ when FiniteElement method is chosen for spatial discretization, so I'll turn to the old good TensorProductGrid. With this method, the b.c. at $$x=0$$ can be imposed straightforwardly, no extra coding is needed.

efab = 0; efbc = -4/10; DA = DB = DC = 10^-5 10^-4; k1 = 0; k2 = 0;
cAbulk = 10^-3 10^3; A = 1 10^-4; F = 9648533/100; FbyRT = F/(8314/1000 298);
ts = 1 ; tmax = 2 ts; ν = -1; e1 = 0.3;
Ε[t_] = Piecewise[{{e1 + ν t, t <= ts}}, e1 + 2 ν ts - ν t];
inf = 2/10 10^-3;
With[{cA = cA[t, x], cB = cB[t, x], cC = cC[t, x]},
react = k1 cB^2 - k2 cA cC;
eq = {
D[cA, t] == DA D[cA, x, x] + react,
D[cB, t] == DB D[cB, x, x] - 2 react,
D[cC, t] == DC D[cC, x, x] + react};
ic = {cA == cAbulk, cB == 0, cC == 0} /. t -> 0;
bcinf = {cA == cAbulk, cB == 0, cC == 0} /. x -> inf;
bc0 = {DA D[cA, x] + DB D[cB, x] + DC D[cC, x] == 0,
cA == Exp[FbyRT (Ε[t] - efab)] cB,
cB == Exp[FbyRT (Ε[t] - efbc)] cC} /. x -> 0;
i = F A (2 DA D[cA, x] + DB D[cB, x]) /. x -> 0];

mol[tf : False | True, sf_ : Automatic] := {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}
sollst = NDSolveValue[{eq, ic, bc0, bcinf}, {cA, cB, cC}, {t, 0, tmax}, {x, 0, inf},
Method -> mol[True, 10^3]];


ibcinc warning will pop up, but don't worry, because it has already been taken into account by adjusting DifferentiateBoundaryConditions sub-option. You may check this post for more details. Let's check if we've succeeded to reproduce the plots in the text book:

time = {0.3, 0.7, 1., 1.3, 1.7, 2.0};
label = CharacterRange["A", "F"];
ilst = i /. Thread[{cA, cB, cC} -> sollst];

ParametricPlot[{-Ε[t], 10^3 ilst}, {t, 0, tmax},
AspectRatio -> 1/GoldenRatio, PlotRange -> All, AxesLabel -> {"-E/V", "I/mA"}]~Show~
ListPlot[Association@Thread[label -> Table[{-Ε[t], 10^3 ilst}, {t, time}]],
PlotStyle -> Black]


GraphicsGrid@Partition[#, 3] &@
Table[Plot[(sollst[time[[i]], 10^-3 x] // Through)/cAbulk// Evaluate, {x, 0, inf 10^3},
PlotRange -> All, PlotLabel -> label[[i]],
AxesLabel -> {"x/mm", "[A]/[\!$$\*SubscriptBox[\(A$$, $$bulk$$]\)]"}], {i, 6}]


As one can see, the results agree well with those in the screenshot.

• Great! That shows what need to be fixed in the equations to get this to work. Well done. – user21 Apr 7 at 7:11
• It is very nice solution (+1) even it is not FEM based solution. But we have message NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. – Alex Trounev Apr 7 at 9:06
• @Alex The ibcinc warning isn't a problem here because we've already taken it into consideration by adjusting DifferentiateBoundaryConditions option. Please refer to the linked post for more details. – xzczd Apr 7 at 10:51
• Nice investigative work. I do have a question though. Numerics aside doesn't the solution $c_{A}[ t,x] =c_{A,bulk} ;c_{B}[ t,x] =c_{C}[ t,x] =0$ satisfy the equations, initial conditions, and boundary conditions? – Tim Laska Apr 7 at 23:06
• @Tim It's against $\frac{[\text A]_{x=0}}{[\text B]_{x=0}}=\exp \left[ \frac{F}{R T}\left\{E-E_f^0(\text A/\text B)\right\}\right]$ :) . – xzczd Apr 8 at 2:14

Thanks to @xzczd answer we can reproduce solution with FEM. First we should note that $$[A]+[B]+[C]=[A]_{bulk}$$ in a case of equal $$D_A=D_B=D_C$$, therefore we can exclude equation for $$[A]$$ and resolve bc at $$x=0$$ as follows

ss = Solve[{bulk - cC - cB ==
Exp[FbyRT (\[CapitalEpsilon][t] - efab)] cB,
cB == Exp[FbyRT (\[CapitalEpsilon][t] - efbc)] cC}, {cB, cC}]

{cB0, cC0} = {cB, cC} /. ss[[1]]


Then we use code similar to xzczd as

efab = 0; efbc = -4/10; DA = DB = DC = 10^-5 10^-4; k1 = 0; k2 = 0;
cAbulk = 10^-3 10^3; A = 1 10^-4; F = 9648533/100; FbyRT =
F/(8314/1000 298);
ts = 1; tmax = 2 ts; \[Nu] = -1; e1 = 0.3;
\[CapitalEpsilon][t_] =
Piecewise[{{e1 + \[Nu] t, t <= ts}}, e1 + 2 \[Nu] ts - \[Nu] t];
inf = 2/10 10^-3;
With[{cA = cAbulk - cB[t, x] - cC[t, x], cB = cB[t, x],

cC = cC[t, x]}, react = k1 cB^2 - k2 cA cC;
eq = {D[cB, t] == DB D[cB, x, x] - 2 react,
D[cC, t] == DC D[cC, x, x] + react};
ic = {cB == 0, cC == 0} /. t -> 0;
bcinf = DirichletCondition[{cB == 0, cC == 0}, x == inf];
bc0 = DirichletCondition[{cB == cB0, cC == cC0}, x == 0];
i = F A (2 DA D[cA, x] + DB D[cB, x]) /. x -> 0];


FEM solution

Needs["NDSolveFEM"]

large = inf; {cBfun, cCfun} =
NDSolveValue[{eq, ic, bc0, bcinf} /. bulk -> cAbulk, {cB, cC}, {t, 0,
tmax}, {x, 0, large},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> large/1000}}}]


Visualization

time = {0.3, 0.7, 1., 1.3, 1.7, 2.0};
label = CharacterRange["A", "F"];
ilst = i /. Thread[{cB, cC} -> {cBfun, cCfun}];

ParametricPlot[{-\[CapitalEpsilon][t], 10^3 ilst}, {t, 0, tmax},
AspectRatio -> 1/GoldenRatio, PlotRange -> All,
AxesLabel -> {"-E/V", "I/mA"}]~Show~
ListPlot[Association@
Table[{-\[CapitalEpsilon][t], 10^3 ilst}, {t, time}]],
PlotStyle -> Black]

Table[Plot[{cAbulk - cBfun[time[[i]], 10^-3 x] -
cCfun[time[[i]], 10^-3 x], cBfun[time[[i]], 10^-3 x],
cCfun[time[[i]], 10^-3 x]}/cAbulk // Evaluate, {x, 0, inf 10^3},
PlotRange -> All, PlotLabel -> label[[i]],
AxesLabel -> {"x/mm", ""}], {i, 6}]


Now we can test Mass Transport model as follows

ClearAll["Global*"]

Needs["NDSolveFEM"]

(*Experimental Parameters*)
efab = 0; efbc = -4/10; DA = DB = DC = 10^-5 10^-4; k1 = 0; k2 = 0;
cAbulk = 10^-3 10^3; A = 1 10^-4; F = 9648533/100; FbyRT =
F/(8314/1000 298);
ts = 1; tmax = 2 ts; \[Nu] = -1; e1 = 0.3;
\[CapitalEpsilon][t_] =
Piecewise[{{e1 + \[Nu] t, t <= ts}}, e1 + 2 \[Nu] ts - \[Nu] t];
inf = 2/10 10^-3; large = inf;
cA[t_, x_] := cAbulk - cB[t, x] - cC[t, x];
ss = Solve[{bulk - cC - cB ==
Exp[FbyRT (\[CapitalEpsilon][t] - efab)] cB,
cB == Exp[FbyRT (\[CapitalEpsilon][t] - efbc)] cC}, {cB, cC}];
{cB0, cC0} = {cB, cC} /. ss[[1]];
vars = {{cB[t, x], cC[t, x]}, t, {x}};
pars = <|"DiffusionCoefficient" -> {{DB, 0}, {0, DC}},
"MassReactionRate" -> {{2 k1 cB[t, x], 0}, {0, k2 cA[t, x]}},
"MassSource" -> {{2 k2 cA[t, x] cC[t, x]}, {k1 cB[t, x]^2}},
"BoundaryConditionConcentration" -> <|
"MassConcentration" -> {cB0, cC0}|>,
"BoundaryConditionInf" -> <|"MassConcentration" -> {0, 0}|>|>;

ops = MassTransportPDEComponent[vars, pars];

ics = {cB[0, x] == 0, cC[0, x] == 0};
\[CapitalGamma]flux =
MassFluxValue[x == 0, vars, pars, "BoundaryConditionMassFlux"];
\[CapitalGamma]cond =
MassConcentrationCondition[x == 0, vars, pars,
"BoundaryConditionConcentration"];
\[CapitalGamma]condinf =
MassConcentrationCondition[x == large, vars, pars,
"BoundaryConditionInf"];

{cBfun, cCfun} =
NDSolveValue[{ops == {0,
0}, \[CapitalGamma]cond, \[CapitalGamma]condinf, ics} /.
bulk -> cAbulk, {cB, cC}, {t, 0, tmax}, {x, 0, large},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> large/1000}}},
MaxStepSize -> 0.001];


Visualization

time = {0.3, 0.7, 1., 1.3, 1.7, 2.0};
label = CharacterRange["A",
"F"]; With[{cA = cAbulk - cB[t, x] - cC[t, x], cB = cB[t, x],
cC = cC[t, x]}, j = F A (2 DA D[cA, x] + DB D[cB, x]) /. x -> 0];
ilst = j /. Thread[{cB, cC} -> {cBfun, cCfun}];

ParametricPlot[{-\[CapitalEpsilon][t], 10^3 ilst}, {t, 0, tmax},
AspectRatio -> 1/GoldenRatio, PlotRange -> All,
AxesLabel -> {"-E/V", "I/mA"}]~Show~
ListPlot[Association@
Table[{-\[CapitalEpsilon][t], 10^3 ilst}, {t, time}]],
PlotStyle -> Black]

Table[Plot[{cAbulk - cBfun[time[[i]], 10^-3 x] -
cCfun[time[[i]], 10^-3 x], cBfun[time[[i]], 10^-3 x],
cCfun[time[[i]], 10^-3 x]}/cAbulk // Evaluate, {x, 0, inf 10^3},
PlotRange -> All, PlotLabel -> label[[i]],
AxesLabel -> {"x/mm", ""}], {i, 6}]


Finally we can test nonlinear Mass Transport Model with k1 = 10^3; k2 = 5 10^2;. This solution is differ from above and also takes time to compute data

• With MaxStepSize -> 0.001 the first plot will be smooth :) . – xzczd Apr 10 at 2:50
• @xzczd Thank you! It is only some preparation for the Mass Transport model testing. ;) – Alex Trounev Apr 10 at 10:25

My suspicion is that the equations from the textbook are misleading for a FEM solution. In my brief reading on cyclic voltammetry (e.g. here or here as I do not have access to the book by Compton), the electron transfers occur at the surface of the electrode and not in the bulk. Therefore, the $$\color{Red}{Red}$$ terms in the following equations should be removed and recast as boundary conditions.

$$\begin{array}{l} \frac{{\partial [A]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[A]}}{{\partial {x^2}}} {\color{Red}{ + {k_1}{[B]^2} - {k_2}[A][C]}}\\ \frac{{\partial [B]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[B]}}{{\partial {x^2}}} {\color{Red}{ - 2{k_1}{[B]^2} + 2{k_2}[A][C]}}\\ \frac{{\partial [C]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[C]}}{{\partial {x^2}}} {\color{Red}{ + {k_1}{[B]^2} - {k_2}[A][C]}} \end{array}$$

I also suspect the surface reactions are incorrect. Species [A], initially, is the only non-zero concentration. Therefore, all the surface reactions would be zero according to this mechanism. Perhaps, the author meant something like shown in $$\color{Green}{Green}$$:

$$\begin{array}{l} \frac{{\partial [A]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[A]}}{{\partial {x^2}}} {\color{Green}{ - {k_1}{[A]^2} + {k_2}[B][C]}}\\ \frac{{\partial [B]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[B]}}{{\partial {x^2}}} {\color{Green}{ + {k_1}{[A]^2} - {k_2}[B][C]}}\\ \frac{{\partial [C]}}{{\partial t}} = {D_A}\frac{{{\partial ^2}[C]}}{{\partial {x^2}}} {\color{Green}{ - {k_2}[B][C]}} \end{array}$$

Now, a reaction could proceed with only [A] as a starting reagent.

This implies that my answer to your previous question on cyclic voltammetry is a better starting point since it has the same domain equations.

For brevity, I will build the recast Neumann values by hand. I do believe, however, that this would be a good opportunity to think about the future enhancements for the finite element method to include surface reaction boundary conditions since they can come in many flavors depending on the limiting conditions that are assumed. Also, as reaction networks become more complex, it is quite easy to include excess reactions that are not linearly independent. So, it is easy to create ill-posed FEM models. This might be considered to be a stoichiometry preprocessing step that needs to be validated before the construction of the FEM model.

As @xzczd pointed out there are many terms that are undefined or never used. The screen capture also shows $$k_1=0$$, which I assume is unintended. So, I am going to modify my previous answer that seems to have the terms better defined. The key is to get the reaction rate constants $$k_1[T]$$ and $$k_2[T]$$ as we did for $$k_c[T]$$ and $$k_a[T]$$ in the previous problem. Since I do not quite know how to do this with the information provided, I will simply reuse the previous problem's definitions.

(*Experimental Parameters*)
cAbulk := 1*10^-3;
k0 := 1;
rtbyf := 25.7 10^-3(*volt*);
dA := 10^-5; dB := 10^-5; dC := 10^-5
α := 0.5; β := 0.5; T = 298; ef0 = 0;
ts := 1; tmax = 2 ts; ν := -1; e1 := 0.5;
large = 6 Sqrt[dA tmax];
e[t_] := Piecewise[{{e1 + ν t,
0 <= t <= ts}, {e1 + 2 ν ts - ν t, ts <= t <= 2 ts}}]
k1[t_] := k0 Exp[-α/rtbyf (e[t] - ef0)]
k2[t_] := k0 Exp[β/rtbyf (e[t] - ef0)]
f := 96485.33; ar := 1;
(*Current*)
i[t_, x_] := f*ar (D[2 dA*cA[t, x] + dB cB[t, x]]) /. x -> 0

(*PDE set up*)
vars = {{cA[t, x], cB[t, x], cC[t, x]}, t, {x}};
pars = <|"DiffusionCoefficient" -> {{dA, 0, 0}, {0, dB, 0}, {0, 0,
dC}}, "BoundaryCondition1" -> <|
"MassConcentration" -> {cAbulk, 0, 0}|>|>;
ops = MassTransportPDEComponent[vars, pars];
Γflux = {NeumannValue[(k1[t] cB[t, x]^2 -
k2[t] cA[t, x] cC[t, x]), x == 0],
NeumannValue[-2 (k1[t] cB[t, x]^2 - k2[t] cA[t, x] cC[t, x]),
x == 0],
NeumannValue[(k1[t] cB[t, x]^2 - k2[t] cA[t, x] cC[t, x]), x == 0]};
Γcondinf =
MassConcentrationCondition[x == large, vars, pars,
"BoundaryCondition1"];
ics = {cA[0, x] == cAbulk, cB[0, x] == 0, cC[0, x] == 0};
(*Solve PDE*)
{cAfun, cBfun, cCfun} =
NDSolveValue[{ops == Γflux, Γcondinf,
ics}, {cA, cB, cC}, {t, 0, tmax}, {x, 0, large},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> large/1000}}}];
Plot3D[{cAfun[t, x], cBfun[t, x], cCfun[t, x]}, {t, 0, tmax}, {x, 0,
large}, PlotRange -> All]


The PDE system solves, but we see that the results are rather uninteresting because of what I stated above that there are no driving forces due to the initial conditions.

We can solve the alternative mechanism in $$\color{Green}{Green}$$. You should note, that we run into stability issues at about 1 s. I would probably review the mechanism carefully before trying to stabilize NDSolve.

Γflux = {NeumannValue[-(k1[t] cA[t, x]^2 -
k2[t] cB[t, x] cC[t, x]), x == 0],
NeumannValue[(k1[t] cA[t, x]^2 - k2[t] cB[t, x] cC[t, x]), x == 0],
NeumannValue[(k2[t] cB[t, x] cC[t, x]), x == 0]};
{cAfun, cBfun, cCfun} =
NDSolveValue[{ops == Γflux, Γcondinf,
ics}, {cA, cB, cC}, {t, 0, tmax}, {x, 0, large},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> large/1000}}},
MaxStepSize -> 0.0001];
Plot3D[{cAfun[t, x], cBfun[t, x], cCfun[t, x]}, {t, 0, tmax}, {x, 0,
large}, PlotRange -> All]


This PDE system solves, so it is a starting point to modify the parameters. I think that I would need the book to do the appropriate mapping so that we can match the figures. We may need to rederive since, based on the current information presented, there appear to be some bugs that need to be ironed out.

• Very nice (+1), but again we have the problem that author don't understand subject. What about mixing boundary condition $D_A\partial _x[A]+D_B\partial _x [B]+ D_C \partial _x [C]=0$ at x=0? Probably we should put 2 k2 [B][C] in the first equation. – Alex Trounev Apr 6 at 19:47
• @AlexTrounev Thank you very much! My PhD was in the surface chemistry of supported metal catalysts. A somewhat related field but that was 30 years ago. You may be correct. Preferably, I would prefer to see the surface reaction stoichiometry explicitly written out because the law of mass action needs to be obeyed. If that is obeyed, I think the mixing boundary condition would not be independent. Unfortunately, it is a bit like archaeology at the moment because the presented pieces do not seem to fit together and more digging into the sources will be required to complete the picture. – Tim Laska Apr 6 at 21:09
• @AlexTrounev another thing I just noticed is that it is impossible to produce any of species C with the current reaction mechanism since initially, it starts at zero. So, with the first reaction scheme, we are unable to produce anything. With the second reaction scheme, we are unable to produce species C. We will need yet another reaction scheme to produce all species. – Tim Laska Apr 6 at 21:21
• @Tim Laska If you read the caption in the figure, it provides the values for k1 and k2. I made k1 and k2 equal to 0 to test out the simpler case in which all the terms that you highlighted were 0. My query was around whether the BCs were implemented correctly because I was getting the error mentioned in the OP. – Walser Apr 7 at 6:59
• @Walser Do the reactions occur at the electrode's surface, or in the bulk, or both? I am not an expert in cyclic voltammetry, and I do not have access to your reference, but the references that have looked at do not have bulk source terms. I would suspect the reactions rates at the surface are significantly faster than the bulk. – Tim Laska Apr 7 at 18:05