Laplace equation in cylindrical coordinates

[alexrao]

Can anyone help with the solution of the Laplace equation in cylindrical coordinates
\frac{\partial^{2} p}{\partial r^{2}} + \frac{1}{r} \frac{\partial p}{\partial r} + \frac{\partial^{2} p}{\partial z^{2}} = 0

with Neumann no-flux boundaries:
\frac{\partial p}{\partial r} \left(0,z\right) = 0
\frac{\partial p}{\partial r} \left(Rmax,z\right) = 0
\frac{\partial p}{\partial z} \left(r,Zmax\right) = 0

and a Dirichlet upper boundary:
p(r,0) = f(r)

 

[arildno]

With a trial solution using separation of variables, p(r,z)=P(r)*Q(z), we get:
\frac{P''(r)+\frac{1}{r}P'(r)}{P(r)}=\frac{1}{Q(z)}Q''(z)

Thus, we get two ordinary diff.eqs,
P''(r)+\frac{1}{r}P'(r)=CP(r) (*)
Q''(z)=CQ(z)(**)

You might try to work with these two, (*) clearly being the most problematic one.

 

[LCKurtz]

With a trial solution using separation of variables, p(r,z)=P(r)*Q(z), we get:
\frac{P''(r)+\frac{1}{r}P'(r)}{P(r)}=\frac{1}{Q(z)}Q''(z)

Thus, we get two ordinary diff.eqs,
P''(r)+\frac{1}{r}P'(r)=CP(r) (*)
Q''(z)=CQ(z)(**)

You might try to work with these two, (*) clearly being the most problematic one.

But maybe not too problematic. Multiply it through by r² and you have a form of Bessel's equation.

 

[alexrao]

Thanks for getting me started, and pardon the ignorance of this geochemist who hasn't taken a pde class. I'm trying to follow you, and can't understand how you got that 2nd ode. Shouldn't (**) be

Q''(z) = -C Q(z) ?

 

[arildno]

You are BOTH right. Mea culpa.

 

[alexrao]

Moving along here...

So the solution to (**):

Q(z) = Acosh(\lambda z)+Bsinh(\lambda z)

As Arildno predicted, the solution to (*) is more problematic, since this is my first experience with Bessel functions. Here's where I am so far...

P(r) = CJ_{0}(\lambda r)+DY_{0}(\lambda r)

And here I set the second term must go to zero for the solution to be bounded.

So the final solution is of the form:
p(r,z) = [Acosh(\lambda z)+Bsinh(\lambda z)]J_{0}(\lambda r)

Using the 2nd Neumann BC:
\frac{\partial p}{\partial r}(r=Rmax=9.5) = 0
This means that
J'_{0}(9.5\lambda) = \frac{\partial J_{0}}{\partial r}= 0
But I'm not sure how this helps me right now, so I'll keep it in mind...

Using the 3rd Neumann BC:
\frac{\partial p}{\partial z}(z=Zmax=20) = 0
This means that
[Acosh(20\lambda)+Bsinh(20\lambda)]=0
Ah, this is more helpful. I learned that A=-B, and I think at the bottom boundary, this whole term drops out.

Using the 1st Neumann BC: (I guess I could have done these in order...)
\frac{\partial p}{\partial r}(r=0) = 0
This means that
J'_{0}(0\lambda) = \frac{\partial J_{0}}{\partial r}= 0
Again I'm not really sure how this helps me right now, so I'll keep it in mind...

So now we have:
p(r,z)=A[cosh(\lambda z)-sinh(\lambda z)]J_{0}(\lambda r)

with
J'_{0}(0\lambda, r=0) = J'_{0}(9.5\lambda, r=9.5) = 0

Now here's where it gets tricky and I get stuck...

Using the fourth Dirichlet upper (z=0) boundary condition, I get:

p(r,0)=f(r)=\sum_{\infty}^{n=1}A_{n}[cosh(0)-sinh(0)]J_{0}(\lambda _{n}r)=\sum_{\infty}^{n=1}A_{n}J_{0}(\lambda _{n}r)

Unfortunately, f(r) is not a constant, but rather an ugly empirical function of r.

So I have written:

f(r)=\sum_{\infty}^{n=1}A_{n}J_{0}(\lambda _{n}r)

And I have no idea where to go from here. I gather this will eventually end up in Matlab or R, but first I need a better understanding of the Bessel function J0 and A. I've read that the Bessel function of the first kind of order n (in my case n=0, right?) can be expressed as:

J_{n}(x) = \sum_{\infty}^{k=1}\frac {(-1)^{k}(x/2)^{n+2k}}{k!\Gamma (n+k+1)}}

Again, any help would be greatly appreciated.

 

[arildno]

Now, assuming that the class of J_0's, with scaled arguments, represents a COMPLETE BASIS for functions on R (i.e, that any function is representable as a linear combination of these J_0's), then the A_n's are simply the required coefficients.

 

[LCKurtz]

I haven't checked all your details, but your P equation and boundary conditions appear to be a Sturm-Liouville system which answers questions about orthogonality of the eigenfunctions, gives formulas for the eigenfunction expansion coefficients and settles convergence. For example, see:

http://www.efunda.com/math/ode/Sturm_liouville.cfm

 

[alexrao]

Thanks again. Arildno, you mentioned that the A_{n}'s are the required coefficients. But if I want to define a pressure distribution in this cylinder, don't I also need to figure out \lambda? And is it true that \lambda in P(r) and \lambda in Q(z) are not the same?

 

[arildno]

No, the \lambda_{n} MUST be the same numbers; otherwise, your diff.eq won't be satisfied. (Remember, it is directly related to your C's!)