Integral $\int_0^1 \int_0^1 \cdots \int_0^1\frac{x_{1}^2+x_{2}^2+\cdots+x_{n}^2}{x_{1}+x_{2}+\cdots+x_{n}}dx_{1}\, dx_{2}\cdots \, dx_{n}=?$

Ah, we may simply integrate by parts!

Denote $h(t)=(2 - e^{-t}( 2 + 2t+t^2))(1-e^{-t})^{n-1}$. Integrating by parts $(n+1)$ times we get $$\int_0^\infty h(t)t^{-n-2}dt=\frac1{(n+1)!}\int_0^\infty h^{(n+1)}(t)t^{-1}dt,$$ all off-integral terms vanish since $h(0)=h'(0)=\dots=h^{(n+1)}(0)=0$ and at infinity we have $h(t)=O(1)$, $h^{(i)}(t)=o(1)$ for $i>0$. We have $h^{(n+1)}(t)=\sum_{k=1}^n a_ke^{-kt}+u(t)$, where $u(t)=t\times \text{polynomial}(e^{-t},t)$. Note that $\sum a_k=h^{(n+1)}(0)=0$, thus $$\int_0^\infty \sum a_ke^{-kt}t^{-1}dt=\int_0^\infty \sum a_k(e^{-kt}-e^{-t})t^{-1}dt=-\sum a_k \log k$$ by Frullani integrals. It is easy to see that $$a_k=(-1)^{k+n+1}\binom{n-1}{k-1}k^{n-1}(n^2+n-2k).$$ It remains to evaluate $\int_0^\infty u(t) t^{-1}dt$. We have $$u(t)=-(2t+t^2)(g(t))^{(n+1)}-2(n+1)t(g(t))^{(n)},\,g(t)=e^{-t}(1-e^{-t})^{n-1}.$$ Therefore $$\int_0^\infty u(t) t^{-1}dt=-2\int g^{(n+1)}(t)dt-2(n+1)\int g^{(n)}(t)dt-\int_0^{\infty}tg^{(n+1)}(t)dt.$$ We get $\int (g^{(n)}+tg^{(n+1)})=tg^{(n)}$, and the definite integral against $(0,\infty)$ equals 0. It remains $\int_0^\infty -2g^{(n+1)}-(2n+1)g^{(n)}=2g^{(n)}(0)+(2n+1)g^{(n-1)}(0)$. We have $$g(t)=(1-t+\dots)t^{n-1}(1-t/2+\dots)^{n-1}=t^{n-1}-\frac{n+1}2t^n+\dots,$$ $g^{(n-1)}(0)=(n-1)!$, $g^{n}(0)=-\frac{(n+1)!}2$, $2g^{(n)}(0)+(2n+1)g^{(n-1)}(0)=-(n^2-n-1)(n-1)!$, confirming the guess of Sylvain JULIEN.

PREVIOUS NON COMPLETE VERSION

Still not a complete answer, but a method to be completed or improved.

Denote $f_n(t)=\frac{2 - e^{-t}\left ( 2 + 2t+t^2 \right )}{t^3}\left ( \frac{1-e^{-t}}{t} \right )^{n-1}$. Assume that we have found the coefficients $c_1,c_2,\dots,c_n$ and polynomials $g_1,\dots,g_n$ such that $$f_n(t)=2t^{-n-2}+\sum_{k=1}^n(g_k(t^{-1})e^{-kt})'+\sum_{k=1}^n c_k \frac{e^{-kt}}t.\,\,(*)$$ Then $\sum c_k=0$ (else $f_n$ would have a non-zero residue at 0, which is absurd). We have $\int_0^\infty \sum_{k=1}^n c_k \frac{e^{-kt}}t dt=-\sum c_k\log k$ by Frullani integrals. The integral of the other part of our sum $(*)$ is minus the limit at zero of the function $-\frac2{n+1}t^{-n-1}+\sum g_k(t^{-1})e^{-kt}$. The limit must exist, since the initial integral converges.

