Bessel Functions / Eigenvalues / Heat Equation

[phioder]

Hello

Trying to calculate and simulate with Matlab the Steady State Temperature in the circular cylinder I came to the book of Dennis G. Zill Differential Equations with Boundary-Value Problems 4th edition pages 521 and 522

The temperature in the cylinder is given in cylindrical coordinates by:

u(r,z)= u_0 [Sum from n=1 to Infinite] of:


sinh( lambda_n*z ) * J_0( lambda_n * r )
________________________________________________

lambda_n * sinh(4 * lambda_n) * J_1(2 *lambda_n}


My problems:
-I don't understand very well the Bessel Function either the Eigenvalues and need a bit of help
-PDE Knowledge and simulations is basic

Information:
With the separation of variables method in cylindrical coordinates and having U as temperature the equations are defined as follows:

Initial Conditions:
u(2,z)=0 0<z<4
u(r,0)=0 0<r<2

Boundary Condition:
u(r,4)=u_0 0<r<2

u=R(r)Z(z)

r*R'' + R' + ((lambda)^2)*r*R = 0 Cauchy-Euler equation

Z'' + 0 - ((lambda)^2) * Z = 0

With solutions:

R = c_1 * J_0 ( lambda * r ) + c_2 * Y_0 (lambda * r)
Z = c_3 * cosh( lambda * z ) + c_4 * sinh(lambda * z)

The book states "the assumption that the function u is bounded at r = 0 demands that c_2 = 0"

The condition u(2,z) = 0 implies that R(2) = 0

The positive eingenvalues lambda_n of the problem are defined by:

J_0(2*lambda)=0

Now I come to my questions:
1.- What is meant by "the function u is bounded at r = 0"?

Is it right to understand that c_2 = 0 because the Bessel Function of the Second Kind of Order Zero (Y_0) tends to minus infinite while aproaching to r=0 from the right side, what is meant by bounded at r=0?


2.- I did some research on the Bessel Functions of the First and Second Kinds, solved the Bessel equation step by step and "more or less" understood it. My problem is that I don't understand neither how to calculate the eigenvalues lambda_n of the steady state temperature in a circular cylinder.


Does the equation J_0(2*lambda) = 0 means that:

2*lambda_{1} = 2.4048
2*lambda_{2} = 5.5201
2*lambda_{3} = 8.6537
.
.
.
lambda_{1} = 2.4048 / 2 ?

Or in words said: The eigenvalues are defined by the division by two of the x value where J_0 is a zero or a root?


3.- If we go back to the final solution there are two terms a J_0(lambda_n*r) and a J_1(2*lambda_n) and my goal is to implement this terms on Matlab to understand better the temperature U.

So my approach would continue trying to define in Matlab a vector for J_0(lambda_n * r), is it right to think that having two vectors of the same size lambda_n and r, being r defined from 0 to 2, find out which is the value of the bessel function J_0 at say J_0( (2.4048/2)*r )?

Unfortunately I can't write my post in a less complex way, hope it is understood, any help, hint or tip would be kindly appreciated.

Best Regards

 

[coomast]

Hello phioder, first of all welcome to this forum.

It took me a couple of hours to read again my notes and look into the book on Bessel functions. I don't use it every day. The work that you are doing seems to be all right for solving the equation. However you stated that you "more or less" understand it, so I will describe the solution in total in order that you might gain more insight. First one needs to state that the solution is indeed (no theta dependence due to symmetry):

$$u(r,z)=R(r)\cdot Z(z)$$

Substituting this into the original differential equation

$$\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+ \frac{\partial^2 u}{\partial z^2} = 0$$

gives:

