In this paper, a Pareto front (PF)-based sampling algorithm, PF-Voronoi sampling method, is proposed to solve the problem of computationally intensive multi-objectives of medium size. The Voronoi diagram is introduced to classify the region containing PF prediction points into Pareto front cells (PFCs). Valid PFCs are screened according to the maximum crowding criterion (MCC), maximum LOO error criterion (MLEC), and maximum mean MSE criterion (MMMSEC). Sampling points are selected among the valid PFCs based on the Euclidean distance. The PF-Voronoi sampling method is applied to the coupled Kringing and NASG-II models and its validity is verified on the ZDT mathematical cases. The results show that the MCC criterion helps to improve the distribution diversity of PF. The MLEC criterion and the MMMSEC criterion reduce the number of training samples by 38.9% and 21.7%, respectively. The computational cost of the algorithm is reduced by more than 44.2%, compared to EHVIMOPSO and SAO-MOEA algorithms. The algorithm can be applied to multidisciplinary, multi-objective, and computationally intensive complex systems.