# Steady state heat equation in a cylinder

1. Dec 7, 2012

### fluidistic

1. The problem statement, all variables and given/known data
I'm unable to solve a problem of heat equation in a cylinder in steady state. The problem is a cylinder of radius a and a height L. The boundary condition are $T(\rho , \theta , 0)=\alpha \sin \theta$, $T(\rho, \theta , L)=0$ and $\frac{\partial T}{\partial \rho} (a, \theta , z)=0$.
I must calculate the temperature inside the cylinder.

2. Relevant equations
The temperature $T(\rho, \theta , z)$ satisfies the Laplace equation inside the cylinder. Namely $\triangle T=0$.

3. The attempt at a solution
I solved the Laplace equation in cylindrical coordinates using separation of variables looking for solutions of the form $T=R(\rho ) \Theta (\theta) Z(z)$.
This produced 3 ODE's, namely:
$$\begin{array} \frac{Z''}{Z}=k^2 \\ \frac{\Theta '' }{\Theta } =-m^2 \\ \rho ^2 R'' + \rho R' + R(k^2 \rho ^2 - m^2) \end{array}$$. Solving the 3 equations, I got the eigenfunctions of the problem, they have the form $T_{m,k}=J_m (k \rho )[C_m \cos m \theta + D_m \sin m \theta ](A_k e^{kz}+B_k e^{-kz})$.

However the function that satisfies the boundary conditions is a linear combination of these eigenfunctions.
Namely $T(\rho , \theta , z)=\sum _{m=0}^\infty \sum _{k=0}^\infty J_m (k \rho ) [C_m \cos (m \theta )+ D_m \sin (m \theta )](A_k e^{kz}+B_k e^{-kz})$.

So far so good. Now I apply the 2nd boundary condition. This gives me that $A_k=-B_ke^{-2kL}$ so that T simplifies to $T(\rho , \theta , z )=T(\rho , \theta , z)=\sum _{m=0}^\infty \sum _{k=0}^\infty J_m (k \rho ) [C_m \cos (m \theta )+ D_m \sin (m \theta )][B_k (e^{-kz}-e^{k(z-2L)})]$.

Now I apply the first boundary condition. This gives me $\sum _{k=0}^\infty \sum _{m=0}^\infty E_k J_m (k\rho )[C_m \cos (m \theta ) +D_m \sin (m\theta )]=\alpha \sin \theta$. Where $E_k=B_k(1-e^{-2kL})$, is a constant for a given k.

So I'm looking for the constant $C_m$ and $B_m$'s. By inspection it's obvious that $C_m=0 \forall m$ and $D_m =0 \forall m \neq 1$.

So I am left with $\sum _{k=0}^\infty E_k J_1 (k \rho ) D_1 \sin \theta = \alpha \sin \theta$. This gives me that $D_1= \frac{\alpha}{\sum _{k=0}^\infty E_k J_1 (k \rho )}$.
Which is impossible since $D_1$ is a constant and MUST NOT depend on rho, while here I reached that it does depend on rho. I've absolutely no idea where I went wrong. Any pointer would be immensily appreciated.

2. Dec 7, 2012

Let k^2 be negative. Then k is a pure imaginary number, and exp(kx) and exp(-kz) become periodic functions.
See where that takes you: the Magical Realm of Bessel Functions!

3. Dec 7, 2012

### pasmith

You have
$$R(\rho) = J_m(k\rho)$$
so your boundary condition requires that $J'_m(ka) = 0$. There are a countable number of such zeroes, so you can label them as $x_{mn}$ for $n \geq 1$.

Thus
$$R(\rho) = J_m(x_{mn}(r/a)).$$

But you're missing a family of eigenfunctions: z is bounded, you can take Z(z) to be cos or sin as well as exp. If you do, you'll find that R will be a Modified Bessel function.

4. Dec 7, 2012

### fluidistic

First of all thanks guys.
I see but how does this notation helps me? I express a bessel function of the m'th kind with an argument depending on the zero of a Bessel function of another kind.
With respect to:
I don't think I can take those. Remember, $T_1$ is the stationary solution with the upper top kept at $T=0$ and the bottom at $T=\alpha \sin \theta$.
If I wait forever, I DO NOT expect the solution to be sinusoidally dependent on z. I'd expect something exponential for sure.

Edit: Also I do not see at all where I went wrong. All my steps seemed justified to me.