$$\frac{R''}{R}+\frac{1}{r} \cdot \frac{R'}{R} = - \frac{Z''}{Z} = -\lambda^2$$

giving then two ordinary differential equations:

$$Z''-\lambda^2Z=0$$

$$r^2R''+rR'+\lambda^2r^2R=0$$

The first one is easily solved as:

$$Z(z) = A cosh(\lambda z)+B sinh(\lambda z)$$

The second one is not a Cauchy-Euler equation, but the one of Bessel with solution:

$$R(r) = C J_0(\lambda r)+D Y_0(\lambda r)$$

Because of the diverging behavior of the Bessel function of the second kind, you must set D equal to 0 because you would end up with a non-physical solution. That is what is meant by "the solution needs to be bounded". The principal solution is now all ready (C is absorbed by A and B):

$$u(r,z)=\left[A cosh(\lambda z)+B sinh(\lambda z)\right] \cdot J_0(\lambda r)$$

With the boundary condition:

$$u(2,z)=0=\left[A cosh(\lambda z)+B sinh(\lambda z)\right] \cdot J_0(2\lambda)$$

we get:

$$J_0(2\lambda)=0$$

and the solutions to this equation are now called:

$$\lambda_n$$

for all n=1,2,3,...

The other boundary condition:

$$u(r,0)=0=A \cdot J_0(\lambda r)$$

From which we obtain A=0. The principal solution is now:

$$u_n(r,z)=B_n \cdot sinh(\lambda_n z) \cdot J_0(\lambda_n r)$$

Because this is valid for all n and superposition is valid, we get for the total solution the following:

$$u(r,z)=\sum_{n=1}^{\infty}B_n \cdot sinh(\lambda_n z) \cdot J_0(\lambda_n r)$$

The only boundary condition left is the following:

$$u(r,4)=f(r)=\sum_{n=1}^{\infty}B_n \cdot sinh(4\lambda_n) \cdot J_0(\lambda_n r)$$

From the theory of orthogonal functions it can be shown that the following is true:

$$B_n \cdot sinh(4\lambda_n) = \frac{2}{J_1^2(2\lambda_n)} \int_{0}^{2}\frac{r}{2}J_0(\lambda_n r)f(r)d\left(\frac{r}{2}\right)$$

With:

$$f(r)=u_o$$

you can derive:

$$B_n=\frac{u_0}{\lambda_n \cdot sinh(4\lambda_n) \cdot J_1(2\lambda_n)}$$

After substitution this, it gives you the solution you were looking for. In order to obtain the lambda values, you need to look for the following:

$$J_0(2\lambda)=0$$

The values were this function is zero, are tabulated, so the first three values are indeed:

$$\lambda_1 = \frac{2.4048}{2}$$

$$\lambda_2 = \frac{5.5201}{2}$$

$$\lambda_3 = \frac{8.6537}{2}$$

The final solution for the problem given is now:

$$u(r,z) = u_0 \cdot \sum_{n=1}^{\infty}\frac{sinh(\lambda_n z) \cdot J_0(\lambda_n r)}{\lambda_n \cdot sinh(4 \lambda_n) \cdot J_1(2 \lambda_n)}$$

I can't help you on the MATLAB problems, I'm not familiar with it. However if I was to create a program, I would make an array of an array of reals (2D) ranging from 0 to 2 and from 0 to 4 for the two coordinates. The number of points is to your liking. This for every term of the sum, as much as you need them. Then it shouldn't be too hard to display the results.

Hope this helped you a bit further.

 

[phioder]

Hello Coomast

Thank you very much for your welcome message and answer. I read some of the forum threads and find the forum
really great and friendly. Thank you very much to review all my post and for writing it in a clear form.

Unfortunatelly I'm still having some problems understanding the bessel fuctions

After substitution this, it gives you the solution you were looking for. In order to obtain the lambda values, you need to look for the following:

$$J_0(2\lambda)=0$$

The values were this function is zero, are tabulated, so the first three values are indeed:

$$\lambda_1 = \frac{2.4048}{2} = 1.2024$$

$$\lambda_2 = \frac{5.5201}{2}$$

$$\lambda_3 = \frac{8.6537}{2}$$

What I understand is that to get the values of each term of the sumatory then the lambda values have to be substituted in the final solution equation

For:

$$\lambda_1 = 1.2024$$

and

$$n = 1$$

The first term of the sumatory would be:

$$\frac{sinh(1.2024 \cdot z) \cdot J_0(1.2024 \cdot r)}{1.2024 \cdot sinh(4 \cdot 1.2024) \cdot J_1(2 \cdot 1.2024)}$$

To calculate this value one needs

$$J_0(1.2024 \cdot r)}$$

and

$$J_1(2 \cdot 1.2024)}$$

Now to the questions for this specific case is it right to think of the values as:

1.- The values of the function [tex]J_1[\tex] at [tex]x=2.4048[\tex]
2.- The values of the function [tex]J_0[\tex] at [tex]x=1.2024 \cdot r [\tex] for all [tex]r[\tex]s between 0 and 2

Best Regards

 

[coomast]

Hello phioder,

