To calculate the surface integral over a surface we must start by finding the surface normal to the surface at each point, ${\rm d}{\bf S}(x,y,z) = {\bf A}(x,y,z){\rm d}x{\rm d}y + {\bf B}(x,y,z){\rm d}x{\rm d}y + {\bf C}(x,y,z){\rm d}y{\rm d}z$. With this in hand we simply take the dot product with the vector field ${\bf V}$, derive the correct integration limits and perform the integrals.

The surface integral over the cube splits into $6$ integrals over each of the faces. The surface normals are here easy to find as they are always aligned with one of the coordinate axis.

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$

For example on the right face at $x=1$ we have ${\rm d}{\bf S} = {\bf e}_x\,{\rm d}y{\rm d}z$. The integration limits on this face is simply $0\leq y,z \leq 1$ so the surface integral over this face becomes

$$\int_{\text{Right face}}{\bf V}\cdot {\rm d}{\bf S} = \int_{\text{Right face}}x^2{\rm d}y{\rm d}z = \int_0^1\int_0^1 1^2{\rm d}y{\rm d}z = 1$$

Repeating this proceedure on all the other faces and summing up will give you $\int_{S}{\bf V}\cdot {\rm d}{\bf S}$.

Using the divergence theorem is often an easier way to compute surface integrals. The divergence theorem says that

$$\int_{S}{\bf V}\cdot {\rm d}{\bf S} = \int_{v}\nabla\cdot {\bf V}\, {\rm d}x{\rm d}y{\rm d}z$$

where the integral is taken over the volume bounded by the surface $S$, i.e. the interior of the cube. For your example $\nabla\cdot{\bf V} = 2(x+y+z)$ and since the cube is defined by $0\leq x,y,z\leq 1$ this gives us

$$\int_{S}{\bf V}\cdot {\rm d}{\bf S} = \int_0^1\int_0^1\int_0^12(x+y+z)\, {\rm d}x{\rm d}y{\rm d}z$$

which are just three simple volume integrals to compute. We can use the symmetry of the integrand to simply this down to $3\int_0^12x{\rm d}x\int_0^1{\rm d}y\int_0^1{\rm d}z$.