This is a follow-up post for this question. I am trying to solve the Poisson-Nernst-Planck equation in the form similar to eqns 9 and 10 in this research paper.
Different than the previous question, I consider two ion species ($i=1,2$ and $z_1=1, z_1=-1$) for $\rho_e$, where $\rho_e$ is defined as
Compared with the previous question, I also made some changes to the boundary conditions/parameter value. I consider dirichlet boundary condition for $\varphi$, ($\varphi(x_L)=-u$ and $\varphi(x_R)=0$) and $c$, ($c_1=c_2=\overline{c}$). Then, $j_1 \neq 0$ and $j_2 = 0$
I managed to partially solve the previous question (only for certain parameter values)
Here is my code:
ClearAll["Global`*"]
(*Parameters*)
xRef = 10^(-8);
rhoRef = 0.1*nA;
j = 0.5; (*normalized conc/(m^2s)*)
rhoFStarVal = 120;(*normalized fixed charge concentration*)
rhoBStarVal = 20;(*normalized bulk concentration*)
phiBC = -0.25;(*in thermal voltage unit*)
kB = 1.380649*10^-23;
absTemp = 298.15;
beta = 1/(kB*absTemp);
eps0 = 8.8541878128*10^(-12);
eps = 78.6;
e = 1.60217663*10^(-19);
nA = 6.02214076*10^(23);
dfsvt = 1.33*10^(-9);
xL = 0;(*in xRef*)
xMemL = 0.75; (*in xRef*)
xMemR = 1.25;(*in xRef*)
xR = 2;(*in xRef*)
stepFunc[x_] = 0.5*(Tanh[100*x] + 1);
filterFunc[x_] = stepFunc[x - xMemL]*stepFunc[-(x - xMemR)];
rhoFStar[x_] = rhoFStarVal* filterFunc[x];
zf = -1;
z1 = 1;(*Counter ion valence*)
z2 = -1;(*Co ion valence*)
deltaU = -1/beta;
vW = 30*10^(-30);
debye2 = (kB*absTemp*eps*eps0)/(e^2*rhoRef);
(*Equations*)
elctrcTrnsprtV2 =
debye2/xRef^2*
phi''[x] == -(zf*rhoFStar[x] + z1*rho1fStar[x] + z2*rho2fStar[x])
ctionTrnsprt =
rho1fStar'[x] == -z1*rho1fStar[x]*phi'[x] - (j*xRef)/dfsvt
coionTrnsprt = rho2fStar'[x] == -z2*rho2fStar[x]*phi'[x]
(*NDSolve*)
totalSolV2 =
NDSolve[{ctionTrnsprt, coionTrnsprt, elctrcTrnsprtV2, phi[xL] == 0,
phi[xR] == phiBC, rho1fStar[xL] == rhoBStarVal,
rho2fStar[xL] == rhoBStarVal}, {phi, rho1fStar, rho2fStar}, {x, xL,
xR},
MaxStepSize->0.0001]
(*Define the concentration function concentration[r,t]*)
phiSol = totalSolV2[[1]][[1]][[2]];
rho1fStarSol = totalSolV2[[1]][[2]][[2]];
rho2fStarSol = totalSolV2[[1]][[3]][[2]];
(*Solution Plotting*)
Plot[phiSol[x], {x, xL, xR}, PlotRange -> All,
PlotLabel -> "phi (KbT/e), steady"]
Plot[rho1fStarSol[x], {x, xL, xR}, PlotRange -> {0, All},
PlotLabel -> "Free Counter Ion (0.1 mM/L), steady"]
Plot[rho2fStarSol[x], {x, xL, xR}, PlotRange -> Full,
PlotLabel -> "Free Co Ion (0.1 mM/L), steady"]
Solving the system at rhoBStarVal=5 gives solution
with acceptable numerical residual (LHS-RHS)
However, when i set rhoBStarVal=20, I encountered this error
To solve for when rhoBStarVal=20, I also tried FEM method (Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}}). It manages to output a set of results. However, the solution for rho2fStar exhibits oscillating behaviour, which looks wrong.
