Aloha :)
Anstatt der kartesichen Koordinaten \((x,y,z)\) sind hier die Koordinaten \((\rho,\theta,\varphi)\) vorgegeben:$$\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}(R+\rho\sin\theta)\cos\varphi\\(R+\rho\sin\theta)\sin\varphi\\\rho\cos\theta\end{pmatrix}\quad;\quad \rho\in[0;r]\;\;\;;\;\;\;\theta,\varphi\in[0;2\pi]$$Beim Übergang zu den neuen Koordinaten wird das Volumenelment \(dV=dx\,dy\,dz\) verzerrt. Über die Stärke dieser Verzerrung gibt die Funktional-Determinante Auskunft:
$$\frac{dx\,dy\,dz}{d\rho\,d\theta\,d\varphi}=\begin{vmatrix}\partial_\rho x & \partial_\theta x & \partial_\varphi x\\\partial_\rho y & \partial_\theta y & \partial_\varphi y\\\partial_\rho z & \partial_\theta z & \partial_\varphi z\end{vmatrix}=\left|\begin{array}{rrr}\sin\theta\cos\varphi & \rho\cos\theta\cos\varphi & -(R+\rho\sin\theta)\sin\varphi\\\sin\theta\sin\varphi & \rho\cos\theta\sin\varphi & (R+\rho\sin\theta)\cos\varphi\\\cos\theta & -\rho\sin\theta & 0\end{array}\right|$$$$=\left|\begin{array}{rrr}\sin\theta\cos\varphi & \rho\cos\theta\cos\varphi & -R\sin\varphi\\\sin\theta\sin\varphi & \rho\cos\theta\sin\varphi & R\cos\varphi\\\cos\theta & -\rho\sin\theta & 0\end{array}\right|+\left|\begin{array}{rrr}\sin\theta\cos\varphi & \rho\cos\theta\cos\varphi & -\rho\sin\theta\sin\varphi\\\sin\theta\sin\varphi & \rho\cos\theta\sin\varphi & \rho\sin\theta\cos\varphi\\\cos\theta & -\rho\sin\theta & 0\end{array}\right|$$Die zeite Determinante ist exakt die Funktional-Determinante beim Übergang zu Kugelkoordinaten. Ihren Wert \(\rho^2\sin\theta\) kennen wir.$$=\rho R\left|\begin{array}{rrr}\sin\theta\cos\varphi &\cos\theta\cos\varphi & -\sin\varphi\\\sin\theta\sin\varphi & \cos\theta\sin\varphi & \cos\varphi\\\cos\theta & -\sin\theta & 0\end{array}\right|+\rho^2\sin\theta$$$$=\rho R\left[-\sin\varphi(-\sin^2\theta\sin\varphi-\cos^2\theta\sin\varphi)-\cos\varphi(-\sin^2\theta\cos\varphi-\cos^2\theta\cos\varphi)\right]+\rho^2\sin\theta$$$$=\rho R[(-\sin\varphi(-\sin\varphi)(\sin^2\theta+\cos^2\theta)-\cos\varphi(-\cos\varphi)(\sin^2\theta+\cos^2\theta)]+\rho^2\sin\theta$$$$=\rho R(\sin^2\varphi+\cos^2\varphi)+\rho^2\sin\theta=\rho R+\rho^2\sin\theta$$
Damit haben wir, was wir brauchen:$$dV=dx\,dy\,dz=(\rho R+\rho^2\sin\theta)\,d\rho\,d\theta\,d\varphi$$und können das Volumen ausrechnen:
$$V=\iiint\limits_VdV=\int\limits_0^rd\rho\int\limits_0^{2\pi}d\varphi\int\limits_0^{2\pi}d\theta(\rho R+\rho^2\sin\theta)$$$$\phantom{V}=R\int\limits_0^r\rho\,d\rho\underbrace{\int\limits_0^{2\pi}d\varphi}_{=2\pi}\underbrace{\int\limits_0^{2\pi}d\theta}_{=2\pi}+\int\limits_0^r\rho^2\,d\rho\int\limits_0^{2\pi}d\varphi\underbrace{\int\limits_0^{2\pi}\sin\theta\,d\theta}_{=0}$$$$\phantom{V}=4\pi^2R\left[\frac{\rho^2}{2}\right]_0^r=4\pi^2R\frac{r^2}{2}=\boxed{2\pi^2r^2R}$$