You have written the first term of the series solution with the value of

$$\lambda_1=1.2024$$

which is correct. Now in order to calculate the following:

$$J_1(2 \cdot 1.2024)=J_1(2.4048)$$

you need a table or a computer code. It is approximately 0.52... For the other function:

$$J_0(1.2024 \cdot r)$$

you need to put in the value of r as you stated and then obtain the function value again with tables or coding. If this is what you meant (post was a bit awkward) it is OK.
Remember that every term of the series has values in the two ranges for r and z. Somehow I suspect you will need a 2D array for every term. Once you have p.e. 2D arrays for the first 10 terms individually, you can display them individually as they are and study the individual behavior of the terms or as a sum for the total solution. If you are only interested in the total solution, you will suffice with one 2D array. In calculating every term, you can just sum them with the previous value in the array.

 

[phioder]

you need to put in the value of r as you stated and then obtain the function value again with tables or coding. If this is what you meant (post was a bit awkward) it is OK.
Remember that every term of the series has values in the two ranges for r and z. Somehow I suspect you will need a 2D array for every term. Once you have p.e. 2D arrays for the first 10 terms individually, you can display them individually as they are and study the individual behavior of the terms or as a sum for the total solution. If you are only interested in the total solution, you will suffice with one 2D array. In calculating every term, you can just sum them with the previous value in the array.

Hello coomast,

Thank you again for your quick answer. Think I'm slowly getting the idea of the Bessel functions. To calculate those tables there is some code at the mathworks webpage

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6794&objectType=FILE

And with the Matlab bessel function it is possible get the values for $$J_0$$ and $$J_1$$The problem now is when evaluating the term $$sinh(4\lambda_n)$$ which values are extremely small and suddenly become very big

At the end I'm expecting to get some similar 3d surfaces as the ones shown here:

http://www.mrl.columbia.edu/ntm/level2/ch03/html/l2c03s05.html

But I think I'm missing something while evaluating the final solution equation or that I don't understand the problem at all, so I will keep trying

Thank you very much and best regards
Phioder

 

[coomast]

Hello phioder,

Indeed, the value for $$sinh(4\lambda_1)=sinh(4\cdot 1.2024)=61.34$$

However if you take the following:

$$\frac{sinh(\lambda_1 z)}{\lambda_1 \cdot sinh(4\lambda_1)}$$

for z=0, this is 0 and for z=4, you have approximately 0.83, giving again a "reasonable" result. It might be useful in calculating this first and then multiply it with the factor of the ratio of the Bessel functions.

It might turn out that after a few terms the solution does not change much anymore.

If anything is uncertain, please post the question, I would like to see a picture of the solution when it is finished.

 

[phioder]

Hello phioder,

Indeed, the value for $$sinh(4\lambda_1)=sinh(4\cdot 1.2024)=61.34$$

However if you take the following:

$$\frac{sinh(\lambda_1 z)}{\lambda_1 \cdot sinh(4\lambda_1)}$$

for z=0, this is 0 and for z=4, you have approximately 0.83, giving again a "reasonable" result. It might be useful in calculating this first and then multiply it with the factor of the ratio of the Bessel functions.

It might turn out that after a few terms the solution does not change much anymore.

If anything is uncertain, please post the question, I would like to see a picture of the solution when it is finished.

Hello Coomast

All the denominator was calculated for n=1,2,3, in the cases for n=2 and n=3 the value is very large, so the sum goes away, the denominator becomes a constant. I also tried to replace the sinh with it's exponential equivalent to reduce the equation. I'm trying to find out how Matlab mesh() works. I'll be back late on Wednesday, need to go now

Best Regards and thank you
phioder

 

[phioder]

The final solution for the problem given is now:

$$u(r,z) = u_0 \cdot \sum_{n=1}^{\infty}\frac{sinh(\lambda_n z) \cdot J_0(\lambda_n r)}{\lambda_n \cdot sinh(4 \lambda_n) \cdot J_1(2 \lambda_n)}$$

Hello

The temperature in the cylinder is given by the above equation, but is dimensionless, what about if one would like to calculate the temperature of cylinders of different materials?

Is there a way to add to that equation, the density, thermal conductivity or other variables that would describe more physically and thermally the cylinder in question? Some type of constant? or the equation needs to be solved completely different?

Any hint or book would be kindly appreciated

Very Best Regards
Phioder

 

