This is by no means a complete answer to your question. This question hasn't received much attention, so I have been jotting down a few ideas.
If $f$ is a normalized distribution, we can find the mean value of $||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||$ by integrating over all possible choices of $\mathbf u$ and $\mathbf v$, weighted by the probability density at those values.
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||}=\iint_{\mathbb R^2}\iint_{\mathbb R^2}||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||f(\mathbf u)f(\mathbf v)d\mathbf vd\mathbf u$$
We can simplify by letting $\mathbf t=\mathbf q-\mathbf p$.
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||}=\iint_{\mathbb R^2}\iint_{\mathbb R^2}||\mathbf u-(\mathbf t+\mathbf v)||f(\mathbf u)f(\mathbf v)d\mathbf vd\mathbf u$$
In order to make use of the symmetry of $f$, the integral can be evaluated in polar coordinates. Some rather tedious geometry can be used to show the following identity. (In order to derive this result, construct the quadrilateral whose sides are $\mathbf u$, $\mathbf v$, and $\mathbf t$, define angles describing their directions, and use law of cosines to find the length of the fourth side.)
$$||\mathbf u-(\mathbf t+\mathbf v)||=\sqrt{u^2+v^2+t^2-2ut\cos\theta+2vt\cos\phi-2uv\cos(\theta+\phi)}=D(u,v,\theta,\phi)$$
Where $u=||\mathbf u||$, $v=||\mathbf v||$, $t=||\mathbf t||$, $\theta$ is the angle between $\mathbf t$ and $\mathbf u$, and $\phi$ is the angle between $\mathbf t$ and $\mathbf v$ (we are essentially letting $\mathbf t$ be the $\theta=\phi=0$ axis). The integral can now be rewritten in polar coordinates, utilizing the fact that $f$ does not depend on $\theta$ or $\phi$.
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||}=\int_0^{2\pi}\int_0^\infty\int_0^{2\pi}\int_0^\infty f(u)f(v)D(u,v,\theta,\phi)v\ dv\ d\phi\ u\ du\ d\theta$$
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||}=\int_0^\infty\int_0^\infty uvf(u)f(v)\left(\int_0^{2\pi}\int_0^{2\pi} D(u,v,\theta,\phi)d\phi d\theta\right)dvdu$$
Unfortunately, the inner integral is not an easy one to evaluate; I am not certain that it has a closed form. Because $f$ is symmetric, I doubt the choice of a particular (nontrivial) $f$ will make an analytic solution easily attainable.
Interestingly, if we instead consider the average squared distance, the integral collapses nicely, thanks to the periodicity of the cosine terms in $D^2$.
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||^2}=\int_0^\infty\int_0^\infty uvf(u)f(v)\left(\int_0^{2\pi}\int_0^{2\pi} D^2(u,v,\theta,\phi)d\phi d\theta\right)dvdu$$
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||^2}=4\pi^2\int_0^\infty\int_0^\infty uv(u^2+v^2+t^2)f(u)f(v)dvdu$$
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||^2}=4\pi^2\int_0^\infty\int_0^\infty (u^3vf(u)f(v)+uv^3f(u)f(v)+uvt^2f(u)f(v))dvdu$$
Each of these three terms is can be separated into a product of single intergals in $u$ and $v$. The result is compactly expressed using the shorthand $f_n=\int_0^\infty r^nf(r)dr$.
$$\overline{||(\mathbf p+\mathbf u)-(\mathbf q+\mathbf v)||^2}=4\pi^2f_1\left(2f_3+||\mathbf p-\mathbf q||^2f_1\right)$$
Taking the square root of this quantity will give the root mean squared distance, which, while not equal to the mean absolute distance, is still useful in some contexts. Also, it should have a closed form for a gaussian $f$.
