All your notation seems wrong. So I am fixing it.
Let us consider the space $\mathrm{X} = \mathscr{B}^1_{\mathbf{R}}(\mathrm{U} \times \mathrm{I})$ of differentiable real valued bounded function with continuity defined on the product of the open set $\mathrm{U}$ of $\mathbf{R}^d$ and the compact interval $\mathrm{I} = [a, b]$ of $\mathbf{R}.$ Endow $\mathrm{X}$ with the structure of complete normed space by setting $\|f\|_\mathrm{X} = \|f\| + \left\| \mathbf{D} f \right\|$ (the sup-norm and the operator norm). (The fact that $\mathrm{X}$ is complete follows from the mean value theorem.)
Define $\varphi:\mathrm{X} \to \mathrm{Y} = \mathscr{C}^b_\mathbf{R}(\mathrm{U}),$ the right hand side is the space of continuous and bounded functions with the sup norm, by
$$\varphi(f) = u_f, \quad u_f(x) = (1 + f(x))^p \int\limits_a^b \mathbf{D}_2f(x, t)\ dt.$$
Now,
$$\begin{align*}
u_{f + h}(x) &= (1 + f(x)+ h(x))^p \int\limits_a^b (\mathbf{D}_2f(x, t) + \mathbf{D}_2h(x, t))\ dt \\
&= \big[(1 + f(x))^p + p(1 + f(x))^{p-1} h(x) + o(1+f(x); h(x))\big]\\
&\times \int\limits_a^b (\mathbf{D}_2f(x, t) + \mathbf{D}_2h(x, t))\ dt \\
&= u_f(x)+p(1+f(x))^{p-1} h(x)\int\limits_a^b \mathbf{D}_2f(x, t)\ dt + (1 + f(x))^p \int\limits_a^b \mathbf{D}_2h(x, t)\ dt \\
&+ o(1+f(x); h(x))\int\limits_a^b (\mathbf{D}_2f(x, t) + \mathbf{D}_2h(x, t))\ dt, \\\\
\end{align*}$$
where $o(a; s)$ is a sum of powers of $s^\alpha,$ where $\alpha$ runs from 2 until p, and the coefficients are powers of $a$ as well. Because $f$ is bounded, there exists a constant $c = c(f)$ such that $$\|o(1 + f; h)\| \leq c(\|h\|^2+\ldots+\|h\|^p) = o(\|h\|).$$
It is easy to check that the function $L_f:h \mapsto L_f(h)$ given according to the rule
$$x\mapsto p(1+f(x))^{p-1} h(x)\int\limits_a^b \mathbf{D}_2f(x, t)\ dt + (1 + f(x))^p \int\limits_a^b \mathbf{D}_2h(x, t)\ dt$$
belongs to $\mathrm{Y}$ (that is, $L_f:\mathrm{X} \to \mathrm{Y}$) and it is linear. The desired derivative of $\varphi$ is therefore $\mathbf{D}\varphi(f) = L_f.$ Q.E.D.
Best Answer
The Frechet derivative $DE$, if it exists, is unique and satisfies
$$E(u+h)=E(u)+DE(h)+r(h),\ $$ where $r(h)$ is $o(h).$ So, if we can find a candidate that satisfies the equation, we are done.
Claim (admittedly with the foreknowledge that the claim is true):
$$DE(h)=\int_{\Omega}\langle \nabla u,\nabla h\rangle$$
The proof is a calculation:
$$E(u+h)-E(u)=\frac{1}{2}\left (\int_{\Omega} | \nabla (u+h)|^2-\int_{\Omega} | \nabla (u)|^2\right )=\frac{1}{2}\left (\int_{\Omega} \langle\nabla (u+h),\nabla (u+h)\rangle-\int_{\Omega} | \nabla (u)|^2\right )=\int_{\Omega}\langle \nabla u,\nabla h\rangle+\frac{1}{2}\int_{\Omega}\langle \nabla h,\nabla h\rangle,$$
from which we see that, setting $r(h)=\frac{1}{2}\int_{\Omega}\langle \nabla h,\nabla h\rangle$ and noting that is is $o(h)$, we have
$$DE(h)=\int_{\Omega}\langle \nabla u,\nabla h\rangle.$$