We can easily find the energy equation for incompressible fluid as a temperature equation:
$$ \rho c \frac{DT}{Dt}=\nabla\cdot(k\ \nabla T) + \tau_{xx}\frac{\partial u}{\partial x}+\tau_{yx}\frac{\partial u}{\partial y} +\tau_{zx}\frac{\partial u}{\partial z}+\tau_{xy}\frac{\partial v}{\partial x}+\tau_{yy}\frac{\partial v}{\partial y}+\tau_{zy}\frac{\partial v}{\partial z}+\tau_{xz}\frac{\partial w}{\partial x}+\tau_{yz}\frac{\partial w}{\partial y}+\tau_{zz}\frac{\partial w}{\partial z}+S_i (1) $$
where $\frac{D()}{Dt}$ is the material derivative, $\nabla\cdot \vec e$ is the divergent operator applied on vector $\vec e$, $\nabla \phi$ is the gradient operator applied on scalar function $\phi$ and $S_i$ are source terms.
Using the Newtonian model for the viscous stresses, we can write it as
$$ \rho c \frac{DT}{Dt}=\nabla\cdot(k\nabla T)+\Phi+S_i (2) $$
and $\Phi$ is the dissipation function and is the multiplication of the viscosity $\mu$ and the rates of deformation.
We can also find another way of writing this equation as: $$ \frac{\partial(\rho T)}{\partial t} + \nabla\cdot(\rho T \vec u)=\nabla\cdot(k\nabla T)+\Phi+S_i (3) $$ which is more suitable for a finite volume analysis.
This equation is similar to the conservative form of the fluid flow equation for a general variable $\phi$, which can be usefully written in the following form:
$$ \frac{\partial(\rho \phi)}{\partial t} + \nabla\cdot(\rho \phi \vec u)=\nabla\cdot(\Gamma\ \nabla T)+S_{\phi} (4) $$
which is a general equation, which can be formally integrated over a volume control and we use the Gauss-Divergence theorem to obtain the discretized equations in a finite volume analysis.
My questions are about these same equations, but written in cylindrical coordinates. The material derivative can be written as
$$ \frac{D\phi}{Dt}=\frac{\partial\phi}{\partial t}+\vec u \nabla(\phi) (5) $$
So, if we use (5) on (2) and the \nabla operator in cylindrical coordinates, we get to the energy equation in cylindrical coordinates (in steady state form, disregarding the time variations):
$$ \rho c \left( u \frac{\partial T}{\partial x} + \frac{v_{\theta}}{r}\frac{\partial T}{\partial \theta}+v_r\frac{\partial T}{\partial r} \right) = k \frac{\partial^2 T}{\partial x^2}+\frac{k}{r}\frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r} \right) + \frac{k}{r^2}\frac{\partial^2 T}{\partial \theta^2} + \mu \Phi (6) $$
which is absolutely right and uses the gradient operator in cylindrical coordinates.
However, if I want to to a finite volume analysis, equation (4) is more suitable, because we want to use the Divergence Theorem to change the volume integral of the advective term in a surface (flux) integral. However, if we use the $\nabla$ operator in cylindrical coordinates, the term relative to the $r$ coordinate in the advective (LHS of the equation) should be
$$ \frac{\rho v_r}{r} \frac{\partial(r T)}{\partial r} (7) $$
instead of the $v_r\frac{\partial T}{\partial r}$ term in eq 6.
What am I doing wrong or using wrong? Is the conservative general form for a general variable $\phi$ only suitable for the cartesian coordinates? How do I treat the advective term in finite volume analysis in cylindrical coordinates?