For a hydrogen-like atom in 3 spatial dimensions, the rewriting of the radial part
$$R(r)~=~\frac{u(r)}{r}$$
is not performed to keep the $u(r)$ part regular, as OP suggests, but usually because the 3D radial equation in terms of the $u$ function has the same form as a 1D Schrödinger equation.
Imagine that the radial wave function goes as a power
$$R(r) ~\sim ~ r^{p} \qquad {\rm for} \qquad r~\to~ 0,
\qquad p~\in~\mathbb{R}.$$
On general grounds, one can impose the following list of consistency conditions, listed with the weakest condition first and the strongest condition last.
Normalizability of the wave function
$$\infty~>~\langle\psi|\psi\rangle~=~\int d^3r~|\psi(\vec{r})|^2 ~\propto~ \int_0^{\infty} r^{2}dr~|R(r)|^2 .$$
Integrability at $r=0$ yields that the power $p>-\frac{3}{2}$. In other words, this normalizability condition does not by itself imply that $R(r)$ or $u(r)$ should be regular at $r=0$, which is also the conclusion of many of the other answers.
The expectation value of the potential energy $V$ should be bounded from below,
$$-\infty~<~\langle\psi| V|\psi\rangle~=~\int d^3r~V(r)|\psi(\vec{r})|^2~\propto~-\int_0^{\infty} rdr~|R(r)|^2. $$
Integrability at $r=0$ yields that the power $p>-1$. In other words, $u(r)$ should be regular for $r\to 0$.
The kinetic energy operator (or equivalently, the Laplacian $\Delta$) should behave self-adjointly for two wave functions $\psi_1(\vec{r})$ and $\psi_2(\vec{r})$,
$$\langle\psi_1| \Delta\psi_2\rangle~=~-\langle\vec{\nabla}\psi_1| \cdot\vec{\nabla}\psi_2\rangle,$$ without picking up pathological contributions at $r=0$. A detailed analysis shows that the powers of the radial parts of $\psi_1(\vec{r})$ and $\psi_2(\vec{r})$ should satisfy $p>-\frac{1}{2}$.
In comparison, the actual bound state solutions have non-negative $p=\ell\in \mathbb{N}_0$, and therefore satisfy these three conditions.