Continuous approximation to binomial distribution Consider an integer variable $k$ that follows a binomial distribution,
$$\binom{N}{k}p^{k}\left(1-p\right)^{N-k}$$
with total draws $N$ and probability of success $p$. I am interested in the fraction of successes, $f = k/N$. The distribution of $f$ can be approximated as a normal distribution with mean $p$ and variance $\sqrt Np(1-p)$. This approximation is good if both $Np$ and $N(1-p)$ are sufficiently large, but it has the problem that there is a finite (though small) probability that $f<0$ or $f>1$, because the normal distribution has infinite support.
Is there an approximation to the distribution of $f$ that has support $f\in[0,1]$? We can assume that $Np$ and $N(1-p)$ are both sufficiently large.
 A: An obvious candidate would be the beta distribution, since this is the conjugate to the binomial distribution and it is on the appropriate support.  To allow for continuity correction and avoid poor approximation at the edges, it is desirable to approximate each discrete $x = 0, 1, ...., N$ by an equal-sized continuous interval.  This means that we approximate the binomial mass function by an integral of the beta density over one of $N+1$ equal-sized intervals on its support:
$$\text{Bin}(x|N,p) \approx \int \limits_{B(x)}^{B(x+1)} \text{Beta}(\theta|\alpha, \beta)d\theta
\quad \quad \quad
B(x) \equiv \frac{x}{N+1}.$$
The approriate parameters $\alpha$ and $\beta$ can be found using the method-of-moments (MOM), which requires us to solve the following two moment equations:
$$\frac{\alpha}{\alpha+\beta} = p \quad \quad \quad \frac{\alpha \beta}{(\alpha+\beta)(\alpha+\beta+1)} = \frac{p(1-p)}{N+1}.$$
Solving for the required parameters yields the values $\alpha = pN$ and $\beta = (1-p)N$ so our approximation to the binomial is:
$$\text{Bin}(x|N,p) \approx \frac{\Gamma(N)}{\Gamma(pN)\Gamma((1-p)N)} \int \limits_{B(x)}^{B(x+1)} \theta^{pN-1} \theta^{(1-p)N-1} d\theta.$$
A: Suppose you see $X$ Successes in $n$ trials. Then the estimate of the binomial
success probability is $\hat p = X/n.$ You cannot know the exact distribution of $X$ or of $\hat p$ because you don't know $p.$ As I showed in my answer below,
there are various methods to get a 'confidence interval' for $p,$ using the fact that $X \sim \mathsf{Norm}(\mu = np, \sigma = \sqrt{np(1-p)}).$  Two Answers have also shown that $X$ is approximately distributed according to a beta distribution. But you say:

No. I am interested in the distribution of the number of successes. In the limit of large number of trials, the fraction of successes should follow an approximately continuous distribution.

This response is puzzling. And you have not responded to @whuber's Comments.
It seems this matter will remain unresolved until you explain what you are doing
and exactly why answers and comments to date are not what you are looking for.

It seems as if you might be trying to get a confidence interval (CI) for $p$ based on seeing $X$ Successes in $n$ trials. 
Wald interval: Then a point estimate for $p$ is $\hat p = X/n.$ A traditional 'Wald' 95% confidence interval is of the form 
$$\hat p \pm 1.96 \sqrt{\frac{\hat p(1-\hat p)}{n}}.$$
It uses two approximations: (a) that $\hat p$ is approximately normally
distributed and (b) that the standard error of $\hat p.$
$$SD(\hat p) = \sqrt{\frac{p(1-p)}{n}} \approx \sqrt{\frac{\hat p(1-\hat p)}{n}}.$$
Because these two approximations are best for large $n$ and $p$ relatively near
$1/2,$ this style of CI does not provide the nominal 95% coverage probability
in many practical situations. You mentioned one difficulty: the CI has length $0$ if $X = 0$ or $X = n.$
Example: If $X = 30$ and $n = 50,$ a 95% Wald CI is $(0.464, 0.736).$
n = 50;  p.hat = 30/50; pm = c(-1,1); wald.ci = p.hat + pm*1.96*sqrt(p.hat*(1-p.hat)/n)
wald.ci
[1] 0.4642072 0.7357928

Agresti-Coull CI: The improved Agresti-Coull style of CI is of the form
$$\tilde p \pm 1.96 \sqrt{\frac{\tilde p(1-\tilde p)}{\tilde n}},$$
where $\tilde n =n+4$ and $\tilde p = (X+2)/\tilde n).$
Example: For the data above, this CI is $(0.456, 0.729).$
n.tilde=54; p.hat=32/54; pm=c(-1,1); agresti.ci=p.hat+pm*1.96*sqrt(p.hat*(1-p.hat)/n)
agresti.ci
[1] 0.4563968 0.7287884

The procedure @Ben (+1) describes can be adapted to give a 95% Bayesian probability interval ('credible interval'), using $\mathsf{Beta}(1,1) \equiv \mathsf{Unif}(0,1)$ as a non-informative prior distribution and a binomial
likelihood function proportional to $p^X(1-p)^{n-X}.$ 
Then the posterior
distribution is $\mathsf{Beta}(X+1, n-X+1)$ and the interval estimate
has quantiles .025 and .975 of that distribution as its endpoints:
$(0.461, 0.724)$ in our example. 
qbeta(c(.025, .975), 31, 21)
[1] 0.461141 0.724157

Note: If you are interested in interval estimation for binomial data,
you may want to look at a similar Q & A and the link at the end.
A: Express factorial functions on the Binomial coefficient as Gamma functions: $\Gamma(x + 1)$.
$${K\choose k} \equiv \frac{K!}{k!(K-k)!} = \frac{\Gamma(K + 1)}{\Gamma(k + 1)\Gamma(K - k + 1)},\ \forall\ \ \mathbb{N} \in \{0, 1, ..., K\},$$
where $K$ is the total number of trials, and $k$ the number of successful trials.
The Binomial distribution can as such be expressed as:
$$\frac{\Gamma(K + 1)}{\Gamma(k + 1)\Gamma(K - k + 1)}p^k(1-p)^{K-k},$$
where $p$ is the probability of success per trial.
This extends the Binomial distribution to the real-number plane. However, the extension allows for non-zero frequencies of occurences less than $0$ and greater than $K$. Hence, it needs to be truncated to the interval $[0, K]$. This range-restriction makes necessary a normalizing constant ensuring that the density on this interval equal $1$. The constant can be obtained by:
$$C = \int_0^K \frac{\Gamma(K + 1)}{\Gamma(k + 1)\Gamma(K - k + 1)}p^k(1-p)^{K-k}\ dk$$
A continuous analogue to the Binomial distribution can as such be expressed as:
$$ P(X; K, k, p) = \begin{cases} \frac{1}{C}\left[\frac{\Gamma(K+1)}{\Gamma(k + 1)\Gamma(K - k + 1)}p^k(1-p)^{K-k}\right],& \text{if}\ 0 \leq k \leq K,\\ 0,& \text{otherwise.} \end{cases} $$
