Skip to main content
added 946 characters in body
Source Link
Alex Trounev
  • 48.8k
  • 3
  • 51
  • 115

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

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

Source Link
Alex Trounev
  • 48.8k
  • 3
  • 51
  • 115

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