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.
 A: Let $n = 20, N = 10.$ then the population from which we can sample (without replacement)
consists of $N$ copies of $1, 2, \dots, n.$ Then the maximum sample size is $200$ and a sample x of size $200$ will exactly match the population. Using R:
pop = rep(1:20, N)
summary(pop);  length(pop);  sd(pop)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    5.75   10.50   10.50   15.25   20.00 
[1] 200       # population size
[1] 5.780751  # population SD (with denom 199)

x = sample(pop, 200)
summary(x);  length(x);  sd(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    5.75   10.50   10.50   15.25   20.00 
[1] 200       # sample size
[1] 5.780751  # sample SD

That part is trivial. So your question must have something to
do with samples of size less than 200.  However, you give no
clue what that question is.
Let's look at a sample of size $50.$
set.seed(2021)
y = sample(pop, 50)
summary(y);  length(y);  sd(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    7.00   11.50   10.92   15.00   20.00 
[1] 50
[1] 5.382663
table(y)
y
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
  2  2  2  1  1  3  3  5  2  2  2  4  3  4  3  3  2  1  2  3 

A: 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))

