Given the densities $\varphi$ and $\psi$ of two independent random variables, the probability that their product is less than $z$ is $ \iint_{xy< z}\varphi(x)\psi(y)\,\mathrm{d}x\,\mathrm{d}y\tag{1} $ Letting $w=xy$ so that $x=w/y$ yields $ \iint_{w< z}\varphi\left(\frac{w}{y}\right)\psi(y)\,\mathrm{d}\frac{w}{y}\,\mathrm{d}y=\iint_{w< z}\varphi\left(\frac{w}{y}\right)\psi(y)\,\mathrm{d}w\,\frac{\mathrm{d}y}{y}\tag{2} $ Taking the derivative of $(2)$ with respect to $z$ gives the density of the product of the random variables to be $ \phi(z)=\int\varphi\left(\frac{z}{y}\right)\psi(y)\,\frac{\mathrm{d}y}{y}\tag{3} $ We can compute the expected value using this distribution as $ \begin{align} \mathrm{E}(Z) &=\int z\phi(z)\,\mathrm{d}z\\ &=\iint z\,\varphi\left(\frac{z}{y}\right)\psi(y)\,\frac{\mathrm{d}y}{y}\,\mathrm{d}z\\ &=\iint xy\,\varphi(x)\psi(y)\,\mathrm{d}y\,\mathrm{d}x\tag{4} \end{align} $ which is exactly what one would expect when computing the expected value of the product directly.
In the same way, we can also compute $ \begin{align} \mathrm{E}(Z^2) &=\int z^2\phi(z)\,\mathrm{d}z\\ &=\iint z^2\,\varphi\left(\frac{z}{y}\right)\psi(y)\,\frac{\mathrm{d}y}{y}\,\mathrm{d}z\\ &=\iint x^2y^2\,\varphi(x)\psi(y)\,\mathrm{d}y\,\mathrm{d}x\tag{5} \end{align} $ again getting the same result as when computing this directly.
The variance is then, as usual, $\mathrm{E}(Z^2)-\mathrm{E}(Z)^2$.