2
$\begingroup$

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. enter image description here

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 enter image description here

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

enter image description here

with acceptable numerical residual (LHS-RHS)

enter image description here

However, when i set rhoBStarVal=20, I encountered this error enter image description here

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.

enter image description here

$\endgroup$

1 Answer 1

2
$\begingroup$

As one of the possibly solutions we can decrease accuracy and use explicit method as follows

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}, AccuracyGoal -> 5, PrecisionGoal -> 5, 
  Method -> "ExplicitEuler"]


(*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"]} 

Figure 1

Second method is the Shooting method

sol2 = 
 NDSolve[{ctionTrnsprt, coionTrnsprt, elctrcTrnsprtV2, phi[xL] == 0, 
   phi[xR] == phiBC, rho1fStar[xL] == rhoBStarVal, 
   rho2fStar[xL] == rhoBStarVal}, {phi, rho1fStar, rho2fStar}, {x, xL,
    xR}, AccuracyGoal -> 5, PrecisionGoal -> 5, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {phi[xL] == 0, phi'[xL] == -.7, 
      rho1fStar[xL] == rhoBStarVal, rho2fStar[xL] == rhoBStarVal}}]

Visualization

{Plot[phi[x] /. sol2[[1]], {x, xL, xR}, PlotRange -> All, 
  PlotLabel -> "phi (KbT/e), steady"],
 Plot[rho1fStar[x] /. sol2[[1]], {x, xL, xR}, PlotRange -> {0, All}, 
  PlotLabel -> "Free Counter Ion (0.1 mM/L), steady"],
 Plot[rho2fStar[x] /. sol2[[1]], {x, xL, xR}, PlotRange -> Full, 
  PlotLabel -> "Free Co Ion (0.1 mM/L), steady"]}

Figure 2

$\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.