Steady state heat equation in a cylinder

[fluidistic]

Homework Statement

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.

Homework Equations

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

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.

 

[tadchem]

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!

 

[pasmith]

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

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

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.

 

[fluidistic]

First of all thanks guys.

You have
<br /> R(\rho) = J_m(k\rho)<br />
so your boundary condition requires that J&#039;_m(ka) = 0. There are a countable number of such zeroes, so you can label them as x_{mn} for n \geq 1.
Thus
<br /> R(\rho) = J_m(x_{mn}(r/a)).<br />

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:

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.

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:

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!

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

 

[pasmith]

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

To add to that: mathematical intuition tells me to expect
<br /> T(r,\theta,z) = \alpha \sin(\theta) \sum_{n=1}^{\infty} A_n J_1(k_nr/a) Z_n(z)<br />
where k_n is the nth positive zero of J_1&#039;(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:
<br /> \sum_{n=1}^{\infty} A_n J_1(k_nx) = 1<br />
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.

 

[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.

 

[fluidistic]

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.

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...

 

[fluidistic]

To add to that: mathematical intuition tells me to expect
<br /> T(r,\theta,z) = \alpha \sin(\theta) \sum_{n=1}^{\infty} A_n J_1(k_nr/a) Z_n(z)<br />
where k_n is the nth positive zero of J_1&#039;(x) and
Z_n(z) = \cosh(k_n z) - \coth(k_n L)\sinh(k_n z)

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.

 

[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.

 

[pasmith]

To add to that: mathematical intuition tells me to expect
<br /> T(r,\theta,z) = \alpha \sin(\theta) \sum_{n=1}^{\infty} A_n J_1(k_nr/a) Z_n(z)<br />
where k_n is the nth positive zero of J_1&#039;(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.

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

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.

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.

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.

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.

 

[pasmith]

I used the relation at http://en.wikipedia.org/wiki/Bessel_function for the orthogonality.

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:

<br /> \int_0^1 xJ_p(k_nx) J_p(k_mx)\,\mathrm{d}x = \delta_{nm}<br /> \frac12(1 - p^2k_n^{-2})J_p^2(k_n)<br />
where k_n and k_m are zeroes of J_p&#039;(x).

 

[fluidistic]

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:

<br /> \int_0^1 xJ_p(k_nx) J_p(k_mx)\,\mathrm{d}x = \delta_{nm}<br /> \frac12(1 - p^2k_n^{-2})J_p^2(k_n)<br />
where k_n and k_m are zeroes of J_p&#039;(x).

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.