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, $T_k\overset D=X_{(k+1)}S_k$, assumingwe 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$.