A thermomechanical model is developed within a large deformation setting in order to simulate the interactions between martensitic phase transformations and crystalline damage growth at the austenitic grain level. Subgrain information is included in the model via the crystallographic theory of martensitic transformations. The damage and transformation characteristics are dependent of the specific martensitic transformation systems activated during a loading process, which makes the model strongly anisotropic. The state of transformation for the individual transformation systems is represented by the corresponding volume fractions. The state of damage in the austenite and in the martensitic transformation systems is reflected by the corresponding damaged volume fractions. The thermodynamical forces energetically conjugated to the rate of volume fraction and the rate of damaged volume fraction are the driving forces for transformation and crystalline damage, respectively. The expressions for these driving forces follow after constructing the specific form of the Helmholtz energy for a phase-changing, damaging material. The model is used to analyze several three-dimensional boundary value problems that are representative of microstructures appearing in multiphase carbon steels containing retained austenite. The analyses show that the incorporation of damage in the model effectively limits the elastic stresses developing in the martensitic product phase, where the maximum value of the stress strongly depends on the toughness of the martensite. Furthermore, in an aggregate of randomly oriented grains of retained austenite embedded in a ferritic matrix the generation of crystalline damage delays the phase transformation process, and may arrest it if the martensitic product phase is sufficiently brittle. The response characteristics computed with the phase-changing damage model are confirmed by experimental results.