Consider making the substitution $k = p/\hbar$ in the first expression, while simultaneously defining
$$
\sqrt{\hbar} \,a(p) = \phi(k)
$$
then the first integral will become
$$
\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \sqrt{\hbar}\, a(p) e^{i\left(\frac{p}{\hbar} x - \frac{\hbar}{2m}\frac{p^2}{\hbar^2}t\right)} \frac{dp}{\hbar}
$$
which is exactly the second integral. It's just a change of notation!
So we measure some definite value of momentum or energy which is an eigenvalue of the momentum or Hamiltonian (since the operators commute for a free particle).
The result of single measurement can be single value, but in case of quantities that have continuous domain we cannot say this is with certainty the actual value of the quantity. With any such measurement we always have uncertainty of the outcome that is greater than zero. This is unavoidable in practice, we do not have means to measure continuous variables with infinite accuracy.
Then we would in principle collapse the wave function to some stationary state $\Psi_k$ but in this case we know that this is not possible (physically realizable).
It is not important here whether such process is physically realizable; this depends on interpretation of the theory. There are interpretations that do not consider collapse as a result of measurement to be physical process at all, irrespective of whether the result is normalizable.
What is important here is that there is no normalizable function that would be eigenfunction of position operator (and there is no one that would be eigenfunction of the momentum operator). Therefore we cannot base our understanding of the theory on such fictive functions. Particle with definite position or momentum with zero uncertainty cannot be represented by normalized $\psi$ function.
So do we measure a particular value of the momentum with some measurement uncertainty, with the uncertainty giving us the spread of values of the observable
Yes, all measurements of position or momentum of particles have finite uncertainty, so the probability that the measured value equals the actual value that was sought is 0. When we look at particle tracks from bubble chamber, the track is thin but finite width, limiting the uncertainty of the particle coordinate to small but finite distance. In practice, I think microns at best.
Or do we not measure a particular value ever but rather a range of values for a given measurement?
When one particle is measured, one value plus uncertainty is usually recorded. If many particles are measured, then many values and uncertainties are recorded. In any case, no result is ever absolutely accurate, there is always some uncertainty.
Best Answer
We usually define $p=\hbar k$. In this post, I'll take $\hbar =1$ (see natural units) to simplify the notation. This means that $k=p$. This is not necessary, but IMHO this makes the arguments more transparent (besides the fact that natural units are indeed more natural once you get used to them). I usually recommend every learning physicist to try to use natural units whenever they can. If you don't like to, you can exchange $k\leftrightarrow p/\hbar$ in my answer.
In general, the expectation value of the momentum is given by $$ \langle k\rangle= \int_{-\infty}^{+\infty} \mathrm dk\ k\ |\phi(k)|^2 \tag{1} $$ which means that $|\phi(k)|^2$ has to decrease faster than $1/k^3$ as $k\to\infty$ for $\langle k\rangle$ to be well defined. In fact, we also want $\langle k^2\rangle$ to be well defined, which means we usually want $\phi(k)$ to be $\mathcal O(k^{-2})$ as $k\to\infty$.
If $\phi(k)$ has a large spread (it doesn't decrease fast enough), then the momentum is not well-defined, because the integrals that define momentum are divergent.
Well, $\phi(k)$ is the wave-function in momentum space, so all the information about the momentum of the system is contained in $\phi(k)$.
Note that $\phi(k)$ is the Fourier Transform of $\Psi$, so that this statement is actually equivalent to the uncertainty principle of the Fourier Transform.
UPDATE
Here we prove that both expressions $(1)$ and $(2)$ are equivalent: $$ \langle k\rangle=-i\int_{-\infty}^{+\infty} \mathrm dx \ \Psi^*\partial_x\Psi \tag{2} $$
The proof is non-trivial, but in the end is just the differentiation property of the Fourier Transform. To prove the equivalence, we'll make use of the known fact $$ \int_{-\infty}^{+\infty} \mathrm dx\ \mathrm e^{iqx}=2\pi\delta(q) \tag{3} $$ for any $q\in\mathbb R$. Here, $\delta(q)$ is the Dirac delta function.
First, note that $$ \Psi=\frac{1}{\sqrt{2\pi}}\int \mathrm dk\ \phi(k)\;\mathrm e^{i(kx-\omega t)} \tag{4} $$ so that $$ \partial_x\Psi=\frac{i}{\sqrt{2\pi}}\int \mathrm dk\ k\;\phi(k)\;\mathrm e^{i(kx-\omega t)} \tag{5} $$
If we make the change of variable $k\to k_1$ in $(4)$ and $k\to k_2$ in $(5)$, and plug these into $(2)$ we get $$ \begin{align} \langle k\rangle&=\frac{1}{2\pi}\int \mathrm dx \int\mathrm dk_1\int\mathrm dk_2\ \overbrace{\phi^*(k_1)\;\mathrm e^{-i(k_1x-\omega t)}}^{\Psi^*}\ \overbrace{k_2\;\phi(k_2)\;\mathrm e^{i(k_2x-\omega t)}}^{\partial_x\Psi}=\\ &=\frac{1}{2\pi}\int\mathrm dk_1\int\mathrm dk_2\ k_2\ \phi^*(k_1)\;\phi(k_2)\int \mathrm dx\ \mathrm e^{i(k_2-k_1)x} \end{align} $$
Next, use $(3)$ to simplify the $x$ integral, where $q\equiv k_2-k_1$: $$ \begin{align} \langle k\rangle&=\int\mathrm dk_1\int\mathrm dk_2\ k_2\ \phi^*(k_1)\;\phi(k_2)\;\delta(k_2-k_1)=\\ &=\int\mathrm dk_1 k_1\ |\phi(k_1)|^2 \end{align} $$ which is just $(1)$ upon the change of variables $k_1\to k$. As you can see, both expressions are equivalent, which means that you can use whichever you like the most.