Edit 2:
That would be the way to go for when I consider T_2, that is the time dependent solution to the PDE. In that case I agree, I will have a periodic function for Z(z) since it must vanish at both z=0 and z=L.
Also notice that I already entered the realm of Bessel function, in case you missed it. :P

Last edited: Dec 7, 2012
5. Dec 8, 2012

### pasmith

Checking a handbook of integrals yields
$$\int x J_n(ax) J_n(bx)\,\mathrm{d}x = \frac{x(aJ_n(bx)J_n'(ax) - bJ_n(ax)J_n'(bx))}{b^2 - a^2},\qquad a \neq b \\ \int x J_n^2(ax)\,\mathrm{d}x = \frac{x^2}{2} (J_n'(ax))^2 + \frac{x^2}{2}\left(1 - \frac{n^2}{a^2x^2}\right)(J_n(ax))^2$$
Taking a and b as zeroes of $J_n'$ and integrating from 0 to 1 gives the following orthogonality relation:
$$\int_0^1 xJ_n(ax) J_n(bx)\,\mathrm{d}x = \left\{\begin{array}{r@{\quad:\quad}l} 0 & a \neq b \\ \frac12(1 - n^2a^{-2})J_n^2(a) & a = b \end{array}\right.$$

To add to that: mathematical intuition tells me to expect
$$T(r,\theta,z) = \alpha \sin(\theta) \sum_{n=1}^{\infty} A_n J_1(k_nr/a) Z_n(z)$$
where $k_n$ is the nth positive zero of $J_1'(x)$ and
$$Z_n(z) = \cosh(k_n z) - \coth(k_n L)\sinh(k_n z)$$
so that $Z_n(0) = 1$ and $Z_n(L) = 0$. On further reflection I don't think eigenfunctions of the form $I_1(kr)\exp(ikz)$ actually add anything you don't already have.

That's basically what you got. But you didn't do the further step of working out the $A_n$.

So on z = 0:
$$\sum_{n=1}^{\infty} A_n J_1(k_nx) = 1$$
where $x = r/a$.

The method for determining the $A_n$ is to multiply both sides by $xJ_1(k_mx)$ and integrate from x = 0 to 1. Then you can apply the orthogonality condition above.

Last edited: Dec 8, 2012
6. Dec 8, 2012

### Dustinsfl

I am a little confused about your R equation. If we are solving the ss heat, i.e. the Laplace equation in polar coordinates, we should have $\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 }{\partial \theta^2} = 0$, correct?

Let's assume the solution is of the form
$$u(r,\theta,t) = R(r)e^{\pm i\omega t}e^{\pm im\theta}.$$
Using the Laplace in polar coordinates, we obtain
$$\frac{1}{r}R'e^{\pm i\omega t}e^{\pm im\theta} + R''e^{\pm i\omega t}e^{\pm im\theta} - \frac{m^2}{r^2}Re^{\pm i\omega t}e^{\pm im\theta} =0\Rightarrow r^2R''+rR'-m^2R = 0$$
In order to get your R equation, we would have to solve the Helmholtz equation. But that arises in the solution to the wave equation of a circular or annular membrane.

7. Dec 8, 2012

### fluidistic

No, it's a cylinder and there's a dependency on z. So it's not polar coordinates that we must use but cylindrical coordinates.

And also it's clearly stated in the problem statement that it's the steady state temperature that I'm looking for, so that there's no dependency on t...

8. Dec 9, 2012

### fluidistic

Hi and thanks for helping me, but how do I justify this?
I had $T(\rho , \theta , z)=\sum _{m=0}^\infty \sum _{k=0}^\infty J_m (k \rho ) [C_m \cos (m \theta )+ D_m \sin (m \theta )][B_k (e^{-kz}-e^{k(z-2L)})]$ where in fact k should start from 1 instead of 0.
So I am at the point of $T(\rho , \theta , z)=\sum _{m=0}^\infty \sum _{k=1}^\infty J_m (k \rho ) [C_m \cos (m \theta )+ D_m \sin (m \theta )][B_k (e^{-kz}-e^{k(z-2L)})]$
I don't understand how you get rid of a constant of separation. Your $k_n$ looks like my k. Hmm I think you collapsed the infinite series of m and only 1 term (m=1) survived. I'll have to think about it.

Edit: I think I see my error THANKS TO YOU!!! k has no reason to be an integer!!! You pointed out that it's worth $x_{mn}/a$. I'll try to advance further.

Last edited: Dec 9, 2012
9. Dec 9, 2012

### fluidistic

I'll try to keep my notation up to the end but using your ideas pasmith.
I denote $k_{mn}=x_{mn}/a$ where $x_{mn}$ is the nth zero of $J_m'$, the derivative of the Bessel function of the first kind of order m.
I reach therefore that $T(\rho ,\phi , z)= \sum _{m=0}^\infty \sum _{n=1}^\infty J_m (k_{mn} \rho ) [C_m \cos m \phi + D_m \sin m \phi ][B_{km}(e^{-k_{mn}z}-e^{k_{mn}(z-2L)})]$.
Now I apply the first boundary condition. Thus $T(\rho ,\phi , 0)= \sum _{m=0}^\infty \sum _{n=1}^\infty J_m (k_{mn} \rho ) [C_m \cos m \phi + D_m \sin m \phi ][\underbrace{B_{km}(1-e^{-2Lk_{mn}})}_{E_{mn}}]=\alpha \sin \phi$.
By inspection $C_m=0$, $D_m=0$ for all $m \neq 1$. This makes that $T(\rho , \theta , z) = \sum _{n=1}^\infty J_1 (k_{1n} \rho ) \sin (\theta ) F_n e^{-k_{1n}z}$ where $F_n=D_1E_{1n}$.
Applying the boundary condition I should get the $F_n$'s, as you said, using the orthogonality of the Bessel function but I'm having a trouble.
I used the relation at http://en.wikipedia.org/wiki/Bessel_function for the orthogonality.
I get that $\sum _{n=1}^\infty F_n \underbrace{\int_0^ a \rho J_1 (k_{1n} \rho ) J_1 (k_{1l}\rho )d \rho }_{\frac{\delta _{nl}}{2}[J_1'(k_{1n})]^2} =\int _0 ^ a \alpha J_1 (k_{1n} \rho ) \rho d\rho$.
So that $\sum _{n=1}^ \infty F_n \frac{1}{2}[J_1' (k_{1,n})]^2=\alpha \int _0^a J_1 (k_{1n} ) \rho d\rho$. I'm stuck here, how do I get the $F_n$'s? I multiplied by rho and the Bessel function and integrated from 0 to a the radius (I did not normalize as you did, but that should still work).
I'd be glad to know how to get the $F_n$'s. Thanks guys so far.

10. Dec 9, 2012

### pasmith

I should of course have
$$Z_n(z) = \cosh(k_n z/a) - \coth(k_n L/a)\sinh(k_n z/a)$$

Yes. The fourier series for $\alpha\sin(\theta)$ was obviously going to be $\alpha\sin(\theta)$, so there was no need keep track of an unnecessary index.

You may like to take note of the orthogonality relation in my previous post; I don't think I've seen an orthogonality relation for eigenvalues based on zeroes of derivatives of Bessel functions in any methods textbook, even though the relation for eigenvalues based on zeroes of functions is standard.

11. Dec 9, 2012

### pasmith

That's the standard relation for eigenvalues based on zeroes of functions, not zeroes of derivatives, which is the one you want and the one I posted earlier.

For clarity:

$$\int_0^1 xJ_p(k_nx) J_p(k_mx)\,\mathrm{d}x = \delta_{nm} \frac12(1 - p^2k_n^{-2})J_p^2(k_n)$$
where $k_n$ and $k_m$ are zeroes of $J_p'(x)$.

12. Dec 9, 2012

### fluidistic

I used $\int_0^1 x J_\alpha(x u_{\alpha,m}) J_\alpha(x u_{\alpha,n}) dx = \frac{\delta_{m,n}}{2} [J_{\alpha+1}(u_{\alpha,m})]^2 = \frac{\delta_{m,n}}{2} [J_{\alpha}'(u_{\alpha,m})]^2,\!$
So I have $\sum _{n=1}^ \infty F_n \frac{1}{2}[J_1' (k_{1,n})]^2=\alpha \int _0^a J_1 (k_{1n} ) \rho d\rho$. Maybe I'm off by a factor of $a^2$ or so on the left side. However I do not know how to isolate the $F_n$'s from the expression.
And once I have the $F_n$ I'm done since $T(\rho , \theta , z) = \sum _{n=1}^\infty J_1 (k_{1n} \rho ) \sin (\theta ) F_n e^{-k_{1n}z}$.
Edit: thanks for the clarifications in post #10.