The vector field in cylindrical coordinates is given by $$ \vec{F} = F_{\rho} \unitv{\rho} + F_{\phi} \unitv{\phi} + F_{z} \unitv{z} $$ We place our tiny cube-like volume in the field \(\vec{F}\) at some point \(\vec{P}=(\rho,\phi,z)\).
First we consider the flux passing through the surfaces normal to the \(z\) direction (the top and bottom faces of the cube). Since the distances involved are small, we can approximate the field passing through these as $$ F_{top} = F_z(\rho,\phi,z + \Delta z) \approx F_z(\rho,\phi,z) + \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ It should be noted that the other components of the field don't play a part here because they are parallel to these surfaces so nothing going in the, say, \(\rho\) direction leaves via the \(z\) direction. Similarly for the bottom surface, $$ F_{bottom} = F_z(\rho,\phi,z - \Delta z) \approx F_z(\rho,\phi,z) - \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ The areas of these surfaces are equal and independent of \(z\). $$ A_z = \rho\Delta\phi\Delta\rho $$ Therefore, the flux through these z-surfaces is (remembering that the flux through the bottom surface is into the volume and so is negative) $$ \Psi_z = F_{top}A_z - F_{bottom}A_z = (F_{top} - F_{bottom})A_z $$ $$ \Psi_z = (F_z + \pard{z}{F_z}\frac{\Delta z}{2} - (F_z - \pard{F_z}{z}\frac{\Delta z}{2}))\rho\Delta\phi\Delta\rho $$ $$ \Psi_z = \pard{F_z}{z}\rho\Delta \phi \Delta\rho \Delta z $$ Recognize that the term \( \rho\Delta \phi \Delta\rho \Delta z \) is the differential volume in cylindrical coordinates, so we can write the previous equation as $$ \begin{equation} \label{eq:fluxz} \Psi_z = \pard{F_z}{z}\Delta V \end{equation} $$ Where \(\Delta V\) is the differential volume element. I should note that if you were to do this derivation in Cartesian coordinates, the x and y component derivations would be identical (except for the subscripts). But in cylindrical coordinates, things get a little more interesting from here on out. Looking now at the \(\phi\) direction we consider the front and back faces of the volume. (I will omit the function arguments to the partial derivatives in the service of brevity with the understanding that we are talking about the point \(P=(\rho,\phi,z)\)). $$ F_{back} = F_\phi(\phi + \Delta \phi) \approx F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ F_{front} = F_\phi(\phi - \Delta \phi) \approx F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ A_\phi = \Delta\rho\Delta z $$ Because the positive \(\phi\) direction is toward the back of the volume, $$ \Psi_{\phi} = (F_{back} - F_{front})A_{\phi} $$ $$ \Psi_{\phi} = (F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} - (F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2}))\Delta\rho\Delta z $$ $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z $$ The term \(\Delta\phi\Delta\rho\Delta z\) is almost the differential volume. And maybe because we know where we are headed, we bend the equation to our will $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z = \pard{F_\phi}{\phi}\frac{\rho}{\rho}\Delta\phi\Delta\rho\Delta z $$ $$ \begin{equation} \label{eq:fluxphi} \Psi_\phi = \frac{1}{\rho}\pard{F_\phi}{\phi}\Delta V \end{equation} $$ Now we try to handle the \(\rho\) component. We need to be careful here. But first things first. The vector field is defined similarly to previous components. $$ F_{right} = F_\rho(\rho + \Delta \rho) \approx F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ $$ F_{left} = F_\phi(\rho - \Delta \rho) \approx F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ Here is where it gets tricky. The area definitions are different because the two surfaces are not equal in size and so must be treated differently. $$ A_{right} = (\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ $$ A_{left} = (\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ First computing the flux on the right side. $$ \Psi_{right} = F_{right}A_{right}=(F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{right} = (F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z + \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{right} = F_\rho\rho\Delta\phi\Delta z+\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Similarly on the left hand side $$ \Psi_{left} = F_{left}A_{left}=(F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{left} = (F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z - \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{left} = F_\rho\rho\Delta\phi\Delta z-\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z - \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Computing \(\Psi_\rho\) (finally, after all that algebra) by subtracting and cancelling terms $$ \Psi_\rho = \Psi_{right} - \Psi_{left} $$ $$ \Psi_\rho = F_\rho\Delta\rho\Delta\phi\Delta z + \pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z $$ We make the same modification to the first term that we did with \(\Psi_\phi\) and substituting \(\Delta V\) $$ \begin{equation} \label{eq:fluxrho} \Psi_\rho = \frac{1}{\rho}F_\rho\Delta V + \pard{F_\rho}{\rho}\Delta V \end{equation} $$ This last equation is correct but we normally don't see it written this way, so we factor out a \(\frac{1}{\rho}\) and invoke the reverse chain rule. $$ \Psi_\rho = \frac{1}{\rho}(F_\rho + \rho\pard{F_\rho}{\rho})\Delta V = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho)\Delta V $$ And so the total flux through our tiny volume is $$ \Psi_{total} = \Psi_\rho + \Psi_\phi + \Psi_z = \left(\frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z}\right)\Delta V $$ The divergence is the ratio of the flux through a closed surface as the volume of that closed surface approaces 0. $$ \nabla\cdot\vec{F} = \lim\limits_{\Delta V \to 0} \frac{\Psi_{total}}{\Delta V} = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z} $$
The derivation of divergence in spherical coordinates procedes along similar lines as the \(\phi\) and \(\rho\) terms. It's a little more complicated because of the added \(\sin\theta\) terms.