$E = \int_{\Omega}\left \| \nabla u(x)\right \|^2 dx$
So, Dirichlet's energy measures the integral of the squared norm of the gradient. Why squared norm? What would we get if we use just a norm? It's still going to be non-negative.
If I calculated a $u(x)$ that minimizes $E = \int_{\Omega}\left \| \nabla u(x)\right \| dx$ (not squared) would it be worse than Dirichlet? My goal is exactly what I specify: the gradient at any point should have minimum length.
Is the reason for the squared norm minimisation equivalent to the role of the square in the least squares fit? (here)
I think the square comes from physical modelling the bending energy.
But then, the square of a Hilbert space norm has several nice properties: it is uniformly convex and twice continuously differentiable with constant second derivative.
The Dirichlet problem $$ \min \int_\Omega \frac12|\nabla u|^2 - uf \ dx $$ is of course not equivalent to $$ \min \int_\Omega |\nabla u| - uf \ dx. $$ The latter problem is much harder to solve and analyze.
For instance, it is much more difficult to prove that the second problem has a solution with integrable gradient.
The minimizer of the Dirichlet energy satisfies $$ -\Delta u = f, $$ whereas the minimizer of the second problem satiesfies $$ -\nabla\cdot\left( \frac{\nabla u}{\|u\|} \right) =f. $$ To have a finite-dimensional analogue, consider the minimization of $$ \sum_{i=1}^n \frac12|a_i|^2 - a^Tc $$ and of $$ \sum_{i=1}^n |b_i| - b^Tc. $$ The first one has solution $a_i=c_i$, the second has a solution only if $|c_i|\le 1$ for all $i$ (which is then $b_i=c_i$).