Physics Stack Exchange is a question and answer site for active researchers, academics and students of physics. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to solve the transient heat equation on a ball. The boundary condition is the same over the hole outer surface. So this should reduce to a one-dimensional problem in radial direction. However I cannot use the one-dim heat equation, since the surface through which the heat flows goes quadratic with the radius.

Is the following approach correct: Take the heat equation, transform it into sperical coordinates and eliminate the derivatives in angular directions. I think the last part is justified, because of the rotational symmetry of the boundary condition. So I end up with this equation:

$$ \frac{\partial T}{\partial t} = \alpha \frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2\frac{\partial T}{\partial r}\right) $$ with thermal diffusivity $\alpha$. Is this the correct reformulation of the heat equation on a ball to a one-dimensional problem?

EDIT: When I solve this with the finite element method, will there be a problem in the centre of the ball? The $1/r^2$ factor might cause a problem.

share|cite|improve this question
I suggest the substitution: $u=rT(r,t)$. Then the equation will simplify:$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2u}{\partial r^2}$ – Martin Gales Apr 28 '11 at 10:29
up vote 3 down vote accepted

That sounds about right to me. Although I've never done this for the heat equation specifically, that's much the same approach we use to solve the Schrödinger equation in spherical coordinates, e.g. for the hydrogen atom. Your full solution will be products of the radial functions $T(r,t)$ with appropriate spherical harmonics.

When determining the radial part of the solution, assuming this goes like the Schrödinger equation, one of the boundary conditions you'll need to use is consistency at the origin: basically you need your solution for $T(r)$ to be finite and differentiable at $r=0$. I can come back and expand on this later, or perhaps someone else will give you a more detailed explanation before I get to it.

share|cite|improve this answer
OK, thanks so far! But since I have full rotational symmetry, will the full solution really involve spherical harmonics? – Till B Apr 28 '11 at 9:26
@Till B: Any short answer to that will be misleading. It will not involve spherical harmonics other than the constant Y_00 one, but that's still a spherical harmonic. – wnoise Apr 28 '11 at 9:31
@TillB Yeah, but your symmetry assumption would anyway kill everything except $Y_0^0$. – mbq Apr 28 '11 at 9:37
Yeah, if your boundary condition is spherically symmetric, you're restricted to the spherically symmetric harmonic. In general, the coefficients of each of the spherical harmonics in your solution will be determined by the external boundary condition. – David Z Apr 28 '11 at 14:24

your calculation is fine but there is no need to employ numerical methods for the given equation at hand. The equation can be solved using separation of variables with a physical ansatz for the radial part. Of course you need not restrict your problem to a radial symmetry of the solution, the full Laplace will lead to spherical harmonics for the angular parts. The following derivation can directly be applied in this case.

Solution using separation of variables

Taking $$T(r,t) = f(t)\cdot B(r)$$

results in

$$f'(t) B(r) = \alpha f(t)\frac{1}{r^2}\partial_r \left(r^2 B'(r) \right)$$

Now we see that we can re-write this as

$$\frac{f'(t)}{f(t)} = \alpha \frac{1}{r^2}\frac{\partial_r \left(r^2 B'(r) \right)}{B(r)}$$

and see that both sides are independent of each other because of the different dependencies, hence they must be a constant, say $-\lambda$, for the left hand side $$\frac{f'(t)}{f(t)} = -\lambda$$ for which we find the solution $$f(t) = f_0 e^{-\lambda t}\ .$$

$\lambda$ can be seen as a dissipation constant defining the loss of heat per time unit. I found a nice illustration of this process on youtube or on wikipedia for another geometry:

wiki heat dissipation

Now to the radial part, $$\alpha \frac{1}{r^2}\partial_r \left(r^2 B'(r) \right) = -\lambda B(r)$$

for which we find

$$B(r) = c_1 \frac{e^{-\mathrm{i}\kappa r}}{r} + c_2\frac{e^{\mathrm{i}\kappa r}}{r} $$ with $$\kappa = \sqrt{\lambda/ \alpha}\ .$$

Q: Can you tell me what kind of solutions we just arrived at?

Application to a spherical body

Now this was the general solution to the problem but we have to apply it to your situation. First of all, wee need some initial thermal distribution: $$T_0(r) = T(r,0)$$ to determine the dynamics from this point. Second of all, we need to know what it means to have some sphere in our system. We don't have a choice but to assign some distribution to $\alpha$, determining different heat coefficients in different regions, e.g. $$\alpha = \begin{cases} \alpha_{i} & r\leq R\\ \alpha_{0} & r>R\end{cases} $$

The only thing left to do is to determine the constants $\lambda$ and $c_i$ from the continuity of $T$ at $r=R$ and meaningful values at $r = \infty$ ($T\rightarrow 0$) and $r = 0$ ($T$ finite).

I will leave this step to you.
Having the solution at hand you can compare it to your simulation results and further extend it, e.g. to a system of shells or special boundary conditions.


share|cite|improve this answer
Hi Robert! Thank you for your long and detailed answer. Actually I know the separation of variables ansatz. However my boundary condition is time-dependent and most important it comes itself out of another numerical computation. So I do not have an analytical expression of my BC and this is why I use a numerical method to solve it. – Till B Apr 28 '11 at 12:41

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.