Processing math: 100%
This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Poisson equation in axisymmetric cylindrical coordinates

+1 vote

I am trying to derive the equation for the heat equation in cylindrical coordinates for an axisymmetric problem. The pde is the following:

\frac{1}{r} \frac{\partial}{\partial r}\left( K r \frac{\partial T}{\partial r} \right) + \frac{\partial}{\partial z} \left( K \frac{\partial T}{\partial z} \right) = 0

By checking online and some trial and error I found the following variational form:

T = TrialFunction(V)
q = TestFunction(V)

r = Expression('x[0]', degree=1)
a = k*(Dx(T,0)*Dx(q,0)*r + Dx(T,1)*Dx(q,1))*dx()

I also assumed k is uniform and can be factored out. I understand where the second term comes from: its from integration by parts or Green's first identity in 1D.

How does one get the first term? Is it correct? I checked it against the analytic solution for heat transfer between two concentric cylinders and I get the same temperature distribution.

If I don't use r and just use x[0] dolfin says that x is not defined. What should I be using?

Thank you

asked Jul 1, 2017 by alexmm FEniCS User (4,240 points)
edited Jul 6, 2017 by alexmm

2 Answers

+1 vote
 
Best answer

The second term needs x[0] factor.

The easiest way to get it is probably transforming the bilinear form in rectangular coordinate to the one with cylindrical coordinate using the chain rule. In vector calculus,

dx dy = r*dr d(theta)
dr/dx = cos(theta)
dr/dy = sin(theta)

If T and q are independent of the theta variable, it is fairly simple because

dT/dx = (dT/dr)(dr/dx) = cos(theta)(dT/dr)
dT/dy = (dT/dr)(dr/dy) = sin(theta)(dT/dr)

and therefore

(dT/dx)(dq/dx) + (dT/dy)(dq/dy) = (dT/dr)*(dq/dr)

Finally, the integration from 0 to 2pi in theta must be averaged by dividing by 2pi because the equation describes the temperature only on one section, not the whole cylinder.

I don't know a way to use x[0] directly in variational forms.

answered Jul 6, 2017 by JeonghunLee FEniCS User (1,050 points)
selected Jul 7, 2017 by alexmm

Thank you for the answer. If I try doing the same for (dT/dz)(dq/dz) I get the same in cylindrical coordinates. However there is an r at the end of the first term in the bilinear form and not in the second term.

That's from the change dx dy --> r dr d theta.

But what about the second term? Shouldn't it have an r at the end as well?

That's right. I fixed it.

What do you mean? If you use the chain rule like you said you get for the all the terms:

r = Expression('x[0]', degree=1)
a = k*(Dx(T,0)*Dx(q,0) + Dx(T,1)*Dx(q,1))*dx()

But I have:

r = Expression('x[0]', degree=1)
a = k*(Dx(T,0)*Dx(q,0)*r + Dx(T,1)*Dx(q,1))*dx()

How does the second term not end up with an r?

I mean that both terms should have r factor as
a = k(Dx(T,0)Dx(q,0) + Dx(T,1)Dx(q,1))r*dx()
This is same as the formula in the other answer.

Ok got it thank you

0 votes

First, there is a typo in your equation, one should reed:
\frac{1}{r}\frac{\partial}{\partial r}\left( K r \frac{\partial T}{\partial r} \right) + \frac{\partial}{\partial z} \left( K \frac{\partial T}{\partial z} \right) = 0
If K is assumed to be constant, you can take K=1.
To obtain a weak formulation, you need to choose an inner product. In cylindrical coordinates with axissymmetry, the natural one is:
(u,v)=\int_{\mathbb{R}^+\times\mathbb{R}}uv r dr dz
Now if you take the inner product between the equation and a test function q, after integrating by parts, you get:
(\partial_r T,\partial_r q) + (\partial_z T,\partial_z q)=0
i.e. in FEniCS

(Dx(T,0)*Dx(q,0) + Dx(T,1)*Dx(q,1))*x[0]*dx()

You can also get the answer directly, since in spherical coordinates the gradient is given by
grad(u) = \partial_r u e_r + \partial_z u e_z
Finally, be aware of the boundary condition at r=0, you probably want a solution smooth that is smooth at the origin!

answered Jul 5, 2017 by neliju FEniCS Novice (210 points)

Thanks for the answer. I fixed the typo. The answer you gave has an extra x[0] for the second term compared to the bilinear form I have.

I'm pretty sure my weak formulation is correct. You also obtain this result if you write the weak formulation in Cartesian coordinates, i.e. (grad(u),grad(w)) and then write the gradient in cylindrical coordinates and the integration dx dy dz=r dr dz.

...