2
$\begingroup$

My original question was here. Thank you to those who took the time to answer. Essentially I was trying to replicate the concentration profiles for a simple cyclic voltammetry experiment. My code seems to match the boundary conditions but the answer doesn't match the one in the book and is obviously incorrect at places. I'll provide my code at the end but first I'd like to provide the complete system (from Compton's Understanding Voltammetry, chapter 4):

Diffusion for species A and B

The boundary conditions are:

enter image description here

The electrochemical rate constants and the potential equations are: enter image description here

enter image description here

The concentration profile should come as follows for various values of potential; and the rest of the parameters are given in the caption:

enter image description here

Now for my code. The parameter large is to model infinity:

ClearAll["Global`*"]

(*Experimental Parameters*)
cAbulk := 1*10^-3;
k0:= 10^-5;
rtbyf:= 25.7 10^-3(*volt*);
dA:= 10^-5; dB:= 10^-5;
\[Alpha]:=0.5; \[Beta]:=0.5; T= 298; ef0= 0;
ts:= 1;\[Nu]:= -1;e1:= 0.5;

large=10;

(*Equations for Reaction and other governing equations*)
eqn1=D[cA[x,t],t]-Div[dA Grad[cA[x,t],{x}],{x}];
eqn2=D[cB[x,t],t]-Div[dB Grad[cB[x,t],{x}],{x}];


e[t_]:= Piecewise[{{e1+ \[Nu] t, 0<= t<= ts}, {e1+ 2\[Nu] ts-\[Nu] t,ts<= t<= 2 ts}}]
kc:= k0 Exp[-\[Alpha]/rtbyf (e[t]- ef0)]
ka:= k0 Exp[\[Beta]/rtbyf (e[t]- ef0)]


(*Solution*)
{ncA,ncB}=NDSolveValue[{eqn1==NeumannValue[kc/dA cA[x,t]-ka/dA cB[x,t],x==0&&(t>=0)],
DirichletCondition[cA[x,t]==cAbulk,t==0],
DirichletCondition[cA[x,t]==cAbulk,(x==large)&&(t>=0)],

eqn2==NeumannValue[-kc /dB cA[x,t]+ka/dB cB[x,t],x==0&&(t>=0)],
DirichletCondition[cB[x,t]==0,t==0],
DirichletCondition[cB[x,t]==0,x==large]},{cA,cB},{x,0,large},{t,0,2ts},
AccuracyGoal->20,PrecisionGoal->20,
InterpolationOrder->All];


