Recreational probability problem: distance between two darts My answer to this question has inspired me to create this recreational probability problem:
A person is throwing darts at a circular board of radius $r$. Nothing is known about this player’s skill. The player throws two darts, and they both land on the board. Let $L$ be the continuous random variable representing the distance between the two darts. Determine


*

*the cumulative distribution function $$P_L(u) = \operatorname{P}(u\le L\le 2r)$$

*the expected value of $L$, $$\operatorname{E}(L)$$


Feel free to diverge from the prescribed notation. I will post my own answer in 48 hours later this week eventually.
Good luck, and have fun!
 A: 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$$

A: Hopefully this will get the ball rolling, and maybe others will join in ¯\_(ツ)_/¯. I did not realize how tricky it would be for me to concrete my thought processes in regard to this problem, but now that I am on autumnal/Thanksgiving break, I have had time to sit down and tackle it head-on!

$\large{\mathbf{\text{Conventions & Symbols}}}$
Let us model the dart board as a two-dimensional set of polar points, $$D=\{(\rho,\theta) : \rho\le r\}$$ where, by convention, the constraints $\rho\ge0$ and $\theta\in[0,2\pi)$ apply. Let $A$ be the point $(\rho_A,\theta_A)$ at which the first dart lands, and let $B$ be the point $(\rho_B,\theta_B)$ at which the second dart lands.
Let $m(S)$ be the measure of the area of the region $S\subseteq D$; for example, $m(D)=\pi r^2$. Since we do not know anything about the player’s skill, we will assume that the probability of the dart landing in any given region is the same as the dart landing in any other given region if the regions are the same size. Thus we will assume the model $$\operatorname{P} (A\in S) = \frac{m(S)}{m(D)} = \frac{ m(S) }{ \pi r^2 }$$
$\large{\mathbf{\text{Area Surrounding $\boldsymbol A$}}}$
In order to define $P_L(u)$, one must find the area of the region $S_A(u)$, which is the intersection of a circle of radius $u$ surrounding $A$ and the dartboard $D$. If $\rho_A+u\le r$, then the circle is completely contained by the dartboard, and the area of $S_A(u)$ is simply $\pi u^2$.

If, however, $\rho_A+u\gt r$, then the region is not a circle.

You can interact with these visuals on Desmos.
This answer includes the formula for calculating the area of the overlapping region
$$X=u^{2}\arccos\left(\frac{\rho_A^{2}+u^{2}-r^{2}}{2\rho_A u}\right)
+ r^{2}\arccos\left(\frac{\rho_A^{2}+r^{2}-u^{2}}{2\rho_A r}\right)
- \frac12 \sqrt{(-\rho_A+u+r)(\rho_A-u+r)(\rho_A+u-r)(\rho_A+u+r)}
$$
Hence $$m\bigl(S_A(u)\bigr) =
\begin{cases}
\pi u^2 & \rho_A\in[0,r-u] \\
X & \rho_A\in(r-u,r) \\
\end{cases}$$
$\large{\mathbf{\text{Hence Determining $\boldsymbol{P_L(u)}$}}}$
Now, consider that $B$ falls outside the region $S_A(u)$ or on its boundary. This would imply that $u\le L\le 2r$. (However, it should be noted that $L$ will only equal $2r$ if $A$ lands on the boundary of $D$ and $B$ lands opposite $A$; nevertheless, the inequality is still technically accurate.)
By the model of calculating probability outlined in $\mathbf{\text{Conventions & Symbols}}$,
$$\operatorname{P}(u\le L\le 2r) = P_L(u) = 1-\frac{m\bigl(S_A(u)\bigr)}{m(D)} = 1-\frac{m\bigl(S_A(u)\bigr)}{\pi r^2}$$
$\large{\mathbf{\text{Hence Computing $\boldsymbol{\operatorname{E}(L)}$}}}$
According to this post, if $X\ge0$ then $\operatorname{E}(X)=\int_0^\infty \operatorname{P}(X>x)\,dx$.  Given the limits for $L$, I believe that this could be logically manipulated to imply
$$\begin{align}
\operatorname{E}(L) &= \int_0^{\infty}\operatorname{P}(L>x)\,dx \\
&= \int_0^\infty \operatorname{P}(x\le L)\, dx \\
&= \int_0^{2r}\operatorname{P}(x\le L \le 2r)\, dx \\
&= \int_0^{2r}\operatorname{P}(u\le L \le 2r)\, du \\
&= \int_0^{2r}P_L(u)\, dx \\
\end{align}$$
I don’t really care to give this in a closed form, and I don’t have access to a graphing calculator with a strong enough computer algebra system to arrive at a numerical value for $\operatorname{E}(L)$.

Well, I look forward to hearing what others think! If you have found any errors, please let me know or feel free to indicate so in the answer section so that I can up-vote it. If you have any alternative solutions (perhaps using vectors), then please contribute that as well!