[coomast]

Hello

The temperature in the cylinder is given by the above equation, but is dimensionless, what about if one would like to calculate the temperature of cylinders of different materials?

Is there a way to add to that equation, the density, thermal conductivity or other variables that would describe more physically and thermally the cylinder in question? Some type of constant? or the equation needs to be solved completely different?

Any hint or book would be kindly appreciated

Very Best Regards
Phioder

The solution is not exactly dimensionless as such, it can very easily made dimensionless by rewriting it as:

$$\frac{u(r,z)}{u_0}= \sum_{n=1}^{\infty}\frac{sinh(\lambda_n z) \cdot J_0(\lambda_n r)}{\lambda_n \cdot sinh(4 \lambda_n) \cdot J_1(2 \lambda_n)}$$

The unit of $$\lambda_n$$ is $$m^{-1}$$. The unit of $$u$$ is $$K$$. Now to answer the question of different materials, it does not change anything in case you are only considering the steady-state solution. This means that no time dependency is introduced into your solution. To see this consider the general heat equation (without generated heat):

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{1}{r^2} \cdot \frac{\partial^2 u}{\partial \varphi ^2}+ \frac{\partial^2 u}{\partial z^2} \right)$$

If you have a time independent solution as in the original post, this is equal to 0. And therefore the material properties vanish. Also the $$\varphi$$ is left out due to symmetry. This can be seen in a somewhat physical way as well. Indeed if we have established a "steady-state" state, the complete temperature profile in the volume has reached it's final value. Also the heat flow is steady and one can "feel" that this is not influenced by the material properties anymore. Everything is now steady-state. A good book, mmm, there are several books on this subject. It depends a bit on what exactly you are looking for phioder. For the mathematics the one I use sometimes is: "Fourier Analysis with applications to boundary value problems" by Murray R. Spiegel, Schaum's theory and problem solver. I know that people or like or dislike this series very much, I have always found it a nice background for additional learning through exercises. I adore doing calculations and exercises :-) For a more theoretical study of partial differential equations, I can't really advice one, my original text is from university. I can recommend if you are also interested in ordinary differential equations (also here opinions vary): "Elementary Differential Equations and Boundary Value Problems" by Boyce and DiPrima, John Wiley & Sons. If you want to know more about heat transfer, I can recommend the following work: www.artikel-software.com/file/ht.pdf I use this one sometimes at work.

Hope this helps you ahead.

 

[phioder]

Dear coomast,

Thank you for the great and clear explanation and the book recommendation. I will try to solve the equation $$\frac{\partial u}{\partial t}$$ living out $$\varphi$$, steady state is left behind and transient is in the horizon. I will give a try to the book of Spiegel. Next step is to have some fluxes at the boundaries but I haven't found a good book that explains and shows the solution of this kind of problems (Neumann). Carslaw and Jaeger is complex but still helps a bit.

Best Regards
phioder

 

[coomast]

Hello phioder,

Thank you for the nice comment. I can't help you with loking into Carslaw and Jaeger, it isn't at my disposal. If you want any help on this new problem (with time), post the exact equation and boundary values. I will see if I can help. I suspect an exponential decay is part of the solution. It is the most likely physical solution. The Neumann conditions can be applied by taking a derivative of some kind, it depends on the exact problem.

hope to hear from you,

Coomast

 

[phioder]

Hello coomast

The first disappointment on the trip was that separation of variables does not work anymore, the equation has three variables and I suppose this method does not function anymore.

The equation is following

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

and I think that the solution is not that strait ahead as in the other cases. On the way I found a similar problem (steady state) solved with a Hankel Transform and in terms of s.

We should stay with the original conditions so that the problem does not get very complicated

Initial Conditions:
$$u(2,z)=0, 0<z<4$$
$$u(r,0)=0, 0<r<2$$

Boundary Condition:
$$u(r,4)=u_0, 0<r<2$$

In the mean time I will be trying to found a solution too
Thank you for all this help :)
Phioder

 

[coomast]

Hello coomast

The first disappointment on the trip was that separation of variables does not work anymore, the equation has three variables and I suppose this method does not function anymore.

The equation is following

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

and I think that the solution is not that strait ahead as in the other cases. On the way I found a similar problem (steady state) solved with a Hankel Transform and in terms of s.

We should stay with the original conditions so that the problem does not get very complicated

