I will prove the stronger result without the subtraction of $\bar{f}$. As we know $\int |f|^2 = \int |f - \bar{f}|^2 + \int |\bar{f}|^2$, the result without subtracting $\bar{f}$ would imply what you want.
Given $r_1, r_2$, let $\rho = \sqrt{r_2r_1}$ and $\rho_2 = r_2 / \rho$. Note that $r_1 / \rho = 1/\rho_2$. Let $g(x) = f(x / \rho)$, then the desired inequality becomes equivalent to
$$ \int |g|^2 \overset{?}{\lesssim} (\rho_2 - \rho_1)^2 \int |\nabla g|^2 $$
for all $g\in C^\infty_c(B_{\rho_2} \setminus B_{1/\rho_2})$.
Case 1:
When $\rho_2 \geq 2$, we have $\rho_2 -\rho_1 \geq \rho_2 / 2$, so our desired result would follow if
$$ \int |g|^2 \overset{?}{\lesssim} \rho_2^2 \int |\nabla g|^2 $$
But this holds for ALL $C^\infty_c(B_{\rho_2})$ functions (by scaling from the $B_1$ case), and so it holds in our case too.
We see that our enemy is when $\rho_2 \searrow 1$.
Case 2: Now suppose $\rho_2 \leq 2$, then we have that every $r\in [1/\rho_2,\rho_2]$ is within a factor of $2$ of $1$. So our integration (in polar coordinates)
$$ \int_{1/\rho_2}^{\rho_2}\int_{\mathbb{S}^{d-1}} \ldots r^{d-1} dr d\theta \approx \int_{1/\rho_2}^{\rho_2} \int_{\mathbb{S}^{d-1}} \ldots dr d\theta $$
(we can remove the $r$ weight coming from the spherical area). But then we have that
$$ \int_{1/\rho_2}^{\rho_2} |g|^2 dr \lesssim (\rho_2 - 1/\rho_2)^2 \int_{1/\rho_2}^{\rho_2} |\partial_r g|^2 dr $$
using the one-dimensional Poincare inequality. Now observing that pointwise $|\nabla g|^2 \geq |\partial_r g|^2$, and that we can add on the outside a spherical integral, we have that the desired result also holds for $\rho_2 \leq 2$.
Conclusion: What you asked about is true.
Incidentally, the "mean-zero" version of the claim is false without the requirement that $f\in C^\infty_c$. To see a counterexample, set $Y$ to be a non-trivial spherical harmonic with eigenvalue $\lambda > 0$, so $\int_{\partial B_1} Y dS = 0$; normalize so that $\int_{\partial B_1} Y^2 dS = 1$. Consider the function given in polar coordinates as
$$ f(r,\theta) = Y(\theta)$$
Note that $f$ no longer has compact support on $B_{r_2}\setminus B_{r_1}$.
Note also that $f$ has also zero mean (inherited from $Y$). Now
$$ \int |f - \bar{f}|^2 = \int |f|^2 = \int_{r_1}^{r_2} r^{d-1} dr $$
while
$$ (r_2 - r_1)^2 \int |\nabla f|^2 = (r_2 - r_1)^2 \int_{r_1}^{r_2} \lambda r^{d-3} dr $$
Now if we fix $r_1$ and write $r_2 = r_1 + \mu$, we have, as $\mu\searrow 0$ that
$$ \int |f - \bar{f}|^2 = r_1^{d-1} \mu + O(\mu^2) $$
while
$$ (r_2 - r_1)^2 \int |\nabla f|^2 = r_1^{d-3} \mu^3 + O(\mu^4) $$
as the latter shrinks to zero faster than the former, an inequality of the form
$$ \int |f - \bar{f}|^2 \overset{?}{\leq} C (r_2 - r_1)^2 \int |\nabla f|^2 $$ cannot be true.
It may be worth asking what breaks in the proof given in part 1 when we remove the condition that $f$ is compactly supported; the reason is subtle.
The analysis of Case 1 is still largely ok, if you make some changes. That is to say, the problem is not with the case where $(r_2 - r_1)$ is large. This is also borne out by the counterexample which only kicks in when $(r_2 - r_1) \ll 1$.
The key is in the analysis of Case 2. If we don't have the compact support assumption, we cannot appeal to the compact support version of the 1 dimensional Poincare inequality. We would like to replace this by the mean zero version. However, we run into a problem: when doing so we are required to subtract the mean of $r\mapsto f(r,\theta)$ for a fixed angle, and for different $\theta$s the various "one dimensional means" may be different. This means you cannot integrate the "one dimensional" inequality over the sphere to obtain the higher dimensional one. This is also seen in our radially symmetric counterexample where for each fixed $\theta$, the function is constant and hence equal to its mean.