We present a new method for modeling discrete constant mean curvature (CMC) surfaces, which arise frequently in nature and are highly demanded in architecture and other engineering applications. Our method is based on a novel use of the CVT (centroidal Voronoi tessellation) optimization framework. We devise a CVTCMC energy function defined as a combination of an extended CVT energy and a volume functional. We show that minimizing the CVT-CMC energy is asymptotically equivalent to minimizing mesh surface area with a fixed volume, thus defining a discrete CMC surface. The CVT term in the energy function ensures high mesh quality throughout the evolution of a CMC surface in an interactive design process for form finding. Our method is capable of modeling CMC surfaces with fixed or free boundaries and is robust with respect to input mesh quality and topology changes. Experiments show that the new method generates discrete CMC surfaces of improved mesh quality over existing methods.