Like you said, it's easiest to consider the cube $[-1,1]^n$. I'll try to generalize the argument made in the link you submitted. We can assume by symmetry that the center of the sphere is the origin. Consider a $k-1$ sphere in $k-1$ spherical coordinates $\theta_i$ for $i=1,k-1$ with the range of $\theta_{k-1}$ being $[0,2\pi)$ and each other having range $[0,\pi]$, and a radius of $r$. For any orthonormal set $\{x^i\}_{i=1}^k$ of vectors in $\Bbb{R}^n$ the sphere can be parametrized as
$$S(\theta_1,\ldots,\theta_{k-1}) = r(\cos(\theta_1)x^1 + \sin(\theta_1)\cos(\theta_2)x^2+\cdots+\sin(\theta_1)\sin(\theta_2)\cdots\cos(\theta_{k-1})x^{k-1}\\
+\sin(\theta_1)\cdots\sin(\theta_{k-1})x^{k})$$
Each component of this must lie in $[-1,1]$. Let $x^j_i = \langle e_i,x^j\rangle$ be the $i$th component of $x^j$. The $i$th component of $S$ is then
$$S_i = r(\cos(\theta_1)x^1_i + \sin(\theta_1)\cos(\theta_2)x^2_i+\cdots+\sin(\theta_1)\sin(\theta_2)\cdots\cos(\theta_{k-1})x^{k-1}_i\\
+\sin(\theta_1)\cdots\sin(\theta_{k-1})x^{k}_i)$$
This can be written as $r\langle x_i^T, \Theta\rangle$ where $\Theta = (\cos\theta_1,\dots, \sin\theta_1\dots\sin\theta_{k-1})$. We have $||\Theta|| = 1$. Thus by the identity $\langle x,y\rangle = ||x|| ||y|| \cos(\theta) $ and the fact that $\Theta$ ranges over the entire sphere (hence $\theta$ can take the value $0$) we get $r||x|| \le 1 $ in order to fit inside the interval $[-1,1]$ and this gives
$$\sum_{j=1}^k (x_i^j)^2 \le \frac{1}{r^2}.$$
Summing over $i$ gives
$$\sum_{i=1}^n\sum_{j=1}^k (x_i^j)^2 \le \frac{n}{r^2}$$
and by changing the order or summation and using orthonormality we get
$$k\le \frac n {r^2}$$
or, as desired
$$r\le \sqrt\frac{n}{k}.$$
Scaling back down to a unit cube gives the bound that you proposed.
EDIT:
I was originally wrong about the way of attaining the bound. Googling the answer I found an article which proved a slightly general result: it's an interesting read.