The equation at the left of the interface is
\begin{equation} \displaystyle\frac{\partial C_i}{\partial t} = D_i \nabla^2 C_i - z_i \frac{D_i}{RT}F \nabla \cdot (C_i \nabla \phi_2) \end{equation}
The equation at the right of the interface is
\begin{equation} \displaystyle\frac{\partial \epsilon C_i}{\partial t} = \nabla \cdot \biggl(\epsilon D_i \nabla C_i + z_i \frac{ D_i}{RT}F\epsilon C_i \nabla \phi_2 \biggr) + S(x) \end{equation}
with an extra source term accounted for pore wall reaction and precipitation (also an extra dependent variable $\epsilon$ for porosity).
At $t=1e-5$ the concentration for $C_1$ looks like
At $t=1e-2$ the concentration for $C_1$ looks like
But since the value of the source term is negative, I suppose it should be a upward parabola at the interface.
From time $t=1e-3$ to $t=1e-2$, the distribution for $\phi_2$ changes dramatically at the left of the interface though.
Some concentrations on the other hand are smooth across the domain though, but $C_1$ and $C_8$, which have the two largest initial condition values, are discontinuous like the picture shown at the interface.
After $t=1e-2$ the RADAU5 solver converges very slowly (I guess it won't converge to a meaningful values anyway). I'm using cell-centered finite volume for it is easier to implement the Neumann boundary conditions. The discontinuous point is at the first cell right at the left of the interface. I'm using $N=100$; for $N=400$ the discontinuous point even show earlier.
What could be the problems? And how should I resolve it?
Thanks!
Note:
It's assumed that concentration and potential are continuous at the interface, which in fact they are, but in a cell-centered grid, I'm thinking I can't just naively do nothing at the interface.