  For $\Omega$ varying among open bounded sets in ${\mathbb R} ^n$, we consider
shape functionals $J (\Omega)$ defined as the infimum over a Sobolev space of
an integral energy of the kind $\int _\Omega[ f (\nabla u) + g (u) ]$, under
Dirichlet or Neumann conditions on $\partial \Omega$. Under fairly weak
assumptions on the integrands $f$ and $g$, we prove that, when a given domain
$\Omega$ is deformed into a one-parameter family of domains $\Omega
_\varepsilon$ through an initial velocity field $V\in W ^ {1, \infty} ({\mathbb
R} ^n, {\mathbb R} ^n)$, the corresponding shape derivative of $J$ at $\Omega$
in the direction of $V$ exists. Under some further regularity assumptions, we
show that the shape derivative can be represented as a boundary integral
depending linearly on the normal component of $V$ on $\partial \Omega$. Our
approach to obtain the shape derivative is new, and it is based on the joint
use of Convex Analysis and Gamma-convergence techniques. It allows to deduce,
as a companion result, optimality conditions in the form of conservation laws.
