Neumann boundary conditions are not satisfied

Question

I am solving a system of 3 ODEs and I want to implement neumann BC (d/dx = 0 @x=0). However, the solution does not seem to obey neumann BC when parameter A1/A2/A3 is non zero, as I notice that the solution curves at x=0 is not flat.

List of constant parameters:

s = 1; (*ratio between Kc and Kb*)
p = 5; (*ratio between Db and Da*)
q = 1; (*ratio between Dc and Da*)
u = 0.5; (*ratio between surface concentration of A and initial \
concentration of active sites*)
TM = 1; (*Thiele modulus*)
Rba = 0.6; (*ratio between surface concentration of species B and A*)
\

Rca = 0.5; (*ratio between surface concentration of species C and A*)
\

Kb = 1.38*10^-23;  (*boltzman constant*)
T = 25 + 273; (*temperature*)
(*coefficients for function of non-ideal chemical potential*)
A1 = 1;
A2 = 1;
A3 = 1;
B1 = 1;
B2 = 1;
B3 = 1;
C1 = 1;
C2 = 1;
C3 = 1;

Some intermediate expressions:

(*Non-ideal chemical potential function of each species*)
chemicalpotential1[x_] = 
 A1/2*(1 - Tanh[B1 (x - C1)])*Kb*T; (*for species A*)
chemicalpotential2[x_] = 
 A2/2*(1 - Tanh[B2 (x - C2)])*Kb*T; (*for species B*)
chemicalpotential3[x_] = 
 A3/2*(1 - Tanh[B3 (x - C3)])*Kb*T; (*for species C*)

(*Expressions that involve non-ideal chemical potential function in \
the first order differential*)
diffintermediate1[x_] = 
 a[x]*x^2*chemicalpotential1'[x];(*for species A*)
diffintermediate2[x_] = 
 b[x]*x^2*chemicalpotential2'[x];(*for species B*)
diffintermediate3[x_] = (phi[x] - 1)*x^2*
  chemicalpotential3'[x];(*for species C*)

The system of ODEs and BCs:

ODE1 = a''[x] + 2/x*a'[x] + 1/(Kb*T)/x^2*diffintermediate1'[x] - 
   TM^2*(1 + s)*phi[x]*a[x] == NeumannValue[0, x == 0]; ODE2  = 
 b''[x] + 2/x*b'[x] + 1/(Kb*T)/x^2*diffintermediate2'[x] + 
   TM^2/p*phi[x]*a[x] == NeumannValue[0, x == 0];
ODE3  = phi''[x] + 2/x*phi'[x] + 1/(Kb*T)/x^2*diffintermediate3'[x] - 
   s*TM^2*u/q*phi[x]*a[x] == NeumannValue[0, x == 0]; ODES = {ODE1, 
  ODE2 , ODE3}; (*Grouping the ODEs for input in NDSolveValue \
function*)

(*Dimensionless dirichlet BCs for species A, B, and C*)
bcs = {  a[1] == 1 , b[1] == Rba, phi[1] == 1 - u*Rca}

Solving and plotting the ODEs,

{asol, bsol, phisol} = 
  NDSolveValue[{ODES, bcs}, {a, b, phi}, {x, 0, 1}, 
   Method -> "FiniteElement"];
csol[x_] = (1 - phisol[x])/u; (*transforming phi into dimensionless c*)

(*Plotting the concentration profiles*)
Plot[{asol[x], bsol[x], phisol[x], csol[x]}, {x, 0, 1}, 
 PlotLegends -> "Expressions"]

Here is the code for testing the gradient at x=0, I realize that gradient@x=0 is more away from 0 when the A1/A2/A3 has larger magnitude, which oppose what I want for neumann BC.

(*Testing the gradient at x=0*)
differentiationdsd[x_] = phisol'[x];
differentiationdsd[0]

Answer 1

Not sure how to explain the non-convergence of phisol'[0]. (Though phisol'[0] == -0.0302716 is close to zero, it won't get closer to zero with smaller MaxCellMeasure. This might be related to the interpolation method of FiniteElement. ) In this answer I just want to point out that

… the solution does not seem to obey neumann BC when parameter A1/A2/A3 is non zero, as I notice that the solution curves at x=0 is not flat

is more of an illusion. The solution found by FiniteElement method is reliable, I'll verify this by solving the problem with finite difference method (FDM).

I'll use pdetoae for generation of finite difference equations.

varlst = {a, b, phi};
newode = Map[# x &, ODES, {2}] /. NeumannValue[__] :> 0 // Expand
bcadd = Table[var'[0] == 0, {var, varlst}];

domain = {0, 1};
points = 25; 
grid = Array[# &, points, domain]; 
difforder = 2;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[varlst[x], grid, difforder];
del = #[[2 ;; -2]] &;

ae = del /@ ptoafunc@newode;
aebc = ptoafunc@{bcadd, bcs};

MapThread[(guess[#][x_] = #2) &, {varlst, bcs[[All, -1]]}];

solrulelst = 
  FindRoot[Flatten@{ae, aebc}, 
   Flatten[#, 1] &@Table[{var[x], guess[var][x]}, {var, varlst}, {x, grid}]];

solmatlst = ArrayReshape[solrulelst[[All, -1]], {3, points}];

solfunclst = 
  ListInterpolation[#, {grid}, InterpolationOrder -> difforder] & /@ solmatlst;

solfunclst[[3]]'[0]
(* -3.9968*10^-15 *)

As we can see, with 25 grid points and 2nd order difference formula, the obtained $\phi'(0)$ is very close to zero, while these solutions are visually the same as those obtained with FiniteElement:

ref = Plot[{asol[x], bsol[x], phisol[x]}, {x, 0, 1}, 
   PlotLegends -> "Expressions"];
Show[ListPlot@solfunclst, ref]

Here I've made use of an undocumented syntax of ListPlot. See this post for more info.