In this work we consider shape optimization of systems, which are governed by external Bernoulli free boundary problems. A pseudo-solid approach for solving discrete free boundary problems is introduced. The solution strategy readily allows us to obtain geometrical sensitivities of the system, which can then be used to solve e.g. inverse design problems. Numerical examples show that the location of the free boundary can, to some extent, be controlled by changing the shape of the other component of the boundary.