We consider the variational approach to prove the existence of solutions of second order stationary Mean Field Games on a bounded domain $\Omega\subseteq \mathbb{R}^{d}$, with Neumann boundary conditions, and with and without density constraints. We consider Hamiltonians which growth as $|\cdot|^{q'}$, where $q'=q/(q-1)$ and $q>d$. Despite this restriction, our approach allows us to prove the existence of solutions in the case of rather general coupling terms. When density constraints are taken into account, our results improve those in \cite{MesSil}. Furthermore, our approach can be used to obtain solutions of systems with multiple populations.