I'm currently trying to code a FEM-program (finite element method) to solve a heat conduction problem including a few tricky things like phase change and a chemical reaction.
In a FEM-book [REDDY - An introduction to nonlinear finite element analysis] I've found a formulation for a problem in cylindrical space (dependent on $r,z$) that is applicable to the problem I'm trying to solve. However I can't fully follow the derivation that has a minor difference to my derivation. I'll try to document my derivation as comprehensive as possible. I don't know whether such detail is necessary - you might just have to take a look at the very last couple of passages I wrote to give me an answer.
I first do the formulation of eq. 1 in cylindrical coordinates:
The governing equation is $$-\operatorname{div}(a\cdot \operatorname{grad}u)=-\nabla\cdot(a\cdot\nabla u)=f \qquad (\text{eq. }1)$$
with: $$a=\begin{bmatrix}a_{rr}&0&0\\0&a_{\phi\phi}&0\\0&0&a_{zz}\end{bmatrix}.$$
To my knowledge in cylindrical coordinates $(r,\phi,z)$ the following statements are valid:
\begin{align*}\operatorname{div}u&=\cfrac{1}{r}\cfrac{\partial}{\partial r}(r\cdot u)+\cfrac{1}{r}\cfrac{\partial u}{\partial \phi}+ \cfrac{\partial u}{\partial z}\\ \operatorname{grad} u&=\cfrac{\partial u}{\partial r}+\cfrac{1}{r}\cfrac{\partial u}{\partial \phi}+\cfrac{\partial u}{\partial z}. \end{align*}
With those statements eq. 1 becomes
$$ -\cfrac{1}{r} \cfrac{\partial}{\partial r}\left(a_{rr}r\cfrac{\partial u}{\partial r}\right) -\cfrac{1}{r^2}\cfrac{\partial}{\partial \phi}\left(a_{\phi \phi}\cfrac{\partial u}{\partial \phi}\right) -\cfrac{\partial}{\partial z}\left(a_{33}\cfrac{\partial u}{\partial z} \right) -f=0\qquad (\text{eq. }1a) $$
Considering that by the definition of the problem $u$ is independent of $\phi$, which means $\cfrac{\partial u}{\partial \phi}=0,$ we get:
$$ -\cfrac{1}{r} \cfrac{\partial}{\partial r}\left(a_{rr}r\cfrac{\partial u}{\partial r}\right) -\cfrac{\partial}{\partial z}\left(a_{zz}\cfrac{\partial u}{\partial z} \right) -f(r,z)=0\qquad (\text{eq. }1b) $$
So far, so good, this is the same equation the book gives me.
The next step in FEM is to discretize the region and develop a weak form of the equation on an element level.
To get the "weak form" of eq. 1b $u$ is approximated (discretization of the domain):
$$u \cong \sum_{j=0}^n u_j\varphi_j.$$
$u_j$ symbolizes the value of $u$ at node $j,$ $\varphi_j$ is the approximation function for node $j$ ($\varphi_j=1$ at node $j,$ and $\varphi_j=0$ at every other node). With this approximation, eq. 1b cannot be exactly zero any more. I'll call eq. 1c "Error" $E(r,z)$ for the moment. In FEM this error gets multiplied with a weight function $w$ and integrated over the domain $\Omega$. This weighted residue statement is equal to zero.
$$\int_\Omega (E(r,z)\cdot w)\, dr\, dz\, d\phi=0, \quad\text{or}$$
$$ \int_\Omega \left(-\cfrac{1}{r} \cfrac{\partial}{\partial r}\left(a_{rr}r\cfrac{\partial u}{\partial r}\right) -\cfrac{\partial}{\partial z}\left(a_{zz}\cfrac{\partial u}{\partial z} \right) -f(r,z)\right)\cdot w\space dr\space d\phi \space dz =0. $$
If we substitute:
$\biggl(a_{rr}r\cfrac{\partial u}{\partial r}\biggr)=F_1$ and $\biggl(a_{zz}\cfrac{\partial u}{\partial z}\biggr)=F_2$ we get:
$$ \int_\Omega \biggl(-\cfrac{1}{r} \cfrac{\partial}{\partial r}F_1 -\cfrac{\partial}{\partial z}F_2 -f(r,z)\biggr)\cdot w\space dr\space d\phi \space dz$$
$$=\int_\Omega \biggl(-\cfrac{1}{r} \cfrac{\partial F_1}{\partial r}w -\cfrac{\partial F_2}{\partial z}w -f(r,z)w\biggr)\space dr\space d\phi \space dz $$
With:
$$\cfrac{\partial}{\partial r}(wF)=\cfrac{\partial w}{\partial r}F+\cfrac{\partial F}{\partial r}w \space \rightarrow \space \cfrac{\partial F}{\partial r}w=\cfrac{\partial}{\partial r}(wF)-\cfrac{\partial w}{\partial r}F$$
we can replace:
$$\int_\Omega \left(-\cfrac{1}{r} \left(\cfrac{\partial}{\partial r}(wF_1)-\cfrac{\partial w}{\partial r}F_1\right) -\left(\cfrac{\partial}{\partial z}(wF_2)-\cfrac{\partial w}{\partial z}F_2\right) -f(r,z)w\right)\, dr\, d\phi \, dz. $$
Further, due to the divergence theorem the following expression can be rewritten:
$$\int_\Omega\cfrac{\partial}{\partial r}(wF)dV=\oint_\Gamma(wF)n_r\,ds$$
and we get:
$$\int_\Omega \left(\cfrac{1}{r} \cfrac{\partial w}{\partial r}F_1+\cfrac{\partial w}{\partial z}F_2-f(r,z)w\right)\, dr\, d\phi \, dz-\oint_\Gamma\left(\cfrac{1}{r}wF_1\right)n_r\,ds-\oint_\Gamma(wF_2)n_z\,ds,$$
with $F_1$ and $F_2$ re-substituted the LHS becomes:
$$\int_\Omega \left(\cfrac{1}{r} \cfrac{\partial w}{\partial r}a_{rr}r\cfrac{\partial u}{\partial r}+\cfrac{\partial w}{\partial z}a_{zz}\cfrac{\partial u}{\partial z}\right)\, dr\, d\phi \, dz =2\pi\int_\Omega \left(a_{rr}\cfrac{\partial w}{\partial r}\cfrac{\partial u}{\partial r}+a_{zz}\cfrac{\partial w}{\partial z}\cfrac{\partial u}{\partial z}\right)\, dr\, dz.$$
In the book however it is (with no derivation at all where does the $r$ come from??):
$$2\pi\int_\Omega \left(a_{rr}\cfrac{\partial w}{\partial r}\cfrac{\partial u}{\partial r}+a_{zz}\cfrac{\partial w}{\partial r}\cfrac{\partial u}{\partial z}\right)\cdot r\, dr\, dz.$$
For the RHS I get:
$$2\pi\int_\Omega wf(r,z)\, dr \, dz +\oint_\Gamma\left[\left(\cfrac{1}{r}wa_{rr}r\cfrac{\partial u}{\partial r}\right)n_r+\left(wa_{zz}\cfrac{\partial u}{\partial z}\right)n_z\right]\,ds $$
$$=2\pi\int_\Omega wf(r,z)\, dr \, dz +\oint_\Gamma\left[\left(wa_{rr}\cfrac{\partial u}{\partial r}\right)n_r+\left(wa_{zz}\cfrac{\partial u}{\partial z}\right)n_z\right]\,ds. $$
Again the formulation in the book is very similar and but has an additional $r$ in each integral - which I can't grasp where it's from:
$$2\pi\int_\Omega wf(r,z)\cdot r\, dr \, dz +\oint_\Gamma\left[\left(wa_{rr}\cfrac{\partial u}{\partial r}\right)n_r+\left(wa_{zz}\cfrac{\partial u}{\partial z}\right)n_z\right]\cdot r \, ds. $$
I'm not quite sure whether it is a valid mathematical operation to put an $r$ in every integral - is it? In this case my solution would be equivalent. Is there any particular reason to where this $r$ comes from?
Okay, I found the answer to my question. During the transformation to cylindrical coordinates or rather when you integrate them over an area / volume, you also have to transform $dA \, or \, dV$. I took a screenshot from a (German) textbook, the graphic should be self-explainatory:
So $dA=r \, d\phi \, dz \, \rightarrow \, dV= r \, dr \, d\phi \, dz$, which explains the additional $r$ in every integral statement that confused me.