Phenomena related to phase change heat transfer are often encountered in engineering. These phenomena are regarded to be complex, since not only phase transition from solid to liquid occurs but also movement of fluid interface has to be taken into consideration. Detailed numerical modeling of these complex systems is essential to better understand them and optimize industrial designs. Lagrangian methods are promising for simulating such complex systems. The Moving Particle Semi-implicit (MPS) method, which is one of the Lagrangian methods, is employed here to simulate the free surface fluid flows involving heat transfer and phase change. On the other hand, the existing MPS method could not apply Neumann boundary condition such as heat flux in the heat transfer simulations. This is because the surface direction could not be readily defined on the surface of the spherical fluid particles in the MPS method. Hence, prescribing the heat fluxes becomes problematic in the existing MPS method. To solve this problem, a new heat flux model is developed, where the divergence operator is applied in the heat transfer simulation. Simple verification tests are performed to demonstrate the heat flux model, where the calculation results are compared against analytically derived solutions. In addition, application of the signed distance function is also investigated in the heat transfer simulation for arbitrary shaped boundary. In simple verification tests, the computation results are shown to agree well with the analytical solutions. Consequently, adequacy of the novel heat transfer model developed here is shown in the Lagrangian fluid dynamics simulation. (C) 2016 Elsevier Ltd. All rights reserved.