Here is another method.

Since the uniform distribution on the unit disk is rotationally symmetric (invariant under orthogonal transformations), the problem can be reduced to a random walk problem in $\mathbb{R}$.

Let $Z_1,Z_2,\ldots$ be independent, and identically uniformly distributed on the unit disk, with $Z_i=(X_i,Y_i)$.
Clearly $X_1,X_2\ldots$ are i.i..d with density $p(x)=\frac{2}{\pi}\sqrt{1-x^2}1_{[-1,1]}(x)$.

For any rotationally symmetric
random vector $V=(V_1,V_2)$ with finite expectation of $|V|:=\sqrt{V_1^2+V_2^2}$  it holds that
 $$\mathbb{E}(|V|)=\frac{\pi}{2}\mathbb{E}(|V_1|)$$
 Applying this to $Z^{(n)}:=Z_1+\ldots+Z_n$ we get that

\begin{align*}
\mathbb{E}\big(\frac{Z^{(n)}}{n}\big)&=\frac{1}{n}\frac{\pi}{2}\mathbb{E}|X_1+\ldots + X_n|\\
                                 &=\frac{1}{n}\big(\frac{2}{\pi}\big)^{n-1}\,\int_{-1}^1\ldots\int_{-1}^1|x_1+\ldots+x_n|\prod_{i=1}^n\sqrt{1-x_i^2}\,dx_1\,\ldots dx_n\\
                                 \end{align*}

 The integrals can be solved explicitly for $n=1$ and $n=2$, but that seems to be as far as it goes. But (using the central limit theorem and uniform integrability) it is easy to see that $$\mathbb{E}\big(\frac{Z^{(n)}}{\sqrt{n}}\big)\longrightarrow \sqrt{\frac{\pi}{8}}$$

**Remarks:**                                                   
(1) a related problem appeared here https://mathoverflow.net/questions/300134/an-interesting-triple-integral                     
(2) for similar considerations see Feller II (1971), p. 30 ff