Initial Conditions:
$$u(2,z)=0, 0<z<4$$
$$u(r,0)=0, 0<r<2$$

Boundary Condition:
$$u(r,4)=u_0, 0<r<2$$

In the mean time I will be trying to found a solution too
Thank you for all this help :)
Phioder

Hello phioder,

I will have a look at it the next days.

 

[phioder]

Hello coomast

If possible I would thank you if we come back to the steady state problem, unfortunately I still don't understand many things and I'm not ripe for transients

The equation would be following

$$0 = \frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2}$$

Initial/Boundary Conditions change to:
$$\frac{\partial u}{\partial z}=0, 0<z<4$$
$$\frac{\partial u}{\partial r}=0, 0<r<2$$
$$u(r,4)=u_0, 0<r<2$$

Be free to set the initial/boundary conditions, my goal is to understand how to solve the problem with conditions defined with fluxes

In the mean time will be trying to understand by my self
Thank you for all your help
Phioder

 

[coomast]

Hello coomast

The first disappointment on the trip was that separation of variables does not work anymore, the equation has three variables and I suppose this method does not function anymore.

The equation is following

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

and I think that the solution is not that strait ahead as in the other cases. On the way I found a similar problem (steady state) solved with a Hankel Transform and in terms of s.

We should stay with the original conditions so that the problem does not get very complicated

Initial Conditions:
$$u(2,z)=0, 0<z<4$$
$$u(r,0)=0, 0<r<2$$

Boundary Condition:
$$u(r,4)=u_0, 0<r<2$$

In the mean time I will be trying to found a solution too
Thank you for all this help :)
Phioder

It can be separated. I will explain it.

Hello coomast

If possible I would thank you if we come back to the steady state problem, unfortunately I still don't understand many things and I'm not ripe for transients

The equation would be following

$$0 = \frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2}$$

Initial/Boundary Conditions change to:
$$\frac{\partial u}{\partial z}=0, 0<z<4$$
$$\frac{\partial u}{\partial r}=0, 0<r<2$$
$$u(r,4)=u_0, 0<r<2$$

Be free to set the initial/boundary conditions, my goal is to understand how to solve the problem with conditions defined with fluxes

In the mean time will be trying to understand by my self
Thank you for all your help
Phioder

No problem, we'll go by it one step at the time. OK, First the separation of the variables. This is based on the following assumption. The solution of the differential equation can be written as a product of functions involving one variable at the time. This means that for the equation:

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

the following solution holds:

$$u(r,z,t)=R(r)\cdot Z(z)\cdot T(t)$$

Substituting this gives then (leaving out the explicit variable naming and calling $$\frac{k}{c_p \rho}=a^2$$):

$$RZT'=a^2\left( R''ZT+\frac{1}{r}R'ZT+RZ''T\right)$$

Now the easiest thing to do with this is to divide it by the solution. This means divide it by RZT, this gives (some like to call this the trick to solving it):

$$\frac{RZT'}{RZT}=a^2 \left(\frac{R''ZT}{RZT}+ \frac{1}{r}\frac{R'ZT}{RZT}+\frac{RZ''T}{RZT}\right)$$

Now we can rewrite it as:

