Edit: The below answer addresses a separate question than that in the updated question. @WillieWong's comments above are most relevant to the question when the function depends on the difference $\mathbf{u}-\mathbf{v}$ instead of the distance.
'The function depends only on the distance between the two vectors, and their dot product.'
In other words the function can be taken to only be a function of the dot product, or distance (since $\|x-y\|^2=2-2\langle x,y \rangle$, for $x,y$ on the sphere).
One can expand the function $f$ into its Gegenbauer expansion $f(\langle x,y \rangle)=\sum\limits_{k=0}^{\infty} a_k C_{k}(\langle x,y \rangle)$ where $C_{k}(t)$ are Gegenbauer polynomials. Assuming the convergence is good enough to allow exchanging the sum and integral, the addition formula gives $$\int_{\mathbb{S}^{d-1}}\int_{\mathbb{S}^{d-1}} \sum\limits_{k=0}^{\infty} a_k C_{k}(\langle x,y \rangle) dx dy=\sum\limits_{k=0}^{\infty}\sum\limits_{m=1}^{\dim V_k}a_k b_k \left(\int_{\mathbb{S}^{d-1}} Y_{m,k}(x)dx \right)^2,$$
where spherical harmonics are denoted $Y_{m,k}$, $V_k$ is the finite-dimensional space of spherical harmonics of degree $k$, and $b_k$ are constants (which depend on how one normalizes the spherical harmonics $Y_{m,k}$). This reduces the double integral to a sum of integrals. It is not clear if $dx$ denotes the surface measure on the sphere. If so, then all terms with $k\geq 1$ are zero, meaning the integral is just the constant term in the Gegenbauer expansion of $f$.