No, we don't generally have $\mathbb E(Z) = \frac{n}{2k+1}$. For instance, for $n\le k$ we have $Z=0$ with probability $1$.
To find the expected median for large $n$ for $k=1$, with $3$ bins, we can expand the multinomial coefficients around $m=\frac n3$ using Stirling's approximation:
\begin{align}
&\binom n{m-x,m-y,m+x+y}\\
={}&\frac{n!}{(m-x)!(m-y)!(m+x+y)!}
\\\propto{}&\exp(-(m-x)\log(m-x)-(m-y)\log(m-y)-(m+x+y)\log(m+x+y)
\\
\propto{}&
\exp\left(-\frac1{2m}\left(x^2+y^2+(x+y)^2\right)-\frac1{6m^2}\left(x^3+y^3-(x+y)^3\right)\right)\;,
\end{align}
where terms independent of $x$ and $y$ and terms of order $\frac{x^2}{m^2}$, $\frac{y^2}{m^2}$, $\frac{x^4}{m^3}$ and $\frac{y^4}{m^3}$ have been dropped.
For $n\to\infty$, we can regard this as a continuous density $f(x,y)$. The probability that two $Y_i$ are equal goes to zero, so we don't have to include symmetry factors, and the expected median is
$$
\mathbb E(m-y)=m-\frac{\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\,f(x,y)}{\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,f(x,y)}\;.
$$
The $\frac1m$-fold denominator is
\begin{align}
&\frac1m\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,\exp\left(-\frac1{2m}\left(x^2+y^2+(x+y)^2\right)-\frac1{6m^2}\left(x^3+y^3-(x+y)^3\right)\right)\\
={}&\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)-\frac1{6m^{\frac12}}\left(x^3+y^3-(x+y)^3\right)\right)\\
={}&\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)\right)+O\left(n^{-\frac12}\right)\\
={}&\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\int_0^\infty r\mathrm dr\,\exp\left(-r^2\left(1+\sin\phi\cos\phi\right)\right)+O\left(n^{-\frac12}\right)\\
={}&\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\,\frac12\frac1{1+\sin\phi\cos\phi}+O\left(n^{-\frac12}\right)\\
={}&\frac\pi{3\sqrt3}+O\left(n^{-\frac12}\right)
\end{align}
(Wolfram|Alpha computation for the angular integral).
The $\frac1m$-fold numerator is
\begin{align}
&\frac1m\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\exp\left(-\frac1{2m}\left(x^2+y^2+(x+y)^2\right)-\frac1{6m^2}\left(x^3+y^3-(x+y)^3\right)\right)\\
={}&m^\frac12\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)-\frac1{6m^{\frac12}}\left(x^3+y^3-(x+y)^3\right)\right)\\
={}&\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\left(m^\frac12-\frac16\left(x^3+y^3-(x+y)^3\right)\right)\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)\right)+O\left(n^{-\frac12}\right)\;.\\
\end{align}
The first term, proportional to $n^\frac12$, vanishes:
\begin{align}
&\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)\right)\\
={}&\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\int_0^\infty r\mathrm dr\,r\sin\phi\exp\left(-r^2\left(1+\sin\phi\cos\phi\right)\right)
\\
={}&\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\sin\phi\frac{\sqrt\pi}{4(1+\sin\phi\cos\phi)^\frac32}\\
={}&0
\end{align}
(Wolfram|Alpha computation for the angular integral).
The second, constant term is finite:
\begin{align}
&-\frac16\int_0^\infty\mathrm dx\int_{-x/2}^x\mathrm dy\,y\left(x^3+y^3-(x+y)^3\right)\exp\left(-\frac12\left(x^2+y^2+(x+y)^2\right)\right)\\
={}&-\frac16\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\int_0^\infty r\mathrm dr\,r^4\sin\phi\left(\cos^3\phi+\sin^3\phi-(\cos\phi+\sin\phi)^3\right)\quad\exp\left(-r^2\left(1+\sin\phi\cos\phi\right)\right)
\\
={}&-\frac16\int_{\arctan\left(-\frac12\right)}^\frac\pi4\mathrm d\phi\frac{\sin\phi\left(\cos^3\phi+\sin^3\phi-(\cos\phi+\sin\phi)^3\right)}{\left(1+\sin\phi\cos\phi\right)^3}
\\
={}&\frac1{18}
\end{align}
(Wolfram|Alpha computation for the angular integral).
You can check that the contribution from the terms dropped in the initial approximation is also $O\left(n^{-\frac12}\right)$, so the expected median is
\begin{align}
&m-\frac{\frac1{18}}{\frac\pi{3\sqrt3}}+O\left(n^{-\frac12}\right)\\
={}&\frac n3-\frac1{2\pi\sqrt3}+O\left(n^{-\frac12}\right)\;.
\end{align}
Here's code that I used to check this result. The numerical results indicate that the shift in the expected median slightly increases with $k$, but not much. It would be interesting to see whether it has a limit for $k\to\infty$, but I don't see how to determine it other than with prohibitive integrations using hyperspherical coordinates in $2k$ dimensions in analogy with the above calculation.
\mid
instead of|
. $\endgroup$