Phononic crystals can be designed to show band gaps—ranges of frequencies whose propagation is strongly attenuated in the material. In essence, their working principle is based on destructive interference of waves reflecting from the periodic arrangement of material interfaces (i.e., Bragg scattering). Consequently, capturing accurately the behavior at material interfaces requires appropriate numerical modeling and computational design techniques. However, the commonly used density-based representation in popular topology optimization methods results in a diffuse staircased boundary. The heavily refined finite element meshes required to compensate for this boundary description results in exceedingly large and expensive optimization problems. In this paper, we demonstrate the adverse effect of the density-based boundary description. Furthermore, we propose a level set-based topology optimization procedure with an enriched finite element method that shows improved performance when compared to the density-based approach.