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

enter image description here

At $t=1e-2$ the concentration for $C_1$ looks like

enter image description here

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. enter image description here enter image description here

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?



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.

  • $\begingroup$ Here is an exact same question but without answer. Unfortunately, I don't have enough reputation to post in the comment section and see if he resolve the problem :( $\endgroup$ – Witty Viper of Hidden Glen Sep 2 '15 at 8:22
  • $\begingroup$ Where is the interface in your plot? How are you discretizing your problem? (Finite elements? Finite differences? Finite volumes?) $\endgroup$ – Geoff Oxberry Sep 2 '15 at 23:59
  • $\begingroup$ The interface is between the jump node and the node on its left. I'm using finite volumes with cell-centered grid. I found that I have a discontinuity on the coefficient between the interface, since the flux on the right of the interface gets a multiple of porosity $\epsilon$, while the flux on the left do not. I simply average the coefficient though. (Regarding this post) Also the potential across the interface is non-smooth. $\endgroup$ – Witty Viper of Hidden Glen Sep 3 '15 at 3:58
  • $\begingroup$ What mathematical conditions must you enforce at the interface? $\endgroup$ – Geoff Oxberry Sep 3 '15 at 5:49
  • $\begingroup$ HEY! I found this paper and replace the variable coefficient to harmonic means (which work for discontinuous coefficient as the author said). There's no discontinuity at the interface! Although he mentioned some other technique, but they're too hard for me though. The solution seems to have steep boundary layer (on the concentration) to tackle next though. $\endgroup$ – Witty Viper of Hidden Glen Sep 3 '15 at 6:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.