Frailty models are essential tools in survival analysis for addressing unobserved heterogeneity and random effects in the data. These models incorporate a random effect, the frailty, which is assumed to impact the hazard rate multiplicatively. In this paper, we introduce a novel class of frailty models in both univariate and multivariate settings, using phase-type distributions as the underlying frailty specification. We investigate the properties of these phase-type frailty models and develop expectation-maximization algorithms for their maximum-likelihood estimation. In particular, we show that the resulting model shares similarities with the Gamma frailty model, has closed-form expressions for its functionals, and can approximate any other frailty model. Through a series of simulated and real-life numerical examples, we demonstrate the effectiveness and versatility of the proposed models in addressing unobserved heterogeneity in survival analysis.