As a follow-up discussion with Adam, I evaluated the integral till $d=4$ for high-dimensional unitary $Z$, with ${Z}_{kl}{=}\exp(\frac{i2\pi k}{d})\delta_{kl}$. For $d{=}4$, I haven't considered the other possibility that Adam mentioned, i.e., $\tilde{U}{=}V_1\sigma_zV_1^\dagger{\oplus}e^{i\alpha}V_2\sigma_zV_2^\dagger$, $\sigma_z$ is Pauli-z.
I referred to this paper to parameterise the normalized, unitary invariant measure $d\psi$ on the manifold of $|\psi\rangle{=}\sum_{j=1}^d \sqrt{r_j}e^{i\theta_j}|j\rangle$ ($r_j{\geq}0$, $\theta_j{\in}[0,2\pi]$). The parameterisation is introduced after Section V, Lemma 4. For completeness, the parameterisation for $d\psi$ is given by the following Dirac-delta representation:
$$
\begin{align}
d{\psi}\equiv\frac{\Gamma(d)}{2\pi^d}\delta\Big(1-\sum_{j=1}^d r_j\Big)\ \prod_{j=1}^d dr_jd\theta_j.
\end{align}
$$
Here $\Gamma(d){=}(d-1)!$ is the Gamma function. First, let us check if the above parameterisation agrees with the $d{=}2$ case. In this case, $|\langle \psi |Z|\psi\rangle|{=}|r_1-r_2|$, and the integral $I_{d=2}$ becomes
$$
\begin{align}
I_{d=2}&= \frac{1}{4\pi^2}\int_{r_1{=}0}^1\int_{r_2{=}0}^1 |r_1-r_2|\delta\Big(1-r_1-r_2\Big)dr_1dr_2 \int_{0}^{2\pi}\int_{0}^{2\pi} d\theta_1d\theta_2 \\
&=\int_{r_1,r_2{=}0}^1 |r_1-r_2|\delta\Big(1-r_1-r_2\Big)dr_1dr_2\\
&{=}\frac{1}{2}.
\end{align}
$$
Then I moved on to $d{=}3$, in which case, $|\langle \psi |Z|\psi\rangle|{=}|r_1+\omega r_2 + \omega^2 r_3|$, with $\omega{=}\exp(\frac{i2\pi}{3})$. With little intention to work it out myself, I fed it to Mathematica and found $I_{d=3}=\frac{1}{3}+\frac{\ln(2+\sqrt{3})}{6\sqrt{3}}{\approx}0.460058$, not conforming to my guess, $\frac{1}{\sqrt{d+2}}$. Similarly, I calculated for $I_{d=4}{\approx}0.405806$ with a clumsy analytical expression as shown in the attached screenshot.

I cross-verified the evaluation with numerical estimations with $10^6$ samples. Here are the list of numerical results for $d{\in}[2,6]$, $[0.4998, 0.4600, 0.4058, 0.3695, 0.3409]$.
So at this stage, I don't have much intuition about the general analytical expression. However the upper bound of $\frac{1}{\sqrt{d+1}}$ is helpful. I will be open to further interesting perspectives on this.