$$\frac{1}{a^2}\frac{T'}{T} = \frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \frac{Z''}{Z}$$

This can now be transformed into two equations by stating that the left hand and right hand side both must be equal to a constant, thus we have:

$$\frac{1}{a^2}\frac{T'}{T} = \frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \frac{Z''}{Z}= -\lambda^2$$

From which:

$$T'+a^2\lambda^2 T=0$$
$$\frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \frac{Z''}{Z}=- \lambda^2$$

Now, the last equation can be rewritten as:

$$\frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \lambda^2= -\frac{Z''}{Z}$$

Again here the left hand and right hand side must be equal to a constant, however not the same as before, we get:

$$\frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \lambda^2= -\frac{Z''}{Z}= -\beta^2$$

From which now:

$$Z''-\beta^2Z=0$$
$$\frac{R''}{R}+ \frac{1}{r}\frac{R'}{R}+ \lambda^2+\beta^2=0$$

The latter one can be transformed into:

$$r^2R''+rR'+(\lambda^2+\beta^2)r^2R=0$$

Now you end up with three ordinary differential equations. See that you understand it completely and apply it to the equation you gave without the time involved. Try to apply it to others as well, p.e. in spherical coordinate systems.

The solution to these equations are given by:

$$T(t)=A\cdot e^{-a^2\lambda^2t}$$
$$Z(z)=Bcosh(\beta z)+Csinh(\beta z)$$
$$R(r)=DJ_0\left(\sqrt{\lambda^2+\beta^2}r\right)+ EY_0\left(\sqrt{\lambda^2+\beta^2}r\right)$$

Try to remember the solutions of some of these elementary equations. The next step is applying the boundary conditions. This is for the next time. Meanwhile, practice on this.

 

[coomast]

Hello phioder,

Let's take Fourier-Bessel one step further. First a few remarks. Solving partial differential equations using this method is not the only one available. There are others which can be used in good practice, p.e. Green's functions, conform transformations to name a few analytical ones apart from the numerical techniques that have evolved in such a manner that it can almost be stated that the classical physics exists in numerical form. So it is very likely that you will run into an equation which you can't solve using Fourier-Bessel series expansion. I have to admit that I'm not very familiar with these other analytic techniques. However the classical diffusion, wave and potential equations can be solved if the boundary conditions are "nice". If a solution is found then it can be shown that it is indeed the only one. There exists a theorem on this. Also the constants we used to separate the equation into smaller ones left out an intermediate step. It is normally written as a constant and this constant can have a positive, zero or negative value. One needs to test these in order to eliminate the trivial solutions to our problem. I did it directly because of the nature of the problem. I assume that the constants can be written like I did because of the physical possibilities of the problem.

Another remark on the solution I gave in the previous post is that the solution we found is "a" solution. This is called a principal solution. There exist many more of this type, all with different parameters. Due to the superposition theorem the sum of all these solutions is also a solution, as you might suspect. It is this sum that we use in applying other boundary conditions to arrive at the final solution to our problem.

So, the next step are the boundary conditions. These need to be physical or you end up with a solution that has nothing to do with reality. There are two types. First there are what I call the boundary conditions who put a value to a boundary of the problem. These need to completely describe the boundary or your problem is not well defined. The second type is what I call initial condition. This involves equations where time comes in the picture. Another condition very often used (this should always be done actually) is stating that the solution must be bounded. This means that "infinity" is not allowed as solution. Again, the physics don't allow this to happen.

So, the boundary conditions need to be complete or your problem is not well defined. In most cases one can almost immediately eliminate part of the principal solution in case of the boundedness condition. However, this is not always so. Mostly the conditions are given as 0. This makes life a lot easier and gives the opportunity of defining some of the parameters. In case these are not 0, one can use some tricks to transform the problem into one with conditions that are 0. Some information can be found in Boyce & DiPrima. Applying all of these conditions leaves us with the final principal solution. The sum of these is then the final complete solution which needs to be "tested" against the remaining conditions.

Let's apply this to the equation with the time involved. The conditions you gave are not good. Consider the first one, it was:

$$u(2,z)=0, 0<z<4$$

This is not well defined, time has been forgotten. It should be:

$$u(2,z,t)=0, 0<z<4$$

This is not to be a difficult person, it needs to be written because it states that this condition does not change in time, which is important to know. The second one was:

$$u(r,0)=0, 0<r<2$$

Same remark, it should be:

$$u(r,0,t)=0, 0<r<2$$

The final one was:

$$u(r,4)=u_0, 0<r<2$$

to become:

$$u(r,4,t)=u_0, 0<r<2$$

These are all boundary conditions, there is no initial condition. This is where the difficulty comes in the picture. It can't be solved because one does not know what the relation in time is. Several cases exist. P.e. We can consider the volume of the cylinder to be with a temperature distribution that is known at time 0. The boundaries determine what the distribution will be at a later time. Very often all boundaries are placed at a 0 value. The cylinder must "cool" down to zero for large time. Another possibility is that the cylinder was initially at 0 and will be "heated" by the walls to it's steady-state distribution. So as you can see there are several cases.

Take a look at the exact problem you want to solve and come back with the right conditions so we can continue with it. If you want to add fluxes at the boundary, you need to state that the derivative of the temperature is some value at the boundary. This comes from the law of Newton, which describes the heat flow at a boundary. If you have gotten the right conditions I will explain what to do with them. They need to be transformed as well, just like the initial partial differential equation.

Are you planning on buying the books I mentioned in a previous post?

Take your time to study this, it is a lot.

 

[phioder]

There are others which can be used in good practice, p.e. Green's functions... So it is very likely that you will run into an equation which you can't solve using Fourier-Bessel series expansion. I have to admit that I'm not very familiar with these other analytic techniques.

Hello coomast

This will be a quick post, thank you again for your answers and the clear solution with separation of variables. Trying to associate physical quantities to the steady state problem the Poisson Eq. crossed the way

http://www.engr.unl.edu/~glibrary/glibcontent/cyl_3d_mar2005b/node4.html

followed by those Green functions you are talking about. Now I will try to understand the meaning of the norms, eigenconditions and kernel functions.

About the books, I took a quick "look" to the online one, got the one of Schaum and the third one is on the way.

Will keep trying to found a way to get the solution, understand more the boundaries conditions to define them right and read carefully your posts.

About the boundaries the best way I can define them is following
$$\frac{\partial u }{\partial z }=0$$ at $$r = 2$$ for $$0<z<4$$
$$\frac{\partial u }{\partial r }=0$$ at $$r = 2$$ for $$0<z<4$$
$$\frac{\partial u }{\partial z }=0$$ at $$z = 0$$ for $$0<r<2$$
$$\frac{\partial u }{\partial r }=0$$ at $$z = 0$$ for $$0<r<2$$
meaning that I want at the side and bottom of the cylinder no change of temperature, outside of the cylinder the is a constant temperature of say 20 degrees celsius

I want to heat all the top with a constant distributed temperature
$$u(r,4,t)=u_0, 0<r<2$$

Thank you very much and best regards
phioder

PS: Actually I suppose the problem is to solve the steady heat conduction throughout a circular cylinder or a rod?

 

[coomast]

Hello phioder, just to keep you up to date. I'm working on some of the boundary issues, but it will take until the end of the week. It's very busy at work.

 

[phioder]

Hello phioder, just to keep you up to date. I'm working on some of the boundary issues, but it will take until the end of the week. It's very busy at work.

Hello coomast

Thank you for all your assistance and help
Does the solution to the problem changes a lot if the heat source is not uniform?

Best Regards
phioder

 

[coomast]

Hello phioder,

I finally found some time to work on the equation.

First some remarks on the boundary conditions you gave. The way to describe a zero heat flow at a boundary is by stating the following for the top and bottom of the cylinder:

$$\frac{\partial u}{\partial z}=0$$

For the bottom this is for z=0, and for the top z=4, for all r. At the circular bondary you can write:

$$\frac{\partial u}{\partial r}=0$$

For r=2, 0<z<4

The equation:

$$\frac{\partial u}{\partial z}=0$$

for r=2, 0<z<4 does not hold here because the heat flow is perpendicular to the surface area. This means that the first and last boundary condition you gave are not applicable. The second and third one describe the bottom and circular surface of the cylinder. The top is:

$$u(r,4,t)=u_0, 0<r<2$$

Now there is nothing here about the time. In case I assume that the temperature initially was zero in the cylinder you have enough conditions to define the problem. To summarise the entire question:

Solve the following equation in a cylinder with radius 2 and height 4:

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

where the cylinder initially was at u=0. The boundary conditons are:

$$\frac{\partial u}{\partial z}=0$$

for z=0,

$$\frac{\partial u}{\partial r}=0$$

for r=2, and the top z=4 is at a temperature

$$u=u_0$$

The solution must be bounded.

This is becoming a complicated problem. Because of the time dependency of the solution and the fact that the top boundary is not at a value 0 together with the fluxes (derivatives) set at 0 for the other boundaries it is not an easy one. Are you sure this is the case you are looking into? It might very well be that there is no solution at all using the Fourier-Bessel technique. I will give it a try if you state that this is indeed the one you need.

The solution will change in case the heat source is not uniform. Assuming that this is what you mean, that the u_0 is not a constant. The problem will become more difficult in case the heat source is phi dependent. However if the heat source is only depending on the radius, and a solution is possible for the constant heat source, it has one in this more general case as well.

If we were to solve this, the first thing to do is to split the problem up into two smaller ones. One which gives the steady-state and the other the transient part of the complete solution. I haven't done this yet on such a complicated model, but if you assure that the question is the one you need I will give it a try.

There is something similar in Spiegel, exercise 6.91, however the dependency of z is not present due to the assumption of an infinite long cylinder. This makes the problem a bit more solvable.

Hope this helps you one step further.

 

[phioder]

Hope this helps you one step further.

Hello coomast

Thank you very much for your answer. You make special comments that make me further understand small details I haven't understand yet.

After doing some research, as much as I understand now, the problem is indeed complicated and not trivial. During all this thread I have tried to split the problem and state basic questions because my knowledge and experience in PDE is very limited.

To reduce the complexity splitting the problem was based on:
1.- 2D steady heat conduction, no transients, no time dependency, no phi dependency, dependencies only on r and z
2.- For a first try the heat source considered as uniform
3.- For a first try not von Neumann boundary conditions, the heat needs to enter the cylinder from the upper part, all other sides are insulated

The equation is following:
$$\frac{-g(r,z)}{k} = \frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2}$$

From this link http://www.engr.unl.edu/~glibrary/glibcontent/cyl_3d_mar2005b/node4.html what I understand is that the function g(r,z) is the heat source. They give a solution with green functions, unfortunately I don't understand them to implement a program.

If the heat souce is constant the equation is following:
$$\frac{u_0}{k} = \frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2}$$

Unfortunately I began to loose focus and my post can't be well formulated, my apologies

Best Regards
phioder

PS: What does it mean that a solution is "bounded" or "not bounded"

 

[coomast]

OK, phioder, now I understand the heat source. It is nothing more than the extension
of the equation we used all along. We had the following:

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)$$

