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!)