# Diffusion inside sphere

Imagine I have the following situation: inside a sphere containing mostly water, solute A is diffusing from the center to the end of the sphere:

At r=0, mass of A is continuously generated through the equation:

$\frac{\partial C(r=0,t)}{\partial t} = \dot{r}_{g}$

where $\dot{r}_{g}$ is a constant giving the mass of A formed per volume, per time.

Also,

$\frac{\partial C(r=0,t)}{\partial r} = 0$

due to symmetry reasons. In my opinion, at the end, this translates into the boundary condition $C(r=0,t) = \dot{r}_{g}\cdot t$

Moreover, at $t = 0$, $C(r,t=0) = 0$, because no A was still generated inside the sphere. So this is my initial condition.

From $r = 0$ to $r = R$, the diffusion equation for A applies:

$\frac{\partial C(r,t)}{\partial t} = D[\frac{\partial^2 C(r,t)}{\partial r^2} + \frac{2}{r}\frac{\partial C(r,t)}{\partial r}]$

What I would like to know is if there is a solution to this partial differential equation and how can we obtain it.

I tried to use separation of variables:

$C(r,t) = R^{*}(r)\cdot T(t)$

converting the original PDE into:

$\frac{1}{T}\frac{\partial T}{\partial t} = D[\frac{1}{R^{*}}\frac{\partial^2 R^{*}}{\partial r^2} + \frac{2}{r\cdot R^{*}}\frac{\partial R^{*}}{\partial r}]$

Although the following equation can be solved:

$\frac{1}{T}\frac{\partial T}{\partial t} = \lambda$

where $\lambda$ is the separation constant, I don't know how to proceed with the other equation:

$r\frac{\partial^2 R^{*}}{\partial r^2} + 2\frac{\partial R^{*}}{\partial r} - \frac{\lambda}{D}rR^{*} = 0$

in order to make the boundary and initial conditions that I specified valid.

• Might Mathematics be better suited for this math question? – Kyle Kanos Jan 10 '17 at 18:39

Drop the unnecessary asterisk for simplicity:

$$C(r,t)=R(r)T(t)$$

After separation, drop the partials (separation separates the PDE into two ODEs), move $D$ to the lefthand side and use separation constant $-\lambda^2$, so you get two ODEs:

