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.$$
In general when $f: \Bbb{R}^n \to \Bbb{R}$ is the euclidean norm, defined by
\begin{equation}
f(x) = |x| = \sqrt{\sum_{i=1}^n x_i^2}
\end{equation}
then away from the origin, $f$ is $C^{\infty}$ and its derivative is given by
\begin{equation}
f'(x)(\eta) = \left\langle \dfrac{x}{|x|}, \eta \right\rangle
\end{equation}
(for all $x \neq 0$, and for all $\eta \in \Bbb{R}^n$).
Now, for convenince, define $g: \Bbb{R} \times \Bbb{R}^n \to \Bbb{R}^n$,
\begin{equation}
g(\varepsilon,x) = \nabla u(x) + \varepsilon \nabla v(x)
\end{equation}
The energy functional's derivative you wish to calculate is given by
\begin{align}
J'(u)(v) &= \dfrac{d}{d \varepsilon} \bigg|_{\varepsilon = 0} J(u+ \varepsilon v) \\
&= \dfrac{d}{d \varepsilon} \bigg|_{\varepsilon = 0} \int_{ x \in\Omega} \dfrac{\left[ (f \circ g)(\varepsilon, x) \right]^p}{p} \\
&= \int_{ x \in\Omega} \dfrac{\partial}{\partial \varepsilon} \bigg|_{\varepsilon = 0} \dfrac{\left[ (f \circ g)(\varepsilon, x) \right]^p}{p} ,
\end{align}
This last equality is by Leibniz's integral rule for differentiation under the integral sign. Now, inside, we just have to make use of the chain rule inside. Doing so yields
\begin{align}
J'(u)(v) &= \int_{ x \in\Omega} \dfrac{p \cdot \left[ f(g(0,x)) \right]^{p-1}}{p} \cdot f'[g(0,x)] \left( \dfrac{\partial g}{\partial \varepsilon}(0,x) \right) \\
&= \int_{ x \in\Omega} |\nabla u(x)|^{p-1} \cdot \left \langle \dfrac{\nabla u(x)}{|\nabla u(x)|} , \nabla v(x) \right \rangle \\
&= \int_{x \in \Omega} |\nabla u(x)|^{p-2} \cdot \left \langle \nabla u(x), \nabla v(x) \right \rangle
\end{align}
Modulo some notation, this is precisely what you wanted to prove.
By the way, my entire answer assumes that all the functions are nice and differentiable, and in particular that if $ p< 2$, then $\nabla u(x)$ doesn't vanish anywhere, so that the norm is differentiable. If this isn't the case, then the argument will necessarily have to be more involved... I'm not too sure how one might proceed in this case though.
If $p\geq 2$, then we can write $|\nabla u|^p = \langle \nabla u, \nabla u\rangle^{p/2}$. In this case, the functions will be differentiable even regardless of whether or not $\nabla u$ vanishes; this is because the function $t \mapsto t^{p/2}$ is everywhere differentiable, and the inner product $\langle \cdot , \cdot \rangle$ is also everywhere differentiable, hence their composition will be as well. Of course, in this case, we will have to slightly modify the above argument, but the final answer is the same
Best Answer
Depending on the context, there can be a few different answers.
In cases outside of Poisson's equation, variational problems often emerge initially as finding the extrema of a functional which sometimes can be intepreted as "energy." For Poisson's equation this is not the case historically (the PDE was derived first, and the variational formulation later). But depending on what you study, it can be the case that more often than not, the PDE comes out of the optimization problem, and not vice versa.
Knowing that a solution (should one exist) must be a minimizer of some energy functional can sometimes be useful for deducing that one must exist and/or be unique. For instance, given the form of the energy functional, one may be able to conclude that it is coercive and weakly lower semi-continuous. This is enough to conclude that a minimizer exists; moreover, such a minimizer can, in principle, be constructed by taking a weak infimizing sequence and upgrading the weak convergence to strong convergence (after possibly modding out the symmetries of the underlying problem: this is referred to sometimes as compensated/concentration compactness). In the case of Poisson's equation, the existence and uniqueness problems are relatively tame and may be possible to resolve with other methods, but this sort of thing becomes valuable for more difficult nonlinear PDEs. Therefore it's a good idea to learn how to solve Poisson's equation using energy methods, since this easy case generalizes to a large class of contexts.