Let's look at a completely general case, where the object whose volume is to be measured can be convex, concave, or a set of smaller objects, in $D$ dimensions. Do note, however, that I am not a mathematician, and am looking at the problem at hand as one to be solved in practice; the math below is treated as a practical tool, instead of as the focus of interest.
The first step is to find the axis-aligned bounding box for the object. This does not need to be the minimum bounding box, although the minimum bounding box is most efficient. Essentially, it is just a matter of finding the coordinate range ($x_{d, min}$ to $x_{d, max}$) along each dimension $d$ that encloses the volume to be measured. The volume of this hyperbox is then
$$\bbox{V_{box} = \prod_{d=1}^D \bigl( x_{d, max} - x_{d, min} \bigr)}$$
In almost all practical cases, you can determine $x_{d, min}$ and $x_{d, max}$ from the nature of the problem. For a $D$-dimensional hypersphere centered at origin, $x_{d, min} = -r$ and $x_{d, max} = r$, where $r$ is the radius of the hypersphere.
In the leftover cases, it is best to calculate and describe the result as the intersection of the object whose volume is to be measured, and the axis-aligned bounding hyperbox you use. There is little sense in trying to construct a method that would find the axis-aligned bounding box for you; any reliable method would be too slow to be practical, even on large cluster supercomputers.
Monte Carlo integration
Generate N uniform random points within the hyperbox, finding out whether each point is inside or outside the volume.
If there are $k$ points inside (and $N-k$ outside), the estimated volume of the object is
$$V_{estimate} = V_{box} \frac{k}{N}$$
The estimate of the error in $V_{estimate}$ is
$$\bbox{\delta V_{estimate} = V_{box}\frac{\sigma_N}{\sqrt{N}}}$$
where $\sigma_N / \sqrt{N}$ is the standard error of the mean. This is not a strict bound, however, and this may be an underestimate, especially if the volume to be measured has lots of small details (that random sampling may not uncover).
Deterministic regular rectangular lattice sampling
A regular rectangular lattice is fit inside the axis-aligned bounding box. If $n_d$ is the number of lattice cells ($n_d + 1$ the number of lattice points) along each dimension $d$, with $k_d = 0, 1, \dots, n_d$, the coordinates of each sampled point are
$$\bbox{x_d = x_{d, min} + \frac{k_d}{n_d} \bigl( x_{d, max} - x_{d, min} \bigr)}$$
and there are a total of
$$\bbox{N_{cells} = \prod_{d = 1}^{D} n_d} , \quad
\bbox{N_{points} = \prod_{d = 1}^{D} (n_d + 1)}$$
cells and lattice points in the lattice, inside the axis-aligned bounding box; with each cuboid cell having $2^D$ vertices or lattice points at its corners.
There are two submethods of estimating the volume using this lattice.
Direct regular rectangular lattice sampling
Let $n$ be the number of lattice points inside the volume. The volume estimate is then
$$\bbox{V_{estimate} = V_{box} \frac{n}{N_{points}}}$$
While this is simple and straightforward, it is difficult to estimate the error in the volume in practice. (There is one, obviously, but I'll omit it, because in practice, as an estimate, it is too often misleading.)
Cell counting
Let $n_{in}$ be the number of lattice cells where all the corner lattice points are inside the volume, and $n_{out}$ the number of lattice cells where all the corner lattice points are outside the volume.
The difference, $N_{cells} - n_{in} - n_{out}$, is the number of lattice cells intersecting the surface of the volume, with at least one corner/lattice point inside, and one outside the volume. If we assume that the detailed features of the volume are roughly the same size or larger than the cells (not smaller), we can use the difference as a practical estimate of error.
By estimating that half the surface cells are within the volume, we get
$$\bbox{V_{estimate} = V_{box} \left( \frac{n_{in}}{N_{cells}} + \frac{N_{cells} - n_{in} - n_{out}}{2 N_{cells}} \right)}$$
and the error estimate in the volume estimate
$$\bbox{\delta V_{estimate} = V_{box} \left( \frac{N_{cells} - n_{in} - n_{out}}{2 N_{cells}} \right) }$$
Note, however, that this method allows an incremental refinement of the estimates, by examining only the surface cells further. (Either by subdividing each into say $2^D$ subcells by "halving" the cell cuboid along each axis, or by Monte Carlo integration, or by any other method.) The key is to start with a cell size as small as the smallest interesting detail in the volume.
Best Answer
The question is which part of the hypercube volume $C$ is also inside the hypersphere, i.e. in $S\cap C$. Since the hypersphere is fully contained in the hypercube, $S\cap C=S$. So according to Wikipedia you want
$$\frac{V(S)}{V(C)} =\frac{A^r\frac{\pi^{r/2}}{\Gamma(\frac r2+1)}}{(2A)^r} =\left(\frac{\sqrt\pi}{2}\right)^r\cdot\frac1{\Gamma(\frac r2+1)} $$
Now $\frac{\sqrt\pi}2\approx0.886<1$ and $\lim_{x\to+\infty}\Gamma(x)\to+\infty$ so yes, this does tend to zero.
Since $\Gamma\left(\frac r2+1\right)=\frac r2\Gamma\left(\frac r2\right)$ the formula for $V(S)$ agrees with what you used, except that for $r=0$ you'd have $\Gamma(0)$ undefined. But for $r\to\infty$ that's irrelevant.
This whole question reminds me of this post of mine about hypershperes and hypercubes…