Close-contact melting refers to the process of a heat source melting its way into a phase-change material. Of special interest is the close-contact melting velocity, or more specifically the relative velocity between the heat source and the phase-change material. In this work, we present a novel numerical approach to simulate quasi-steady, heat flux driven close-contact melting. It extends existing approaches found in the literature, and, for the first time, allows to study the impact of a spatially varying heat flux distribution. We will start by deriving the governing equations in a Lagrangian reference frame fixed to the heat source. Exploiting the narrowness of the melt film enables us to reduce the momentum balance to the Reynolds equation, which is coupled to the energy balance via the velocity field. We particularize our derivation for two simple, yet technically relevant geometries, namely a 3d circular disc and a 2d planar heat source. An iterative solution procedure for the coupled system is described in detail and discussed on the basis of a convergence study. Furthermore, we present an extension to allow for rotational melting modes. Various test cases demonstrate the proficiency of our method. In particular, we will utilize the method to assess the efficiency of the close-contact melting process and to quantify the model error introduced if convective losses are neglected. Finally, we will draw conclusions and present an outlook to future work.