The work below is in many ways similar to what you've done, but I derive some closed formulas and numerical approximations.
For the sake of convention, I also work with the usual cumulative distribution function, meaning $F_X(x)=\mathbb{P}(X\leq x)$.
Your version reverses the inequality, but that's simply $1-F_X(x)$.
Let $D_R(p)\subset \mathbb{R}^2$ be the disk of radius $R$ centered on $p\in\mathbb{R}^2$.
Let $X,Y$ be uniform random variables on $\mathcal{D}=D_r(0)$.
We wish to calculate $\mathbb{P}(\lVert X-Y \rVert\leq u)$.
Observe that
\begin{align}
\mathbb{P}(\lVert X-Y \rVert\leq u)
&= \mathbb{P}( X - Y \in D_u(0))\\
&= \mathbb{P}\left( X \in D_u(Y)\right)\\
&=\int_{\mathbb{R}^2}\mathbb{P}\left( X \in D_u(y)\right) \cdot f_Y(y)\, dA(y)\\
&=\int_{\mathbb{R}^2} \left(\int_{D_u(y)}f_X(x)\,dA(x)\right) f_Y(y)\, dA(y),
\end{align}
where $dA$ indicates the area element.
Since $X$ and $Y$ are uniform we have that $f_X(x)=\frac{1}{\mu(\mathcal{D})}\cdot\chi_{\mathcal{D}}(x)$ and similarly for $f_Y$, where $\chi_S$ denotes the indicator function of $S$.
It follows that
\begin{align}
\mathbb{P}(\lVert X-Y \rVert\leq u)
&=\frac1{\mu(\mathcal{D})^2}\,\int_{\mathbb{R}^2} \left(\int_{D_u(y)}\chi_{\mathcal{D}}(x)\,dA(x)\right) \chi_{\mathcal{D}}(y)\, dA(y)\\
&=\frac1{\mu(\mathcal{D})^2}\,\int_{\mathbb{R}^2} \mu\left(D_u(y)\cap\mathcal{D}\right) \cdot \chi_{\mathcal{D}}(y)\, dA(y)\\
&=\frac1{\mu(\mathcal{D})^2}\,\int_{\mathcal{D}} \mu\left(D_u(y)\cap\mathcal{D}\right)\, dA(y)
\end{align}
Now, is there a pretty formula for this?
Let's try polar coordinates.
We'll get:
$$\mathbb{P}(\lVert X-Y \rVert\leq u)=
\frac1{\mu(\mathcal{D})^2}\,
\int_0^{2\pi}\,\int_0^r \rho\cdot \mu\left(D_u(\rho\,\cos\theta, \rho\,\sin\theta)\cap\mathcal{D}\right)\,d\rho \,d\theta\tag{$*$}$$
We will write $D_u$ as shorthand for $D_u(\rho\,\cos\theta, \rho\,\sin\theta)$.
For any combination of $r$, $u$ and $\rho$, there are three cases.
$\qquad(1)$: $u \leq r-\rho$
In this case, $D_u$ lies entirely within $\mathcal{D}$, so that $\mu\left(D_u\cap\mathcal{D}\right)=$ $\mu\left(D_u\right)=\pi u^2$.
$\qquad(2)$: $r-\rho < u < r+\rho$
In this case, $\partial D_u$ intersects $\partial\mathcal{D}$ at two distinct points. The asymmetric lens formed by $D_u\cap \mathcal{D}$ has area
\begin{align}
\mu(D_u\cap\mathcal{D})= \mathcal{A}(\rho,r,u)=\,
&r^2\arccos\left(\frac{\rho^2+r^2-u^2}{2\rho r}\right)
+u^2\arccos\left(\frac{\rho^2+u^2-r^2}{2\rho u}\right)\\
-\,&\frac12 \,\sqrt{(-\rho+r+u)(\rho+r-u)(\rho-r+u)(\rho+r+u)}
\end{align}
$\qquad(3)$: $u \geq r+\rho$
In this case, $\mathcal{D}$ lies entirely within $D_u$, so that $\mu(D_u\cap\mathcal{D})=\mu(\mathcal{D})=\pi r^2$.
Since the integral $(*)$ is over $\rho$, we would do well to break the previous cases in terms of $\rho$.
$\qquad(\text{I})$: $u\leq r$
Then the cases are $\rho \leq r -u$, which corresponds to $(1)$ and $\rho > r-u$, which corresponds to $(2)$.
Hence, for $0\leq u \leq r$ we have that $(*)$ becomes
\begin{align}
&\frac2{\pi r^4}\left(
\int_0^{r-u} \rho\cdot \pi u^2\,d\rho
+\int_{r-u}^r\rho\cdot\mathcal{A}(\rho,r,u) \,d\rho
\right)\\
=\,&\frac{u^2\,(r-u)^2}{r^4} +
\frac2{\pi r^4}\left(
\int_{r-u}^r\rho\cdot\mathcal{A}(\rho,r,u) \,d\rho
\right)
\end{align}
Now, Maple does find an antiderivative (I can provide if needed) for $\rho\cdot\mathcal{A}(\rho,r,u)$, and after some manipulation we get
\begin{align}
\mathbb{P}(\lVert X-Y \rVert\leq u; \,u\leq r)=\,&
\frac{u^2\,(r-u)^2}{r^4}
-\frac{3u^2}{\pi r^2}\cdot\arcsin\left(\frac{u}{2r}\right)
-\frac1\pi\cdot\arcsin\left(1-\frac{u^2}{2r^2}\right)\\
-\,&\left(\frac{u}{2\pi r^2}+\frac{u^3}{4\pi r^4}\right)\sqrt{4r^2-u^2}
+\frac{u^2}{\pi r^2}\cdot \arctan\left(\frac{u}{\sqrt{4r^2-u^2}}\right)\\
+\,&\frac12
+2\frac{u^3}{r^3}
-\frac{u^4}{r^4}
\end{align}
As a sanity check, this expression is $0$ when $u=0$, and is increasing and $<1$ for $u\in [0,r]$.
Moreover, it's easy to see that this expression is actually a function of $\frac{u}r$, so the problem is scale invariant, as expected.
$\qquad(\text{II})$: $u > r$
Then the cases are $\rho \leq u-r$, which corresponds to $(3)$ and $\rho > u-r$, which corresponds to $(2)$.
Hence, for $u > r$ we have that $(*)$ becomes
\begin{align}
&\frac2{\pi r^4}\left(
\int_0^{u-r} \rho\cdot \pi r^2\,d\rho
+\int_{u-r}^r\rho\cdot\mathcal{A}(\rho,r,u) \,d\rho
\right)\\
=\,&\frac{(u-r)^2}{r^2} +
\frac2{\pi r^4}\left(
\int_{u-r}^r\rho\cdot\mathcal{A}(\rho,r,u) \,d\rho
\right)
\end{align}
Like before, we can derive a closed form.
We'll get:
\begin{align}
\mathbb{P}(\lVert X-Y \rVert\leq u;\, u>r)=\,&
\frac{(u-r)^2}{r^2}
-\frac{3u^2}{\pi r^2}\cdot\arcsin\left(\frac{u}{2r}\right)
-\frac1\pi\cdot\arcsin\left(1-\frac{u^2}{2r^2}\right)\\
-\,&\left(\frac{u}{2\pi r^2}+\frac{u^3}{4\pi r^4}\right)\sqrt{4r^2-u^2}
+\frac{u^2}{\pi r^2}\cdot \arctan\left(\frac{u}{\sqrt{4r^2-u^2}}\right)\\
-\,&\frac12
+2\frac{u}{r}
\end{align}
As a sanity check, this expression is once again increasing and tends to $1$ as $u\to {2r}^-$.
We have that both closed forms agree when $u=r$, with value $1-\frac{3\sqrt{3}}{4\pi}\simeq 0.5865$.
Moreover, this expression is also actually a function of $\frac{u}r$.
Let $L=\lVert X-Y\rVert$ and $F_L(u)=\mathbb{P}(L\leq u)$ be its cumulative distribution function.
Since our random variable is non-negative, we can can calculate its expected value via
$$\mathbb{E}(L)=\int_0^\infty 1- F_L(u) \,du$$
Now, $F_L(u)=1$ whenever $u\geq 2r$, so the ingreal reduces to
\begin{align}
\mathbb{E}(L)
&=\int_0^{2r} 1- F_L(u) \,du\\
&=2r-\int_0^{2r}F_L(u)\,du
\end{align}
One can directly compute the integral to obtain that $\int_0^{2r}F_L(u)\,du=2r - \frac{128r}{45\pi}$.
It follows that:
$$\mathbb{E}(L)=\frac{128r}{45\pi}\simeq 0.9054\,r$$
