NowFor a clearer context, I am basically solving for a transient Poisson-Nernst-Planck (PNP) system. It consists of three major equations: one transport PDE for c1f (co-ion), one transport PDE for c2f (counter-ion), and one poisson equation (the ODE) that governs the electrical potential phi induced by the salt ion and membrane fixed charge group. Additionally, I also consider the binding of the counter-ion at the membrane fixed charge, termed as c2b or $C_{2\text{b}}^* $.
The computational domain consists of the salt solution (0<xStar<3) and the membrane (3<xStar<6). Here, I consider one property difference between the salt solution and the membrane, which is the existence of fixed-charge cF[xStar] in the membrane, while none in the salt solution. This is a piecewise function that is approximated by xStar-dependent hyperbolic tangent function. To note, it is expected to have a steep change in mathematical solution value at the solution-membrane interface, and I believe this might be one of the possible origins for the solving issue.
I believe that a logically-sound numerical solution should exist for my system. This is because I managed to obtain the solution at tStar=0 and steady state.
My code to obtain solution at tStar=0:
(*Phi Solution for Initial Condition*)
ClearAll["Global`*"]
(*Parameters*)
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);
xRef = 10*10^(-6);
tRef = 10*10^(-2);
xStarMem0 = 3;
xStarMemMid = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);
(*BCs*)
phiDirichBC = phi[0] == 0;
phiNuemBC = NeumannValue[0, xStar == xStarMemMid];
(*IC Equations*)
c1fIC[xStar_] = -z2/z1*cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2fIC[xStar_] = cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2bIC[xStar_] =
cF[xStar]*(nA*c2fIC[xStar]*vW)/(nA*c2fIC[xStar]*vW +
Exp[beta*deltaU]);
elctrcTrnsprtIC =
D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2fIC[xStar] + c2bIC[xStar]) +
z1*c1fIC[xStar]) + phiNuemBC;
(*NDSolve*)
totalSol =
NDSolve[{elctrcTrnsprtIC, phiDirichBC}, {phi}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" ->
0.0001}}}];
(*Define the concentration function concentration[r,t]*)
phiSolIC[xStar_] = phi[xStar] /. totalSol[[1]];
(*Solution Plotting*)
Plot[Evaluate@phiSolIC[xStar], {xStar, 0, xStarMemMid},
PlotRange -> All, FrameLabel -> {"xStar", "Elec Pot (V)"},
PlotLabel -> "Electrical Potential, Initial"]
Plot[c1fIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Co Ion, Initial"]
Plot[c2fIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Ct Ion, Initial"]
Plot[c2bIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Binded Ct Ion, Initial"]
My code to obtain steady-state solution:
ClearAll["Global`*"]
(*Parameters*)
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);
xRef = 10*10^(-6);
tRef = 10*10^(-2);
xStarMem0 = 3;
xStarMemMid = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);
(*BCs*)
phiDirichBC = phi[0] == 0;
phiNuemBC = NeumannValue[0, xStar == xStarMemMid];
(*Equations*)
c1f[xStar_] = -z2/z1*cSaltBar*Exp[-beta*z1*e*phi[xStar]] ;
c2f[xStar_] = cSaltBar*Exp[-beta*z2*e*phi[xStar]] ;
c2b[xStar_] =
cF[xStar]*(nA*c2f[xStar]*vW)/(nA*c2f[xStar]*vW + Exp[beta*deltaU]);
elctrcTrnsprt =
D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2f[xStar] + c2b[xStar]) +
z1*c1f[xStar]) + phiNuemBC;
(*Input Parameter Plotting*)
Plot[cF[x], {x, 0, xStarMemMid}, PlotRange -> All,
PlotLabel -> "Fixed Charge Distribution (mol/m3)"]
(*NDSolve*)
totalSol =
NDSolve[{elctrcTrnsprt, phiDirichBC}, {phi}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];
(*Define the concentration function concentration[r,t]*)
phiSol[xStar_] = phi[xStar] /. totalSol[[1]];
(*Solution Plotting*)
Plot[Evaluate@phiSol[xStar], {xStar, 0, xStarMemMid},
PlotRange -> All, FrameLabel -> {"xStar", "Elec Pot (V)"},
PlotLabel -> "Electrical Potential, steady"]
Plot[c1f[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Co Ion, steady"]
Plot[c2f[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Ct Ion, steady"]
Plot[c2b[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Binded Ct Ion, steady"]
To solve the Poisson-Nernst-Planck in transient state, I have adopted the Method of Lines approach from question 78493, which is commonly used to solve for integro-differential equation and PDE-ODE-coupled system. To discretize my system, I am using the pdetoode function from the accepted answer of question 127980. Not only that, I also differentiated both sides of my ODE wrt timeStar, making it a PDE. Subsequently, I included the initial condition for phi, which is the variable to be solved from my ODE.
