Let $f$ be in $W^{2,p}(\mathbb{R}^n)$ for $n\geq 3$ and $p>n/2$, with $f=0$ at the origin. I want to show that the integral $$\int_{B(0,r)} (f |x|^{-2})^p dV <\infty$$ for some small $r>0$. A quick application of Sobolev embedding gives $f\leq C |x|^{2-n/p}$, which then immediately shows $$\int_{B(0,r)} (f|x|^{-2+\epsilon})^p dV<\infty$$ for any $\epsilon>0$.

A Hardy-type inequality seems like it should work in this situation. For instance, from this (though this is hardly the original reference), if $p>n$, Theorem 2 holds for this instance. Let $\Omega$ be an open domain, $d(x)$ be the distance to $\partial \Omega$, and $\Omega_r = \{x\in \Omega: d(x)<r\}$. Then we get $$\int_{\Omega_r} (f d(x)^{-2})^p dV \leq c \int_{\Omega_r} |\nabla^2 f|^p dx$$ for any open domain and $r>0$ small enough. Since all we care about is near the origin, taking $\Omega = B(0,r_0)\setminus\{0\}$ gives the inequality I want.

My question is essentially if that last inequality holds more generally for $p\in (n/2,n]$. Obviously it doesn't hold for some domains, but I only care about my specific circumstance. And really, I don't care if it is exactly of that form. For instance, in Evans' PDE book, 5.8.4 Theorem 7, he proves for $n\geq 3$ and $f\in H^1(B(0,r))$ that $$\int_{B(0,r)} f^2 |x|^{-2} dV \leq C \int_{B(0,r)} |\nabla f|^2 + f^2r^{-2} dV.$$ Something like that would be more than sufficient.