The solution to a free boundary problem of Bernoulli type, also known as Alt-Caferelli problem, is studied via shape optimization techniques. In particular, a novel energy- gap cost functional approach with a state constraint consisting of a Robin condition is proposed as a shape optimization reformulation of the problem. Accordingly, the shape derivative of the cost is explicitly determined, and using the gradient information, a Lagrangian-like method is used to formulate an efficient boundary variation algorithm to numerically solve the minimization problem. The second order shape derivative of the cost is also computed, and it is shown how the Hessian information can be utilized to formulate a second-order optimization procedure based on a Newton-type method for resolving shape optimization problems. The feasibility of the newly proposed method and its comparison with the classical energy-gap type cost functional approach is then presented through various numerical results.