Based on the Gurson-Tvergaard-Needleman (GTN) model and Hill’s quadratic anisotropic yield criterion, a combined experimental-numerical study on fracture initiation in the process of thermal stamping of Mg alloy AZ31 sheets was carried out. The aim is to predict the formability of thermal stamping of the Mg alloy sheets at different temperatures. The presented theoretical framework was implemented into a VUMAT subroutine for ABAQUS/EXPLICIT. Internal damage evolution due to void growth and coalescence developed at different temperatures in the Mg alloy sheets was observed by scanning electron microscopy (SEM). Moreover, the thermal effects on the void growth, coalescence, and fracture behavior of the Mg alloy sheets were analyzed by the extended GTN model and forming limit diagrams (FLD). Parameters employed in the GTN model were determined from tensile tests and numerical iterative computation. The distribution of major and minor principal strains in the specimens was determined from the numerical results. Therefore, the corresponding forming limit diagrams at different stress levels and temperatures were drawn. The comparison between the predicted forming limits and the experimental data shows a good agreement.