Mass Transport Model

Question

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["NDSolve`FEM`"]
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];
TableForm[%] // TraditionalForm;

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.

Answer 1

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.

Answer 2

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["NDSolve`FEM`"]

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@
   Thread[label -> 
     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["NDSolve`FEM`"]

(*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];
TableForm[%] // TraditionalForm;

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@
   Thread[label -> 
     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

Answer 3

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.