Which can be generalized for the case with a heat source to become:

$$\frac{\partial u}{\partial t} = \frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)+ \frac{E}{c_p \rho}$$

In which E the heat source is. The units of E are W/m^3. And for the steady-state this is then:

$$0=\frac{k}{c_p \rho} \cdot \left(\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2} \right)+ \frac{E}{c_p \rho}$$

Or:

$$\frac{-E}{k}=\frac{\partial^2 u}{\partial r^2}+\frac{1}{r} \cdot \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial z^2}$$

This is indeed not trivial, I am not sure whether this can be solved by the method of Fourier-Bessel series expansion. I will have a look in my books and see what I can find. Nevertheless the separation of the variables and the use of the series together with the study of the very interesting Bessel functions has not been a waste of time.

[Edit: name changed]

A solution that is bounded is simply stated one that is finite. It is one that does not diverge to infinity for some point(s). Consider a graphic of the Bessel function of the second type of order 0, you will see that for x going towards 0, the function goes to -infinity. This is not realistic in our problem. If we state that the solution must be bounded, we simply eliminate these solutions from the set of possible solutions. This does not mean that it is a bad solution, it is a perfect valid one for the math, but not for the physics, it is not allowed here. Otherwise a solution that is unbounded is then one that can diverge to infinity.

 

