Sampling of a continuous function: Kronecker's or Dirac's delta?
Modeling of the sampling process via multiplication of a continuous-time signal by a train of Dirac impulses is the most common interpretation in my experience. If you dig into it deeply enough, you'll find some disagreement about the mathematical precision of this approach*, but I wouldn't worry about it; it's just a convenient model for the process. There are no impulse generators inside your cell phone's ADC generating periodic bolts of lightning that multiply their analog inputs.

As you noted, you can't calculate the continuous-time Fourier transform of the Kronecker delta function, as its domain isn't continuous (it is limited to the integers). The Dirac delta function, in contrast, has a simple Fourier transform, and the effect of multiplying a signal by a train of Dirac impulses is easy to show due to its sifting property.

*: As an example, if you're going to be mathematically precise, you would say that the Dirac delta isn't a function at all, but a distribution instead. But at an engineering level, these issues are really just semantics.

Edit: I'll address the comment below. You gave your mental model of the sampling process as:

$$
f_s(t) = \sum_{k=1}^N \int_{t_k-\epsilon_k}^{t_k+\epsilon_k} f(t) \delta(t-t_k) dt.
$$

The problem with this interpretation is that the typical ideal sampling model does not have that integration built in. Instead, it's a pure multiplication of the input signal by a Dirac impulse train. If you look more closely at the equation you showed for $f_s(t)$, you'll see that the right hand side actually has no independent variable; $t$ is the dummy variable of integration. For any $\epsilon_k &gt; 0$ in the above, according to the Dirac impulse's sifting property, you would get: 

$$
f_s(t)=\sum_{k=1}^N f(t_k),
$$ 

which isn't correct. Instead, the model for the sampled signal is:

$$
f_s(t)=\sum_{k=-\infty}^{\infty} f(t) \delta(t-kT)
$$

Which is very similar to the above, except generalizing for an infinitely long impulse train along the time axis and assuming that the data is uniformly sampled at time instants $t_k = kT$. The Fourier transform of the resulting signal is:

$$
\begin{align}
F_s(\omega) &amp;= \int_{-\infty}^{\infty}f_s(t) e^{-j \omega t} dt \\
&amp;= \int_{-\infty}^{\infty} \sum_{k=-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&amp;= \sum_{k=-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&amp;= \sum_{k=-\infty}^{\infty} f(kT) e^{-j \omega kT}
\end{align}
$$

If we define the discrete, sampled version of the signal $f(t)$ to be $x[n] = f(nT)$, then you are left with:

$$
F_s(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n}
$$

which is exactly the definition of the discrete-time Fourier transform.