How $Var[e^{\frac{-1}{X+a}}]$ varies with $n$ where $X \sim Bin(n,p)$? I have a binomial random variable $X \sim Bin(n,p)$. I am interested in the variance of a function $f(X)$ given by :
$f(X)=e^{\frac{-1}{X+a}}$. Here $a>0$.
Specifically, I would like to know how $Var[f(X)]$ varies with $n$?
I am unable to come with a way to proceed. Any help is appreciated.

To supplement the answer given by @Ben, consider the following plot

which is produced with the following R code.
# Moment function
m <- function(k, n, p, a){
  res <- rep(NA, length(n))
  for(i in seq_along(n)){
    xx <- 0:n[i]
    ff <- exp(-1/(xx+a))
    xf <- k*log(ff) + dbinom(xx, n[i], p, log=TRUE) #For computational stability
    res[i] <- exp(matrixStats::logSumExp(xf))
  }
  return(res)
}

# Variance function
v <- function(n, p, a){
  m(2, n, p, a) - m(1, n, p, a)^2
}

# Make plot
par(mfrow=c(2,2))
nn <- 0:30
aa <- c(0.01, 0.1, 1, 10)
pp <- c(0.25, 0.5, 0.75)

for(i in 1:4){
  plot(NULL, xlim=c(0, 30), ylim=c(0, max(v(nn, pp[1], aa[i]))),
             xlab="sample size", ylab="Var(f(X))", 
             cex.lab=1.4, main=paste0("a = ", aa[i]), cex.main=1.4)
  if(i == 1){
    legend('topright', c("p=0.25", "p=0.50", "p=0.75"), 
           lty=1:3, lwd=2, bty='n', cex=1.4)
  }
  for(j in 1:3){
    lines(nn, v(nn, pp[j], aa[i]), col=i, lty=j, lwd=2)
  }
}

