This surely allows to normalize the wavefunction, and to carry on calculations, but isn't equivalent to saying that we are sure to find the electron inside a particular unit cell, at every time?
This doesn't normalize the wavefunction; $\int_{-\infty}^\infty |\psi_{nk}(x)|^2 \mathrm dx$ is infinite, as always. However, $\int_0^a |u_{nk}(x)|^2\mathrm dx$ is always going to be finite, and sometimes it happens to be convenient to re-scale $u_{nk}$ to make that quantity equal to $1$. Alternatively, sometimes it is helpful to choose the normalization $\int_0^a |u_{nk}(x)|^2 \mathrm dx = a$.
For example, the "orthogonality" of differing eigenstates suggests that $\int \psi^*_{nk}(x) \psi_{mq}(x) \mathrm dx \propto \delta(k-q)\delta_{nm}$. Indeed, if $u_{nk}$ is periodic with period $a$, then it can be expanded as a Fourier series:
$$u_{nk}(x) = \frac{1}{\sqrt{a}}\sum_{r\in \mathbb Z} c^{(nk)}_r e^{-2\pi i rx/a} \iff c^{(nk)}_r = \frac{1}{\sqrt{a}}\int_0^a u_{nk} (x) e^{2\pi irx/a} \mathrm dx$$
Therefore,
$$\psi_{nk}(x) = e^{ikx} u_{nk}(x) = \frac{1}{\sqrt{a}}\sum_{r\in \mathbb Z}c^{(nk)}_r e^{i(k-2\pi r/a)x}$$
As a result,
$$\int \mathrm dx \ \psi^*_{nk}(x) \psi_{mq}(x) = \frac{1}{a}\sum_{r,s\in \mathbb Z} c_r^{(nk)*} c_{s}^{(mq)} \int \mathrm dx\ e^{-i[k-q + 2\pi (r-s)/a]x}$$
$$= \frac{2\pi}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*} c_{s}^{(mq)} \delta\big(k-q + \frac{2\pi}{a} (r-s)\big)$$
However, $k$ and $q$ are both taken from the first Brillouin zone, and so they differ by less than $2\pi/a$. As a result, if $r-s\neq 0$ then the delta function will never "fire" and the result is zero. Therefore, we can replace it with $\delta(k-q) \delta_{rs}$. From there, we obtain
$$\int \mathrm dx \psi^*_{nk}(x) \psi_{mq}(x) = \frac{2\pi}{a} \sum_{r\in \mathbb Z} c^{(nk)*}_r c^{(mk)}_r \delta(k-q)$$
We now notice that
$$\int_0^a u_{nk}^*(x) u_{mk}(x) \mathrm dx = \frac{1}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*}c_s^{(mk)} \int_0^a \mathrm dx \ e^{2\pi i(r-s)x/a}$$
$$= \frac{1}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*}c_s^{(mk)} a \delta_{rs} = \sum_{r\in \mathbb Z}c^{(nk)*}_rc^{(mk)}_r$$
Defining the Bloch function inner product
$$\langle u_{nk}, u_{mk}\rangle \equiv \int_0^a u_{nk}^*(x) u_{mk}(x) \mathrm dx$$
yields
$$\int \psi^*_{nk}(x)\psi_{mk}(x) \mathrm dx = \frac{2\pi}{a} \langle u_{nk},u_{mk}\rangle\delta(k-q)$$
Finally, we note that $u_{nk}$ and $u_{mk}$ are eigenfunctions of the Bloch Hamiltonian $H(k)$ with eigenvalues $E_n(k)$ and $E_{m}(k)$ respectively; since $H(k)$ is self-adjoint, it follows that we may choose $u_{nk}$ and $u_{mk}$ to be orthogonal (of course, if $E_n(k) \neq E_m(k)$ then they are automatically orthogonal). Choosing them to be normalized such that $\langle u_{nk},u_{mk}\rangle = a\delta_{nm}$ yields the final result
$$\int_{-\infty}^\infty \psi^*_{nk}(x) \psi_{mq}(x) \mathrm dx = 2\pi \delta_{nm} \delta(k-q)$$
as expected by analogy with the plane waves of a free particle. On the other hand, choosing the OP's normalization $\langle u_{nk},u_{mk}\rangle = \delta_{nm}$ yields an extra factor of $a$ in the above; however, if you only care about the $u_{nk}$'s (which is common when examining the topology of the band structure, for example) then this may be more convenient.