We consider the evolution of a multi-phase system where the motion of the interfaces is driven by anisotropic curvature and some of the phases are subject to volume constraints. The dynamics of the phase boundaries is modeled by a system of Allen-Cahn type equations for phase field variables resulting from a gradient flow of an appropriate Ginzburg-Landau type energy. Several ideas are presented in order to guarantee the additional volume constraints. Numerical algorithms based on explicit finite difference methods are developed, and simulations are performed in order to study local minima of the system energy. Wulff shapes can be recovered, i.e. energy minimizing forms for anisotropic surface energies enclosing a given volume. Further applications range from foam structures or bubble clusters to tessellation problems in two and three space dimensions.