Throughout the last 25–30 years, optimal shape design has become more and more important in engineering applications. Many problems that arise in structural mechanics, fluid dynamics and electromagnetics lead to the minimization of functionals defined over a class of admissible domains, governed by the solutions of boundary value problems or initial boundary value problems. Most of the related optimization methods are based on shape gradients which has two major drawbacks. Firstly, the shape gradients need to be analytically computed which is a nontrivial task. Secondly, the shape gradients need also to be implemented. We thus intend to develop gradient-free optimization methods for the numerical solution of shape optimization problems where the focus is on free boundary problems.