$$\frac{1}{DT}\frac{T'}{T}=-\lambda^2\tag{1}$$ $$rR''+2R'+\lambda^2rR=0\tag{2}$$

(Note: a sign error has been corrected here)

The latter is a Sturm-Liouville equation and it solves to:

$$R(r)=c_1\frac{\sin \lambda r}{r}+c_2\frac{\cos \lambda r}{r}\tag{3}$$

This offers the possibility of a fairly simple solution, at least for certain types of boundary conditions.

For example we can make it a bound problem by specifying: $$C(0,t)=C_0$$

Simplify a lot by a simple substitution:

$$u(r,t)=C(r,t)-C_0$$

So that the first boundary condition becomes homogeneous:

$$u(0,t)=0$$

As second boundary condition we can choose:

$$C_r(R,t)=0\implies u_r(R,t)=0\implies R'(R)=0,$$

where $R$ is the radius of the sphere.

Using the first BC in $(3)$, we can see that necessarily:

$$c_2=0,$$

or else $R(0)\to +\infty$.

So we have:

$$R(r)=c_1\frac{\sin \lambda r}{r}$$

Applying the second BC gives us:

$$c_1\Big(-\frac{\sin \lambda R}{R^2}+\lambda\frac{\cos \lambda R}{R}\Big)=0$$

Which, assuming $c_1\neq 0$ and $R\neq0$, reworks to:

$$\boxed{\lambda=\frac {\tan \lambda R}{R}}$$

This transcendental equation has no analytical solutions but does have an infinite number of numerical solutions: $\lambda_1$, $\lambda_2$, $\lambda_3$... $\lambda_n$. These are the eigenvalues of the originl PDE.

A particular solution of $R$ is thus:

$$R_n(r)=c_1\frac{\sin \lambda_n r}{r}$$

Equation $(1)$ simply solves to:

$$T_n(t)=c_3e^{-D\lambda_n^2 t}$$

The particular solution to $u(r,t)$ is thus:

$$u_n(r,t)=A_n\frac{\sin \lambda_n r}{r}\times e^{-D\lambda_n^2 t}$$

Due to the Superposition Principle the overal solution becomes:

$$u(r,t)=\displaystyle\sum_{n=1}^{+\infty}\Big(A_n\frac{\sin \lambda_n r}{r}\times e^{-D\lambda_n^2 t}\Big)$$

The coefficients $A_n$ are found from the initial condition for $r>0$ and using the substitution higher up:

$$C(r,0)=0\implies u(r,0)=-C_0$$

$$\implies -C_0=\displaystyle\sum_{n=1}^{+\infty}A_n\frac{\sin \lambda_n r}{r}$$

Using Fourier we then get the coefficients $A_n$ from:

$$A_n=\frac2R\int_0^R\Big(-C_0\frac{\sin \lambda_n r}{r}\Big)dr$$

Finally, back-substituting we get:

$$C(r,t)=C_0+u(r,t)$$

Note: a different and probably simpler solution caould be obtained by using the second BC:

$$C(R,t)=0$$

Physically this corresponds to a sphere constantly being washed with fresh ($C=0$) solvent.

in order to make the boundary and initial conditions that I specified valid.

That's putting the horse before the cart: you need to choose the right boundary conditions first. Your ODE will then solve correctly.

And I'm not sure you've chosen the right BCs:

due to symmetry reasons. In my opinion, at the end, this translates into the boundary condition

$C(r=0,t) = \dot{r}_{g}\cdot t$

Not necessarily. You seem to have simply integrated:

$\frac{\partial C(0,t)}{\partial t} = \dot{r}_{g}$

But have omitted an integration constant.

I think you also lack a BC for $r=R$. Often in diffusion problems we set:

$$\Big(\frac{\partial C(r,t)}{\partial r}\Big)_{r=R}=0$$

Physically this means the solute cannot diffuse beyond $r=R$.

• Thank you Gert for the answer. What I am thinking now, is that probably I will never be able to get a finite concentration at r = 0. Since this is a sphere and at the middle point the area is zero, concentration has to be infinite in order to have some flux: J = A.dC/dr – Marcelo Jan 6 '17 at 18:40
• One way to avoid that is to specify a small central area ($r \leq r_0$) where $C=C_0$ (constant). Then change to a new variable: $u=C-C_0$, this way the boundary condition is homogeneous. – Gert Jan 6 '17 at 19:37
• Your initial condition becomes $u(r,0)=-C_0$. When all is done and dusted your solution will be: $C(r,t)=C_0+u(r,t)$. – Gert Jan 6 '17 at 20:02
• I was just looking around on the internet and what I found closest to the solution I was looking for is the Green's function that satisfies: $\frac{\partial C}{\partial t} - k(\frac{\partial^{2} C}{\partial x^{2}}+\frac{\partial^{2} C}{\partial y^{2}}+\frac{\partial^{2} C}{\partial z^{2}}) = \delta(r-r_{0})$. The solution to this equation involves the gaussian function (en.wikipedia.org/wiki/Green's_function). However, that solution is for a non-bounded domain and for a point source that happens at t = 0 – Marcelo Jan 10 '17 at 10:19
• What I don't know is how to specify that $\frac{\partial C}{\partial r} = 0$ for $r = R$ (because diffusion does not proceed further that $r=R$) and for $r = 0$ (if function is C1 and for symmetry reasons). I believe this is done with the method of images. And also how to specify that point source is continuous. – Marcelo Jan 10 '17 at 10:26