\begin{equation} \frac{\partial C_i}{\partial t} = D_i \nabla^2 C_i - \frac{I \cdot \nabla t_i}{z_i F} - \sum_{i'} \frac{z_{i'}}{z_i} D_{i'}\nabla \cdot (t_i\nabla C_{i'}) \end{equation}
\begin{equation} \nabla \eta = \frac{I}{\kappa} + \frac{F}{\kappa} \sum_i z_i D_i \nabla C_i \end{equation} The dependent variables are $C_i$ and $\eta$ ($t_i$ is a function of $C_i$). In the 1D case.
I have some thought on difference scheme (TRBDF2) to try on the first equation. The second equation looks simple, but I'm afraid a simple centered difference will not work. I've try something like
\begin{equation} \left. D_i\frac{\partial C_i}{\partial x}\right \vert_{k-1/2} = D_{i, k-1/2}\frac{C_{i,k}-C_{i,k-1}}{\Delta x} \approx \frac{D_{i, k}+D_{i,k-1}}{2}\frac{C_{i,k}-C_{i,k-1}}{\Delta x} \end{equation}
but with no luck when implementing in the nonlinear solver MINPACK. What will be a good scheme on differencing this two equations? Also, should I run analysis on stability issue every time I came up with a difference scheme before I compute them?
Another set of equation look like this, where the dependent variables are $C_i$, $\eta$, $i_2$, $\epsilon$ and $\epsilon_k$. I've try difference scheme and put on some efforts but in vain. Is there any general guideline as to the way of discretization? I'm using the method of line by the way.
\begin{equation} \frac{\partial \epsilon C_i}{\partial t} = \nabla \cdot \biggl( \epsilon D_i \nabla C_i \biggr) t-\frac{i_2\cdot \nabla t_i}{z_i F} - \sum_j ai_j\biggl( \frac{t_i}{z_i F}+\frac{s_{i,j}}{n_j F} \biggr) - \sum_{i'} \frac{z_{i'}}{z_i}\epsilon D_{i'}\nabla \cdot (t_i\nabla C_{i'}) - R_i \end{equation}
\begin{equation} \nabla \eta = -\biggl( \frac{I-i_2}{\sigma} \biggr) + \frac{i_2}{\epsilon \kappa} + \frac{F}{\kappa} \sum_i z_i D_i \nabla C_i \end{equation}
\begin{equation} \nabla \cdot i_2 = a \sum_j i_j \end{equation}
\begin{equation} \frac{\partial \epsilon}{\partial t} = -\sum_k \overset{\sim}{V_k} k_k\epsilon_k\biggl(\prod_i C_i^{\gamma_i,k}-K_{sp,k}\biggr) \end{equation}
\begin{equation} \frac{\partial \epsilon_k}{\partial t} = \overset{\sim}{V_k} k_k\epsilon_k\biggl(\prod_i C_i^{\gamma_i,k}-K_{sp,k}\biggr) \end{equation}
EDIT To be more specific I'm now trying to solve the initial condition for $\eta$ and $i_2$ which is not given. Since there is no concentration gradient, in the separator (of the battery)
\begin{equation} \nabla \eta = \frac{I}{\kappa} \end{equation} \begin{equation} i_2=I \end{equation}
In the porous cathode
\begin{equation} \nabla \eta = -\biggl( \frac{I-i_2}{\sigma} \biggr) + \frac{i_2}{\epsilon \kappa} \end{equation}
\begin{equation} \nabla \cdot i_2 = a \sum_j i_j \end{equation}
$i_j$ is Butler-Volmer equation, which in the Newton's method tends to explode when initial guess is bad.
\begin{equation} i_j = i_{o,j}\biggl[ \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{p_{i,j}}\exp\biggl(\frac{\alpha_{aj} F}{RT}(\eta_j-U_j)\biggr) - \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{q_{i,j}}\exp\biggl(-\frac{\alpha_{cj} F}{RT}(\eta_j-U_j)\biggr) \biggr] \end{equation}
I discretize the spatial derivative with second order difference, and the output shows a jump every other mesh point (I did not encounter this when I wrote my own code with Newton's method and analytic Jacobian, but I am afraid later on I will make another error when writing down analytic Jacobian since there are so many variables. So I try the MINPACK with numerical Jacobian but with no luck).
I felt like something is wrong in the boundary. Since $\eta$ is not continuous in the separator/porous cathode interface, so I imposed $\nabla \eta_s = \nabla \eta_c$. This is correct because $\eta=\phi_1-\phi_2$, and I assumed $\phi_2$ is continuous in the interface, $\phi_1$ does not exist in the separator but $\nabla \phi_1=0$ in the interface. Particularly in the interface mesh points I wrote
\begin{equation} \frac{-3\eta_{s,n}+4\eta_{s,n-1}-\eta_{s,n-2}}{2\Delta x} = \frac{3\eta_{c,1}-4\eta_{c,2}+\eta_{c,3}}{2\Delta x} \end{equation}
Thanks!
