Starting with the case where $n = m + 1$, if I understand correctly you can take a point sampled from the uniform distribution on the surface of the unit sphere in $p$-space and obtain $d(x,S)^2$ by ignoring the contribution from the first $m$ coordinates. So using the representation of the uniform distribution on the sphere as a vector of $p$ standard normal random variables: $$x = (X_1, X_2, \dots , X_p)/\sqrt{\sum\limits_{i=1}^n {X_i^2}}$$ we get $$d(x,S)^2 = \sum\limits_{i=m+1}^n{X_i^2}/\sum\limits_{i=1}^n{X_i^2}$$ which is I think, after looking at Wikipedia for $\chi^2$ distributions, a Beta distribution with parameters determined by $m$ and $p$. It looks like the distribution for $\sum\limits_i d(x_i,S)^2$ should then be a sum of iid Beta distributions.
Added in edit. If this is purely about the random subspace I don't think there is much that can be said without more info about $X$. If all the columns are orthogonal then $\sum\limits_i d(x_i,S)^2 = n-m$ while for any $\epsilon>0$ there are configurations where it is an upper bound for the sum for any simplex. That means that $P(\sum\limits_i d(x_i,S)^2 > n-m) = 0$ is the best that can be said without further qualification.
Something might be possible by using the covariance matrix or volumes and Cayley Menger determinants of subsimplices of the simplex made of the points in $X$ and the origin to control the lengths but it won't be straightforward because the geometry of simplices is a complicating factor. If the object of interest were a hypercuboid instead of a simplex, the problem would become purely combinatorical whereas in the simplex the relevant lengths vary between the pairwise distances and each point's opposite face. However even in the rectilinear case the distributions of sums of squared lengths depend strongly on the distribution of individual lengths and can be very complex even if the size of the sample of subsets of X is huge.