I am working on a project involving computer simulation of temperature-controlled release of hydrogen from spherical pellets. My goal is to make a model that estimates the temperature curve needed to maintain a constant mass flux of hydrogen out of the system. I have been told to neglect heat transfer due to the high thermal conductivity of the substance being used, and thus can assume that temperature changes are instant and uniform for this preliminary model. Additionally, I assume that, for this system, the chemical potential gradient simplifies to the concentration gradient. However, even with these big simplifications, I am having trouble finding a function for temperature that works for the model.
The main problem I am having is solving Fick’s Second Law for the system. I have not taken a partial differential equations course, and the only solution tactic I know how to use personally is separation of variables, which this model has proven resistant to (or I’ve just been erroneous in my calculations).
Since spatial temperature gradients are being neglected, diffusion takes place in the radial direction only; thus, Fick’s second law is simplified to
$$\frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial C}{\partial r}\right) = \frac{\partial C}{\partial t}$$
Where $C$ is concentration, $r$ is radius, and $t$ is time.
In this case, it is desired for the molar flux leaving the system to be constant and, in all cases of symmetrical diffusion, the flux at the midpoint is zero. Additionally, at the beginning, the concentration is assumed to be uniform and known. Therefore, the boundary conditions become
$$ \begin{align} @\,r &= R &\quad N &= -D \frac{dC}{dr} = f\\ @\,r &= 0 &\quad N &= -D \frac{dC}{dr} = 0\\ @\,t &= 0 &\quad C(r) &= C_i \end{align} $$
Where $R$ is the radius of the pellet, $C_i$ is the initial concentration, $N$ is the molar flux, and $D$ is the diffusivity. For the system in question, hydrogen diffuses along the inner surface of pores in the pellet, and thus the diffusivity is defined as
$$D = \frac{va^2}{z}\exp{\left(\frac{-E}{k_B T(t)}\right)}$$
Where $v$ is the vibrational frequency of hydrogen, a is the bond hopping distance, $z$ is a coordination number, $E$ is the potential energy barrier to diffusion, $k_B$ is the Boltzmann constant, and $T(t)$ is temperature as a function of time. This tracer diffusivity is then adjusted using experimentally determined values to account for the fact that diffusion takes place on the pore surface, not throughout the entirety of the sphere.
Like I said earlier, I have tried using separation of variables on this problem but it has proven ineffective. Mostly, I’m lost on trying to figure out what way to approach this problem. There are so many methods available for solving PDEs and I don’t have the time to learn all of them and hope one works. If someone could point me in the direction of a solution strategy it would be much appreciated.
Also the only computational tool I have available to me is MATLAB, so I don’t have access to a good computational tool for solving PDEs. I am aware of the finite element method, but it seems that coding a finite element solver in MATLAB from scratch would take a lot of time even if I was already familiar with its intricacies, and thus I would only consider it as a last resort if no other approach would work. However, MATLAB is useful when it comes to solving convergent infinite summations and numeric computation of integrals, so solution methods that involve those are fine.
I’m new to this site so if this doesn’t belong here or if I’m missing something from my question statement please tell me and I’ll do what I can to fix it.