First of all, the integration region looks nasty:

Visualization of integration region as stated

As you suggest, let $x=u^3+v^3$, $y=u^3+w^3$ and $z=v^3+w^3$. Then

$$ \mathrm{d} u \,\, \mathrm{d} v \,\, \mathrm{d} w = \frac{4 \mathrm{d} x \,\, \mathrm{d} y \,\, \mathrm{d} z }{ ( (x+y-z)^2 (x+z-y)^2 (y+z-x)^2 )^{1/3}} $$ The region becomes $\vert x y z \vert < 1$, and integral becomes

$$ \int_{( x \cdot y \cdot z )^2 < 1} \frac{4 \mathrm{d} x \,\, \mathrm{d} y \,\, \mathrm{d} z }{ ( 4(x+y+z)^2 (x+y-z)^2 (x+z-y)^2 (y+z-x)^2 )^{1/3}} $$

From this form it seems like the integral would diverge at the origin, since the integrand scales as as $t^{-8/3}$ if all coordinates are rescaled by scale factor $t$.


Added: The integrand, as well as the region of integration, is symmetric with respect to each of $x\to-x$, $y \to -y$, $z \to -z$ as well as any permutations of $(x, y, z)$. Hence we can restrict the integration to $0 < x < y < z$ for the sake of determining the convergence. Consider $$ \int_0^\infty \mathrm{d} x \int_x^\infty \mathrm{d} y \int_y^\infty \mathrm{d} z \frac{\chi_{x y z < 1}}{{ ( (x+y+z)^2 (x+y-z)^2 (x+z-y)^2 (y+z-x)^2 )^{1/3}} } $$ In this region $x y z < 1$ implies $0<x<1 \land x < y < \frac{1}{\sqrt{x}} \land y < z < \frac{1}{x y}$:

$$ \int_0^1 \mathrm{d} x \int_x^{\frac{1}{\sqrt{x}}} \mathrm{d} y \int_y^{\frac{1}{x y}} \mathrm{d} z \frac{1}{{ ( (x+y+z)^2 (x+y-z)^2 (x+z-y)^2 (y+z-x)^2 )^{1/3}} } $$ The only factor in the denominator which could vanish here is $(x+y-z)^2$. This occurs in the following sub-region: $$ 0 < x < \frac{1}{\sqrt[3]{2}} \land x < y < -\frac{x}{2} + \sqrt{\frac{1}{x} + \frac{x^2}{4}} $$

Numerically, the integral does seem to diverge, but the proof escape me.