By dyadic decomposition, it suffices to obtain bounds on the quantities $$ \int_{|x-a_i| \sim R_i \forall i} \left[ \frac{V(x,a_1,\dots,a_n)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}}\ dx$$ which, when summed over dyadic powers of two $R_1,\dots,R_n > 0$, is bounded by $O(1)$, uniformly over $a_1,\dots,a_n$.
By permuting we may assume that $R_1 \leq \dots \leq R_n$. The most important geometric feature of the set $\{a_1,\dots,a_n\}$ turns out to be its diameter; by rescaling we can normalize this diameter to be $1$. The triangle inequality then lets us assume $R_n \gtrsim 1$. In fact the triangle inequality leaves us with just two scenarios: the "far field" scenario $$ R_1 \sim \dots \sim R_n \gtrsim 1 \quad (1)$$ and the "near field" scenario $$ R_1 \leq \dots \leq R_n \sim 1 \quad (2).$$
One could refine the range of possibilities of the near field further by inspecting the finer metric geometry of the $a_1,\dots,a_n$ beyond just using the diameter normalization, but fortunately we do not have to do so for this problem.
As implicitly observed by the OP, the volume $V(x,a_1,\dots,a_n)$ is of size at most $O(R_1 \cdots R_n)$ by multiplying all the lengths emenating from $x$. On the other hand, since $a_1,\dots,a_n$ lie in a diameter one set and $x$ is at distance $O(R_1)$ from $a_1$, we also have the bound $O(R_1)$ by the base-times-height formula. [Note that this is a significantly superior bound in the far field case, which was identified by the OP as the case where the previous estimates were poor.] Finally, the condition $|x-a_1| \sim R_1$ restricts $x$ to a set of measure $O(R_1^n)$. Thus the above integral can be bounded by $$ \lesssim \left[\frac{\min(R_1 \cdots R_n, R_1)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n.$$
In the far field case (1), this simplifies to $O(R_n^{-2})$ which is summable to $O(1)$. In the near field case (2), we bound the above by $$ \left[\frac{(R_1 \cdots R_n)^2}{(R_1 \cdots R_n)^{n+1}}\right]^{\frac{1}{n-1}} R_1^n = \prod_{j=1}^n \frac{R_1}{R_j} \leq \frac{R_1}{R_n} \sim R_1$$ which is summable to $O(1)$ since $R_1 \lesssim 1$ (and there are only $O(\log^{n-2} (1/R_1))$ choices of the remaining scales $R_2,\dots,R_n$ for any fixed $R_1$).
EDIT: one can refactor the above proof to avoid dyadic decomposition as follows.
We again normalize $\mathrm{diam}(a_1,\dots,a_n)=1$, and by translation we also normalize $a_1=0$. We now split into the far field $|x| \geq 2$ and the near field $|x| < 2$. In the far field we have $V(x,a_1,\dots,a_n) \lesssim |x|$ by the base-times-height formula and $|x-a_i| \sim |x|$ by the triangle inequality, so the contribution of the far field is $$ \lesssim \int_{|x| \geq 2} \left[ \frac{|x|^2}{|x|^{n(n+1)}}\right]^{\frac{1}{n-1}}\ dx = \int_{|x| \geq 2} \frac{1}{|x|^{n+2}}\ dx = O(1).$$ In the near field we argue as in the OP to bound things by $$ \lesssim \int_{|x| \leq 2} \frac{dx}{\prod_{i=1}^n |x-a_i|}.$$ By the triangle inequality, at least one of the $|x-a_i|$ is comparable to one. Deleting this factor and then applying the AM-GM inequality, one can bound this by $$ \lesssim \int_{|x| \leq 2} \sum_{i=1}^n \frac{dx}{|x-a_i|^{n-1}} = O(1)$$ and the claim follows.