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.
NDSolveValue
should return a list with three elements." No,NDSolveValue
already fails after the first warning, what's returned is an unevaluatedNDSolveValue[…]
. Just execute theNDSolveValue[…]
separately and observe. 2. As indicated by the description ofNDSolveValue::fembcdepderiv
, you can not have derivatives inNeumannValue
. (Please observe what's insideΓflux
.) Do noticeMassTransportPDEComponent
, etc. are no more than generator of PDEs and b.c.s, in other words, ifNeumannValue
isn't able to do something,MassTransportPDEComponent
, etc. won't help either. $\endgroup$TensorProductGrid
instead ofFiniteElement
. $\endgroup$