In Shell's Thermodynamics and Statistical Mechanics, he has the following statement:
I understand the presence of the $N!$ in (2.17) - if two particles Bob and Carl exchanged their $\mathbf{r},\mathbf{ p}$ coordinates, it wouldn't make a difference in the ensemble.
I don't understand the presence of the $h^{3N}$ factor. Why is that necessary?
This explanation is due to David Tong's lecture notes (see lecture 2), which I am paraphrasing: The single particle classical partition function is normalized with a factor of $h^3$, where $h$ is some constant with units $Js$. This is to make the the partition function dimensionless (to cancel the units of phase space volume element). For classical stat mech, the value of the constant is not important as partition function is only a counting tool. Physical observables will depend on log derivatives of the partition function and so the constant $h$ will always cancel out. However the actual value can be obtained by deriving the classical partition function as a limit of the quantum partition function by taking the classical limit. You can look at David Tong's derivation of this for further details but the essence is that the fundamental scale which has units of $Js$ is the Planck constant and it is also defines the smallest cell of volume in phase space (you cannot resolve higher due to uncertainty principle).
In computing the partition function, equation 2.17 is an integral over $d\bf{r}$ and $d\bf{p}$ -- over all of the phase space volume that is defined by the particle's coordinates $q$ and $p$. The quantity $d\bf{r}$ $d\bf{p}$ is a differential volume in phase space. This describes 3 coordinates for each of the $N$ particles, so you have $3N$ dimensions to integrate over, for both $q$ and $p$.
Another way of writing the differential phase space volume in the partition function would be
$$ d\mathbf{q}d\mathbf{p} = d^{3N}q\, d^{3N}p$$
We can use this to integrate over and find our partition function, but the caveat is that your result is not dimensionless, despite the fact that the partition function should have no units (e.g., the Boltzman factor in the canonical ensemble -- the integrand -- would be unitless). Therefore, we instead integrate over
$$ \frac{d^{3N}q\, d^{3N}p}{h^{3N}}$$
since the Planck constant has the equivalent unit of phase space volume. Simply, it is a correction for units, ensuring that our answer for the partition function makes good physical sense (by not having units!).
