Solving Fick's second law of diffusion in 1D sphere using FiPy

I'm trying to solve the second law of diffusion for spheres PDE using fipy. I've checked the documentation but there's no example for such a case, so I was wondering if it is actually possible to do it, as I was not successful reaching for the adequate equation definition structure. I consider azimuthal and zenith angular symmetries, so the equation I need to solve results in the following.

Of course, boundary conditions are fixed at r=0 and r=R at fixed values and the initial concentration is also known. I also tried to follow the ideas presented in here but didn't get any clear result for it. Any ideas would be welcomed.

The code I'm using at the moment is the following:

from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable
from builtins import range

nr = 100 #number of cells on the mesh
r_ca = 8.5e-6 #sphere radius

Iapp = -1*1.656 #Current [A] discharge
F = 96485.33212331001 #Faraday constant
S_ca = 1.1167 #Surface of cathode
D_ca = 1e-14 #Diffusion coef.

Xinit_ca = 0.4952 #[p.u.]
Xinit_an = 0.7522

BoundaryR0 = 0 #Fixed flux at r=0
BoundaryR1_ca = -Iapp/(S_ca*F) #Fixed flux at r=r_ca

mesh = SphericalGrid1D(nr=nr, Lr=r_ca) 
X_ca = CellVariable(mesh=mesh, name='Concentration cathode', value=Xinit_ca) 
X_ca.faceGrad.constrain(BoundaryR1_ca, mesh.facesRight) #Fixed flux boundary condition
X_ca.faceGrad.constrain(BoundaryR0, mesh.facesLeft)

eq_ca = TransientTerm() == DiffusionTerm(coeff=D_ca) # dX/dt = D/r^2 * d/dr(r^2*dX/dr)

tstep = 1 #s
Nstep = 1000

for step in range(Nstep):
    eq_ca.solve(var=X_ca, dt=tstep)

viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)

Solution

There are a few things going on:

Some solvers don't like the spherical mesh, probably because of the enormous range in cell volumes. The SciPy LinearLUSolver seems to work. Judicious preconditioning might help other solvers.
Eq. (45) in the paper you linked below defines the flux, but you are constraining the gradient. You need to divide by the diffusivity.
X_ca is in units of [stoichiometry], but BoundaryR1_ca is in units of mol/(m**2 * s) (or mol/m**4 after dividing by D_ca. You need to divide by C_ca_max to get [stoichiometry]/m, as you're solving something halfway between Eq. (43) and Eq. (52).
No-flux is the natural boundary condition for FiPy, so there's no need to constrain at mesh.facesLeft.
The gradient at mesh.facesRight should be constrained to a vector (even in 1D).

With the changes below, I get behavior that looks consistent with Fig. 7 in Mayur et al..

diff --git a/spherical.py b/spherical.py
index e6c2969..1eee346 100644
--- a/spherical.py
+++ b/spherical.py
@@ -1,4 +1,4 @@
-from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable
+from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable, LinearLUSolver
 from builtins import range
 
 nr = 100 #number of cells on the mesh
@@ -13,19 +13,22 @@ Xinit_ca = 0.4952 #[p.u.]
 Xinit_an = 0.7522
 
 BoundaryR0 = 0 #Fixed flux at r=0
-BoundaryR1_ca = -Iapp/(S_ca*F) #Fixed flux at r=r_ca
+BoundaryR1_ca = -Iapp/(51410*D_ca*S_ca*F) #Fixed flux at r=r_ca
 
 mesh = SphericalGrid1D(nr=nr, Lr=r_ca) 
 X_ca = CellVariable(mesh=mesh, name='Concentration cathode', value=Xinit_ca) 
-X_ca.faceGrad.constrain(BoundaryR1_ca, mesh.facesRight) #Fixed flux boundary condition
-X_ca.faceGrad.constrain(BoundaryR0, mesh.facesLeft)
+X_ca.faceGrad.constrain([BoundaryR1_ca], mesh.facesRight) #Fixed flux boundary condition
 
 eq_ca = TransientTerm() == DiffusionTerm(coeff=D_ca) # dX/dt = D/r^2 * d/dr(r^2*dX/dr)
 
 tstep = 1 #s
 Nstep = 1000
+solver = LinearLUSolver()
+
+viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)
 
 for step in range(Nstep):
-    eq_ca.solve(var=X_ca, dt=tstep)
+    eq_ca.solve(var=X_ca, dt=tstep, solver=solver)
+    if step % 100 == 0:
+        viewer.plot()
 
-viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)