We consider shape functionals obtained as minima on Sobolev spaces of classical integrals having smooth and convex densities, under mixed Dirichlet-Neumann boundary conditions. We propose a new approach for the computation of the second order shape derivative of such functionals, yielding a general existence and representation theorem. In particular, we consider the p-torsional rigidity functional for p >= 2.