Many important industrial processes, such as additive manufacturing, involve rapid mass, flow and heat transport between gas, liquid and solid phases. Various associated challenges, such as the large density ratio between gas and condensed phases, make accurate, robust thermal multi-phase flow simulations of these processes very difficult. In order to address some of the associated challenges, a computational framework for thermal multi-phase flows is developed based on the finite element method (FEM). A unified model for thermal multi-phase flows similar to the models widely used in the manufacturing community is adopted. The combination of the level-set method and residual-based variational multi-scale formulation (RBVMS) is used to solve the governing equations of thermal multi-phase flows. Phase transitions between solid and liquid phases, i.e., melting and solidification, are considered. Interfacial forces, including surface tension and Marangoni stress, are taken into account and handled by a density-scaled continuum surface force model. A robust fully coupled solution strategy is adopted to handle various numerical difficulties associated with thermal multi-phase flow simulations, and implemented by means of a matrix-free technique using Flexible GMRES. The mathematical formulation and its algorithmic implementation are described in detail. Four numerical test cases are presented to demonstrate the capability of the proposed formulation. The first case is a benchmark example of solidification of aluminum in a graphite mold, the second case is a thermo-capillary droplet migration problem, the third case is a spot laser melting problem, and the fourth case is the melting of metal with an interior gas bubble. The computational results are compared with analytical, experimental and simulation data from other researchers, with good agreement in cases where such data is available.