Now I cheat a bit. Note that if we write $f_n(t)=\sum_{k=0}^n q_k(1/t) e^{-kt}$ for polynomials $q_k$, then $c_k$ equals to the residue of the $q_k(1/t) e^{-kt}$, which is pretty computable. If I am not mistaken, $$q_k(t^{-1})=(-1)^kt^{-n-2}\left(2\binom{n-1}k+\binom{n-1}{k-1}(2+2t+t^2)\right)=\\=(-1)^kt^{-n-2}\left(2\binom{n}k+\binom{n-1}{k-1}(2t+t^2)\right).$$ Thus $$c_k=2(-1)^{k+n+1}\binom{n}k\frac{k^{n+1}}{(n+1)!}+2(-1)^{k+n}\binom{n-1}{k-1}\frac{k^{n}}{n!}+(-1)^{k+n+1}\binom{n-1}{k-1}\frac{k^{n-1}}{(n-1)!}=\\ =(-1)^{k+n+1}\frac{k^{n-1}(n^2+n-2k)}{n(n+1)(n-k)!(k-1)!}.$$ This matches a coefficient of $\log 3$ for $n=4$ from Shahrooz Janbaz's answer, you may check others for be sure. It remains to prove Sylvain JULIEN's guess for the rational part.

Here is another approach, which also gives the rational term.

(I) To see how it works let $n\geq 2$ and consider first the simpler case \begin{align*} \mathbb{E}\bigg(\frac{1}{X_1+\ldots+X_n}\bigg)=\int_0^\infty \bigg(\frac{1-e^{-t}}{t}\bigg)^n\,dt \end{align*} Using $\frac{1}{t^n}=\int_0^\infty \frac{z^{n-1}}{(n-1)!} e^{-zt}\,dz$ we find \begin{align*} \int_0^\infty \bigg(\frac{1-e^{-t}}{t}\bigg)^n\,dt&=\int_0^\infty\int_0^\infty \frac{z^{n-1}}{(n-1)!} e^{-zt}(1-e^{-t})^n\, dz\,dt\\ &=\int_0^\infty\int_0^\infty \frac{z^{n-1}}{(n-1)!} e^{-zt}(1-e^{-t})^n\, dt\,dz\\ &=\int_0^\infty \frac{z^{n-1}}{(n-1)!}\, \mathrm{Beta}(z,n+1)\,dz\\ &=n\,\int_0^\infty \frac{z^{n-1}}{z(z+1)\cdots(z+n)}\,dz \end{align*} The following observation will be the key:

Lemma: let $x_0,\ldots,x_n$ be distinct positive numbers and $k\leq n-1$. Then $$I_{k,n}(x_0,\ldots,x_n):=\int_0^\infty \frac{z^{k}}{(z+x_0)(z+x_1)\cdots(z+x_n)}\,dz=(-1)^{n+k+1}\Delta^n(x^k\log(x);x_0,\ldots,x_n)$$ where (for a real function $f$) $\Delta^n(f;x_0,\ldots,x_n)$ denotes the divided difference of $f$ corresponding to $x_0,\ldots,x_n$.

Recall that (Newton-interpolation)

(1) the divided differences are for $f$ and mutually distinct $x_0,\ldots,x_n$ are defined recursively by $\Delta^0(f;x_0)=f(x_0)$, $\Delta^n(f;x_0,\ldots,x_n)=\frac{\Delta^{n-1}(f; x_1,\ldots,x_n)-\Delta^{n-1}(f; x_0,\ldots,x_{n-1})}{x_n-x_0}$

(2) they are explicitly given by $$\Delta^n(f;x_0,\ldots,x_n)=\sum_{i=0}^n \frac{f(x_i)}{\prod_{j\neq i} (x_i-x_j)}\;\;(**)$$ (3) $$\Delta^n(f;x,x+1,\ldots,x+n)=\frac{1}{n!} \sum_{i=0}^n {n \choose i} (-1)^{n-i} f(x+i)$$

Proof of the lemma: For $k=0, n=1$ we have $$\int_0^\infty \frac{1}{(z+x_0)(z+x_1)}=\frac{\log(x_1)-\log(x_0)}{x_1-x_0}=\Delta^1(\log(x);,x_0,x_1)$$ For $k=0,n>1$ the repeated use of $\frac{1}{(z+a)(z+b)}=\frac{-1}{b-a}\left(\frac{1}{z+b}-\frac{1}{z+a}\right)$ shows that $(-1)^{n+1}I_{0,n}(x)=\Delta^n(\log(x),x)$. The validity for $k>0,n=k+1$ follows from $\frac{z}{z+b}=1-\frac{b}{z+b}$ and $(**)$. End of proof.

