Bound on probabilities of the sum of uniform order statistics Let $X_1,...,X_n$ be i.i.d. random variable with a uniform distribution on [0,1]. Denote by $X_{(1)}\leq X_{(2)} \leq \ldots \leq X_{(n)}$ their order statistics.
Given $k\geq 1$ and $u\in[0,k]$, I want a simple formula for 
$$
p_k(u):=\mathbb{P}[X_{(1)}+\ldots + X_{(k)}\leq u],
$$
or at least a simple lower bound on $p_k(u)$. (By simple, I mean that I don't want a k-dimensional integral...)
Note:
Of course, the case $k=1$ is easy:
$$p_1(u) = 1-(1-u)^n.$$
And, by using the representation 
$$U_{(1)},U_{(2)} = 1-Y^{1/n},1-Y^{1/n}Z^{1/(n-1)},$$
where $X,Y$ are iid $\sim U([0,1])$,
I was able to compute $p_2(u)$ (with the help of Maple):
$$p_2(u) = 1-2(1-\frac{u}{2})^n+ \big(\max(1-u,0)\big)^n$$
 A: As pointed out by user63957, the conditional distribution of $T_k:=X_{(1)}+\ldots + X_{(k)}$ given $X_{(k+1)}=x$ is that of $xS_k$, where $S_k:=U_1+\dots+U_k$ and the $U_i$'s are iid random variables (r.v.'s) uniformly distributed on $[0,1]$. So, we have the key relation 
$$T_k\overset D=X_{(k+1)}S_k,$$
assuming that the $U_i$'s and hence $S_k$ are independent of $X_{(k+1)}$, where $\overset D=$ denotes the equality in distribution. So, the distribution of $T_k$ is the product-convolution of the distributions of $X_{(k+1)}$ and $S_k$. 
The cdf of $S_k$ is given by the Irwin--Hall formula:
\begin{equation}
 P(S_k\le x)=\frac1{k!}\sum_{j=0}^k(-1)^j\binom kj(x-j)_+^k
\end{equation}
for real $x$, where $z_+:=0\vee z$ and $z_+^k:=(z_+)^k$. 
The distribution of $X_{(k+1)}$ is the beta one with parameters $k+1,n-k$. So, for $t>0$ we have 
\begin{align*}
 P(T_k\le t)&=\int_0^1 P(S_k\le t/y)P(X_{(k+1)}\in dy) \\ 
 &=\frac1{k!}\sum_{j=0}^k(-1)^j\binom kj \int_0^1 (t/y-j)_+^k\frac{n!}{k!(n-k-1)!}y^k(1-y)^{n-k-1}dy \\ 
 &=\frac{n!}{(k!)^2(n-k-1)!}\sum_{j=0}^k(-1)^j\binom kj J_{n,k,j}\big(t,(1-t/j)_+\big), \tag{1}
\end{align*} 
where $(1-t/0)_+:=0$ and (using the substitution $y=1-x$) 
\begin{align*}
 J_{n,k,j}(t,u)&:=\int_u^1 (t-j+jx)^k x^{n-k-1}dx
 =\sum_{r=0}^k\binom kr(t-j)^{k-r}j^r \int_u^1 x^{r+n-k-1}dx \\ 
 &=\sum_{r=0}^k\binom kr(t-j)^{k-r}j^r \frac{1-u^{r+n-k}}{r+n-k}.\tag{2}
\end{align*}
Thus, by formulas (1)--(2), the cdf of $T_k$ is expressed as a double sum of products of powers and binomial coefficients. In particular, for $k=1,2$ this yields the expressions obtained by the OP. 
One can also obtain some asymptotics, as follows. 
Suppose that $k\wedge(n-k)\to\infty$. Then both $X_{(k+1)}$ and $S_k$ are asymptotically normal, and we have 
\begin{align*}
T_k\overset D=X_{(k+1)}S_k&=\Big(\frac kn+(Z_1+o_P(1))\sqrt{\frac{k(n-k)}{n^3}}\Big) 
\Big(\frac k2+(Z_2+o_P(1))\sqrt{\frac{k}{12}}\Big) \\ 
&=\frac{k^2}{2n}+(Z_1+o_P(1))\frac k2\sqrt{\frac{k(n-k)}{n^3}} 
+(Z_2+o_P(1))\frac kn\sqrt{\frac{k}{12}} \\ 
&+(Z_1+o_P(1))(Z_2+o_P(1))O\Big(\frac kn\Big) \\
&\overset D=\frac{k^2}{2n}+(Z+o_P(1))\frac{k\sqrt k}{2n}\sqrt{\frac13+\frac{n-k}n}, 
\end{align*}
where $Z_1,Z_2,Z$ are iid standard normal r.v.'s, and $o_P(1)$ denotes any r.v.'s $Y_{n,k}$ that go to $0$ in probability as $k\wedge(n-k)\to\infty$. 
We see that the distribution of $T_k=X_{(1)}+\ldots + X_{(k)}$ is concentrated near $\frac{k^2}{2n}$, which was easy to predict. 
One can refine this asymptotics by using bounds on the errors of the normal approximations for $X_{(k+1)}$ and $S_k$. 
A: If you condition on the value of X(k+1), then X(1)...X(k) with the labels "removed" or randomized are IID Uniform(0,X(k+1)). So you can use IID methods to bound the conditional probability, then integrate over the Beta distribution of X(k+1).
A: (More a comment that is too long for the comment box) ...  Given a sample of size $n$ from a standard Uniform parent, one can show, by induction or otherwise, that the joint pdf of the first $k$ order statistics, is:
$$f(x_{(1)}, x_{(2)}, \dots, x_{(k)}) = \frac{n!}{(n-k)!}\left(1-x_{(k)}\right)^{n-k} \quad  \quad \text{for }\quad  0 < x_{(1)} < x_{(2)} < \dots <x_{(k)} < 1$$
... which has a neat functional form.
Your problem then reduces to:

Given joint pdf $f(x_{(1)}, x_{(2)}, \dots, x_{(k)})$, find the cdf of the sum $X_{(1)} + X_{(2)} + \dots + X_{(k)}$  

... which unfortunately will result in a $k$-part piecewise solution with kinks at $(1, 2, \dots, k)$. 
Using the mathStatica package for  Mathematica, I obtained the same solution the OP provided for $k=2$, and was also able to derive a general exact solution for the sum in the $k = 3$ case (in 3 parts). However, an exact general solution seems to get messy quite rapidly.
To illustrate the behaviour of the cdf, here is a plot of the cdf of the sum in the $k =2$ case (as $n$ increases from 2 to 10, and then with $n =30$):

The intuition, of course, is that when $n$ becomes large, the first two order statistics become very small, and so the cdf of their sum shifts to the left.
Here is the $k=3$ case:

