This paper proposes a novel approach for multiobjective optimization when the criteria of interest rely on a functional output from an expensive-to-evaluate numerical simulator. More specifically, the proposed method is developed in the frame of an automotive application. The aim of this application is to design the shape of an intake port in order to maximize the mass flow (denoted by Q) and the tumble (denoted by T ), which both depend on a 3D velocity field obtained by numerical flow simulation. Since the considered flow simulator is time-consuming, using regular multi-objective genetic algorithms (MOGA) directly on integral quantities depending on the simulator output is prohibitive. Three different Reduced Order Models (ROMs) are presented. The first one consists in directly Kriging the integral quantities Q and T on the basis of the outputs computed at an initial design of experiments, and basing the optimization search on the sequentially obtained couples of response surfaces. The other methods explored in the present work consist in building a parametrized representation of the whole velocity field by different variants of the Proper Orthogonal Decomposition (POD). Instead of directly Kriging Q and T at un-sampled locations, the proposed technique is hence to proceed in two steps: first approximate the functional outcome by Kriging the POD coefficients, and then compute the integral quantities Q and T associated with the approximate 3D field. However, such an approach induces new difficulties since the truncated POD does not preserve the global (integrated) quantities, and that surrogate-based MOGA with this kind of POD are therefore likely to fail locating the (Q, T )-Pareto front accurately. This is what motivates to propose an original constrained POD method (called CPOD) meant to overcome the bias created by the truncation made in regular POD. More precisely, this means modifying the way of calculating the POD coefficients by imposing the integral quantities Q and T based on the truncated POD to match with the actual Q and T values obtained by flow simulation at the design of experiments. A detailed comparison of the Pareto sets obtained from the three ROMs demonstrates the interest of the CPOD approach.