Poisson equation in axisymmetric cylindrical coordinates

Question

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

Answer 1

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.

Comments

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

Answer 2

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!

Comments

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