The first thing to do is to simplify the problem by using spherical coordinates. Because you're really only interested in points on the unit sphere, you can write $x=\cos\phi$, $y=\sin\theta\sin\phi$, and $z=\cos\theta\sin\phi$, so that $x^2+y^2+z^2$ is automatically $1$. (Note that, while the polar axis is conventionally taken to be the $z$ axis, I have taken it to be the $x$ axis, since that is the "special" one in your original definition of $D$). Then you have

\begin{align*}
D&=\sqrt{(\cos\phi-1)^2+\sin^2\theta\sin^2\phi+\cos^2\theta\cos^2\phi}\\
&=\sqrt{\cos^2\phi-2\cos\phi+1+\sin^2\phi}\\
&=\sqrt{2-2\cos^2\phi}=2\sqrt{\frac{1-\cos^2\phi}{2}}.
\end{align*}

We can use a half-angle formula to simplify this even further, finding that
\begin{equation*}
D=2\sin\frac{\phi}{2}.
\end{equation*}

To find the probability density function for $D$, start by finding the CDF:

\begin{align*}
P(D\le b)=P\left(2\sin\frac{\phi}{2}\le b\right)=P\left(\frac{\phi}{2}\le \arcsin\frac{b}{2}\right)=P\left(\phi\le 2\arcsin\frac{b}{2}\right).
\end{align*}

Since the point $(x,y,z)$ is to be uniformly distributed on the unit sphere, the probability density function $f(\theta,\phi)$ should be proportional to the surface area of an infinitesimal patch of the sphere of width $d\theta$ and height $d\phi$:

\begin{equation*}
f(\theta,\phi)\ d\theta\ d\phi=\frac{1}{4\pi}\sin\phi\ d\theta\ d\phi
\end{equation*}

The CDF for $D$ is then given by

\begin{align*}
P(D\le b)&=P\left(\phi\le 2\arcsin\frac{b}{2}\right)=\int_0^{2\arcsin b/2}\int_0^{2\pi} f(\theta,\phi)\ d\theta\ d\phi\\
&=\int_0^{2\arcsin b/2}\int_0^{2\pi} \frac{1}{4\pi}\sin\phi\ d\theta\ d\phi=\frac{1}{2}\int_0^{2\arcsin b/2} \sin\phi\ d\phi\\
&=-\frac{1}{2}\cos\left(2\arcsin\frac{b}{2}\right)+\frac{1}{2}\cos(0)=\frac{1}{2}-\frac{1}{2}\cos\left(2\arcsin\frac{b}{2}\right).
\end{align*}

We can find the PDF for $D$ by differentiating this expression:

\begin{align*}
f(b) &= \frac{d}{db}\left(\frac{1}{2}-\frac{1}{2}\cos\left(2\arcsin\frac{b}{2}\right)\right) = \frac{1}{2}\sin\left(2\arcsin\frac{b}{2}\right)\left(\frac{1}{\sqrt{1-\left(\frac{b}{2}\right)^2}}\right)\\
&=\left(\frac{b}{2}\right)\left(\sqrt{1-\left(\frac{b}{2}\right)^2}\right)\left(\frac{1}{\sqrt{1-\left(\frac{b}{2}\right)^2}}\right)=\frac{b}{2}.
\end{align*}

This is a surprisingly simple result. In fact, it means that $D^2$ is uniformly distributed on $[0,4]$.

There's another interesting solution to the same problem here: http://godplaysdice.blogspot.com/2011/12/solution-to-distance-between-random.html.
This uses a result known as Archimedes' Hat-Box Theorem to deduce that $D^2$ is uniformly distributed on $[0,4]$, from which the distribution of $D$ follows.