One approach is to first calculate the moment generating function (mgf) of $Y_n$ defined by $Y_n = U_1^2 + \dotsm + U_n^2$ where the $U_i, i=1,\dotsc, n$ is independent and identically distributed standard uniform random variables.
When we have that, we can see that
$$ \DeclareMathOperator{\E}{\mathbb{E}}
\E \sqrt{Y_n}
$$
is the fractional moment of $Y_n$ of order $\alpha=1/2$. Then we can use results from the paper Noel Cressie and Marinus Borkent: "The Moment Generating Function has its Moments", Journal of Statistical Planning and Inference 13 (1986) 337-344, which gives fractional moments via fractional differentiation of the moment generating function.
First the moment generating function of $U_1^2$, which we write $M_1(t)$.
$$
M_1(t) = \E e^{t U_1^2} = \int_0^1 \frac{e^{tx}}{2\sqrt{x}}\; dx
$$
and I evaluated that (with help of Maple and Wolphram Alpha) to give
$$ \DeclareMathOperator{\erf}{erf}
M_1(t)= \frac{\erf(\sqrt{-t})\sqrt{\pi}}{2\sqrt{-t}}
$$ where $i=\sqrt{-1}$ is the imaginary unit.
(Wolphram Alpha gives a similar answer, but in terms of the Dawson integral.) It turns out we will mostly need the case for $t<0$. Now it is easy to find the mgf of $Y_n$:
$$
M_n(t) = M_1(t)^n
$$
Then for the results from the cited paper. For $\mu>0$ they define the $\mu$th order integral of the function $f$ as
$$
I^\mu f(t) \equiv \Gamma(\mu)^{-1} \int_{-\infty}^t (t-z)^{\mu-1} f(z)\; dz
$$
Then, for $\alpha>0$ and nonintegral, $n$ a positive integer, and $0<\lambda<1$ such that $\alpha=n-\lambda$. Then the derivative of $f$ of order $\alpha$ is defined as
$$
D^\alpha f(t) \equiv \Gamma(\lambda)^{-1}\int_{-\infty}^t (t-z)^{\lambda-1} \frac{d^n f(z)}{d z^n}\; dz.
$$
Then they state (and prove) the following result, for a positive random variable $X$: Suppose $M_X$ (mgf) is defined. Then, for $\alpha>0$,
$$
D^\alpha M_X(0) = \E X^\alpha < \infty
$$
Now we can try to apply these results to $Y_n$. With $\alpha=1/2$ we find
$$
\E Y_n^{1/2} = D^{1/2} M_n (0) = \Gamma(1/2)^{-1}\int_{-\infty}^0 |z|^{-1/2} M_n'(z) \; dz
$$
where the prime denotes the derivative. Maple gives the following solution:
$$
\int_{-\infty}^0 \frac{n\cdot\left(\erf(\sqrt{-z})\sqrt{\pi}-2e^z\sqrt{-z} \right)e^{\frac{n(-2\ln 2 +2 \ln(\erf(\sqrt{-z}))-\ln(-z) +\ln(\pi))}{2}}}{2\pi(-z)^{3/2}\erf(\sqrt{-z})} \; dz
$$
I will show a plot of this expectation, made in maple using numerical integration, together with the approximate solution $A(n)=\sqrt{n/3-1/15}$ from some comment (and discussed in the answer by @Henry). They are remarkably close:
As a complement, a plot of the percentage error:
Above about $n=20$ the approximation is close to exact. Below the maple code used:
int( exp(t*x)/(2*sqrt(x)), x=0..1 ) assuming t>0;
int( exp(t*x)/(2*sqrt(x)), x=0..1 ) assuming t<0;
M := t -> erf(sqrt(-t))*sqrt(Pi)/(2*sqrt(-t))
Mn := (t, n) -> exp(n*log(M(t)))
A := n -> sqrt(n/3 - 1/15)
Ex := n -> int( diff(Mn(z, n), z)/(sqrt(abs(z))*GAMMA(1/2) ),
z=-infinity..0 , numeric=true)
plot([Ex(n), A(n)], n=1..100, color=[blue, red], legend=
[exact, approx], labels=[n, expectation],
title="expectation of sum of squared uniforms")
plot([((A(n)-Ex(n))/Ex(n))*100], n=1..100, color=
[blue], labels=[n, "% error"],
title="Percentage error of approximation")
