Is the Numerical Dissolution Time Calculation for a Sphere Accurate?

[henxan]

The problem statement
Based on an analytical solution for the concentration profile of a dissolving sphere, I am supposed to use a numerical method to find the time at which the sphere has fully dissolved. This is not so much a question about specific values - but about the technique by which I attain an answer, thus I have not given a lot of constants or source code. My question is more:
- Are there invalid logic in my thinking?

Relevant equations
Analytical solution to the concentration profile is given by:
$$
\begin{equation}
C(r, t) = C_s\cdot C_{int} \frac{R_0}{r^{\text{*}}} \left( 1 - \textbf{erf}\left( \frac{r^{\text{*}} - R_0}{2\sqrt{Dt}} \right) \right),
\end{equation}
$$
where $r^{\text{*}} = R_0 + r.$ The start volume of the sphere:
$$
\begin{equation}
V_0 = \frac{4}{3}\pi R_0^3
\end{equation}
$$
The concentration in this sphere is $C_s = 1$, thus giving an initial amount of $C_{tot}^0=V_0\cdot C_s = 1$. The inter-facial concentration is $C_{int} = 0.0025$, $R_0 = 2.5\cdot10^{-8}$ and $D=8.5542\cdot10^{-15}$.

An attempt at a solution
I make a matrix in Matlab with $t\times C(r)$. I have also solved the concentration profile using Fick's second law, so I am pretty sure i get the correct diffusion profiles at different times. Since I now have a matrix of concentration profiles at different times, i use the following formula to element-wise solve the sum concentration at all my timesteps between 10 and 15 seconds:
$$
\begin{equation}
C_{tot} = \frac{4}{3}\pi \sum_{i=1}^m \left( r_{i+1}^3 - r_i^3 \right)\frac{\left(C_{i+1}+C_i\right)}{2},
\end{equation}
$$
where $r = 0,1\cdot10^{-9},\ldots, 1\cdot10^{-6}$. My thinking is that when $C_{tot}^0 = C_{tot}$, the particle is totally dissolved. Using this method i found a time for total dissolution of the particle at 10.643 s, using a ridiculous small $\Delta r$. So, what is the problem with this? Well, the answer should lie between 14 and 14.5 seconds, thus I am very far off the mark.

 

[mfb]

If you have the analytic solution, why do you need a numerical simulation?
Anyway, you can check if your numerical results are close to the analytic prediction.

How do you apply Fick's law at the boundary of C_s=1?

 

[henxan]

Hi!

Thanks for the reply!

My row elements are from $r=0$ with step $\Delta r$ to max size of $r$. So, at column one:
$r_1=0$ and $C_i = C_{int} = 0.0025$ for all elements.

Hope this was answer to what you were wondering about.

 

[mfb]

The concentration in this sphere is C_s=1

How does that fit to a constant concentration everywhere? If the concentration is the same everywhere, you should not see any time-evolution of the system as that is an equilibrium.

 

[henxan]

The problem we are dealing with is that of diffusion controlled dissolution of a solid. The solid has a concentration (typically set to 1) while the concentration at the liquid solid interface is a fraction of that This is illustrated in the following image (From On The Kinteics of Precipitate Dissolution by M. J. Whelan in Metal Science Journal):

Thus, at the solid-liquid phase boundary, the concentration is constant. While I use $C$ the image above uses $\rho$. C and S corresponds to my S and inf respectively.

 

[mfb]

Thus, at the solid-liquid phase boundary, the concentration is constant.

You get a constant concentration value for a position that depends on time. This has nothing to do with the type of constants you mentioned above.

 

[henxan]

Well, as I have understood it, the analytical solution above, gives the concentration profile from the sphere surface. It does not tell me what the radius of the sphere is, but it tells me what the diffusion profile is at the surfae and at distances from the sphere surface. Am I wrong - and is it here I am at error? oo)

 

[mfb]

The analytic solution has to have some smallest r where it is valid, this should be the surface of your sphere, otherwise the analytic solution is weird.

 

[Chestermiller]

I wouldn't do this using spatial finite differences. I would start out by using the analytic solution to find the concentration gradient at the surface of the sphere as a function of time. I would multiply this by the diffusion coefficient to get the dissolution flux out of the sphere. This would tell me how fast the surface recedes.

Chet

 

[henxan]

Thanks for your replies!

I was aware of the method you describe Chestmiller, though I wanted to see if I could make the "shell method work". However, it seemed like this method required such enormous amount of cells to calculate correct, so I went for the solution you described Chet. Thus, because I eventually were going to use my method with a numerical solution to Fick's second law:
$$
\frac{\partial C}{\partial t} = D\left( \frac{\partial^2 C}{\partial r^2} + \frac{2}{r^{\text{*}}} \frac{\partial C}{\partial r} \right)
$$
Where $r^{\text{*}} = R(t) + r$, I used the following method to recalculate $R$ (In accordance with Chet's suggestion):
$$
\Delta R = \frac{D}{(C_p-C_i)}\frac{\Delta C_1}{\Delta r_1}\Delta t
$$
Giving a new $R(t) = R(t-\Delta t) + \Delta R$. The reasoning by this, for those not familiar with the solution, I tried to illustrate with the following illustration:

The shell method I tried using gave correct estimates for a 1D problem, but very bad results for a 3D problem. I think I have to conclude that the shell method does work, but that it requires an extremely large amount of distinctive shells, thus leading to very high computing power requirements.

Thanks, H.