Abstract: To improve the robustness of quadrilateral mesh generation and the quality of boundary elements on complex combined surfaces, a boundary-first meshing algorithm based on Delaunay-advancing layer method is presented. Firstly, with the coarse constrained Delaunay background grid in surface domain, the normal vector of boundary point is calculated iteratively in parametric domain under the guidance of the position deviation from physical domain. Next, anisotropic or isotropic orthogonal meshes are generated around geometry features by advancing-layer strategy. Finally, Coring method is used to accelerate the creation of inner elements, and quadrilateral meshes are obtained by merging triangles. The meshing results of several complex plane regions and combined surface models show that the proposed method can generate meshes with better equiangular distortion and aspect ratio than popular commercial software. It took only 38.68 s to discretize an engine model with 21772 surfaces using OpenMP framework on a 12-thread PC platform.