[phioder]

This is indeed not trivial, I am not sure whether this can be solved by the method of Fourier-Bessel series expansion. I will have a look in my books and see what I can find. Nevertheless the separation of the variables and the use of the series together with the study of the very interesting Bessel functions has not been a waste of time.

Hello coomast

Absolutely not a waste of time, your answers give me a great feedback to understand better the equation and also included a lot of motivation. IMPO there are many small things that are not explained in books that save paper and ink and on the other side engineers are lost.

This forum is a great place to share knowledge, ideas and questions
phioder

 

[coomast]

This is indeed not trivial, I am not sure whether this can be solved by the method of Fourier-Bessel series expansion. I will have a look in my books and see what I can find.

phioder, I looked trough some of my books in order to find some more information on the solution of the problem you stated. The partial differential equation for conductive heat transport with a heat source seems to be not used in books. At least not with what I could find. The only example I found was one given in my course on heat transfer. It was one describing the heat transfer through an insulated electrical wire. The current causes the wire to heat up. It was however for an infinite long wire, thus without phi and z dependency, giving rise to an ordinary differential equation which had nothing to do anymore with Bessel.

I assume that is not easy to do with the use of the separation of the variables and then to use the Fourier-Bessel expansion. Otherwise I definitely would have bumped into an exercise or an example somewhere. I think the other methods are more applicable. But here I can't help you. It has been too long since I read about it. I never actually had to use them. I would suggest to start a new post with the exact description of the problem you would like to see solved. There are certainly people who are willing to help you with the Green's functions.

Thank you very much for your kind words along the way and I'm very happy that I was of some assistance in getting you to appreciate (amongst a large set with many more interesting) Bessel functions.