# Solving an IVBP for a metal rod dipped from cold to warm water

I have to solve the diffusion equation for a solid rod, which is in a bath of 100 degrees. At $$t=0$$ it is moved to a second bath, where the water temperature is 0 degrees.

I prepare the diffusion equation:

$$$$u_t-\alpha\Delta=0.$$$$

Boundary and initial conditions:

$$$$u(r,\theta,t)=0 \ \ \ \ 0<\phi<\pi, t>0 \\ u(r,\phi,0)=100 \ \ \ \ 0

Since the function $$u$$ is of three variables $$u=R\Theta T$$, I prepare the given partial differential equation with the differential operator on a circle (the cross section of the rod):

$$$$T_tR\Theta-\alpha\bigg[R_{rr}\Theta T+\frac{1}{r}R_r\Theta T+\frac{1}{r^2}\Theta_{\theta\theta}RT\bigg]=0$$$$

Divide by $$R\Theta T$$ and obtain

$$$$\frac{T_t}{T}-\alpha\bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=0$$$$

Forming one ODE and one PDE:

$$$$\frac{T_t}{T}=-\lambda^2 \\ \bigg[\frac{R_{rr}}{R}+\frac{1}{r}\frac{R_r}{R}+\frac{1}{r^2}\frac{\Theta_{\theta\theta}}{\Theta}\bigg]=\frac{\lambda_n^2}{\alpha}$$$$

The first gives $$T(t)=Ce^{-\lambda^2t}$$, the second partial differential equation gives two ordinary differential equations:

$$$$r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}+\frac{\Theta_{\theta\theta}}{\Theta}=\frac{\lambda_n^2}{\alpha}r^2$$$$

which can be further split into two ODEs which must equal some constant $$\mu^2$$, which must be negative to give Bessel solutions:

$$$$r^2\frac{R_{rr}}{R}+r\frac{R_r}{R}-\frac{\lambda_n^2}{\alpha}r^2=\frac{\Theta_{\theta\theta}}{\Theta}=-\mu_k^2$$$$

this gives the Bessel equation:

$$$$R_{rr}+\frac{1}{r}R_r+\bigg(\frac{\lambda_n^2}{\alpha}-\frac{\mu_k^2}{r^2}\bigg)R=0$$$$

and the second-order ODE:

$$$$\Theta_{\theta\theta}+\mu_k^2\Theta=0$$$$

The latter has the solution $$\Theta(\theta)=A\sin\mu_k\theta+B\cos\mu_k\theta$$.

The former is the Bessel equation which gives the Bessel functions:

$$$$R(r)=\bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_n r}{\sqrt{\alpha}}\bigg)+bY_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)$$$$

So if this procedure is right, the general solution would be (ignoring the Bessel function of the second kind:

$$$$u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\bigg(\big(A\sin\mu_k\theta+B\cos\mu_k\theta\big)\big(e^{-\lambda_n^2t})\bigg)$$$$

Then we set the angular part equal to some constant, which is merged with the Bessel constant $$a$$ and get

$$$$u(r,\theta,t)=\sum_{n=1}^\infty \bigg(aJ_{\mu_k}\bigg(\frac{i\lambda_nr}{\sqrt{\alpha}}\bigg)\bigg)\big(e^{-\lambda_n^2t})\bigg)$$$$

Using ICs to find the constants, the zeros of the Bessel are given by $$\alpha_{n,k}$$ and since the core will still be at 100 degrees at t=0, while the surface will be 0, then $$R(0)=100$$ and $$R(R)=0$$, and with the radial function given by:

$$$$R(r)=aJ_{\lambda n}(i\sqrt{\mu}r)$$$$

we solve: $$$$\alpha=aJ_{\lambda n}(i\sqrt{\mu}R)$$$$

and

$$$$100=aJ_{\lambda n}(i\sqrt{\mu}\alpha)$$$$

which gives:

$$$$u(r,t)=\frac{\alpha J_{\lambda_n}(i\sqrt{\mu}r)}{J_{\lambda_n}\big(\frac{100R}{a\alpha}\big)}e^{-\lambda^2t}$$$$

The plot of the time-dependent evolution looks like this (with x=r) and x=0=R (surface) and t goes from 0 to 100. (ignore minus signs)

Does this entire procedure appear reasonable? Thanks

• Comments are not for extended discussion; this conversation has been moved to chat. Commented May 23, 2022 at 22:17

First let $$u=X(\mathbf{x}) T(t)$$. The equation reads $$T_t X - \alpha \Delta(X) T = 0$$. Take the cylindrical Laplacian and kill the $$z$$ and $$\theta$$ derivatives from symmetry considerations, so that $$X(\mathbf{x})=R(\| \mathbf{x} \|)$$; let $$r=\| \mathbf{x} \|$$.

Now you have

$$T_t R - \alpha T \left ( R_{rr} + r^{-1} R_r \right ) = 0.$$

Doing the usual separation of variables trick, you conclude that $$\frac{T_t}{T}=\alpha \frac{R_{rr} + r^{-1} R_r}{R}=\lambda_n$$ for the same constant $$\lambda_n$$. It follows that up to a constant factor that we'll sort out later, $$T=e^{\lambda_n t}$$.

Now in the radial part we have

$$R_{rr} + r^{-1} R_r - \frac{\lambda_n}{\alpha} R = 0.$$

This is actually an adjusted Bessel equation itself. Multiply through by $$r^2$$:

$$r^2 R_{rr} + r R_r + \left ( -\frac{\lambda_n}{\alpha} r^2-0^2 \right ) R = 0.$$

That $$-\frac{\lambda_n}{\alpha}$$ wouldn't be there in the ordinary Bessel equation. The fix comes because the first two terms are invariant under a linear change of variable, so we can replace $$s=\sqrt{-\lambda_n/\alpha} r$$ (dropping the $$n$$ as an abuse of notation) and then we have

$$s^2 R_{ss} + s R_s + \left ( s^2 - 0^2 \right ) R = 0.$$

This gives $$R=c_{1,n} J_0(s) + c_{2,n} Y_0(s)$$, Bessel functions of order zero. For nonsingularity at the origin we force $$c_{2,n}=0$$.

So now the question is about what $$c_{1,n}$$ is. The big catch is that in order for $$J_0(\sqrt{-\lambda_n/\alpha} r)$$ to have zeros, we must have $$\sqrt{-\lambda_n/\alpha}$$ real. (This also makes intuitive sense as it implies exponential decay.)

As a result you choose $$\lambda_n$$ so that $$\sqrt{-\lambda_n/\alpha} R$$ is a zero of $$J_0$$ (where $$R$$ is now the radius of the rod). We can pick the zero of interest to be $$j_{0,n}$$ which is where $$n$$ enters into the picture.

Now we solve that for the eigenvalue, and get $$\lambda_n=-\frac{\alpha j_{0,n}^2}{R^2}$$ which gives

$$u=\sum_n c_n J_0 \left ( j_{0,n} \frac{r}{R} \right ) e^{-\frac{\alpha j_{0,n}^2}{R^2} t}.$$

Now it remains to find the $$c_n$$ to match the initial profile.

• Thanks Ian, I will study this on Monday! Have a good weekend. PS: Upscore my post! Commented May 21, 2022 at 15:41
• @Luthier415Hz I fixed the mistake, it was partially my fault; you can in fact change the radial eigenvalue while keeping the angular eigenvalue frozen. You were moving them in lockstep and that confused me when it shouldn't have.
– Ian
Commented May 21, 2022 at 15:48
• No problem. I think that we have reached some knowledge here anyhow. Commented May 21, 2022 at 16:17