(*Presentation of Results*)
Manipulate[
{Show[
Plot3D[{ ncA[x,t]/cAbulk, ncB[x,t]/cAbulk},{t,0,2ts},{x,0,0.1  large},PlotRange->All, Boxed-> False, Axes -> True, AxesLabel->Automatic, 
MeshFunctions->{#1&}, MeshStyle->Gray, LabelStyle->Directive[Bold,Blue, 12], PlotStyle->Opacity[0.5], WorkingPrecision->20, PlotLegends-> {cA,cB}],
Plot3D[{ ncA[x,t]/cAbulk},{t,dt,dt+0.03},{x,0,0.1  large}, PlotStyle-> Red, PlotRange->All, Boxed-> False, Axes -> True, AxesLabel->Automatic, MeshFunctions->{#1&}, LabelStyle->Directive[Bold,Blue, 12], WorkingPrecision->20]
],

Plot[{ncA[x,dt]/cAbulk, ncB[x,dt]/cAbulk}, {x,0,0.1 large}, PlotRange-> {0,2}, PlotLegends -> {cA, cB}]}, 
{dt,0,2ts}, ControlPlacement->Top]

The output below however has the bulk normalized concentration of specie A exceeding 1 at places; so there is obviously something wrong. Moreover the 2D plot shows the concentration dropping below 0 at intervals of 0.2mm. The finite element solution provided by Tim Laska doesn't do either of those things.

enter image description here

So where did I go wrong?

Thanks in advance.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ 1. You should not use FEM in $t$ direction, for more info check this post: mathematica.stackexchange.com/a/222373/1871 2. You haven't set NeumannValue correctly, please check the document of NeumannValue carefully, especially the Details section. 3. As mentioned under your previous question, if the domain is always regular, consider using the old good TensorProductGrid instead of FiniteElement. 4. As to AccuracyGoal and PrecisionGoal, you may want to check this: mathematica.stackexchange.com/q/118249/1871 $\endgroup$
    – xzczd
    Commented Feb 22, 2021 at 6:24
  • $\begingroup$ 5. large is too large, either decrease it, or set a smaller MaxCellMeasure. $\endgroup$
    – xzczd
    Commented Feb 22, 2021 at 6:44

1 Answer 1

Reset to default
5
$\begingroup$

One can define an "infinite" domain-based solely on $t_{max}$ for a purely diffusive problem. The key disadvantage of using this approach versus anisotropic meshing is that the mesh will change as a function of $t_{max}$ and $\mathcal{D}$.

The analytical solution for one-dimensional diffusion is:

$$\frac{{{C_A}\left( {t,x} \right)}}{{{C_{A0}}}} = Erfc\left( {\frac{x}{{\sqrt {4{\cal D}t} }}} \right) = Erfc\left( \eta \right);\eta = \frac{x}{{\sqrt {4{\cal D}t} }}$$

At $\eta=3$, the complementary function has decayed close to zero as shown below:

Erfc[3] // N
(* 0.0000220905 *)

Thus, a reasonable definition for the infinite domain, ${x_\infty }$ is:

$${x_\infty } = 6\sqrt {{\cal D}{t_{\max }}}$$

Here is a workflow that will put 100 elements across the "infinite domain." The results are close to the textbook solution, but some conversion factors will need to be applied to match the amperage.

ClearAll["Global`*"]
(*Experimental Parameters*)
cAbulk := 1*10^-3;
cAbulk := 1;
k0 := 10^-5;
rtbyf := 25.7 10^-3(*volt*);
dA := 10^-5; dB := 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}}]
kc[t_] := k0 Exp[-α/rtbyf (e[t] - ef0)]
ka[t_] := k0 Exp[β/rtbyf (e[t] - ef0)]

vars = {{cA[t, x], cB[t, x]}, t, {x}};
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}|>|>;
ops = MassTransportPDEComponent[vars, pars];
Γflux = 
  MassTransferValue[x == 0, vars, pars, "BoundaryCondition1"];
Γcondinf = 
  MassConcentrationCondition[x == large, vars, pars, 
   "BoundaryCondition2"];
ics = {cA[0, x] == cAbulk, cB[0, x] == 0};
{cAfun, cBfun} = 
  NDSolveValue[{ops == Γflux, Γcondinf, 
    ics}, {cA, cB}, {t, 0, tmax}, {x, 0, large}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> large/1000}}}];
(*Presentation of Results*)
Manipulate[
 Row[{Row[{ParametricPlot[{-(e[t] - ef0), 
       kc[t] (cAfun[t, 0] - (ka [t] cBfun[t, 0])/kc[t])}, {t, 0, 
       tmax}, ColorFunction -> (ColorData["DarkRainbow"][#3] &), 
      PlotPoints -> 200, AspectRatio -> 1, 
      PlotLegends -> 
       Placed[BarLegend[{"DarkRainbow", {0, tmax}}, 
         LegendFunction -> "Frame", LegendLabel -> "t"], {After, 
         Top}], ImageSize -> 350, 
      Epilog -> {PointSize[Large], 
        Point[{-(e[dt] - ef0), 
          kc[dt] (cAfun[dt, 0] - (ka[dt] cBfun[dt, 0])/kc[dt])}]}]}, 
    ImageMargins -> {{100, 100}, {0, 0}}], 
   Row[{Show[
      Plot3D[{cAfun[t, x], cBfun[t, x]}, {t, 0, tmax}, {x, 0, 
        0.8 large}, PlotRange -> All, Boxed -> False, Axes -> True, 
       AxesLabel -> Automatic, MeshFunctions -> {#1 &}, 
       MeshStyle -> Gray, LabelStyle -> Directive[Bold, Blue, 12], 
       PlotStyle -> Opacity[0.5], WorkingPrecision -> 20, 
       PlotLegends -> {cA, cB}], 
      Plot3D[{cAfun[t, x]}, {t, dt, dt + 0.03}, {x, 0, 0.8 large}, 
       PlotStyle -> Red, PlotRange -> All, Boxed -> False, 
       Axes -> True, AxesLabel -> Automatic, MeshFunctions -> {#1 &}, 
       LabelStyle -> Directive[Bold, Blue, 12], 
       WorkingPrecision -> 20], ImageSize -> 400], 
     Plot[{cAfun[dt, x], cBfun[dt, x]}, {x, 0, 0.8 large}, 
      PlotRange -> {0.01, 1.01}, PlotLegends -> {cA, cB}, 
      PlotLabel -> StringTemplate["E=`` mV"][1000 (e[dt] - ef0)], 
      AspectRatio -> 1.2, ImageSize -> 250]}]}], {dt, 0, tmax, 
  tmax/200}, ControlPlacement -> Top]

Manipulate at -500 mV

Manipulate at 410 mV

The following manipulate reproduces points (A) through (E) from the textbook example quite well.

(* Find times that reproduce textbook cases *)
inv[f_, s_, s0_] := Function[{t}, s /. FindRoot[f - t, {s, s0}]]
einv = inv[1000 e[t], t, 0.9];
times = einv /@ {-300, -412, -500};
einv = inv[1000 e[t], t, 1.1];
times = Join[times[[1 ;; 3]], einv /@ {0, 406, 500}];
(*Presentation of Results*)
Manipulate[
 Row[{Row[{ParametricPlot[{-(e[t] - ef0), 
       kc[t] (cAfun[t, 0] - (ka [t] cBfun[t, 0])/kc[t])}, {t, 0, 
       tmax}, ColorFunction -> (ColorData["DarkRainbow"][#3] &), 
      PlotPoints -> 200, AspectRatio -> 1, 
      PlotLegends -> 
       Placed[BarLegend[{"DarkRainbow", {0, tmax}}, 
         LegendFunction -> "Frame", LegendLabel -> "t"], {After, 
         Top}], ImageSize -> 350, 
      Epilog -> {PointSize[Large], 
        Point[{-(e[dt] - ef0), 
          kc[dt] (cAfun[dt, 0] - (ka[dt] cBfun[dt, 0])/kc[dt])}]}]}, 
    ImageMargins -> {{200, 200}, {0, 0}}], 
   Row[{Show[
      Plot3D[{cAfun[t, x], cBfun[t, x]}, {t, 0, tmax}, {x, 0, 
        0.8 large}, PlotRange -> All, Boxed -> False, Axes -> True, 
       AxesLabel -> Automatic, MeshFunctions -> {#1 &}, 
       MeshStyle -> Gray, LabelStyle -> Directive[Bold, Blue, 12], 
       PlotStyle -> Opacity[0.5], WorkingPrecision -> 20, 
       PlotLegends -> {cA, cB}], 
      Plot3D[{cAfun[t, x]}, {t, dt, dt + 0.03}, {x, 0, 0.8 large}, 
       PlotStyle -> Red, PlotRange -> All, Boxed -> False, 
       Axes -> True, AxesLabel -> Automatic, MeshFunctions -> {#1 &}, 
       LabelStyle -> Directive[Bold, Blue, 12], 
       WorkingPrecision -> 20], ImageSize -> 400], 
     Plot[{cAfun[dt, x], cBfun[dt, x]}, {x, 0, 0.8 large}, 
      PlotRange -> {0.01, 1.01}, PlotLegends -> {cA, cB}, 
      PlotLabel -> StringTemplate["E=`` mV"][1000 (e[dt] - ef0)], 
      AspectRatio -> 1.2, ImageSize -> 250]}]}], {dt, times}, 
 ControlPlacement -> Top]

Textbook cases

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.