Reduced Order Models (ROMs) represent a powerful tool to capture the most important features of a flow field by using a small number of degrees of freedom. Recently, the interest in the use of ROMs for design and optimisation purposes has increased in several different field of engineering. These applications introduce important challenges related to the training of the model and the estimation of the error in the predicted field. In particular, the different ROM procedures share the need of a training stage in which several high-fidelity simulations are performed in order to get a set of snapshots and to build a reference database. The sampling in the space of the design parameters is a critical issue since it influences directly the ability of the model to predict a wide range of configurations. In order to optimise the sampling a possible approach is represented by the use of a recursive Voronoi algorithm which explores the design space and focuses the attention on the regions which require further investigation [1]. When a ROM is used to predict a field which corresponds to a set of parameters not included in the database particular care must be taken. Indeed, the non linear nature of the fluid dynamics equations and the high sensitivity of the flow field to some design parameters (geometry of the body, Reynolds and Mach number,…) make the direct use of ROMs particularly difficult for general purpose and industrial applications. Furthermore, the dynamics equations related to ROMs are often characterised by an unstable behaviour which requires the introduction of ad-hoc dissipation terms. In order to avoid these shortcomings, an hybrid approach can be efficiently used to deal with complex flows [2]. In particular, the computational domain can be split in two regions: a region close to the body in which the effects of the body are directly taken into account by CFD and a far-field region in which the flow is described by the ROM. The coupling between the CFD solver and the ROM requires the definition of an overlapping region on which the coefficients of the ROM are computed in order to fit the CFD solution. As a result, the size of the domain on which the expensive CFD simulation has to be performed is strongly reduced. In this work, the previously described hybrid ROM-CFD approach is used to solve a shape optimisation problem. The POD database is trained on a set of full CFD simulations for several values of the design parameters which have been chosen according to the previously described Voronoi sampling technique. The optimisation process is driven by an evolutionary algorithm which calls the hybrid ROM-CFD simulation in order to evaluate the goal function. [1] Bergmann, M., Colin, T., Iollo, A., Lombardi, D., Saut, O., & Telib, H. (2013). Reduced Order Models at work. Modeling, Simulation and Applications, 9. [2] Buffoni, M., Telib, H., & Iollo, A. (2009). Boundary conditions by low-order modelling. In Computational Fluid Dynamics 2006 (pp. 747-752). Springer Berlin Heidelberg.