This work focused on the modelling of latent heat thermal energy storage systems. The mathematical modelling of a melting and solidification process has time-dependent boundary conditions because the interface between solid and liquid phases is a moving boundary. The heat transfer analysis needs the interface position over time to predict the temperature inside the liquid and the solid regions. This work started by solving the classical two-phase (one-dimensional) Stefan problem through a Matlab implementation of the analytical model. The same physical problem was numerically simulated using ANSYS FLUENT, and the good match of analytical and numerical results validated the numerical model, which was used for a more interesting problem: comparing three different latent heat TES configurations during the discharging process to evaluate the most efficient in terms of maximum average discharging power. The three axial heat conduction structures changed only for the fin shape (rectangular, trapezoidal and fractal), keeping constant the volume fractions of steel, aluminium and PCM to perform a proper comparison. Results showed that the trapezoidal fin profile performs better than the rectangular one, and the fractal fin profile geometry was revealed as the best for faster thermal exchange when the solidifying frontier moves away from the steel ring. In conclusion, the average discharging power for the three configurations was evaluated for a time corresponding to a reference value (10%) of the liquid fraction: the rectangular fin profile provided 950.8 W, the trapezoidal fin profile 979.4 W and the fractal fin profile 1136.6 W, confirming its higher performance compared with the other two geometries.