A two-level domain decomposition method is introduced for general shape optimization problems constrained by nonlinear partial differential equations. The problem is discretized with a finite element method on unstructured moving meshes and then solved by a parallel two-level one-shot Lagrange-Newton-Krylov-Schwarz algorithm. Due to the pollution effects of the coarse to fine interpolation, direct extensions of the one-level method to two-level do not work. To fix the pollution problem, a pollution removing coarse to fine interpolation scheme is introduced in this paper. As applications, we consider the shape optimization of a cannula problem and an artery bypass problem in 2D. Numerical experiments show that our algorithm performs well on a supercomputer with over one thousand processors for problems with millions of unknowns.