Elementary question concerning flux integral
Solution 1:
The volume integral of divergence can be evaluated as
$$\int_0^1\int_0^{2\pi}\int_r^1 (2+4z)r dzd\theta dr = \frac{4 \pi }{3 }$$
Also surface can parametrized as $S(r, \theta) = (r \cos (\theta ),r \sin (\theta ),r)$ and integral can be evaluated as
$$\int_0^1 \int_0^{2\pi }\langle r \cos \theta, r\sin\theta,r^{4}\rangle \cdot (-1)\left( \frac{\partial S}{\partial r} \times\frac{\partial S}{\partial \theta} \right) d\theta dr = \frac \pi 3 $$ $(-)$ because $\displaystyle \frac{\partial S}{\partial r} $ is radially tangent and $\displaystyle \frac{\partial S}{\partial \theta}$ is circularly tangent anticlockwise giving the flux through inner surface.
Solution 2:
Your error lies in your evaluation of the volume integral.
Assuming that you have made a typo in the order of the integration variables, so that it was meant to read
$$ \int^{2\pi}_{0}\int^{1}_{0}\int^{r}_{0} \left(2+4z^{3} \right)dz\,r\,dr\,d\theta $$
then you are describing the area between the cone and the $x$-$y$ plane, rather than between the cone and the plane $z=1$. Instead, you want the first integral to run from $r$ to $1$. If you had the order of integration variables right, then your limits were even more incorrect.
With this correction, you get \begin{align} \int^{2\pi}_{0}\int^{1}_{0}\int^{1}_{r} \left(2+4z^{3} \right)dz\,r\,dr\,d\theta &= \int_0^{2\pi}\int_0^1\left(3r-2r^2-r^5\right)dr\,d\theta\\ &= 2\pi \left[\frac{3r^2}2-\frac{2r^3}3-\frac{r^6}6\right]_0^1 = \frac{4\pi}3 \end{align}
Note that your "correct" answer is also incorrect. Simple geometric considerations (vector, for $z<1$, will point outwards compared with surface of cone) indicate that the flux must be positive. I haven't checked, but I assume that orientation wasn't accounted for.