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$$



Bound on probabilities of the sum of uniform order statistics
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](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution):
\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$. 


