Diffusion System with Boundary Conditions

Question

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):

The boundary conditions are:

The electrochemical rate constants and the potential equations are:

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

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

ClearAll["Global`*"]
(Experimental Parameters)
cAbulk := 110^-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.

So where did I go wrong?

Thanks in advance.

Answer 1

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]

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]