The lemma and (3) now give that \begin{align*} \int_0^\infty \bigg(\frac{1-e^{-t}}{t}\bigg)^n\,dt &=n \Delta^{n-1}(x^{n-2}\log(x);x+1,x+2,\ldots,x+n)\\ &=\frac{n}{(n-1)!}\sum_{i=0}^{n-1}{ n-1 \choose i} (-1)^{n-1-i} (i+1)^{n-2}\log(i+1) \end{align*}

(II) Now let $n\geq 1$ and consider $$ Q_{n+1}:=\mathbb{E}\bigg(\frac{X_1^2+\ldots + X_{n+1}^2}{X_1+\ldots+X_{n+1}}\bigg)$$ Write $$Q_{n+1}=(n+1)\int_0^1 u^2 \bigg(\int_0^\infty e^{-ut}\bigg(\frac{1-e^{-t}}{t}\bigg)^n\,dt\bigg)\,du$$ Proceeding as above shows that for $u>0$ \begin{align*} \int_0^\infty e^{-ut}\bigg(\frac{1-e^{-t}}{t}\bigg)^n\,dt &=\frac{n}{n!} \sum_{i=0}^n {n \choose i} (-1)^{n-i} (u+i)^{n-1}\log(u+i) \end{align*} Thus $$q_{n+1}:=\frac{Q_{n+1}}{n(n+1)}=\frac{1}{n!} \sum_{i=0}^n {n \choose i} (-1)^{n-i} \int_0^1 u^2(u+i)^{n-1}\log(u+i)\,du$$ Now expand $u^2=(u+i-i)^2$ and integrate partial to find that $q_{n+1}=L-R$ where $L=\frac{1}{n!}\sum_{i=0}^n{n\choose i}(-1)^{n-i} \ell(i),\;R=\frac{1}{n!}\sum_{i=0}^n {n \choose i}(-1)^{n-i} r(i)$ with \begin{align*} \ell(i)=&\frac{1}{n+2}\left((i+1)^{n+2}\log(i+1)-i^{n+2}\log(i)\right)\\ &-\frac{2i}{n+1}\left((i+1)^{n+1}\log(i+1)-i^{n+1}\log(i)\right)\\ &+\frac{i^2}{n}\left((i+1)^{n}\log(i+1)-i^n\log(i)\right)\\[0.2cm] r(i)=&\frac{1}{(n+2)^2}\left((i+1)^{n+2}-i^{n+2}\right)\\ &-\frac{2i}{(n+1)^2}\left((i+1)^{n+1}-i^{n+1}\right)\\ &+\frac{i^2}{n^2}\left((i+1)^{n}-i^n\right) \end{align*} Collecting terms in the logarithmic part $L$ shows that the coefficient of $\log(i)$ in $q_{n+1}$ is $c_{i,n+1} =(-1)^{n-i} \frac{i^n}{n}{ n \choose i-1}\left(\frac{n^2+3n+2-2i}{(n+2)!}\right)$, matching Fedor Petrov's answer.

The rational part $R$ can readily be summed (note that only the contributions of the powers $n$ and $n+1$ (with coefficients $-\frac{1}{n(n+1)(n+2)}$ and $\frac{2}{n(n+1)(n+2)}$) need to be considered) to give $$R=\frac{1}{n+2}-\frac{1}{n(n+1)(n+2)}\;\;,$$ so that the rational term of $Q_{n+1}$ is $-\frac{n(n+1)-1}{n+2}=-\frac{(n+2)(n-1)+1}{n+2}$, confirming Sylvain JULIEN's guess.

A little long as a comment:

Thanks to Maple, some first few terms are as follows:

$n=1:$ $I=\frac{1}{2}$

$n=2:$ $I=\frac{2^2\ln(2)}{3}-\frac{1}{3}$

$n=3:$ $I=-2^3\ln(2)+\frac{3^3\ln(3)}{4}-\frac{5}{4}$

$n=4:$ $I={2^6\ln(2)}-\frac{189\ln(3)}{5}-\frac{11}{5}$

$n=5:$ $I=-\frac{5\times 584\ln(2)}{9}+3^4\ln(3)+\frac{5^4\ln(5)}{36}-\frac{19}{6}$

$n=6:$ $I=\frac{6\times 11696\ln(2)}{63}+\frac{6\times 405\ln(3)}{14}-\frac{6\times 6250\ln(5)}{63}-\frac{29}{7}$

and so on. It seems there is some patterns which may help to obtain closed formula for the integral.