Suppose that $X_1,\dotsc , X_n$ are independent random variables uniformly distributed in $[0,1]$. Denote by $Y_i=X_{(i)}$ its order statistics.

The random vector $(Y_1,\dotsc, Y_n)$ has distribution

$$
p(y_1,\dotsc, y_n)=\begin{cases}
n!, & 0\leq y_1\leq \cdots \leq y_n\leq 1, \\
0, & {\rm otherwise}.
\end{cases}
$$

The random variable $Y_k$ has distribution

$$
p_k(y_k)=n!\int_{\substack{0\leq y_1\cdots \leq y_k\\ y_k\leq y_{k+1}\leq \cdots \leq y_n\leq 1}}dy_1\cdots dy_{k-1}dy_{k+1}\cdots dy_{n}
$$
$$
= n!\left(\int_{0\leq y_1\leq \cdots \leq y_{k-1} \leq y_k} dy_1\cdots dy_{k-1}\right)\left(\int_{y_k\leq y_{k+1}\leq \cdots \leq y_{n} \leq 1} dy_{k+1}\cdots dy_{n}\right)
$$

$$
= \frac{n!}{(k-1)!(n-k)!} y_k^{k-1}(1-y_k)^{n-k}
$$

Thus $p_k$ is a $B(k,n+1-k)$-distribution. $\newcommand{\bE}{\mathbb{E}}$ We have

$$
\bE[Y_k]=\frac{n!}{(k-1)!(n-k)!}\int_0^1 y^k(1-y)^{n-k} dy=\frac{n!}{(k-1)!(n-k)!}\frac{\Gamma(k+1)\Gamma(n-k+1)}{\Gamma(n+2)}
$$

$$
= \frac{n!}{(k-1)!(n-k)!}\frac{k!(n-k!)}{(n+1)!}=\frac{k}{n+1}.
$$

$$
\bE[Y_k^2]=\frac{n!}{(k-1)!(n-k)!}\int_0^1 y^{k+1}(1-y)^{n-k} dy= \frac{n!}{(k-1)!(n-k)!}\frac{\Gamma(k+2)\Gamma(n-k+1)}{\Gamma(n+3)}
$$

$$
=\frac{n!}{(k-1)!(n-k)!}\frac{(k+1)!(n-k)!}{(n+2)!}=\frac{k(k+1)}{(n+1)(n+2)}.
$$

$\DeclareMathOperator{\Var}{Var}$. Hence

$$
\Var[Y_k]= \frac{k(k+1)}{(n+1)(n+2)}-\left(\frac{k}{n+1}\right)^2=\frac{k}{n+1}\left(\frac{k+1}{n+2}-\frac{k}{n+1}\right)=\frac{k(n+1-k)}{(n+1)^2(n+2)}.
$$

We deduce

$$
\bE[\Delta_n]=\sum_{k=1}^n \Var[Y_k]=\frac{1}{(n+1)^2(n+2)}\sum_{k=1}^n k(n+1-k).
$$

We have

$$
\sum_{k=1}^nk(n+1-k)= (n+1)\sum_{k=1}^n kn-\sum_{k=1}^n k(k-1)
$$

$$
= \frac{n^2(n+1)}{2}-\sum_{k=1}^nk(k-1).
$$

Now we write

$$
k(k-1)=\frac{1}{3}\Big(\; k^3-(k-1)^3-1\;\Big)
$$

so

$$
\sum_{k=1}^n k(k-1)=\frac{1}{3}\Big(\; n^3 -n\;\Big)=\frac{n(n+1)(n-1)}{3}.
$$

We deduce

$$
\bE[\Delta_n]= \frac{1}{(n+1)^2(n+2)}\left(\;\frac{n^2(n+1)}{2}-\frac{n(n+1)(n-1)}{3}\;\right)=\frac{n(n+1)(n+2)}{6(n+1)^2(n+2)}=\frac{n}{6(n+1)}.
$$
Suppose now that $X_1,\dotsc, X_n$ are i.i.d. with common distribution $p(x)dx$ To compute $\bE[\Delta_n^]$ the following alternate description of $\Delta_n$ could help.

$$\Delta_n= \sum_{k=1}^n Y_k^2-2\sum_{k=1}^n \bE[Y_k]Y_k+\sum_{k=1}^N\bE[Y_k]^2 $$
$$
=\sum_{k=1}^n X_k^2-\sum_{k=1}^n \bE[Y_k]Y_k+\sum_{k=1}^N\bE[Y_k]^2.
$$

Hence
$$
\bE[\Delta_n]=\bE\left(\sum_{k=1}^n X_k^2\right)-\sum_{k=1}^n \bE[Y_k]^2
=n\bE[X_1^2]-\sum_{k=1}^n \bE[Y_k]^2.
$$