Sum of Discrete Uniforms, but each value can be picked no more than N times? Suppose there are i.i.d. variables $X_{1,..n}$ with discrete uniform distribution with the support $[1, n]$. What would be the distribution of such a sum if we introduce the condition that each value cannot appear more than $N$ times in the sum?
E.g., if $N = 1$, we have sampling without replacement, if $N = 2$, any number can only appear twice in the sum. $N = \infty$ is equivalent to sampling with replacement.

I get a reasonable normal distribution approximation with the following code. It is a bit heuristic and I have not yet a strong argument for this method. Below is the reasoning behind it:
I compute the variance by considering the sum as a sort of random walk process of the partial sum $S_k = \sum_{j=1}^{k} X_k$. With this random walk there is an attractive force towards $(n+1)/2$ because the mean of the variables that are not yet chosen will be higher/lower depending on the sum $S_k$. So if $S_k/k$ is larger/lower than $(n+1)/2$ then the expectation of $X_{k+1}$ will be the opposite and with a relative magnitude $\frac{k}{Nn-k}$.
This makes that the variance of $S_k$ is not simply $k \sigma^2$ with $\sigma^2 = \frac{n^2-1}{12}$ the variance of the discrete uniform distributed numbers from which the sample is taken. Instead, we approximate it with:
$$\text{Var}(S_{k+1}) = \frac{1}{(Nn-k)^2} \text{Var}(S_{k}) + \sigma^2$$

### function to simulate sampling with limited numbers
### and return the sum of the sample 
smp = function(n,N) {
  x = rep(1:n, times=N)
  y = sample(x,n,replace = 0)
  return(sum(y))
}


n = 30 ### size of population to sample 
nr = 10^6 ### number of simulation

### compute ranges for plotting
summu = n*(n+1)/2
sumvar = mean((c(1:n)-(n+1)/2)^2)*n

x_lims = summu + c(-3,3)*sumvar^0.5
y_lim = dnorm(summu, summu, sumvar^0.5)*2

plot(-100,-100, 
     xlim = x_lims, ylim = c(0,y_lim), 
     xlab = "sum", ylab = "frequency")


### compute for different values of N
for (i in 1:4) {
  N = c(2,3,4,n)[i]
  col = c(2,3,4,1)[i]
  
  #### repeat sampling and plotting
  z = replicate(nr,smp(n,N))
  h = hist(z, breaks = seq(-0.5,max(z+0.5),1), plot=0)
  points(h$mids, h$counts/nr, col=col, pch=21, bg=col, cex=0.5)

  ### normal approximation estimate
  mu = n*(n+1)/2
  ### compute variance 
  vr1 = mean((c(1:n)-(n+1)/2)^2)
  vr = vr1
  for (j in 1:(n-1)) {
    vr = vr * (1-1/(N*n-j))^2 + vr1
  }

  ### add line for normal approximation
  x = seq(min(x_lims),max(x_lims),0.1)
  lines(x,dnorm(x,mu,vr^0.5), col = col)
}

title("simulation and estimate for n = 30")   
legend(min(x_lims),y_lim, c("N = 2", "N = 3", "N = 4", "N = n"),
       col = c(2,3,4,1), pch = rep(21,4), lty = rep(1,4), pt.bg = c(2,3,4,1))

