GARCH model using high frequency price return

Scaling the length of the discrete time step in GARCH models is, from memory, not at all straightforward. For example, you can't just multiply things by the square root of time, like we do for some other, simpler, processes.

For the case where $delta rightarrow 0$, the convergence is derived in Nelson (1990) "ARCH Models as Diffusion Approximations". I think there has also been some work done on scaling from one discrete step length to another, but the reference escapes me now (but this would effectively provide the academic answer to your question if you can find it).

Regarding your specific questions, replacing $u^2$ with $sum_i r_i^2$ is equivalent to replacing the squared return with the realized variance estimator. The result is not a standard GARCH model (I don't think...), although I suspect it would perform quite well empirically, so just because it doesn't fit nicely into the standard framework is not a reason to abandon it per se.

As to whether you would need to use something like Realized Kernels, this would depend entirely on the the length $delta$. At lower frequencies, e.g. 5-minutes, you can just sum squared returns fairly safely without worrying too much about microstructure noise (although of course this all depends on the liquidity of the underlying asset - the more frequently traded the asset, the higher the frequency you can go to without having to worry about microstructure noise). At higher frequencies, e.g. 5-seconds, you would probably need to use an estimator that corrects for microstructure noise, such as Realized Kernels.

One final point of interest, if $c = 0$, and $a + b = 1$, and if we're prepared to make some additional heroic assumptions, like the realized variance converging to the true variance over the interval $(t - k delta, t)$, then I think the model as you've written it converges, as $delta rightarrow 0$, to standard simple exponential smoothing, e.g.: begin{equation} hat{y}_t = alpha y_{t-1} + (1 - alpha) hat{y}_{t-1} end{equation} Note, in this equation $hat{y}_t$ is the variance forecast and $y_{t}$ is the true variance. The "heroic assumptions" are those necessary for the realized variance to converge to the true variance.

ADDITIONAL: You seem to be very focused on whether your model is "correct" or "legitimate". It is important to understand that those words don't really have any meaning here. All models are wrong. And all estimation techniques are inefficient, except under ideal theoretical assumptions that are rarely (never?) satisfied in the real world. So I'm not really sure how to respond to that aspect of the question. Perhaps this:

What matters is whether the model helps you to reach your end goal. For example, one possible end goal might be find a set of volatility forecasts that are useful predictors for implied volatility (i.e. option prices) or the VIX so you can design a profitable trading rule. In this case, perhaps your metric for "better" and "worse" might be something like: "find a volatility forecast model with smaller out-of-sample squared error than the standard GARCH(1,1)". In my experience, in this context, replacing squared returns with realized variances constructed from higher frequency returns spanning the same interval almost always reduces estimation error and will help you towards that goal.

In regards to methods of estimation, we can apply the same logic. Sure, theory might tell you to construct a maximum likelihood estimator with Gaussian innovations, and as far as I can tell what you've written looks good (but I make no promises here - I'm not going to grind through the maths myself to check your working :-). But that might not give you the best out-of-sample forecast. Indeed, in my experience, robust methods of estimation (i.e. those that reduce/eliminate the impact of observations in the tail of the distribution) tend to provide better estimates when working with financial data.