The phase-field method is widely established for modeling crack propagation in material science. It shows good agreement with analytical solutions and is able to describe a complex fracture behavior. Nevertheless, the models are mostly introduced for homogeneous systems, and there are only a few approaches for heterogeneous systems. So models, that are able to describe the crack propagation in such systems, are highly desirable. Based on a classic crack propagation phase-field approach, existing models are discussed, and a new model is derived. The first model, which is based on an already existing approach, uses homogenized properties and a single-crack order parameter, while the second model, a new approach, introduces multiple crack order parameters, each of which only tracks the damage of a corresponding phase. Furthermore, the issues of the single-crack order parameter model are demonstrated, such as the ability to reproduce the quantitative surface energies for an interfacial crack and the non-physical behavior of a crack propagating along a sloped material interface. In contrast, the novel multi-crack order parameter model prevents distortions for an interfacial crack and propagates along the interface in a more physical manner. In comparison with an analytical solution, based on linear elastic fracture mechanics, the novel model shows a good agreement, even for strongly sloped interfaces, where the single-crack order parameter model fails to reproduce the analytical solution. In a subsequent application to fiber-reinforced polymers, the new model has proven to be able to predict fractures in complex systems, including crack nucleation, branching, and merging. Finally, the applicability to a 3D system is illustrated.