I need to calculate the flux of the vector field $\vec{F}$ through the surface $D$, where $$\vec{F} = \left<z, \, y \sqrt{x^2 + z^2}, \, -x \right> \\ D = \{x^2+6x+z^2\le 0 \,| -1\le y \le 0\}.$$
So there should be a cylinder (height on $y$ axis) shifted by 3 units (center $x=-3$), with cut on $ -1\le y \le 0$.
Cylindrical parametrization:
$$\begin{cases} x=r\cos\theta -3\\ y=y\\ z=r\sin\theta \end{cases} $$
I think I can do this problem in two ways: The first one by calculating the flux for each of the 3 surfaces (1 cylinder, 2 disks), and the second one by using the divergence theorem.
$$\operatorname{div}F = \sqrt{x^2+y^2} \overset{\text{cylindrical}}{\underset{\text{coordinates}}{=}} \sqrt{r^2-9-6r\cos\theta},$$
so with the divergence theorem, $$\int_{0}^{2\pi}\int_{0}^{3\sqrt{3}}\int_{-1}^{0}{r\sqrt{r^2-9-6r\cos\theta} \, dy \, dr \, d\theta}.$$
Is this correct?
I'm not exactly sure where the $3\sqrt{3}$ comes from in your result, but there is indeed more than one way to evaluate this problem.
(1) Direct method
Here is some technical information about this method from MIT's open notes, and some visualization for what the flux of a vector field through a surface is. Let the flux of a vector field $\vec{\mathbf{V}}$ through a surface $\Sigma$ be denoted $\Phi$ and defined $$ \Phi := \iint_\Sigma \mathbf{\vec{V}} \cdot \mathbf{\hat{n}} \, d\sigma. $$ The vector $\mathbf{\hat{n}}$ is the unit outward normal to the surface $\Sigma$. Suppose $\Sigma$ is given by $z = f(x,y).$ Let $\mathbf{\vec{r}}(x,y)$ trace $\Sigma$ such that $$ \mathbf{\vec{r}}(x,y) = \begin{pmatrix} x \\ y \\ f(x,y) \end{pmatrix}. $$ Then the unit normal $\mathbf{\vec{n}}$ is given by $$ \mathbf{\vec{n}} = \frac{\mathbf{\vec{r}}_x \times \mathbf{\vec{r}}_y}{|| \mathbf{\vec{r}}_x \times \mathbf{\vec{r}}_y ||} = \frac{1}{\sqrt{f_x^{\,2} + f_y^{\,2} + 1}} \begin{pmatrix} -f_x \\ -f_y \\ 1 \end{pmatrix}. $$ So given that $\mathbf{\vec{V}} = u(x,y,z) \mathbf{\hat{i}} + v(x,y,z) \mathbf{\hat{j}} + w(x,y,z) \mathbf{\hat{k}}$, the corresponding flux of $\mathbf{\vec{V}}$ through $\Sigma$ is $$ \Phi = \iint_\Sigma \frac{-uf_x - vf_y + w}{\sqrt{f_x^{\,2} + f_y^{\,2} + 1}} \, d\sigma. $$
For the given field, we have $$\mathbf{\vec{V}} = z \mathbf{\hat{i}} + y \sqrt{x^2 + z^2} \mathbf{\hat{j}} - x \mathbf{\hat{k}}, $$ and the surface $\Sigma$ is given such that $(x + 3)^2 + z^2 = 9\ \forall y \in (-1,0).$ Thus we choose to trace the surface of the cylinder with $$ \mathbf{\vec{r}}(x,z) = \begin{pmatrix} x \\ (x + 3)^2 + z^2 - 9 \\ z \end{pmatrix}, $$ where the unit outward normal on the cylinder is $$\mathbf{\hat{n}} = \frac{1}{\sqrt{4(x+3)^2 + 4z^2 + 1}} \begin{pmatrix} -2(x+3) \\ 1 \\ -2z \end{pmatrix} = \frac{1}{\sqrt{37}} \begin{pmatrix} -2(x+3) \\ 1 \\ -2z \end{pmatrix} . $$ The flux is then \begin{align} \Phi &= \underbrace{\iint_{\Sigma_1} \mathbf{\vec{V}} \cdot \begin{pmatrix} 0 \\ -1 \\ 0 \end{pmatrix} \, d\sigma}_{y = -1} + \overbrace{\iint_{\Sigma_2} \mathbf{\vec{V}} \cdot \mathbf{\hat{n}} \, d\sigma}^{\mathbf{\vec{V}}_2\text{ cannot contribute}} + \underbrace{\iint_{\Sigma_3} \mathbf{\vec{V}} \cdot \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \, d\sigma}_{\text{nothing since }y = 0} \\ &= \iint_{\Sigma_1} \sqrt{x^2 + z^2} \, d\sigma + \frac{1}{\sqrt{37}}\iint_{\Sigma_2} -2z(x+3) + y\sqrt{x^2 + z^2} + 2xz \, d\sigma \\ &= \int_0^{2\pi}\int_0^3 r\sqrt{(r \cos\theta - 3)^2 + (r \sin\theta)^2} \, dr \, d\theta - \frac{6}{\sqrt{37}} \int_0^3 r^2 \, dr \, \underbrace{\int_0^{2\pi} \sin\theta \, d\theta}_{0} \\ &= 96. \end{align}
(2) Divergence Theorem
Technical information about this method can also be found in MIT's open notes, and no visualization from the site this time, but the divergence theorem here is explained in casual language. Given everything is nice, the flux of the field through the surface is $$ \iint_\Sigma \mathbf{\vec{V}} \cdot \mathbf{\hat{n}} \, d\sigma = \iiint_M \nabla \cdot \mathbf{\vec{V}} \, dV, $$ where $M$ is the bounded region contained within $\Sigma$.
Applying it to this problem, the divergence theorem takes us straight to the end result of the direct approach. \begin{align} \Phi &= \iint_\Sigma \nabla \cdot \mathbf{\vec{V}} \, dV \\ &= \iiint_M \sqrt{x^2 + z^2} \, dV \\ &= \int_{-1}^0 \int_0^{2\pi} \int_0^3 r\sqrt{(r \cos\theta - 3)^2 + (r \sin\theta)^2} \, dr \, d\theta \, dy \\ &= 96. \end{align}
From the comments, $18\pi$ falls out as the solution if $x$ is not properly shifted over from the origin. Also, Wikipedia is a fairly good source for this material as well.