Three-dimensional case (n=3)
The Wikipedia article on area and volume element states that the area element of the unit sphere is $\sin\phi\,\mathrm d\phi\,\mathrm d\theta$ where I use $\phi$ to denote inclination and $\theta$ to denote azimuth. If you consider your fixed vector as the zenith direction, then my $ \phi$ will be the same as yours.
To turn that area element into a probability density function, you have to divide by the total area of the unit sphere, which is $4\pi$. So you get the probability density function
$$f(\phi)=\frac1{4\pi}\int_0^{2\pi}\sin\phi\,\mathrm d\theta=\frac12\sin\phi$$
You can verify that it's a probability density function by computing
$$\int_0^\pi f(\phi)\,\mathrm d\phi=1$$
The function you gave in your question only integrates to $\frac12$ in this check. Did you perhaps overlook the fact that the angle between the vectors can never exceed $\pi$? In any case, the uniform distribution looks wrong. Angles close to $\frac\pi2$ should be far more common since there is a lot more sphere area in that range than along the axis of the fixed vector.
As an interesting cross reference, equation (5) of the MathWorld article on Sphere Point Picking gives the same PDF I gave above.
Arbitrary dimension
Now let's generalize this to arbitrary dimensions. Parametrizing the $n$-sphere in a similar way can become complicated. But take another look at
$$\int_0^{2\pi}\sin\phi\,\mathrm d\theta$$
and interpret that geometrically. It's the perimeter of a circle of radius $\sin\phi$. This should generalize to higher dimensions in the obvious way: take the surface of the $n-2$-sphere (the sphere in $n-1$ dimensions but with a $n-2$-dimensional surface) of radius $\sin\phi$. Again you have to divide by the surface of the $n-1$-sphere. Looking at Wikipedia we can formulate
$$ f(\phi)=\frac{S_{n-2}}{S_{n-1}}\sin^{n-2}\phi
=\frac{(n-1)\pi^{(n-1)/2}\Gamma(\frac n2+1)}{n\pi^{n/2}\Gamma(\frac{n-1}2+1)}
\sin^{n-2}\phi $$
Wolfram Alpha isn't able to simplify the coefficient in this form, or at least not much. It could be that by using an alternate formulation with factorials, things become easier, but you pay for this with a case distinction between even and odd dimensions. At the end of the day, all that factor does is ensuring that the PDF integrates to $1$. You might as well define
$$ f(\phi)=\frac{\sin^{n-2}\phi}{\int_0^\pi\sin^{n-2}\theta\,\mathrm d\theta} $$
Wolfram Alpha can't simplify that denominator either, but for a concrete $n\in\mathbb N$, computing that integral is fairly easy. On the other hand, computing the relevant $S_i$ using the recurrence formulas given on Wikipedia is even easier, so choose our aproach however you like.
Here are the first few dimensions:
$$\begin{array}{c|ccccccccccc}
n & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 \\\hline
S_{n-2}/S_{n-1} & \frac{1}{\pi} & \frac{1}{2} & \frac{2}{\pi} & \frac{3}{4} & \frac{8}{3\pi} & \frac{15}{16} & \frac{16}{5\pi} & \frac{35}{32} & \frac{128}{35\pi} & \frac{315}{256} & \frac{256}{63\pi}
\end{array}$$
With a sharp eye you will spot many powers of two occurring in there. Namely in the numerators for even dimensions and in the denominators for odd dimensions. If you pick the even dimensions, and compute $\frac{2^{n-2}S_{n-1}}{\pi S_{n-2}}$ you end up with A000984. If you pick the odd dimensions and compute $\frac{2^{n-2}S_{n-2}}{S_{n-1}}$ you end up with A002457. I'm not sure whether using these will be any easier than computing the $S_i$ using recurrence, but I found it worth investigating.
Remarks
Two notes at the end, for others who read this but are confronted with a slightly different setup. Firstly, it doesn't really matter whether you fix one vector and have one random direction, or whether both directions are randomly and uniformly chosen. And secondly, since you only care about the angle, unit length is not a requirement. As long as the direction is uniformly distributed, you could have any distribution at all for the lengths.
