Ductile fracture in metal alloys frequently entails the growth and coalescence of internal voids that originate from second-phase particles. The creation of ultra-tough ductile materials is made possible by a strong hardening capability [1] that constrains void growth and coalescence by the formation of a hardened layer at the boundary of the voids. A large number of homogenized models are available in the literature to describe porous materials. However, the majority of these rely on Gurson's heuristic approach to account for strain hardening, which is only relevant for moderate and/or saturating hardening behaviour. The present study proposes a homogenized model for the description of porous materials exhibiting strong and/or non-saturating hardening behaviour. Limit analysis is used to derive a yield criterion for spherical voids in a matrix material with an inhomogeneous yield stress. A three-parameter yield stress profile is considered, which allows recovering cases of interest as particular cases: constant yield stress and a hardened layer at the void's surface. Three distinct yield criteria are derived, corresponding to void growth/low stress triaxiality, void growth/large stress triaxiality, and void coalescence. These are then combined using a multi-surface plasticity framework. Furthermore, an evolution law is proposed for the yield stress profile, as a function of the hardening behaviour. The model is evaluated against a comprehensive database of porous unit cell FFT simulations. These simulations are conducted under both small and large strain settings, allowing for the assessment of the yield criterion and the evaluation of the model's capability to reproduce stress-strain curves and porosity evolutions. The results demonstrate a high degree of agreement, thereby validating the model. Finally, the model is employed to perform numerical simulations of ductile tearing.