I am looking into the problem where $X$ is a random variable defined on $\mathbb{R}^d$ with $X\sim N(0, I_d)$ and I want to find, $$E\left\{\frac{\langle X-u_1, X-u_2\rangle}{\|X-u_1\|\|X-u_2\|}\right\},$$
interms of $u_1$ and $u_2$, where $\langle \cdot,\cdot\rangle$ and $\|\cdot\|$ denotes the usual inner product and norm on $\mathbb{R}^d$ respectively and $u_1,u_2\in\mathbb{R}^d$ are fixed points.
I do not know whether there is any existing similar problem that has been discussed in this platform or anywhere else; I could not find any. Any help would be greatly appreciated!
Not sure about a closed form, but you can turn the $d$ dimensional integral into a 2 dimensional one by using various integral representations of ${1\over |x|}$, for example: $${1\over \|X-u\| } \sim \int d\lambda e^{-\lambda^2 \|X-u\|^2}.$$ Working from
$$ \mathbb{E}\int d\lambda d\mu (X-u_1)(X-u_2) e^{-\lambda^2 \|X-u_1\|^2 -\mu^2\|X-u_2\|^2}$$
The d-dimensional $dX$ integral now becomes doable, at the expense of the two new variables. I tried pushing a bit, but it becomes rather tedious. The good thing about it is it may allow you to better identify the $u_1^2,u_2^2,u_1\cdot u_2$ dependencies.